Fortran Wiki
c

! ************************************************************
	module drive    
	  integer:: ntype ! =0 for quick rel.,-1 saddle, +1 min, 2 spect
        real*8 :: step,dsig    
        real*8 ::  visc(4)   !  viscosity for marked atoms
 
    
	end module drive
! ************************************************************
	module param  ! Параметры HCNO потенциала
	  real*8 fi0(4,4),betta(4,4),R0(4,4)
	  real*8 K0(4,4),delta(4,4)
	  real*8 H(4,2),ksi(4,2)
	end module param
!c ************************************************************
	module common_var
	  real*8, parameter:: pi=3.14159265d0
	  real*8,allocatable:: x(:),y(:),z(:) ! Координаты атомов
	  real*8,allocatable:: x0(:),y0(:),z0(:) ! Координаты атомов
	  real*8,allocatable:: fx(:),fy(:),fz(:) ! Force components
	  real*8,allocatable:: fx0(:),fy0(:),fz0(:) ! Force components
	  real*8,allocatable:: felx(:),fely(:),felz(:) ! electron Force components
	  real*8,allocatable:: frepx(:),frepy(:),frepz(:) ! Repuls. Force components
      	  real*8,allocatable:: amass(:) ! masses of atoms
 	  real*8,allocatable:: Hess(:,:) ! Hessian
 	  real*8,allocatable:: Hesel(:,:),hesrep(:,:) ! Hessian
 	  real*8,allocatable:: delta3(:) ! minsad displacement
! 	  real*8,allocatable:: ak(:),at(:)  
 	  real*8,allocatable:: ddd(:),e1(:) ! spectrum attributes
	  integer,allocatable:: at_type(:) ! Типы атомов: H=1, C=2, N=3, O=4
	  integer H_count,C_count,N_count,O_count ! Кол-во атомов каждого типа
	  integer at_count ! Общее количество атомов
	  integer Nat3 ! Полное количество обобщенных координат = 3*at_count
	  integer el_count,wfunc_count ! Кол-во электронов и орбиталей
           integer:: indper
           integer:: na1,na2 
          real*8 ::   Ff12
           integer:: na3,na4
          real*8 ::   Ff34
           integer:: na5,na6
          real*8 ::   Ff56
           integer:: na7,na8
          real*8 ::   Ff78
          real*8 :: perx(3),pery(3),perz(3)



	end module common_var
!c ************************************************************
subroutine main(masx,masy,masz,mast,masn,callback)
  !DEC$ ATTRIBUTES DLLEXPORT::main
	use common_var
	use drive
        implicit none
	integer masn
	real*4:: ay,ttt
	real*8:: e,fix,om,sss,e000,e11,eold,eee,delt,ddel,dsw,sway,df
	real*8:: dx,dy,dz,r,ferr
	real*8:: masx(masn),masy(masn),masz(masn)
    integer flag,mast(masn)
        integer :: i,i11,j,ierr,n,nm,moden,k,l 
	
	interface
	subroutine callback(masx,masy,masz,masn)
	  integer masn
	  real*8:: masx(masn),masy(masn),masz(masn)
	end subroutine
	end interface


	  call read_param  
	  call read1_coord(masx,masy,masz,mast,masn)
	             
  
	  call force(e)  
	   
      i=0; flag=0

	    do while (flag.eq.0)
		 i=i+1
            call force(e)
            x=x+step*fx
            y=y+step*fy
            z=z+step*fz
           ferr=sum(dabs(fx))+sum(dabs(fy))+sum(dabs(fz))
          

 
		masx(1:masn)=x(1:masn); masy(1:masn)=y(1:masn); masz(1:masn)=z(1:masn)
			     

         if(ferr/dfloat(at_count)<1.d-2) then
		    flag=1
		  endif
		if (i.gt.10000) then
		    flag=1
		  endif	
		  call callback(masx,masy,masz,masn)  	    
      end do  

 	  
      deallocate (x,y,z,fx,fy,fz,fx0,fy0,fz0,at_type,hess,amass,ddd,e1)
      deallocate (x0,y0,z0,delta3)
      deallocate (felx,fely,felz,frepx,frepy,frepz,hesel,hesrep)
         end subroutine main

!c ****************************************************************


!c ***********************************************************************
   	subroutine ener_rep(e_rep)
!c	вычисляем классическую энергию отталкивания
	use param
	use common_var
	integer i,j
	real*8 e_rep,rij,help_var,dx,dy,dz
	e_rep=0.d0
	do i=1,at_count
	do j=i+1,at_count

        call dxdydz(dx,dy,dz,rij,i,j) 

!	 call dist(i,j,rij) ! rij - расстояние между атомами номер i и j
	 help_var=-betta(at_type(i),at_type(j))
	 help_var=help_var*(rij-R0(at_type(i),at_type(j)))
	 help_var=dexp(help_var)
	 e_rep=e_rep+fi0(at_type(i),at_type(j))*help_var
	enddo
	enddo
	end subroutine ener_rep
!c ******************************************************************
	subroutine force(e)
	use param
	use common_var
	real*8 e,erep,erepn,FFtot,dhtot
	real*8 t,k,p,p1,p2,ks,R,arg,dx,dy,dz
	real*8 hlp,hlp1,hlp2,help_var,Fxx,Fyy,Fzz,buf(3),anor
	real*8 l,m,n
	real*8 S(wfunc_count,wfunc_count)   ! Матрица перекрытия
	real*8 Ham(wfunc_count,wfunc_count) ! Гамильтониан
	real*8 Sold(wfunc_count,wfunc_count),dS(wfunc_count,wfunc_count)   
	real*8 Hamold(wfunc_count,wfunc_count),dHam(wfunc_count,wfunc_count) 
	real*8 ener(wfunc_count)            ! Собственные уровни энергии
	real*8 eifun(wfunc_count,wfunc_count) ! Eigen functions
	real*8 fmatH(wfunc_count,wfunc_count) ! Sum(eifun*eifun) over levels
	real*8 fmatS(wfunc_count,wfunc_count) ! Sum(eifun*eifun*ener) over levels
	integer i,j,i1,j1,k1,k2,el_remain,nxyz
	i1=0; j1=0

        S=0.d0;   k1=-1
!          do lll=1,2        

	i1=0; j1=0
!     Формируем матрицу S
	do i=1,at_count
	j1=0
	if (at_type(i).eq.1) i1=i1+1
	if (at_type(i).gt.1) i1=i1+4 
	do j=1,at_count
!	print*,i,j
!	print*,i1,j1
!	read*
	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  
                  if(i>j) jind=0
                  
	           if(jind==1) then     !  Blocade
	if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,i,j) 
	  l=dx/R; m=dy/R; n=dz/R 

!	  R=(x(i)-x(j))**2+(y(i)-y(j))**2+(z(i)-z(j))**2; R=dsqrt(R)
!	  l=(x(j)-x(i))/R; m=(y(j)-y(i))/R; n=(z(j)-z(i))/R 

	  if (at_type(i).eq.1) then
	    if (at_type(j).eq.1) then
	      j1=j1+1
	      
	      call In1S1S(R,S(i1,j1))
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(S(i1,j1)-Sold(i1,j1))/dhtot
!************************************************
	    endif
	    if (at_type(j).gt.1) then
	      j1=j1+4
	      call In1S2S(R,at_type(j),S(i1,j1-3))
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(S(i1,j1-3)-Sold(i1,j1-3))/dhtot
!************************************************
	      call In1S2P(R,at_type(j),hlp)
		  S(i1,j1-2)=hlp*l; S(i1,j1-1)=hlp*m; S(i1,j1)=hlp*n 		   
	    endif
	   endif
	   if (at_type(i).gt.1) then
	     if (at_type(j).eq.1) then
	       j1=j1+1
	       call In1S2S(R,at_type(i),S(i1-3,j1))
	       call In1S2P(R,at_type(i),hlp); hlp=-hlp
	       S(i1-2,j1)=hlp*l; S(i1-1,j1)=hlp*m; S(i1,j1)=hlp*n
	     endif
	     if (at_type(j).gt.1) then
	       j1=j1+4
	       call In2S2S(R,at_type(i),at_type(j), S(i1-3,j1-3))
	       call In2S2P(R,at_type(i),at_type(j),hlp)
		   S(i1-3,j1-2)=hlp*l; S(i1-3,j1-1)=hlp*m; S(i1-3,j1)=hlp*n 
	       call In2S2P(R,at_type(j),at_type(i),hlp); hlp=-hlp
		   S(i1-2,j1-3)=hlp*l; S(i1-1,j1-3)=hlp*m; S(i1,j1-3)=hlp*n 
		   call In2P2Pq(R,at_type(i),at_type(j),hlp1)
		   call In2P2Ppi(R,at_type(i),at_type(j),hlp2)
		   S(i1-2,j1-2)=hlp1*l*l+hlp2*(1.d0-l*l)
		   S(i1-2,j1-1)=hlp1*l*m-hlp2*l*m
		   S(i1-2,j1)=hlp1*l*n-hlp2*l*n
		   S(i1-1,j1-2)=hlp1*m*l-hlp2*m*l
		   S(i1-1,j1-1)=hlp1*m*m+hlp2*(1.d0-m*m)
		   S(i1-1,j1)=hlp1*m*n-hlp2*m*n
		   S(i1,j1-2)=hlp1*n*l-hlp2*n*l
		   S(i1,j1-1)=hlp1*n*m-hlp2*n*m
		   S(i1,j1)=hlp1*n*n+hlp2*(1.d0-n*n)		   		   		   	   
	     endif	      
	   endif
      endif
	if (i.eq.j) then
	  if (at_type(i).eq.1) then
	     j1=j1+1; S(i1,i1)=1.d0
	  endif
	  if (at_type(i).ne.1) then
	     j1=j1+4
		 S(i1-3,i1-3)=1.d0; S(i1-3,i1-2)=0.d0
		 S(i1-3,i1-1)=0.d0; S(i1-3,i1)=0.d0
	     S(i1-2,i1-3)=0.d0; S(i1-2,i1-2)=1.d0
		 S(i1-2,i1-1)=0.d0; S(i1-2,i1)=0.d0
	     S(i1-1,i1-3)=0.d0; S(i1-1,i1-2)=0.d0
		 S(i1-1,i1-1)=1.d0; S(i1-1,i1)=0.d0
	     S(i1,i1-3)=0.d0; S(i1,i1-2)=0.d0
		 S(i1,i1-1)=0.d0; S(i1,i1)=1.d0
	  endif
	endif
	
	
	         else     !  Blocade


	if (i.ne.j) then
 
	  if (at_type(i).eq.1) then
	    if (at_type(j).eq.1) then
	      j1=j1+1
	      
	    endif
	    if (at_type(j).gt.1) then
	      j1=j1+4
 		   
	    endif
	   endif
	   if (at_type(i).gt.1) then
	     if (at_type(j).eq.1) then
	       j1=j1+1

	     endif
	     if (at_type(j).gt.1) then
	       j1=j1+4
	     endif	      
	   endif
      endif
	if (i.eq.j) then
	  if (at_type(i).eq.1) then
	     j1=j1+1

	  endif
	  if (at_type(i).ne.1) then
	     j1=j1+4

	  endif
	endif
	         
	         
	         
	         end if    !  Blocade
	
	
     	enddo
	enddo

	
	do i1=1,wfunc_count
	do j1=i1,wfunc_count
   s(j1,i1)=s(i1,j1)   ! Symm
	end do
	end do
	
	
	
	
! 	Формируем матрицу гамильтониана Ham
 
  i1=0; j1=0
	do i=1,at_count
	  if (at_type(i).eq.1) then
	   i1=i1+1 
	   do j=1,at_count
	   
	   
	   	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  if(i>j) jind=0
                  
                  
	           if(jind==1) then     !  Blocade
	   
	   
	   
!     Водород с водородом
	      if (at_type(j).eq.1) then
	        j1=j1+1; Ham(i1,j1)=S(i1,j1)*H(1,1)

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	        R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	        R=R+(z(i)-z(j))**2; R=dsqrt(R)
	        hlp1=dexp(-delta(1,1)*(R-R0(1,1)))
	        hlp1=K0(1,1)*hlp1
	        if (i.ne.j) then
                Ham(i1,j1)=Ham(i1,j1)*hlp1
                end if 
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(Ham(i1,j1)-Hamold(i1,j1))/dhtot
!************************************************



	      endif
!     Водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4
			hlp1=(H(1,1)+H(at_type(j),1))/2.d0
			Ham(i1,j1-3)=S(i1,j1-3)*hlp1
  			hlp1=(H(1,1)+H(at_type(j),2))/2.d0
			Ham(i1,j1-2)=S(i1,j1-2)*hlp1
			Ham(i1,j1-1)=S(i1,j1-1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
	        if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(1,at_type(j))
			  hlp1=hlp1*(R-R0(1,at_type(j)))
	          hlp1=dexp(-hlp1)
	          hlp1=hlp1*K0(1,at_type(j))
			  Ham(i1,j1-3)=Ham(i1,j1-3)*hlp1
			  Ham(i1,j1-2)=Ham(i1,j1-2)*hlp1
			  Ham(i1,j1-1)=Ham(i1,j1-1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
			endif 
	      endif
	      
	      
	      else     !  Blocade
	      
	      if (at_type(j).eq.1) then
	        j1=j1+1; 

	      endif
!     Водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4

	        if (i.ne.j) then

			endif 
	      endif
	      
	      	      
	      
	      
	      
	      end if     !  Blocade
	      
	      
	      
	      
	   enddo
	  endif
	  if (at_type(i).gt.1) then
	   i1=i1+4 
	   do j=1,at_count
	   
	   	   
	   	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  if(i>j) jind=0
	           if(jind==1) then     !  Blocade
	   

	   
	   
!     Не водород с водородом
	      if (at_type(j).eq.1) then
	        j1=j1+1
	        hlp1=(H(1,1)+H(at_type(i),1))/2.d0
			Ham(i1-3,j1)=S(i1-3,j1)*hlp1
	        hlp1=(H(1,1)+H(at_type(i),2))/2.d0
			Ham(i1-2,j1)=S(i1-2,j1)*hlp1
			Ham(i1-1,j1)=S(i1-1,j1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
			if (i.ne.j) then	        

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(at_type(i),1)
	          hlp1=hlp1*(R-R0(at_type(i),1))
	          hlp1=dexp(-hlp1)
	          hlp1=K0(at_type(i),1)*hlp1
	          Ham(i1-3,j1)=Ham(i1-3,j1)*hlp1
   			  Ham(i1-2,j1)=Ham(i1-2,j1)*hlp1
			  Ham(i1-1,j1)=Ham(i1-1,j1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
	        endif
	      endif
!	Не водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4
			hlp1=(H(at_type(i),1)+H(at_type(j),1))/2.d0
			Ham(i1-3,j1-3)=S(i1-3,j1-3)*hlp1
  			hlp1=(H(at_type(i),1)+H(at_type(j),2))/2.d0
			Ham(i1-3,j1-2)=S(i1-3,j1-2)*hlp1
			Ham(i1-3,j1-1)=S(i1-3,j1-1)*hlp1
			Ham(i1-3,j1)=S(i1-3,j1)*hlp1
			hlp1=(H(at_type(i),2)+H(at_type(j),1))/2.d0
			Ham(i1-2,j1-3)=S(i1-2,j1-3)*hlp1
			Ham(i1-1,j1-3)=S(i1-1,j1-3)*hlp1
			Ham(i1,j1-3)=S(i1,j1-3)*hlp1
			hlp1=(H(at_type(i),2)+H(at_type(j),2))/2.d0
			Ham(i1-2,j1-2)=S(i1-2,j1-2)*hlp1
			Ham(i1-2,j1-1)=S(i1-2,j1-1)*hlp1
			Ham(i1-2,j1)=S(i1-2,j1)*hlp1
			Ham(i1-1,j1-2)=S(i1-1,j1-2)*hlp1
			Ham(i1-1,j1-1)=S(i1-1,j1-1)*hlp1
			Ham(i1-1,j1)=S(i1-1,j1)*hlp1
			Ham(i1,j1-2)=S(i1,j1-2)*hlp1
			Ham(i1,j1-1)=S(i1,j1-1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
	     	if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(at_type(i),at_type(j))
			  hlp1=hlp1*(R-R0(at_type(i),at_type(j)))
	          hlp1=dexp(-hlp1)
	          hlp1=hlp1*K0(at_type(i),at_type(j))
			  Ham(i1-3,j1-3)=Ham(i1-3,j1-3)*hlp1
			  Ham(i1-3,j1-2)=Ham(i1-3,j1-2)*hlp1
			  Ham(i1-3,j1-1)=Ham(i1-3,j1-1)*hlp1
			  Ham(i1-3,j1)=Ham(i1-3,j1)*hlp1
			  Ham(i1-2,j1-3)=Ham(i1-2,j1-3)*hlp1
			  Ham(i1-2,j1-2)=Ham(i1-2,j1-2)*hlp1
		 	  Ham(i1-2,j1-1)=Ham(i1-2,j1-1)*hlp1
			  Ham(i1-2,j1)=Ham(i1-2,j1)*hlp1
			  Ham(i1-1,j1-3)=Ham(i1-1,j1-3)*hlp1
			  Ham(i1-1,j1-2)=Ham(i1-1,j1-2)*hlp1
			  Ham(i1-1,j1-1)=Ham(i1-1,j1-1)*hlp1
			  Ham(i1-1,j1)=Ham(i1-1,j1)*hlp1
			  Ham(i1,j1-3)=Ham(i1,j1-3)*hlp1
			  Ham(i1,j1-2)=Ham(i1,j1-2)*hlp1
			  Ham(i1,j1-1)=Ham(i1,j1-1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
			endif 
	      endif
	      
	      
	      
	      	        else     !  Blocade
	      
        	        
	      if (at_type(j).eq.1) then
	        j1=j1+1

			if (i.ne.j) then	        

	        endif
	      endif
!	Не водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4

	     	if (i.ne.j) then

			endif 
	      endif
	      	      	        
	      	        
	      	        
	      	        end if   !  Blocade
	      
	      
	      
	      
	   enddo
	  endif
	  j1=0
	  enddo
	  
 	do i1=1,wfunc_count
	do j1=i1,wfunc_count
   Ham(j1,i1)=Ham(i1,j1)   !Symm
	end do
	end do
	  
!          end do             
       sold=s; hamold=ham

!c	Находим собственные значения энергии без собственных векторов:
!c     Ham - гамильтониан, S - матрица перекрытия, wfunc_count - их размерность
!c     ener - получающийся массив собственных чисел,   eifun - eigenvector

         call enonort(ener,ham,s,eifun)
!         call enonort(ener,ham,s,eifun)
          
!c     Расселяем электроны по уровням

!          do lll=1,2
	call ener_rep(erep) ! вычислили классическую энергию отталкивания
!	        end do
  
           	e=erep
           	fmatH=0.d0
           	fmatS=0.d0
	el_remain=el_count; i=0
	do 

       if (el_remain.eq.0) exit
	  i=i+1
	  if (el_remain.eq.1) then
          e=e+ener(i)
          do k1=1,wfunc_count
          do k2=1,wfunc_count
          fmatH(k1,k2)=fmatH(k1,k2)+eifun(k1,i)*eifun(k2,i)
          fmatS(k1,k2)=fmatS(k1,k2)+eifun(k1,i)*eifun(k2,i)*ener(i)
          end do
          end do
          exit 
	  endif
	  e=e+2.d0*ener(i)
          do k1=1,wfunc_count
          do k2=1,wfunc_count
          fmatH(k1,k2)=fmatH(k1,k2)+2.d0*eifun(k1,i)*eifun(k2,i)
          fmatS(k1,k2)=fmatS(k1,k2)+2.d0*eifun(k1,i)*eifun(k2,i)*ener(i)
          end do
          end do
    el_remain=el_remain-2
	end do
  e=e-dfloat(H_count)*H(1,1)                    ! минус водород
	e=e-dfloat(C_count)*(2.d0*H(2,1)+2.d0*H(2,2)) ! минус углерод
	e=e-dfloat(N_count)*(2.d0*H(3,1)+3.d0*H(3,2)) ! минус азот
	e=e-dfloat(O_count)*(2.d0*H(4,1)+4.d0*H(4,2)) ! минус кислород
	
!***************  Force calculation  *******************
         dhtot=1.d-7
         frepx=0.d0;          frepy=0.d0;         frepz=0.d0;
         do k1=1,at_count

         S=Sold; Ham=Hamold

         !********* Fx determination **********
         x(k1)=x(k1)+dhtot
         

	i1=0; j1=0
!     Формируем матрицу S
	do i=1,at_count
	j1=0
	if (at_type(i).eq.1) i1=i1+1
	if (at_type(i).gt.1) i1=i1+4 
	do j=1,at_count
!	print*,i,j
!	print*,i1,j1
!	read*
	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  
                  if(i>j) jind=0
                  
	           if(jind==1) then     !  Blocade
	if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,i,j) 
	  l=dx/R; m=dy/R; n=dz/R 

!	  R=(x(i)-x(j))**2+(y(i)-y(j))**2+(z(i)-z(j))**2; R=dsqrt(R)
!	  l=(x(j)-x(i))/R; m=(y(j)-y(i))/R; n=(z(j)-z(i))/R 

	  if (at_type(i).eq.1) then
	    if (at_type(j).eq.1) then
	      j1=j1+1
	      
	      call In1S1S(R,S(i1,j1))
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(S(i1,j1)-Sold(i1,j1))/dhtot
!************************************************
	    endif
	    if (at_type(j).gt.1) then
	      j1=j1+4
	      call In1S2S(R,at_type(j),S(i1,j1-3))
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(S(i1,j1-3)-Sold(i1,j1-3))/dhtot
!************************************************
	      call In1S2P(R,at_type(j),hlp)
		  S(i1,j1-2)=hlp*l; S(i1,j1-1)=hlp*m; S(i1,j1)=hlp*n 		   
	    endif
	   endif
	   if (at_type(i).gt.1) then
	     if (at_type(j).eq.1) then
	       j1=j1+1
	       call In1S2S(R,at_type(i),S(i1-3,j1))
	       call In1S2P(R,at_type(i),hlp); hlp=-hlp
	       S(i1-2,j1)=hlp*l; S(i1-1,j1)=hlp*m; S(i1,j1)=hlp*n
	     endif
	     if (at_type(j).gt.1) then
	       j1=j1+4
	       call In2S2S(R,at_type(i),at_type(j), S(i1-3,j1-3))
	       call In2S2P(R,at_type(i),at_type(j),hlp)
		   S(i1-3,j1-2)=hlp*l; S(i1-3,j1-1)=hlp*m; S(i1-3,j1)=hlp*n 
	       call In2S2P(R,at_type(j),at_type(i),hlp); hlp=-hlp
		   S(i1-2,j1-3)=hlp*l; S(i1-1,j1-3)=hlp*m; S(i1,j1-3)=hlp*n 
		   call In2P2Pq(R,at_type(i),at_type(j),hlp1)
		   call In2P2Ppi(R,at_type(i),at_type(j),hlp2)
		   S(i1-2,j1-2)=hlp1*l*l+hlp2*(1.d0-l*l)
		   S(i1-2,j1-1)=hlp1*l*m-hlp2*l*m
		   S(i1-2,j1)=hlp1*l*n-hlp2*l*n
		   S(i1-1,j1-2)=hlp1*m*l-hlp2*m*l
		   S(i1-1,j1-1)=hlp1*m*m+hlp2*(1.d0-m*m)
		   S(i1-1,j1)=hlp1*m*n-hlp2*m*n
		   S(i1,j1-2)=hlp1*n*l-hlp2*n*l
		   S(i1,j1-1)=hlp1*n*m-hlp2*n*m
		   S(i1,j1)=hlp1*n*n+hlp2*(1.d0-n*n)		   		   		   	   
	     endif	      
	   endif
      endif
	if (i.eq.j) then
	  if (at_type(i).eq.1) then
	     j1=j1+1; S(i1,i1)=1.d0
	  endif
	  if (at_type(i).ne.1) then
	     j1=j1+4
		 S(i1-3,i1-3)=1.d0; S(i1-3,i1-2)=0.d0
		 S(i1-3,i1-1)=0.d0; S(i1-3,i1)=0.d0
	     S(i1-2,i1-3)=0.d0; S(i1-2,i1-2)=1.d0
		 S(i1-2,i1-1)=0.d0; S(i1-2,i1)=0.d0
	     S(i1-1,i1-3)=0.d0; S(i1-1,i1-2)=0.d0
		 S(i1-1,i1-1)=1.d0; S(i1-1,i1)=0.d0
	     S(i1,i1-3)=0.d0; S(i1,i1-2)=0.d0
		 S(i1,i1-1)=0.d0; S(i1,i1)=1.d0
	  endif
	endif
	
	
	         else     !  Blocade


	if (i.ne.j) then
 
	  if (at_type(i).eq.1) then
	    if (at_type(j).eq.1) then
	      j1=j1+1
	      
	    endif
	    if (at_type(j).gt.1) then
	      j1=j1+4
 		   
	    endif
	   endif
	   if (at_type(i).gt.1) then
	     if (at_type(j).eq.1) then
	       j1=j1+1

	     endif
	     if (at_type(j).gt.1) then
	       j1=j1+4
	     endif	      
	   endif
      endif
	if (i.eq.j) then
	  if (at_type(i).eq.1) then
	     j1=j1+1

	  endif
	  if (at_type(i).ne.1) then
	     j1=j1+4

	  endif
	endif
	         
	         
	         
	         end if    !  Blocade
	
	
     	enddo
	enddo

	
	do i1=1,wfunc_count
	do j1=i1,wfunc_count
   s(j1,i1)=s(i1,j1)   ! Symm
	end do
	end do
	
	
	
	
! 	Формируем матрицу гамильтониана Ham
 
  i1=0; j1=0
	do i=1,at_count
	  if (at_type(i).eq.1) then
	   i1=i1+1 
	   do j=1,at_count
	   
	   
	   	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  if(i>j) jind=0
                  
                  
	           if(jind==1) then     !  Blocade
	   
	   
	   
!     Водород с водородом
	      if (at_type(j).eq.1) then
	        j1=j1+1; Ham(i1,j1)=S(i1,j1)*H(1,1)

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	        R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	        R=R+(z(i)-z(j))**2; R=dsqrt(R)
	        hlp1=dexp(-delta(1,1)*(R-R0(1,1)))
	        hlp1=K0(1,1)*hlp1
	        if (i.ne.j) then
                Ham(i1,j1)=Ham(i1,j1)*hlp1
                end if 
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(Ham(i1,j1)-Hamold(i1,j1))/dhtot
!************************************************



	      endif
!     Водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4
			hlp1=(H(1,1)+H(at_type(j),1))/2.d0
			Ham(i1,j1-3)=S(i1,j1-3)*hlp1
  			hlp1=(H(1,1)+H(at_type(j),2))/2.d0
			Ham(i1,j1-2)=S(i1,j1-2)*hlp1
			Ham(i1,j1-1)=S(i1,j1-1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
	        if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(1,at_type(j))
			  hlp1=hlp1*(R-R0(1,at_type(j)))
	          hlp1=dexp(-hlp1)
	          hlp1=hlp1*K0(1,at_type(j))
			  Ham(i1,j1-3)=Ham(i1,j1-3)*hlp1
			  Ham(i1,j1-2)=Ham(i1,j1-2)*hlp1
			  Ham(i1,j1-1)=Ham(i1,j1-1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
			endif 
	      endif
	      
	      
	      else     !  Blocade
	      
	      if (at_type(j).eq.1) then
	        j1=j1+1; 

	      endif
!     Водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4

	        if (i.ne.j) then

			endif 
	      endif
	      
	      	      
	      
	      
	      
	      end if     !  Blocade
	      
	      
	      
	      
	   enddo
	  endif
	  if (at_type(i).gt.1) then
	   i1=i1+4 
	   do j=1,at_count
	   
	   	   
	   	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  if(i>j) jind=0
	           if(jind==1) then     !  Blocade
	   

	   
	   
!     Не водород с водородом
	      if (at_type(j).eq.1) then
	        j1=j1+1
	        hlp1=(H(1,1)+H(at_type(i),1))/2.d0
			Ham(i1-3,j1)=S(i1-3,j1)*hlp1
	        hlp1=(H(1,1)+H(at_type(i),2))/2.d0
			Ham(i1-2,j1)=S(i1-2,j1)*hlp1
			Ham(i1-1,j1)=S(i1-1,j1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
			if (i.ne.j) then	        

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(at_type(i),1)
	          hlp1=hlp1*(R-R0(at_type(i),1))
	          hlp1=dexp(-hlp1)
	          hlp1=K0(at_type(i),1)*hlp1
	          Ham(i1-3,j1)=Ham(i1-3,j1)*hlp1
   			  Ham(i1-2,j1)=Ham(i1-2,j1)*hlp1
			  Ham(i1-1,j1)=Ham(i1-1,j1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
	        endif
	      endif
!	Не водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4
			hlp1=(H(at_type(i),1)+H(at_type(j),1))/2.d0
			Ham(i1-3,j1-3)=S(i1-3,j1-3)*hlp1
  			hlp1=(H(at_type(i),1)+H(at_type(j),2))/2.d0
			Ham(i1-3,j1-2)=S(i1-3,j1-2)*hlp1
			Ham(i1-3,j1-1)=S(i1-3,j1-1)*hlp1
			Ham(i1-3,j1)=S(i1-3,j1)*hlp1
			hlp1=(H(at_type(i),2)+H(at_type(j),1))/2.d0
			Ham(i1-2,j1-3)=S(i1-2,j1-3)*hlp1
			Ham(i1-1,j1-3)=S(i1-1,j1-3)*hlp1
			Ham(i1,j1-3)=S(i1,j1-3)*hlp1
			hlp1=(H(at_type(i),2)+H(at_type(j),2))/2.d0
			Ham(i1-2,j1-2)=S(i1-2,j1-2)*hlp1
			Ham(i1-2,j1-1)=S(i1-2,j1-1)*hlp1
			Ham(i1-2,j1)=S(i1-2,j1)*hlp1
			Ham(i1-1,j1-2)=S(i1-1,j1-2)*hlp1
			Ham(i1-1,j1-1)=S(i1-1,j1-1)*hlp1
			Ham(i1-1,j1)=S(i1-1,j1)*hlp1
			Ham(i1,j1-2)=S(i1,j1-2)*hlp1
			Ham(i1,j1-1)=S(i1,j1-1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
	     	if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(at_type(i),at_type(j))
			  hlp1=hlp1*(R-R0(at_type(i),at_type(j)))
	          hlp1=dexp(-hlp1)
	          hlp1=hlp1*K0(at_type(i),at_type(j))
			  Ham(i1-3,j1-3)=Ham(i1-3,j1-3)*hlp1
			  Ham(i1-3,j1-2)=Ham(i1-3,j1-2)*hlp1
			  Ham(i1-3,j1-1)=Ham(i1-3,j1-1)*hlp1
			  Ham(i1-3,j1)=Ham(i1-3,j1)*hlp1
			  Ham(i1-2,j1-3)=Ham(i1-2,j1-3)*hlp1
			  Ham(i1-2,j1-2)=Ham(i1-2,j1-2)*hlp1
		 	  Ham(i1-2,j1-1)=Ham(i1-2,j1-1)*hlp1
			  Ham(i1-2,j1)=Ham(i1-2,j1)*hlp1
			  Ham(i1-1,j1-3)=Ham(i1-1,j1-3)*hlp1
			  Ham(i1-1,j1-2)=Ham(i1-1,j1-2)*hlp1
			  Ham(i1-1,j1-1)=Ham(i1-1,j1-1)*hlp1
			  Ham(i1-1,j1)=Ham(i1-1,j1)*hlp1
			  Ham(i1,j1-3)=Ham(i1,j1-3)*hlp1
			  Ham(i1,j1-2)=Ham(i1,j1-2)*hlp1
			  Ham(i1,j1-1)=Ham(i1,j1-1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
			endif 
	      endif
	      
	      
	      
	      	        else     !  Blocade
	      
        	        
	      if (at_type(j).eq.1) then
	        j1=j1+1

			if (i.ne.j) then	        

	        endif
	      endif
!	Не водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4

	     	if (i.ne.j) then

			endif 
	      endif
	      	      	        
	      	        
	      	        
	      	        end if   !  Blocade
	      
	      
	      
	      
	   enddo
	  endif
	  j1=0
	  enddo
	  
 	do i1=1,wfunc_count
	do j1=i1,wfunc_count
   Ham(j1,i1)=Ham(i1,j1)   !Symm
	end do
	end do
	  
               ds=(s-sold)/dhtot; dham=(ham-hamold)/dhtot
               nxyz=1
	       FFtot=-sum(dHam*fmatH)+sum(dS*fmatS) !+erep-erepn

         felx(k1)=FFtot
         x(k1)=x(k1)-dhtot
         !********* End Fx determination **********

         !********* Fy determination **********
         y(k1)=y(k1)+dhtot
	                

	i1=0; j1=0
!     Формируем матрицу S
	do i=1,at_count
	j1=0
	if (at_type(i).eq.1) i1=i1+1
	if (at_type(i).gt.1) i1=i1+4 
	do j=1,at_count
!	print*,i,j
!	print*,i1,j1
!	read*
	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  
                  if(i>j) jind=0
                  
	           if(jind==1) then     !  Blocade
	if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,i,j) 
	  l=dx/R; m=dy/R; n=dz/R 

!	  R=(x(i)-x(j))**2+(y(i)-y(j))**2+(z(i)-z(j))**2; R=dsqrt(R)
!	  l=(x(j)-x(i))/R; m=(y(j)-y(i))/R; n=(z(j)-z(i))/R 

	  if (at_type(i).eq.1) then
	    if (at_type(j).eq.1) then
	      j1=j1+1
	      
	      call In1S1S(R,S(i1,j1))
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(S(i1,j1)-Sold(i1,j1))/dhtot
!************************************************
	    endif
	    if (at_type(j).gt.1) then
	      j1=j1+4
	      call In1S2S(R,at_type(j),S(i1,j1-3))
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(S(i1,j1-3)-Sold(i1,j1-3))/dhtot
!************************************************
	      call In1S2P(R,at_type(j),hlp)
		  S(i1,j1-2)=hlp*l; S(i1,j1-1)=hlp*m; S(i1,j1)=hlp*n 		   
	    endif
	   endif
	   if (at_type(i).gt.1) then
	     if (at_type(j).eq.1) then
	       j1=j1+1
	       call In1S2S(R,at_type(i),S(i1-3,j1))
	       call In1S2P(R,at_type(i),hlp); hlp=-hlp
	       S(i1-2,j1)=hlp*l; S(i1-1,j1)=hlp*m; S(i1,j1)=hlp*n
	     endif
	     if (at_type(j).gt.1) then
	       j1=j1+4
	       call In2S2S(R,at_type(i),at_type(j), S(i1-3,j1-3))
	       call In2S2P(R,at_type(i),at_type(j),hlp)
		   S(i1-3,j1-2)=hlp*l; S(i1-3,j1-1)=hlp*m; S(i1-3,j1)=hlp*n 
	       call In2S2P(R,at_type(j),at_type(i),hlp); hlp=-hlp
		   S(i1-2,j1-3)=hlp*l; S(i1-1,j1-3)=hlp*m; S(i1,j1-3)=hlp*n 
		   call In2P2Pq(R,at_type(i),at_type(j),hlp1)
		   call In2P2Ppi(R,at_type(i),at_type(j),hlp2)
		   S(i1-2,j1-2)=hlp1*l*l+hlp2*(1.d0-l*l)
		   S(i1-2,j1-1)=hlp1*l*m-hlp2*l*m
		   S(i1-2,j1)=hlp1*l*n-hlp2*l*n
		   S(i1-1,j1-2)=hlp1*m*l-hlp2*m*l
		   S(i1-1,j1-1)=hlp1*m*m+hlp2*(1.d0-m*m)
		   S(i1-1,j1)=hlp1*m*n-hlp2*m*n
		   S(i1,j1-2)=hlp1*n*l-hlp2*n*l
		   S(i1,j1-1)=hlp1*n*m-hlp2*n*m
		   S(i1,j1)=hlp1*n*n+hlp2*(1.d0-n*n)		   		   		   	   
	     endif	      
	   endif
      endif
	if (i.eq.j) then
	  if (at_type(i).eq.1) then
	     j1=j1+1; S(i1,i1)=1.d0
	  endif
	  if (at_type(i).ne.1) then
	     j1=j1+4
		 S(i1-3,i1-3)=1.d0; S(i1-3,i1-2)=0.d0
		 S(i1-3,i1-1)=0.d0; S(i1-3,i1)=0.d0
	     S(i1-2,i1-3)=0.d0; S(i1-2,i1-2)=1.d0
		 S(i1-2,i1-1)=0.d0; S(i1-2,i1)=0.d0
	     S(i1-1,i1-3)=0.d0; S(i1-1,i1-2)=0.d0
		 S(i1-1,i1-1)=1.d0; S(i1-1,i1)=0.d0
	     S(i1,i1-3)=0.d0; S(i1,i1-2)=0.d0
		 S(i1,i1-1)=0.d0; S(i1,i1)=1.d0
	  endif
	endif
	
	
	         else     !  Blocade


	if (i.ne.j) then
 
	  if (at_type(i).eq.1) then
	    if (at_type(j).eq.1) then
	      j1=j1+1
	      
	    endif
	    if (at_type(j).gt.1) then
	      j1=j1+4
 		   
	    endif
	   endif
	   if (at_type(i).gt.1) then
	     if (at_type(j).eq.1) then
	       j1=j1+1

	     endif
	     if (at_type(j).gt.1) then
	       j1=j1+4
	     endif	      
	   endif
      endif
	if (i.eq.j) then
	  if (at_type(i).eq.1) then
	     j1=j1+1

	  endif
	  if (at_type(i).ne.1) then
	     j1=j1+4

	  endif
	endif
	         
	         
	         
	         end if    !  Blocade
	
	
     	enddo
	enddo

	
	do i1=1,wfunc_count
	do j1=i1,wfunc_count
   s(j1,i1)=s(i1,j1)   ! Symm
	end do
	end do
	
	
	
	
! 	Формируем матрицу гамильтониана Ham
 
  i1=0; j1=0
	do i=1,at_count
	  if (at_type(i).eq.1) then
	   i1=i1+1 
	   do j=1,at_count
	   
	   
	   	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  if(i>j) jind=0
                  
                  
	           if(jind==1) then     !  Blocade
	   
	   
	   
!     Водород с водородом
	      if (at_type(j).eq.1) then
	        j1=j1+1; Ham(i1,j1)=S(i1,j1)*H(1,1)

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	        R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	        R=R+(z(i)-z(j))**2; R=dsqrt(R)
	        hlp1=dexp(-delta(1,1)*(R-R0(1,1)))
	        hlp1=K0(1,1)*hlp1
	        if (i.ne.j) then
                Ham(i1,j1)=Ham(i1,j1)*hlp1
                end if 
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(Ham(i1,j1)-Hamold(i1,j1))/dhtot
!************************************************



	      endif
!     Водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4
			hlp1=(H(1,1)+H(at_type(j),1))/2.d0
			Ham(i1,j1-3)=S(i1,j1-3)*hlp1
  			hlp1=(H(1,1)+H(at_type(j),2))/2.d0
			Ham(i1,j1-2)=S(i1,j1-2)*hlp1
			Ham(i1,j1-1)=S(i1,j1-1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
	        if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(1,at_type(j))
			  hlp1=hlp1*(R-R0(1,at_type(j)))
	          hlp1=dexp(-hlp1)
	          hlp1=hlp1*K0(1,at_type(j))
			  Ham(i1,j1-3)=Ham(i1,j1-3)*hlp1
			  Ham(i1,j1-2)=Ham(i1,j1-2)*hlp1
			  Ham(i1,j1-1)=Ham(i1,j1-1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
			endif 
	      endif
	      
	      
	      else     !  Blocade
	      
	      if (at_type(j).eq.1) then
	        j1=j1+1; 

	      endif
!     Водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4

	        if (i.ne.j) then

			endif 
	      endif
	      
	      	      
	      
	      
	      
	      end if     !  Blocade
	      
	      
	      
	      
	   enddo
	  endif
	  if (at_type(i).gt.1) then
	   i1=i1+4 
	   do j=1,at_count
	   
	   	   
	   	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  if(i>j) jind=0
	           if(jind==1) then     !  Blocade
	   

	   
	   
!     Не водород с водородом
	      if (at_type(j).eq.1) then
	        j1=j1+1
	        hlp1=(H(1,1)+H(at_type(i),1))/2.d0
			Ham(i1-3,j1)=S(i1-3,j1)*hlp1
	        hlp1=(H(1,1)+H(at_type(i),2))/2.d0
			Ham(i1-2,j1)=S(i1-2,j1)*hlp1
			Ham(i1-1,j1)=S(i1-1,j1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
			if (i.ne.j) then	        

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(at_type(i),1)
	          hlp1=hlp1*(R-R0(at_type(i),1))
	          hlp1=dexp(-hlp1)
	          hlp1=K0(at_type(i),1)*hlp1
	          Ham(i1-3,j1)=Ham(i1-3,j1)*hlp1
   			  Ham(i1-2,j1)=Ham(i1-2,j1)*hlp1
			  Ham(i1-1,j1)=Ham(i1-1,j1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
	        endif
	      endif
!	Не водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4
			hlp1=(H(at_type(i),1)+H(at_type(j),1))/2.d0
			Ham(i1-3,j1-3)=S(i1-3,j1-3)*hlp1
  			hlp1=(H(at_type(i),1)+H(at_type(j),2))/2.d0
			Ham(i1-3,j1-2)=S(i1-3,j1-2)*hlp1
			Ham(i1-3,j1-1)=S(i1-3,j1-1)*hlp1
			Ham(i1-3,j1)=S(i1-3,j1)*hlp1
			hlp1=(H(at_type(i),2)+H(at_type(j),1))/2.d0
			Ham(i1-2,j1-3)=S(i1-2,j1-3)*hlp1
			Ham(i1-1,j1-3)=S(i1-1,j1-3)*hlp1
			Ham(i1,j1-3)=S(i1,j1-3)*hlp1
			hlp1=(H(at_type(i),2)+H(at_type(j),2))/2.d0
			Ham(i1-2,j1-2)=S(i1-2,j1-2)*hlp1
			Ham(i1-2,j1-1)=S(i1-2,j1-1)*hlp1
			Ham(i1-2,j1)=S(i1-2,j1)*hlp1
			Ham(i1-1,j1-2)=S(i1-1,j1-2)*hlp1
			Ham(i1-1,j1-1)=S(i1-1,j1-1)*hlp1
			Ham(i1-1,j1)=S(i1-1,j1)*hlp1
			Ham(i1,j1-2)=S(i1,j1-2)*hlp1
			Ham(i1,j1-1)=S(i1,j1-1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
	     	if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(at_type(i),at_type(j))
			  hlp1=hlp1*(R-R0(at_type(i),at_type(j)))
	          hlp1=dexp(-hlp1)
	          hlp1=hlp1*K0(at_type(i),at_type(j))
			  Ham(i1-3,j1-3)=Ham(i1-3,j1-3)*hlp1
			  Ham(i1-3,j1-2)=Ham(i1-3,j1-2)*hlp1
			  Ham(i1-3,j1-1)=Ham(i1-3,j1-1)*hlp1
			  Ham(i1-3,j1)=Ham(i1-3,j1)*hlp1
			  Ham(i1-2,j1-3)=Ham(i1-2,j1-3)*hlp1
			  Ham(i1-2,j1-2)=Ham(i1-2,j1-2)*hlp1
		 	  Ham(i1-2,j1-1)=Ham(i1-2,j1-1)*hlp1
			  Ham(i1-2,j1)=Ham(i1-2,j1)*hlp1
			  Ham(i1-1,j1-3)=Ham(i1-1,j1-3)*hlp1
			  Ham(i1-1,j1-2)=Ham(i1-1,j1-2)*hlp1
			  Ham(i1-1,j1-1)=Ham(i1-1,j1-1)*hlp1
			  Ham(i1-1,j1)=Ham(i1-1,j1)*hlp1
			  Ham(i1,j1-3)=Ham(i1,j1-3)*hlp1
			  Ham(i1,j1-2)=Ham(i1,j1-2)*hlp1
			  Ham(i1,j1-1)=Ham(i1,j1-1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
			endif 
	      endif
	      
	      
	      
	      	        else     !  Blocade
	      
        	        
	      if (at_type(j).eq.1) then
	        j1=j1+1

			if (i.ne.j) then	        

	        endif
	      endif
!	Не водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4

	     	if (i.ne.j) then

			endif 
	      endif
	      	      	        
	      	        
	      	        
	      	        end if   !  Blocade
	      
	      
	      
	      
	   enddo
	  endif
	  j1=0
	  enddo
	  
 	do i1=1,wfunc_count
	do j1=i1,wfunc_count
   Ham(j1,i1)=Ham(i1,j1)   !Symm
	end do
	end do
	  
               ds=(s-sold)/dhtot; dham=(ham-hamold)/dhtot
               nxyz=2
           	       

	       FFtot=-sum(dHam*fmatH)+sum(dS*fmatS) !+erep-erepn
         fely(k1)=FFtot
         y(k1)=y(k1)-dhtot
         !********* End Fy determination **********
         
         !********* Fz determination **********
         z(k1)=z(k1)+dhtot
	               

	i1=0; j1=0
!     Формируем матрицу S
	do i=1,at_count
	j1=0
	if (at_type(i).eq.1) i1=i1+1
	if (at_type(i).gt.1) i1=i1+4 
	do j=1,at_count
!	print*,i,j
!	print*,i1,j1
!	read*
	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  
                  if(i>j) jind=0
                  
	           if(jind==1) then     !  Blocade
	if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,i,j) 
	  l=dx/R; m=dy/R; n=dz/R 

!	  R=(x(i)-x(j))**2+(y(i)-y(j))**2+(z(i)-z(j))**2; R=dsqrt(R)
!	  l=(x(j)-x(i))/R; m=(y(j)-y(i))/R; n=(z(j)-z(i))/R 

	  if (at_type(i).eq.1) then
	    if (at_type(j).eq.1) then
	      j1=j1+1
	      
	      call In1S1S(R,S(i1,j1))
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(S(i1,j1)-Sold(i1,j1))/dhtot
!************************************************
	    endif
	    if (at_type(j).gt.1) then
	      j1=j1+4
	      call In1S2S(R,at_type(j),S(i1,j1-3))
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(S(i1,j1-3)-Sold(i1,j1-3))/dhtot
!************************************************
	      call In1S2P(R,at_type(j),hlp)
		  S(i1,j1-2)=hlp*l; S(i1,j1-1)=hlp*m; S(i1,j1)=hlp*n 		   
	    endif
	   endif
	   if (at_type(i).gt.1) then
	     if (at_type(j).eq.1) then
	       j1=j1+1
	       call In1S2S(R,at_type(i),S(i1-3,j1))
	       call In1S2P(R,at_type(i),hlp); hlp=-hlp
	       S(i1-2,j1)=hlp*l; S(i1-1,j1)=hlp*m; S(i1,j1)=hlp*n
	     endif
	     if (at_type(j).gt.1) then
	       j1=j1+4
	       call In2S2S(R,at_type(i),at_type(j), S(i1-3,j1-3))
	       call In2S2P(R,at_type(i),at_type(j),hlp)
		   S(i1-3,j1-2)=hlp*l; S(i1-3,j1-1)=hlp*m; S(i1-3,j1)=hlp*n 
	       call In2S2P(R,at_type(j),at_type(i),hlp); hlp=-hlp
		   S(i1-2,j1-3)=hlp*l; S(i1-1,j1-3)=hlp*m; S(i1,j1-3)=hlp*n 
		   call In2P2Pq(R,at_type(i),at_type(j),hlp1)
		   call In2P2Ppi(R,at_type(i),at_type(j),hlp2)
		   S(i1-2,j1-2)=hlp1*l*l+hlp2*(1.d0-l*l)
		   S(i1-2,j1-1)=hlp1*l*m-hlp2*l*m
		   S(i1-2,j1)=hlp1*l*n-hlp2*l*n
		   S(i1-1,j1-2)=hlp1*m*l-hlp2*m*l
		   S(i1-1,j1-1)=hlp1*m*m+hlp2*(1.d0-m*m)
		   S(i1-1,j1)=hlp1*m*n-hlp2*m*n
		   S(i1,j1-2)=hlp1*n*l-hlp2*n*l
		   S(i1,j1-1)=hlp1*n*m-hlp2*n*m
		   S(i1,j1)=hlp1*n*n+hlp2*(1.d0-n*n)		   		   		   	   
	     endif	      
	   endif
      endif
	if (i.eq.j) then
	  if (at_type(i).eq.1) then
	     j1=j1+1; S(i1,i1)=1.d0
	  endif
	  if (at_type(i).ne.1) then
	     j1=j1+4
		 S(i1-3,i1-3)=1.d0; S(i1-3,i1-2)=0.d0
		 S(i1-3,i1-1)=0.d0; S(i1-3,i1)=0.d0
	     S(i1-2,i1-3)=0.d0; S(i1-2,i1-2)=1.d0
		 S(i1-2,i1-1)=0.d0; S(i1-2,i1)=0.d0
	     S(i1-1,i1-3)=0.d0; S(i1-1,i1-2)=0.d0
		 S(i1-1,i1-1)=1.d0; S(i1-1,i1)=0.d0
	     S(i1,i1-3)=0.d0; S(i1,i1-2)=0.d0
		 S(i1,i1-1)=0.d0; S(i1,i1)=1.d0
	  endif
	endif
	
	
	         else     !  Blocade


	if (i.ne.j) then
 
	  if (at_type(i).eq.1) then
	    if (at_type(j).eq.1) then
	      j1=j1+1
	      
	    endif
	    if (at_type(j).gt.1) then
	      j1=j1+4
 		   
	    endif
	   endif
	   if (at_type(i).gt.1) then
	     if (at_type(j).eq.1) then
	       j1=j1+1

	     endif
	     if (at_type(j).gt.1) then
	       j1=j1+4
	     endif	      
	   endif
      endif
	if (i.eq.j) then
	  if (at_type(i).eq.1) then
	     j1=j1+1

	  endif
	  if (at_type(i).ne.1) then
	     j1=j1+4

	  endif
	endif
	         
	         
	         
	         end if    !  Blocade
	
	
     	enddo
	enddo

	
	do i1=1,wfunc_count
	do j1=i1,wfunc_count
   s(j1,i1)=s(i1,j1)   ! Symm
	end do
	end do
	
	
	
	
! 	Формируем матрицу гамильтониана Ham
 
  i1=0; j1=0
	do i=1,at_count
	  if (at_type(i).eq.1) then
	   i1=i1+1 
	   do j=1,at_count
	   
	   
	   	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  if(i>j) jind=0
                  
                  
	           if(jind==1) then     !  Blocade
	   
	   
	   
!     Водород с водородом
	      if (at_type(j).eq.1) then
	        j1=j1+1; Ham(i1,j1)=S(i1,j1)*H(1,1)

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	        R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	        R=R+(z(i)-z(j))**2; R=dsqrt(R)
	        hlp1=dexp(-delta(1,1)*(R-R0(1,1)))
	        hlp1=K0(1,1)*hlp1
	        if (i.ne.j) then
                Ham(i1,j1)=Ham(i1,j1)*hlp1
                end if 
!*********************   TEST   *****************
!              print*,'test1 ',i1,j1,(Ham(i1,j1)-Hamold(i1,j1))/dhtot
!************************************************



	      endif
!     Водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4
			hlp1=(H(1,1)+H(at_type(j),1))/2.d0
			Ham(i1,j1-3)=S(i1,j1-3)*hlp1
  			hlp1=(H(1,1)+H(at_type(j),2))/2.d0
			Ham(i1,j1-2)=S(i1,j1-2)*hlp1
			Ham(i1,j1-1)=S(i1,j1-1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
	        if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(1,at_type(j))
			  hlp1=hlp1*(R-R0(1,at_type(j)))
	          hlp1=dexp(-hlp1)
	          hlp1=hlp1*K0(1,at_type(j))
			  Ham(i1,j1-3)=Ham(i1,j1-3)*hlp1
			  Ham(i1,j1-2)=Ham(i1,j1-2)*hlp1
			  Ham(i1,j1-1)=Ham(i1,j1-1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
			endif 
	      endif
	      
	      
	      else     !  Blocade
	      
	      if (at_type(j).eq.1) then
	        j1=j1+1; 

	      endif
!     Водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4

	        if (i.ne.j) then

			endif 
	      endif
	      
	      	      
	      
	      
	      
	      end if     !  Blocade
	      
	      
	      
	      
	   enddo
	  endif
	  if (at_type(i).gt.1) then
	   i1=i1+4 
	   do j=1,at_count
	   
	   	   
	   	                jind=0
	                if(k1<0) then 
                  jind=1
                  else
                  if(k1==i.or.k1==j) jind=1
                  end if
                  if(i>j) jind=0
	           if(jind==1) then     !  Blocade
	   

	   
	   
!     Не водород с водородом
	      if (at_type(j).eq.1) then
	        j1=j1+1
	        hlp1=(H(1,1)+H(at_type(i),1))/2.d0
			Ham(i1-3,j1)=S(i1-3,j1)*hlp1
	        hlp1=(H(1,1)+H(at_type(i),2))/2.d0
			Ham(i1-2,j1)=S(i1-2,j1)*hlp1
			Ham(i1-1,j1)=S(i1-1,j1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
			if (i.ne.j) then	        

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(at_type(i),1)
	          hlp1=hlp1*(R-R0(at_type(i),1))
	          hlp1=dexp(-hlp1)
	          hlp1=K0(at_type(i),1)*hlp1
	          Ham(i1-3,j1)=Ham(i1-3,j1)*hlp1
   			  Ham(i1-2,j1)=Ham(i1-2,j1)*hlp1
			  Ham(i1-1,j1)=Ham(i1-1,j1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
	        endif
	      endif
!	Не водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4
			hlp1=(H(at_type(i),1)+H(at_type(j),1))/2.d0
			Ham(i1-3,j1-3)=S(i1-3,j1-3)*hlp1
  			hlp1=(H(at_type(i),1)+H(at_type(j),2))/2.d0
			Ham(i1-3,j1-2)=S(i1-3,j1-2)*hlp1
			Ham(i1-3,j1-1)=S(i1-3,j1-1)*hlp1
			Ham(i1-3,j1)=S(i1-3,j1)*hlp1
			hlp1=(H(at_type(i),2)+H(at_type(j),1))/2.d0
			Ham(i1-2,j1-3)=S(i1-2,j1-3)*hlp1
			Ham(i1-1,j1-3)=S(i1-1,j1-3)*hlp1
			Ham(i1,j1-3)=S(i1,j1-3)*hlp1
			hlp1=(H(at_type(i),2)+H(at_type(j),2))/2.d0
			Ham(i1-2,j1-2)=S(i1-2,j1-2)*hlp1
			Ham(i1-2,j1-1)=S(i1-2,j1-1)*hlp1
			Ham(i1-2,j1)=S(i1-2,j1)*hlp1
			Ham(i1-1,j1-2)=S(i1-1,j1-2)*hlp1
			Ham(i1-1,j1-1)=S(i1-1,j1-1)*hlp1
			Ham(i1-1,j1)=S(i1-1,j1)*hlp1
			Ham(i1,j1-2)=S(i1,j1-2)*hlp1
			Ham(i1,j1-1)=S(i1,j1-1)*hlp1
			Ham(i1,j1)=S(i1,j1)*hlp1
	     	if (i.ne.j) then

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 

!	          R=(x(i)-x(j))**2+(y(i)-y(j))**2
!	          R=R+(z(i)-z(j))**2; R=dsqrt(R)
	          hlp1=delta(at_type(i),at_type(j))
			  hlp1=hlp1*(R-R0(at_type(i),at_type(j)))
	          hlp1=dexp(-hlp1)
	          hlp1=hlp1*K0(at_type(i),at_type(j))
			  Ham(i1-3,j1-3)=Ham(i1-3,j1-3)*hlp1
			  Ham(i1-3,j1-2)=Ham(i1-3,j1-2)*hlp1
			  Ham(i1-3,j1-1)=Ham(i1-3,j1-1)*hlp1
			  Ham(i1-3,j1)=Ham(i1-3,j1)*hlp1
			  Ham(i1-2,j1-3)=Ham(i1-2,j1-3)*hlp1
			  Ham(i1-2,j1-2)=Ham(i1-2,j1-2)*hlp1
		 	  Ham(i1-2,j1-1)=Ham(i1-2,j1-1)*hlp1
			  Ham(i1-2,j1)=Ham(i1-2,j1)*hlp1
			  Ham(i1-1,j1-3)=Ham(i1-1,j1-3)*hlp1
			  Ham(i1-1,j1-2)=Ham(i1-1,j1-2)*hlp1
			  Ham(i1-1,j1-1)=Ham(i1-1,j1-1)*hlp1
			  Ham(i1-1,j1)=Ham(i1-1,j1)*hlp1
			  Ham(i1,j1-3)=Ham(i1,j1-3)*hlp1
			  Ham(i1,j1-2)=Ham(i1,j1-2)*hlp1
			  Ham(i1,j1-1)=Ham(i1,j1-1)*hlp1
			  Ham(i1,j1)=Ham(i1,j1)*hlp1
			endif 
	      endif
	      
	      
	      
	      	        else     !  Blocade
	      
        	        
	      if (at_type(j).eq.1) then
	        j1=j1+1

			if (i.ne.j) then	        

	        endif
	      endif
!	Не водород с не водородом
	      if (at_type(j).gt.1) then
	        j1=j1+4

	     	if (i.ne.j) then

			endif 
	      endif
	      	      	        
	      	        
	      	        
	      	        end if   !  Blocade
	      
	      
	      
	      
	   enddo
	  endif
	  j1=0
	  enddo
	  
 	do i1=1,wfunc_count
	do j1=i1,wfunc_count
   Ham(j1,i1)=Ham(i1,j1)   !Symm
	end do
	end do
	  
               ds=(s-sold)/dhtot; dham=(ham-hamold)/dhtot
               nxyz=3
	        

	       FFtot=-sum(dHam*fmatH)+sum(dS*fmatS) !+erep-erepn
         felz(k1)=FFtot
         z(k1)=z(k1)-dhtot
         !********* End Fz determination **********

!******************  Repuls part of force ************

   Fxx=0.d0;  Fyy=0.d0;  Fzz=0.d0
	do i=1,at_count
	 if(i==k1) cycle

        call dxdydz(dx,dy,dz,R,i,k1) 
!	  l=dx/R; m=dy/R; n=dz/R 

! 	  r=(x(i)-x(k1))**2
!    r=r+(y(i)-y(k1))**2
!	  r=r+(z(i)-z(k1))**2
!    r=dsqrt(r)
	 help_var=-betta(at_type(i),at_type(k1))
	 arg=help_var*(r-R0(at_type(i),at_type(k1)))
	 expon=help_var*dexp(arg)
	 help_var=fi0(at_type(i),at_type(k1))*expon/r
	 fxx=fxx-dx*help_var
	 fyy=fyy-dy*help_var
	 fzz=fzz-dz*help_var
!	 fxx=fxx-(x(k1)-x(i))*help_var
!	 fyy=fyy-(y(k1)-y(i))*help_var
!	 fzz=fzz-(z(k1)-z(i))*help_var

	enddo
     frepx(k1)=frepx(k1)+Fxx; frepy(k1)=frepy(k1)+Fyy; frepz(k1)=frepz(k1)+Fzz
!****************** End of Repuls part of force ************

         end do
         fx=felx+frepx;         fy=fely+frepy;         fz=felz+frepz


	end subroutine force
! **********************************************************************
          
       

!***************   Global eigenproblem h*f=diag*S*f
!***  diag - eigenvalues,  f-eigenvectors, h - Hermit. matrix, S -  symm.+positive?? overlap  matrix
        SUBROUTINE enonort(diag,h,s,f)
        use common_var
        implicit real*8(a-h,o-z)
       	integer nm,n
        dimension h(wfunc_count,wfunc_count),h0(wfunc_count,wfunc_count),hh(wfunc_count,wfunc_count),S(wfunc_count,wfunc_count)
        dimension ee1(wfunc_count),diag(wfunc_count),S0(wfunc_count,wfunc_count),F(wfunc_count,wfunc_count)
        dimension bufer(wfunc_count,wfunc_count)
          
          Norb=wfunc_count

         do 62 i=1,Norb
         do 62 j=1,Norb
  62     S0(i,j)=S(i,j)

         do 63 i=1,Norb
         do 63 j=1,Norb
  63     h0(i,j)=h(i,j)

!C       Diagonalization of matrix S(i,j)

        nm=Norb
        n=Norb

        call eisrs1(nm,n,diag,s,ierr,ee1)
        do j=1,n
        if(diag(j)<0.d0)then 
!        print* ,' Error S<0.',j,diag(j)
!        read*
	        diag(j)=-diag(j)
        end if
        
        end do

!C       diag(n) are eigenvalues Sn of matrix S(i,j)
!C       From now on S(i,n) are coefficients in decomposition
!C       of eigenvectors of matrix S(i,j) over atomic functions:
!C       FFn = SUM_i{S(i,n)*Fi} where Fi are atomic orbitals

!C       !!! NOTE: Functions FFn are not normalized ???

!C       'Normalization' of FFFn

         do 60 i=1,Norb
         if(diag(i)<0.d0)then
!         print*,' KARRRAAU U U U LLL diag<0',i,diag(i)
         end if
         do 60 j=1,Norb
  60     S(i,j)=S(i,j)/dsqrt(dabs(diag(j)))


!C Checking ortho-normalization of FFn

         do 64 n=1,Norb
         do 64 m=1,Norb

 64     F(n,m)=0.d0


        do 776 i1=1,Norb
        do 776 i2=1,Norb
 776    bufer(i1,i2)=0.d0
         do  m=1,Norb
         do  i=1,Norb
         do  j=1,Norb
        bufer(m,i)=bufer(m,i)+S(j,m)*S0(i,j)
        end do
        end do
        end do

         do 66 n=1,Norb                      !
         do 66 m=n,Norb                      !
         do 66 i=1,Norb                      !  expansive block
 66     F(n,m)=F(n,m)+S(i,n)*bufer(m,i)  !
        do  n=1,norb
        do  m=n,norb
        F(m,n)=F(n,m)
        end do
        end do

!C       Creation of Matrix HH(n,m)
!C       This is the Hamiltonian matrix over FFn basis set

        do 10 n=1,Norb
        do 10 m=1,Norb
 10     HH(n,m)=0.d0

        do 777 i1=1,Norb
        do 777 i2=1,Norb
 777    bufer(i1,i2)=0.d0
        do  m=1,Norb
        do  i=1,Norb
        do  j=1,Norb
        bufer(i,m)=bufer(i,m)+S(j,m)*H(i,j)
        end do
        end do
        end do
                                !

        do 20 n=1,Norb                         !
        do 20 m=n,Norb                         !  expansive block
        do 20 i=1,Norb                         !

 20     HH(n,m)=HH(n,m)+S(i,n)*bufer(i,m)      !
         
        do  n=1,Norb                         !
        do  m=n,Norb                         !  expansive block
        HH(m,n)=HH(n,m)
        end do
        end do

!C       Diagonalization of matrix HH(n,m)

        nm=Norb
        n=Norb
           
        call eisrs1(nm,n,diag,HH,ierr,ee1)
!C       Calculation of coefficients in decomposition of
!C       eigenvectors FFFn of HH(n,m) over atomic ofrbitals Fi:
!C       FFFn=SUM_m{HH(m,n)*SUM_i{S(i,m)*Fi}}
!C       Let those coefficients form matrix h(i,n),
!C       i.e. h(i,n)=SUM_m{HH(m,n)*S(i,m)}, i.e., the columns
!C       of matrix h(i,n) are coefficients of eigenvectors in terms
!C       of atomic orbitals Fi !!!

        do 40 i=1,Norb
        do 40 n=1,Norb
 40     h(i,n)=0.d0

        do 50 i=1,Norb
        do 50 n=1,Norb
        do 50 m=1,Norb
 50     h(i,n)=h(i,n)+HH(m,n)*S(i,m)


!C       Saving non-normalized self-vectors as F(i,k)
        do 358 i=1,Norb
        do 358 k=1,Norb
 358    f(i,k)=h(i,k)

!C       Normalization of self-vectors
        do 360 k=1,Norb
        SS=0.d0
        do 61 i=1,Norb
 61     SS=SS+h(i,k)*h(i,k)
        do 362 i=1,Norb
 362    h(i,k)=h(i,k)/dsqrt(SS)
 360    continue


!c       CHEKING6 (fulfilment of initial nonorthogonal equation)
        SSS=0.d0
        do 440 n=1,Norb
        SSS1=0.d0
        do 442 i=1,Norb
        SSS2=0.d0
        do 444 k=1,Norb
 444    SSS2=SSS2+h0(i,k)*h(k,n)-s0(i,k)*h(k,n)*diag(n)
 442    SSS1=SSS1+SSS2*SSS2
 440    SSS=SSS+SSS1

        H=H0;  S=S0  


        RETURN
        END

!c ************************************************************
	subroutine read_coord
	use common_var
	integer i,ind
	character(1) s




        alper=1.d6
       	perx(3)=0.d0; pery(3)=0.d0; perz(3)=alper

      	perx(1)=12.52879d0; pery(1)=0.d0; perz(1)=0.d0 !    Graphene C60
        perx(2)=.0d0; pery(2)=13.0203d0; perz(2)=0.d0         	

!      	perx(1)=25.055847d0; pery(1)=0.d0; perz(1)=0.d0 !    Graphene C240
!        perx(2)=.0d0; pery(2)=26.0388d0; perz(2)=0.d0         	

 

	i=1; ind=0
	H_count=0; C_count=0
	N_count=0; O_count=0

	open(1,file='coord.ini')
	do
	read(1,*,end=100) xi,yi,zi,s
!c     1=водород, 2=углерод, 3=азот, 4=кислород
	if (s.eq.'H') then
         H_count=H_count+1; ind=1
	endif
	if (s.eq.'C') then
          C_count=C_count+1; ind=1
	endif
	if (s.eq.'N') then
           N_count=N_count+1; ind=1
	endif
	if (s.eq.'O') then
           O_count=O_count+1; ind=1
	endif
	if (ind.eq.0) then
!	   i=i-1; print*, 'WARNING: Bad chemical element: ',s; read*
	endif
	i=i+1
	enddo
100	close(1)

	at_count=H_count+C_count+N_count+O_count
!        print*,' Nat= ',at_count 
	if ((i-1).ne.at_count) then
!	print*, 'Error in file coord'
	read*; goto 101
	endif


      	if (na1>at_count) then
!	print*, 'Error in drive.ini na1=',na1
        read*
        stop  
	endif


      	if (na2>at_count) then
!	print*, 'Error in drive.ini na2=',na2
        read*
        stop  
	endif


      	if (na3>at_count) then
!	print*, 'Error in drive.ini na3=',na3
        read*
        stop  
	endif


      	if (na4>at_count) then
!	print*, 'Error in drive.ini na4=',na4
        read*
        stop  
	endif


      	if (na5>at_count) then
!	print*, 'Error in drive.ini na5=',na5
        read*
        stop  
	endif


      	if (na6>at_count) then
!	print*, 'Error in drive.ini na6=',na6
        read*
        stop  
	endif


      	if (na7>at_count) then
!	print*, 'Error in drive.ini na7=',na7
        read*
        stop  
	endif

      	if (na8>at_count) then
!	print*, 'Error in drive.ini na8=',na8
        read*
        stop  
	endif


   allocate (x(at_count),y(at_count),z(at_count))
   allocate (x0(at_count),y0(at_count),z0(at_count))
   allocate (fx(at_count),fy(at_count),fz(at_count))
   allocate (fx0(at_count),fy0(at_count),fz0(at_count))
   allocate (felx(at_count),fely(at_count),felz(at_count))
   allocate (frepx(at_count),frepy(at_count),frepz(at_count))
   Nat3=3*at_count
   allocate (amass(Nat3))  
   allocate (Hess(Nat3,Nat3)) 
   allocate (Hesel(Nat3,Nat3),hesrep(Nat3,Nat3)) 
!   allocate (ddd(Nat3),e1(Nat3),delta3(Nat3),ak(Nat3),at(Nat3)) 
   allocate (ddd(Nat3),e1(Nat3),delta3(Nat3)) 
   allocate (at_type(at_count)) 	
  i=1; ind=0
	H_count=0; C_count=0
	N_count=0; O_count=0
		open(1,file='coord.ini')
	do jjj=1,at_count
	
	     read(1,*) x(i),y(i),z(i),s
!c     1=водород, 2=углерод, 3=азот, 4=кислород
	if (s.eq.'H') then
	  at_type(i)=1; H_count=H_count+1; ind=1 
	endif
	if (s.eq.'C') then
	  at_type(i)=2; C_count=C_count+1; ind=1 
	endif
	if (s.eq.'N') then
	  at_type(i)=3; N_count=N_count+1; ind=1 
	endif
	if (s.eq.'O') then
	  at_type(i)=4; O_count=O_count+1; ind=1 
	endif
	if (ind.eq.0) then
	   i=i-1; !print*, 'WARNING: Bad chemical element: ',s; read*
	endif
	i=i+1
	enddo
	close(1)

        do j=1,at_count

        if(at_type(j)==1) then
        amass(1+3*(j-1))=1.d0 
        amass(2+3*(j-1))=1.d0 
        amass(3+3*(j-1))=1.d0 
        end if

        if(at_type(j)==2) then
        amass(1+3*(j-1))= 12.d0
        amass(2+3*(j-1))= 12.d0
        amass(3+3*(j-1))= 12.d0
        end if

        if(at_type(j)==3) then
        amass(1+3*(j-1))=14.d0 
        amass(2+3*(j-1))=14.d0 
        amass(3+3*(j-1))=14.d0 
        end if

        if(at_type(j)==4) then
        amass(1+3*(j-1))=16.d0 
        amass(2+3*(j-1))=16.d0 
        amass(3+3*(j-1))=16.d0 
        end if

        end do
	
	
	el_count=H_count+4*C_count+5*N_count+6*O_count
	wfunc_count=H_count+4*(C_count+N_count+O_count)
!	print*, 'Coordinates have read successfully'
!      print*, 'H: ',H_count,'       C: ',C_count
!	print*, 'N: ',N_count,'       O: ',O_count
101	continue
	end subroutine read_coord
!c ************************************************************
!c ************************************************************
	subroutine read1_coord(masx,masy,masz,mast,masn)
	use common_var
	integer masn,i,ind,mast(masn)
	real*8 masx(masn),masy(masn),masz(masn)
	character(1) s

	i=1; ind=0
	H_count=0; C_count=0
	N_count=0; O_count=0

    do i=1,masn
		if (mast(i).eq.1) H_count=H_count+1
		if (mast(i).eq.2) C_count=C_count+1
		if (mast(i).eq.3) N_count=N_count+1
		if (mast(i).eq.4) O_count=O_count+1 
	enddo

	at_count=H_count+C_count+N_count+O_count
    !print*,' Nat= ',at_count 
   allocate (x(at_count),y(at_count),z(at_count))
   allocate (x0(at_count),y0(at_count),z0(at_count))
   allocate (fx(at_count),fy(at_count),fz(at_count))
   allocate (fx0(at_count),fy0(at_count),fz0(at_count))
   allocate (felx(at_count),fely(at_count),felz(at_count))
   allocate (frepx(at_count),frepy(at_count),frepz(at_count))
   Nat3=3*at_count
   allocate (amass(Nat3))  
   allocate (Hess(Nat3,Nat3)) 
   allocate (Hesel(Nat3,Nat3),hesrep(Nat3,Nat3)) 
!   allocate (ddd(Nat3),e1(Nat3),delta3(Nat3),ak(Nat3),at(Nat3)) 
   allocate (ddd(Nat3),e1(Nat3),delta3(Nat3)) 
   allocate (at_type(at_count)) 	
 
 x(1:at_count)=masx(1:at_count); y(1:at_count)=masy(1:at_count);
 z(1:at_count)=masz(1:at_count); at_type(1:at_count)=mast(1:at_count)

        do j=1,at_count

        if(at_type(j)==1) then
        amass(1+3*(j-1))=1.d0 
        amass(2+3*(j-1))=1.d0 
        amass(3+3*(j-1))=1.d0 
        end if

        if(at_type(j)==2) then
        amass(1+3*(j-1))= 12.d0
        amass(2+3*(j-1))= 12.d0
        amass(3+3*(j-1))= 12.d0
        end if

        if(at_type(j)==3) then
        amass(1+3*(j-1))=14.d0 
        amass(2+3*(j-1))=14.d0 
        amass(3+3*(j-1))=14.d0 
        end if

        if(at_type(j)==4) then
        amass(1+3*(j-1))=16.d0 
        amass(2+3*(j-1))=16.d0 
        amass(3+3*(j-1))=16.d0 
        end if

        end do
	
	
	el_count=H_count+4*C_count+5*N_count+6*O_count
	wfunc_count=H_count+4*(C_count+N_count+O_count)
!	print*, 'Coordinates have read successfully'
!      print*, 'H: ',H_count,'       C: ',C_count
!	print*, 'N: ',N_count,'       O: ',O_count
101	continue
	end subroutine read1_coord
!c ************************************************************
	subroutine read_param
	use param
	use drive
	use common_var

fi0(1,1)= 0.78d0                
fi0(1,2)= 0.561102d0          
fi0(1,3)= 0.543169218862232d0 
fi0(1,4)= 0.977083258568997d0 
fi0(2,2)= 0.943505d0            
fi0(2,3)= 0.98833969237026d0  
fi0(2,4)= 3.75939979834437d0  
fi0(3,3)= 1.0300069814527388d0
fi0(3,4)= 1.58235628534691d0    
fi0(4,4)= 2.534701203510635d0  

	fi0(2,1)=fi0(1,2); 
  fi0(3,1)=fi0(1,3)
	fi0(4,1)=fi0(1,4); 
  fi0(3,2)=fi0(2,3)
	fi0(4,2)=fi0(2,4); 
  fi0(4,3)=fi0(3,4)

betta(1,1)=6.84d0     
betta(1,2)=9.433587d0
betta(1,3)=9.74435453678239d0   
betta(1,4)=7.29873931487849d0  
betta(2,2)=4.912617d0 
betta(2,3)=4.9160194858728d0
betta(2,4)=4.75915316283834d0   
betta(3,3)=2.842712743792253d0 
betta(3,4)=11.6082068698282d0   
betta(4,4)=2.625773945160157d0 

	betta(2,1)=betta(1,2); betta(3,1)=betta(1,3)
	betta(4,1)=betta(1,4); betta(3,2)=betta(2,3)
	betta(4,2)=betta(2,4); betta(4,3)=betta(3,4)	
	!read(1,*)
!c	параметры отталкивающего потенциала R0_ij
 R0(1,1)=  0.75d0            
 R0(1,2)=  1.04512d0      
 R0(1,3)=  0.973807830628576d0   
 R0(1,4)=  1.05109423147364d0
 R0(2,2)=  1.582565d0        
 R0(2,3)=  1.53590349862476d0
 R0(2,4)=  1.22732578158281d0   
 R0(3,3)=  2.2357416600411d0
 R0(3,4)=  1.15862917929351d0    
 R0(4,4)=  1.7912484062831358d0  
  
	R0(2,1)=R0(1,2); R0(3,1)=R0(1,3)
	R0(4,1)=R0(1,4); R0(3,2)=R0(2,3)
	R0(4,2)=R0(2,4); R0(4,3)=R0(3,4)
	!read(1,*); read(1,*)
!c	ПАРАМЕТРЫ КВАНТОВОГО ЭЛЕКТРОННОГО ПОТЕНЦИАЛА
!c	параметры электронного потенциала K0_ij
K0(1,1)=  1.68d0         
K0(1,2)=  1.763801d0     
K0(1,3)=  1.79670074631674d0
K0(1,4)=  1.98485271191831d0
K0(2,2)=  2.06029d0      
K0(2,3)=  2.11783599747328d0
K0(2,4)=  2.31614488021489d0
K0(3,3)=  2.691627179735403d0
K0(3,4)=  2.11413781529936d0
K0(4,4)=  2.6873751888495883d0  

	K0(2,1)=K0(1,2); K0(3,1)=K0(1,3)
	K0(4,1)=K0(1,4); K0(3,2)=K0(2,3)
	K0(4,2)=K0(2,4); K0(4,3)=K0(3,4)
	!read(1,*)
!c	параметры электронного потенциала delta_ij
delta(1,1)= 0.13d0        
delta(1,2)= 0.01435d0     
delta(1,3)= 0.01457306082829996d0   
delta(1,4)= 0.487728960368312d0
delta(2,2)= 0.164262d0    
delta(2,3)= 0.162132198056033d0
delta(2,4)= 0.153835842271449d0 
delta(3,3)= 0.17721992892424443d0  
delta(3,4)= 0.2107205790853d0 
delta(4,4)= 0.2817058969905267d0

	delta(2,1)=delta(1,2); delta(3,1)=delta(1,3)
	delta(4,1)=delta(1,4); delta(3,2)=delta(2,3)
	delta(4,2)=delta(2,4); delta(4,3)=delta(3,4)
	!read(1,*)
!c	параметры гамильтониана H_ialfa
	H(1,1)=-10.7d0  
  H(2,1)=-16.157972d0
	H(2,2)=-10.078261d0  
  H(3,1)=-23.753608060128023d0
	H(3,2)=-10.961350662538958d0 
  H(4,1)=-22.092416999000033d0
	H(4,2)=-11.255111473843405d0 
  H(1,2)=0.d0; 
!c	параметры волновых функций ksi_ialfa
 ksi(1,1)=  2.456644d0  
 ksi(2,1)=  2.991164d0  
 ksi(2,2)=  3.857861d0  
 ksi(3,1)=  4.360867316226992d0  
 ksi(3,2)=  4.31405344045801d0
 ksi(4,1)=  5.980314169054799d0  
 ksi(4,2)=  4.461321601735512d0  
  ksi(1,2)=0.d0;
!c	print*, 'Parameters have read successfully'
	
 step =0.005d0
	
	end subroutine read_param
   !c     ПРОЦЕДУРЫ ДЛЯ ВЫЧИСЛЕНИЯ АНАЛИТИЧЕСКИХ ИНТЕГРАЛОВ ПЕРЕКРЫТИЯ 1S1S, 1S2S,1S2P,2S2S, 2S2P, 2P2Pq, 2p2Ppi
!c	R-расстояние, i,j - типы атомов (1=H; 2=C; 3=N; 4=O; In - выходной параметр - получившийся интеграл)
!c **********************************************************************
	subroutine In1S1S(R,In)
	use param
	real*8 R,p,ks,In
	   ks=ksi(1,1); p=ks*R
	   In=(1.d0+p+p*p/3.d0)*dexp(-p)
	end subroutine In1S1S
!c **********************************************************************
	subroutine In1S2S(R,i,In)
	use param
	real*8 R,p,ks,k,t,p1,p2,hlp1,hlp2,hlp3,In
	integer i
	    ks=(ksi(1,1)+ksi(i,1))/2.d0
	    t=(ksi(1,1)-ksi(i,1))/(2.d0*ks)
	    k=(t+1.d0/t)/2.d0
	    p1=R*ksi(1,1); p2=R*ksi(i,1); p=ks*R 
	    hlp1=dsqrt(1-t*t)/(dsqrt(3.d0)*t*p)
	    hlp2=2.d0*(1.d0+k)*(2.d0-3.d0*k)
	    hlp2=hlp2+(1.d0-2.d0*k)*p1
	    hlp2=-(1.d0-k)*hlp2*dexp(-p1)
	    hlp3=2.d0*(1.d0-k)*(2.d0-3.d0*k)
	    hlp3=hlp3+4.d0*(1.d0-k)*p2+p2*p2
	    hlp3=(1.d0+k)*hlp3*dexp(-p2)
	    In=hlp1*(hlp2+hlp3)
	end subroutine In1S2S
!c **********************************************************************
	subroutine In1S2P(R,i,In)
	use param
	real*8 R,p,ks,k,t,p1,p2,hlp1,hlp2,hlp3,In
	integer i
	  ks=(ksi(1,1)+ksi(i,2))/2.d0
	  t=(ksi(1,1)-ksi(i,2))/(2.d0*ks)
	  k=(t+1.d0/t)/2.d0
	  p1=R*ksi(1,1); p2=R*ksi(i,2); p=ks*R 
	  hlp1=dsqrt((1.d0+t)/(1.d0-t))/(t*p*p)
	  hlp2=6.d0*(1.d0+k)*(1.d0+p1)+2.d0*p1*p1
	  hlp2=-(1.d0-k)*(1.d0-k)*hlp2*dexp(-p1)
	  hlp3=6.d0*(1.d0-k)*(1.d0-k)*(1+p2)
	  hlp3=hlp3+4.d0*(1.d0-k)*p2*p2+p2*p2*p2
	  hlp3=(1.d0+k)*hlp3*dexp(-p2)
        In=-hlp1*(hlp2+hlp3)  ! коррекция знака
	end subroutine In1S2P
!c ***********************************************************************
	subroutine In2S2S(R,i,j,In)
	use param
	real*8 R,p,ks,k,t,p1,p2,hlp1,hlp2,hlp3,In
	integer i,j
!c	    if (i.eq.j) then
		if (dabs(ksi(i,1)-ksi(j,1)).lt.0.1d0) then
	      ks=(ksi(i,1)+ksi(j,1))/2.d0; p=ks*R
	      hlp1=1.d0+p+4.d0*p*p/9.d0+p*p*p/9.d0
	      hlp1=(hlp1+p*p*p*p/45.d0)*dexp(-p)
	      In=hlp1
	    endif
	    if (i.ne.j) then
            ks=(ksi(i,1)+ksi(j,1))/2.d0
		  t=(ksi(i,1)-ksi(j,1))/(2.d0*ks)
		  k=(t+1.d0/t)/2.d0
		  p1=R*ksi(i,1); p2=R*ksi(j,1); p=ks*R 
	      hlp1=dsqrt(1-t*t)/(3.d0*t*p)
	      hlp2=2.d0*(1.d0+k)*(7.d0-12.d0*k*k)
	      hlp2=hlp2+4.d0*(1.d0+k)*(2.d0-3.d0*k)*p1
	      hlp2=hlp2+(1.d0-2.d0*k)*p1*p1;
	      hlp2=-(1.d0-k)*hlp2*dexp(-p1)
	      hlp3=2.d0*(1.d0-k)*(7.d0-12.d0*k*k)
	      hlp3=hlp3+4.d0*(1.d0-k)*(2.d0+3.d0*k)*p2
	      hlp3=hlp3+(1.d0+2.d0*k)*p2*p2;
	      hlp3=(1.d0+k)*hlp3*dexp(-p2)
	      In=hlp1*(hlp2+hlp3)
	    endif
	end subroutine In2S2S
!c ***********************************************************************
	subroutine In2S2P(R,i,j,In)
	use param
	real*8 R,p,ks,k,t,p1,p2,hlp1,hlp2,hlp3,In
	integer i,j
	    if (i.eq.j) then
	      ks=(ksi(i,1)+ksi(j,2))/2.d0
		  t=(ksi(i,1)-ksi(j,2))/(2.d0*ks)
		  k=(t+1.d0/t)/2.d0
		  p1=R*ksi(i,1); p2=R*ksi(j,2); p=ks*R 
	      hlp1=dsqrt((1.d0+t)/(1.d0-t))
		  hlp1=hlp1/(dsqrt(3.d0)*t*p*p)
	      hlp2=6.d0*(1.d0+k)*(3.d0+4.d0*k)*(1.d0+p1)
	      hlp2=hlp2+2.d0*(5.d0+6.d0*k)*p1*p1
	      hlp2=hlp2+2.d0*p1*p1*p1
	      hlp2=-(1.d0-k)*(1.d0-k)*hlp2*dexp(-p1)
	      hlp3=6.d0*(1.d0-k)*(1.d0-k)*(3.d0+4.d0*k)
	      hlp3=hlp3*(1.d0+p2)
	      hlp3=hlp3+4.d0*(1.d0-k)*(2.d0+3.d0*k)*p2*p2
	      hlp3=hlp3+(1.d0+2.d0*k)*p2*p2*p2
	      hlp3=(1.d0+k)*hlp3*dexp(-p2)
	      In=-hlp1*(hlp2+hlp3) ! коррекция знака
	    endif
	    if (i.ne.j) then
	      ks=(ksi(i,1)+ksi(j,2))/2.d0
		  t=(ksi(i,1)-ksi(j,2))/(2.d0*ks)
		  k=(t+1.d0/t)/2.d0
		  p1=R*ksi(i,1); p2=R*ksi(j,2); p=ks*R 
	      hlp1=dsqrt((1.d0+t)/(1.d0-t))
		  hlp1=hlp1/(dsqrt(3.d0)*t*p*p)
	      hlp2=6.d0*(1.d0+k)*(3.d0+4.d0*k)*(1.d0+p1)
	      hlp2=hlp2+2.d0*(5.d0+6.d0*k)*p1*p1
	      hlp2=hlp2+2.d0*p1*p1*p1
	      hlp2=-(1.d0-k)*(1.d0-k)*hlp2*dexp(-p1)
	      hlp3=6.d0*(1.d0-k)*(1.d0-k)*(3.d0+4.d0*k)
	      hlp3=hlp3*(1.d0+p2)
	      hlp3=hlp3+4.d0*(1.d0-k)*(2.d0+3.d0*k)*p2*p2
	      hlp3=hlp3+(1.d0+2.d0*k)*p2*p2*p2
	      hlp3=(1.d0+k)*hlp3*dexp(-p2)
	      In=-hlp1*(hlp2+hlp3)  ! коррекция знака
	    endif
	    IF (dabs(ksi(i,1)-ksi(j,2)).lt.0.2d0) THEN
		  ks=(ksi(i,1)+ksi(j,2))/2.d0; p=ks*R 
	      hlp1=p/(2.d0*dsqrt(3.d0))
	      In=hlp1*(1.d0+p+(7.d0*p*p+2.d0*p*p*p)/15.d0)*dexp(-p)
	      In=-In  ! коррекция знака
	    ENDIF
	end subroutine In2S2P
!c ***********************************************************************
	subroutine In2P2Pq(R,i,j,In)
	use param
	real*8 R,p,ks,k,t,p1,p2,hlp1,hlp2,hlp3,In
	integer i,j
!c	    if (i.eq.j) then
		if (dabs(ksi(i,2)-ksi(j,2)).lt.0.1d0) then
            ks=(ksi(i,2)+ksi(j,2))/2.d0; p=ks*R
      	  hlp1=(-1.d0-p-p*p/5.d0)
	      hlp1=hlp1+(2.d0*p*p*p+p*p*p*p)/15.d0
 	      In=-hlp1*dexp(-p)   ! коррекция знака
	    endif
	    if (i.ne.j) then
	      ks=(ksi(i,2)+ksi(j,2))/2.d0
		  t=(ksi(i,2)-ksi(j,2))/(2.d0*ks)
		  k=(t+1.d0/t)/2.d0
		  p1=R*ksi(i,2); p2=R*ksi(j,2); p=ks*R 
	      hlp1=1.d0/(dsqrt(1.d0-t*t)*t*p*p*p)
	      hlp2=48.d0*(1.d0+k)*(1.d0+k)
	      hlp2=hlp2*(1.d0+p1+(p1*p1/2.d0))
	      hlp2=hlp2+2.d0*(5.d0+6.d0*k)*p1*p1*p1
	      hlp2=hlp2+2.d0*p1*p1*p1*p1
	      hlp2=-(1.d0-k)*(1.d0-k)*hlp2*dexp(-p1)
	      hlp3=48.d0*(1.d0-k)*(1.d0-k)
	      hlp3=hlp3*(1.d0+p2+(p2*p2/2.d0))
	      hlp3=hlp3+2.d0*(5.d0-6.d0*k)*p2*p2*p2
	      hlp3=hlp3+2.d0*p2*p2*p2*p2
	      hlp3=(1.d0+k)*(1.d0+k)*hlp3*dexp(-p2)
	      In=-hlp1*(hlp2+hlp3) ! коррекция знака    
	    endif
	end subroutine In2P2Pq
!c ***********************************************************************
	subroutine In2P2Ppi(R,i,j,In)
	use param
	real*8 R,p,ks,k,t,p1,p2,hlp1,hlp2,hlp3,In
	integer i,j
!c	    if (i.eq.j) then
		if (dabs(ksi(i,2)-ksi(j,2)).lt.0.1d0) then
            ks=(ksi(i,2)+ksi(j,2))/2.d0; p=ks*R
      	  hlp2=1.d0+p+(2.d0*p*p/5.d0)+p*p*p/15.d0
	      In=hlp2*dexp(-p)
	    endif
	    if (i.ne.j) then
	      ks=(ksi(i,2)+ksi(j,2))/2.d0
		  t=(ksi(i,2)-ksi(j,2))/(2.d0*ks)
		  k=(t+1.d0/t)/2.d0
		  p1=R*ksi(i,2); p2=R*ksi(j,2); p=ks*R 
	      hlp1=1.d0/(dsqrt(1.d0-t*t)*t*p*p*p)
	      hlp2=24.d0*(1.d0+k)*(1.d0+k)*(1.d0+p1)
	      hlp2=hlp2+12.d0*(1.d0+k)*p1*p1
	      hlp2=hlp2+2*p1*p1*p1
	      hlp2=-(1.d0-k)*(1.d0-k)*hlp2*dexp(-p1)
	      hlp3=24.d0*(1.d0-k)*(1.d0-k)*(1.d0+p2)
	      hlp3=hlp3+12.d0*(1.d0-k)*p2*p2
	      hlp3=hlp3+2*p2*p2*p2
	      hlp3=(1.d0+k)*(1.d0+k)*hlp3*dexp(-p2)
	      In=hlp1*(hlp2+hlp3)
	    endif
	end subroutine In2P2Ppi
!c ***********************************************************************
!	subroutine dist(i,j,r)
!	use common_var
!	integer i,j
!	real*8 r
!	  if (i.gt.at_count) then
!		print*, 'ERROR in subroutine dist'; read*
!	  endif
!	  if (j.gt.at_count) then
!		print*, 'ERROR in subroutine dist'; read*
!	  endif
!	  r=(x(i)-x(j))*(x(i)-x(j))+(y(i)-y(j))*(y(i)-y(j))
!	  r=r+(z(i)-z(j))*(z(i)-z(j)); r=dsqrt(r)
!	end subroutine dist
!**********************************************************************
        subroutine wri_crd
  	use  common_var

        open(2,file='coord.rel')
        do j=1,at_count
        if(at_type(j)==1) write(2,1) x(j),y(j),z(j)
        if(at_type(j)==2) write(2,2) x(j),y(j),z(j)
        if(at_type(j)==3) write(2,3) x(j),y(j),z(j)
        if(at_type(j)==4) write(2,4) x(j),y(j),z(j)
        end do        
        close(2)
1       format(1x, 3e20.12,3x,'H')
2       format(1x, 3e20.12,3x,'C')
3       format(1x, 3e20.12,3x,'N')
4       format(1x, 3e20.12,3x,'O')


        open(2,file='coord.xyz')
        do j=1,at_count
        if(at_type(j)==1) write(2,5) x(j),y(j),z(j)
        if(at_type(j)==2) write(2,6) x(j),y(j),z(j)
        if(at_type(j)==3) write(2,7) x(j),y(j),z(j)
        if(at_type(j)==4) write(2,8) x(j),y(j),z(j)
        end do        
        close(2)
5       format(1x,'1', 3e20.12,3x)
6       format(1x,'6', 3e20.12,3x)
7       format(1x,'7', 3e20.12,3x)
8       format(1x,'8', 3e20.12,3x)


        end  
!*********************************************************
!******************************************************************
        subroutine EISRS1(NM,N,D,ZR,IERR,E)
        !use a111
        use common_var
        implicit real*8(a-h,o-z)
        INTEGER NM,N,IERR,i
        dimension E(nm),D(nm),ZR(nm,nm)
!C       all eigenvalues and corresponding eigenvectors
!C	of reai symmetric matrix

       	call TRED2(NM,N,D,E,ZR)

!        print *,'matrix obtained in threedigonal type'
        call TQL2(NM,N,D,E,ZR,IERR)
        return
	      end
!*****************************************
      SUBROUTINE TRED2(NM,N,D,E,ZR)
      !use a111
      use common_var
        implicit real*8(a-h,o-z)
      INTEGER I,J,K,L,N,II,NM,JP1
      dimension D(nm),E(nm),ZR(nm,nm)
      
          
      IF (N .EQ. 1) GO TO 320
      DO  II = 2, N
      
         I = N + 2 - II
         L = I - 1
         H = 0.d0
         SCALE = 0.d0
         IF (L .LT. 2) GO TO 130

         DO  K = 1, L
 
         SCALE = SCALE + DABS(ZR(I,K))
         end do



         IF (SCALE .NE. 0.d0) GO TO 140
  130    E(I) = ZR(I,L)
         GO TO 290



  140    continue
  
            
            DO  K = 1, L
            ZR(I,K) = ZR(I,K) / SCALE
            H = H + ZR(I,K) * ZR(I,K)
           end do

         F = ZR(I,L)
         G = -DSIGN(DSQRT(H),F)
         E(I) = SCALE * G
         H = H - F * G
         ZR(I,L) = F - G
         F = 0.d0
         DO  J = 1, L
            ZR(J,I) = ZR(I,J) / (SCALE * H)
            G = 0.d0
            DO  K = 1, J
            G = G + ZR(J,K) * ZR(I,K)
            end do

            JP1 = J + 1
            IF (L .LT. JP1) then
            else
                  DO  K = JP1, L
                  G = G + ZR(K,J) * ZR(I,K)
                  end do
            end if
            E(J) = G / H
            F = F + E(J) * ZR(I,J)
         end do

         HH = F / (H + H)
         DO  J = 1, L
            F = ZR(I,J)
            G = E(J) - HH * F
            E(J) = G
            DO  K = 1, J
               ZR(J,K) = ZR(J,K) - F * E(K) - G * ZR(I,K)
         end do
         end do
         DO  K = 1, L
         ZR(I,K) = SCALE * ZR(I,K)
         end do
  290    D(I) = H


        end do
  

  320 D(1) = 0.d0
      E(1) = 0.d0
      DO  I = 1, N
         L = I - 1
         IF (D(I) .EQ. 0.d0) GO TO 380

         DO  J = 1, L
            G = 0.d0
            DO  K = 1, L
            G = G + ZR(I,K) * ZR(K,J)
            end do
            DO  K = 1, L
               ZR(K,J) = ZR(K,J) - G * ZR(K,I)
          end do
          end do

  380    D(I) = ZR(I,I)
         ZR(I,I) = 1.d0
         IF (L .LT. 1) cycle

         DO  J = 1, L
            ZR(I,J) = 0.d0
            ZR(J,I) = 0.d0
            end do
              end do
  
      RETURN
      END
!*****************************************************
      SUBROUTINE TQL2(NM,N,D,E,ZR,IERR)
      !use a111
            use common_var
        implicit real*8(a-h,o-z)
        real*8::MACHEP
      INTEGER I,J,K,L,M,N,II,NM,MML,IERR
      dimension D(nm),E(nm),ZR(nm,nm)
      MACHEP=2.d0**(-52)
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
      DO  I = 2, N
      E(I-1) = E(I)
      end do
      F = 0.d0
      B = 0.d0
      E(N) = 0.d0

      DO 240 L = 1, N

         J = 0
         H = MACHEP * (DABS(D(L)) + DABS(E(L)))
         IF (B .LT. H) B = H
         DO  M = L, N
            IF (DABS(E(M)) .LE. B) exit
          end do
         IF (M .EQ. L) GO TO 220

  130    IF (J .EQ. 50) GO TO 1000
         J = J + 1
         P = (D(L+1) - D(L)) / (2.d0 * E(L))
         R = DSQRT(P*P+1.d0)
         H = D(L) - E(L) / (P + DSIGN(R,P))
         DO  I = L, N
         D(I) = D(I) - H
         end do
         F = F + H
         P = D(M)
         C = 1.d0
         S = 0.d0
         MML = M - L
         DO 200 II = 1, MML
            I = M - II
            G = C * E(I)
            H = C * P
            IF (DABS(P) .LT. DABS(E(I))) GO TO 150
            C = E(I) / P
            R = DSQRT(C*C+1.d0)
            E(I+1) = S * P * R
            S = C / R
            C = 1.d0 / R
            GO TO 160
  150       C = P / E(I)
            R = DSQRT(C*C+1.d0)
            E(I+1) = S * E(I) * R
            S = 1.d0 / R
            C = C * S
  160       P = C * D(I) - S * G
            D(I+1) = H + S * (C * G + S * D(I))
            DO  K = 1, N
               H = ZR(K,I+1)
               ZR(K,I+1) = S * ZR(K,I) + C * H
               ZR(K,I) = C * ZR(K,I) - S * H
            end do
  200    CONTINUE
         E(L) = S * P
         D(L) = C * P

         IF (DABS(E(L)) .GT. B) GO TO 130
         
  220    D(L) = D(L) + F
  240 CONTINUE
  
      DO  II = 2, N
         I = II - 1
         K = I
         P = D(I)
         DO  J = II, N
            IF (D(J) .GE. P) cycle
            K = J
            P = D(J)
           end do
         IF (K .EQ. I) cycle
         D(K) = D(I)
         D(I) = P
         DO  J = 1, N
            P = ZR(J,I)
            ZR(J,I) = ZR(J,K)
            ZR(J,K) = P
         end do
          end do


      GO TO 1001
 1000 IERR = L
 1001   continue
 
	RETURN
      END
!*************************************************************
        SUBROUTINE dxdydz(dx,dy,dz,romin,i11,i22)
        use common_var
        real*8:: romin,dx0,dy0,dz0,ro,x1,x2,y1,y2,z1,z2 
        real*8::dx,dy,dz,dxmin,dymin,dzmin  
        integer:: i11,i22,i1,i2,i3 
        
        x1=x(i11)
        y1=y(i11)
        z1=z(i11)

        x2=x(i22)
        y2=y(i22)
        z2=z(i22)

        dx=x2-x1
        dy=y2-y1
        dz=z2-z1
        romin=dsqrt(dx**2+dy**2+dz**2)
        return  ! without periodic cond.
        end
!*************************************************************


           SUBROUTINE Hess_cre
	use param
	use common_var

	integer i,ind,j,k1
        real*8:: dd,eee
        real*8:: dx,dy,dz,r,r2,arg,expon,bet,l,m,n 
        real*8::  sxx,syy,szz,sxy 

        x0=x;  y0=y;  z0=z;

        dd=1.e-4
         
!*************  Electron part ***************************
        do i=1,at_count
        x(i)=x(i)-dd
        call force(eee)

        do j=1,at_count
        fx0(j)=felx(j)
        fy0(j)=fely(j)
        fz0(j)=felz(j)
        end do

        x(i)=x(i)+2.*dd 
        call force(eee)

        do j=1,at_count
        hesel(1+3*(i-1),1+3*(j-1))=-0.5*(felx(j)-fx0(j))/dd
        hesel(1+3*(i-1),2+3*(j-1))=-0.5*(fely(j)-fy0(j))/dd
        hesel(1+3*(i-1),3+3*(j-1))=-0.5*(felz(j)-fz0(j))/dd
        end do
        x(i)=x(i)-dd 


        y(i)=y(i)-dd 
        call force(eee)

        do j=1,at_count
        fx0(j)=felx(j)
        fy0(j)=fely(j)
        fz0(j)=felz(j)
        end do

        y(i)=y(i)+2.*dd 
        call force(eee)

        do j=1,at_count
        hesel(2+3*(i-1),1+3*(j-1))=-0.5*(felx(j)-fx0(j))/dd
        hesel(2+3*(i-1),2+3*(j-1))=-0.5*(fely(j)-fy0(j))/dd
        hesel(2+3*(i-1),3+3*(j-1))=-0.5*(felz(j)-fz0(j))/dd
        end do
        y(i)=y(i)-dd 


        z(i)=z(i)-dd
        call force(eee)

        do j=1,at_count
        fx0(j)=felx(j)
        fy0(j)=fely(j)
        fz0(j)=felz(j)
        end do

        z(i)=z(i)+2.*dd
        call force(eee)

        do j=1,at_count
        hesel(3+3*(i-1),1+3*(j-1))=-0.5*(felx(j)-fx0(j))/dd
        hesel(3+3*(i-1),2+3*(j-1))=-0.5*(fely(j)-fy0(j))/dd
        hesel(3+3*(i-1),3+3*(j-1))=-0.5*(felz(j)-fz0(j))/dd
        end do
        z(i)=z(i)-dd

        end do
!************* End Electron part ***************************


!*************  Repuls part ***************************
        do i=1,at_count


        do j=1,at_count  !  (X,XYZ)

         Sxx=0.d0; Sxy=0.d0; Sxz=0.d0
                   Syy=0.d0; Syz=0.d0
                             Szz=0.d0
             if(i==j) then  !!!!!!!!!!!!!!!
         do k1=1,at_count
	 if(i==k1) cycle
        call dxdydz(dx,dy,dz,R,k1,i) 
!	  l=dx/R; m=dy/R; n=dz/R 
          r2=r**2

!         dx=x0(i)-x0(k1)
!         dy=y0(i)-y0(k1)
!         dz=z0(i)-z0(k1)
! 	  r2=dx**2+dy**2+dz**2
!          r=dsqrt(r2)


	 bet=-betta(at_type(i),at_type(k1))
	 arg=bet*(r-R0(at_type(i),at_type(k1)))
	 expon=bet*dexp(arg)
	 expon=fi0(at_type(i),at_type(k1))*expon/r2
	 Sxx=Sxx+expon*(dx**2*bet+(dy**2+dz**2)/r ) 
	 Sxy=Sxy+expon*(dx*dy*bet-dx*dy/r )
 	 Sxz=Sxz+expon*(dx*dz*bet-dx*dz/r ) 
	 Syy=Syy+expon*(dy**2*bet+(dx**2+dz**2)/r )
	 Syz=Syz+expon*(dy*dz*bet-dy*dz/r )
	 Szz=Szz+expon*(dz**2*bet+(dx**2+dy**2)/r )
	 end do

             else

        call dxdydz(dx,dy,dz,R,j,i) 
!	  l=dx/R; m=dy/R; n=dz/R 
          r2=r**2

!         dx=x0(i)-x0(j)
!         dy=y0(i)-y0(j)
!         dz=z0(i)-z0(j)
! 	  r2=dx**2+dy**2+dz**2
!          r=dsqrt(r2)

	 bet=-betta(at_type(i),at_type(j))
	 arg=bet*(r-R0(at_type(i),at_type(j)))
	 expon=bet*dexp(arg)
	 expon=fi0(at_type(i),at_type(j))*expon/r2
	 Sxx=Sxx-expon*(dx**2*bet+(dy**2+dz**2)/r ) 
	 Sxy=Sxy-expon*(dx*dy*bet-dx*dy/r )
	 Sxz=Sxz-expon*(dx*dz*bet-dx*dz/r )
   Syy=Syy-expon*(dy**2*bet+(dx**2+dz**2)/r )
   Syz=Syz-expon*(dy*dz*bet-dy*dz/r )
   Szz=Szz-expon*(dz**2*bet+(dx**2+dy**2)/r )

             end if   !!!!!!!!!!!!!!!!!!!!
    hesrep(1+3*(i-1),1+3*(j-1))=Sxx   !  XX
    hesrep(1+3*(i-1),2+3*(j-1))=Sxy   !  XY
    hesrep(1+3*(i-1),3+3*(j-1))=Sxz   !  XZ
    hesrep(2+3*(i-1),1+3*(j-1))=Sxy   !  YX
    hesrep(2+3*(i-1),2+3*(j-1))=Syy   !  YY
    hesrep(2+3*(i-1),3+3*(j-1))=Syz   !  YZ
    hesrep(3+3*(i-1),1+3*(j-1))=Sxz   !  ZX
    hesrep(3+3*(i-1),2+3*(j-1))=Syz   !  ZY
    hesrep(3+3*(i-1),3+3*(j-1))=Szz   !  ZZ
    
        end do

        end do  !  End  i - cycle
!************* End  Repuls part ***************************
        hess=hesel+hesrep  


!****   Symmetrization  *********
       do i=1,at_count
       do j=i,at_count
        dd=(hess(i,j)+hess(j,i))/2.d0
        hess(j,i)=dd 
        hess(i,j)=dd 
        end do
        end do
!********************************


        RETURN       

        END
 ~~~~