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