
module drive
integer:: ntype
real*8 :: step,dsig
real*8 :: visc(4)
end module drive
module param
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
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(:)
real*8,allocatable:: fx0(:),fy0(:),fz0(:)
real*8,allocatable:: felx(:),fely(:),felz(:)
real*8,allocatable:: frepx(:),frepy(:),frepz(:)
real*8,allocatable:: amass(:)
real*8,allocatable:: Hess(:,:)
real*8,allocatable:: Hesel(:,:),hesrep(:,:)
real*8,allocatable:: delta3(:)
real*8,allocatable:: ddd(:),e1(:)
integer,allocatable:: at_type(:)
integer H_count,C_count,N_count,O_count
integer at_count
integer Nat3
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
subroutine main(masx,masy,masz,mast,masn,callback)
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
subroutine ener_rep(e_rep)
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)
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
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)
real*8 fmatH(wfunc_count,wfunc_count)
real*8 fmatS(wfunc_count,wfunc_count)
integer i,j,i1,j1,k1,k2,el_remain,nxyz
i1=0; j1=0
S=0.d0; k1=-1
i1=0; j1=0
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
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
if (i.ne.j) then
call dxdydz(dx,dy,dz,R,i,j)
l=dx/R; m=dy/R; n=dz/R
if (at_type(i).eq.1) then
if (at_type(j).eq.1) then
j1=j1+1
call In1S1S(R,S(i1,j1))
endif
if (at_type(j).gt.1) then
j1=j1+4
call In1S2S(R,at_type(j),S(i1,j1-3))
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
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
enddo
enddo
do i1=1,wfunc_count
do j1=i1,wfunc_count
s(j1,i1)=s(i1,j1)
end do
end do
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
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)
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
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)
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
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
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
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)
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)
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
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
enddo
endif
j1=0
enddo
do i1=1,wfunc_count
do j1=i1,wfunc_count
Ham(j1,i1)=Ham(i1,j1)
end do
end do
sold=s; hamold=ham
call enonort(ener,ham,s,eifun)
call ener_rep(erep)
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))
dhtot=1.d-7
frepx=0.d0; frepy=0.d0; frepz=0.d0;
do k1=1,at_count
S=Sold; Ham=Hamold
x(k1)=x(k1)+dhtot
i1=0; j1=0
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
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
if (i.ne.j) then
call dxdydz(dx,dy,dz,R,i,j)
l=dx/R; m=dy/R; n=dz/R
if (at_type(i).eq.1) then
if (at_type(j).eq.1) then
j1=j1+1
call In1S1S(R,S(i1,j1))
endif
if (at_type(j).gt.1) then
j1=j1+4
call In1S2S(R,at_type(j),S(i1,j1-3))
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
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
enddo
enddo
do i1=1,wfunc_count
do j1=i1,wfunc_count
s(j1,i1)=s(i1,j1)
end do
end do
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
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)
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
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)
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
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
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
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)
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)
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
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
enddo
endif
j1=0
enddo
do i1=1,wfunc_count
do j1=i1,wfunc_count
Ham(j1,i1)=Ham(i1,j1)
end do
end do
ds=(s-sold)/dhtot; dham=(ham-hamold)/dhtot
nxyz=1
FFtot=-sum(dHam*fmatH)+sum(dS*fmatS)
felx(k1)=FFtot
x(k1)=x(k1)-dhtot
y(k1)=y(k1)+dhtot
i1=0; j1=0
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
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
if (i.ne.j) then
call dxdydz(dx,dy,dz,R,i,j)
l=dx/R; m=dy/R; n=dz/R
if (at_type(i).eq.1) then
if (at_type(j).eq.1) then
j1=j1+1
call In1S1S(R,S(i1,j1))
endif
if (at_type(j).gt.1) then
j1=j1+4
call In1S2S(R,at_type(j),S(i1,j1-3))
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
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
enddo
enddo
do i1=1,wfunc_count
do j1=i1,wfunc_count
s(j1,i1)=s(i1,j1)
end do
end do
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
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)
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
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)
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
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
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
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)
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)
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
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
enddo
endif
j1=0
enddo
do i1=1,wfunc_count
do j1=i1,wfunc_count
Ham(j1,i1)=Ham(i1,j1)
end do
end do
ds=(s-sold)/dhtot; dham=(ham-hamold)/dhtot
nxyz=2
FFtot=-sum(dHam*fmatH)+sum(dS*fmatS)
fely(k1)=FFtot
y(k1)=y(k1)-dhtot
z(k1)=z(k1)+dhtot
i1=0; j1=0
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
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
if (i.ne.j) then
call dxdydz(dx,dy,dz,R,i,j)
l=dx/R; m=dy/R; n=dz/R
if (at_type(i).eq.1) then
if (at_type(j).eq.1) then
j1=j1+1
call In1S1S(R,S(i1,j1))
endif
if (at_type(j).gt.1) then
j1=j1+4
call In1S2S(R,at_type(j),S(i1,j1-3))
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
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
enddo
enddo
do i1=1,wfunc_count
do j1=i1,wfunc_count
s(j1,i1)=s(i1,j1)
end do
end do
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
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)
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
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)
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
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
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
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)
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)
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
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
enddo
endif
j1=0
enddo
do i1=1,wfunc_count
do j1=i1,wfunc_count
Ham(j1,i1)=Ham(i1,j1)
end do
end do
ds=(s-sold)/dhtot; dham=(ham-hamold)/dhtot
nxyz=3
FFtot=-sum(dHam*fmatH)+sum(dS*fmatS)
felz(k1)=FFtot
z(k1)=z(k1)-dhtot
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)
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
enddo
frepx(k1)=frepx(k1)+Fxx; frepy(k1)=frepy(k1)+Fyy; frepz(k1)=frepz(k1)+Fzz
end do
fx=felx+frepx; fy=fely+frepy; fz=felz+frepz
end subroutine force
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)
nm=Norb
n=Norb
call eisrs1(nm,n,diag,s,ierr,ee1)
do j=1,n
if(diag(j)<0.d0)then
diag(j)=-diag(j)
end if
end do
do 60 i=1,Norb
if(diag(i)<0.d0)then
end if
do 60 j=1,Norb
60 S(i,j)=S(i,j)/dsqrt(dabs(diag(j)))
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
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
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
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
HH(m,n)=HH(n,m)
end do
end do
nm=Norb
n=Norb
call eisrs1(nm,n,diag,HH,ierr,ee1)
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)
do 358 i=1,Norb
do 358 k=1,Norb
358 f(i,k)=h(i,k)
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
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
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
perx(2)=.0d0; pery(2)=13.0203d0; 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
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
endif
i=i+1
enddo
100 close(1)
at_count=H_count+C_count+N_count+O_count
if ((i-1).ne.at_count) then
read*; goto 101
endif
if (na1>at_count) then
read*
stop
endif
if (na2>at_count) then
read*
stop
endif
if (na3>at_count) then
read*
stop
endif
if (na4>at_count) then
read*
stop
endif
if (na5>at_count) then
read*
stop
endif
if (na6>at_count) then
read*
stop
endif
if (na7>at_count) then
read*
stop
endif
if (na8>at_count) then
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))
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
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;
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)
101 continue
end subroutine read_coord
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
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))
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)
101 continue
end subroutine read1_coord
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)
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)
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)
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)
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;
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;
step =0.005d0
end subroutine read_param
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
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
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
subroutine In2S2S(R,i,j,In)
use param
real*8 R,p,ks,k,t,p1,p2,hlp1,hlp2,hlp3,In
integer i,j
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
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
subroutine In2P2Pq(R,i,j,In)
use param
real*8 R,p,ks,k,t,p1,p2,hlp1,hlp2,hlp3,In
integer i,j
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
subroutine In2P2Ppi(R,i,j,In)
use param
real*8 R,p,ks,k,t,p1,p2,hlp1,hlp2,hlp3,In
integer i,j
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
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 common_var
implicit real*8(a-h,o-z)
INTEGER NM,N,IERR,i
dimension E(nm),D(nm),ZR(nm,nm)
call TRED2(NM,N,D,E,ZR)
call TQL2(NM,N,D,E,ZR,IERR)
return
end
SUBROUTINE TRED2(NM,N,D,E,ZR)
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 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
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
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
do i=1,at_count
do j=1,at_count
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)
r2=r**2
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)
r2=r**2
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
hesrep(1+3*(i-1),2+3*(j-1))=Sxy
hesrep(1+3*(i-1),3+3*(j-1))=Sxz
hesrep(2+3*(i-1),1+3*(j-1))=Sxy
hesrep(2+3*(i-1),2+3*(j-1))=Syy
hesrep(2+3*(i-1),3+3*(j-1))=Syz
hesrep(3+3*(i-1),1+3*(j-1))=Sxz
hesrep(3+3*(i-1),2+3*(j-1))=Syz
hesrep(3+3*(i-1),3+3*(j-1))=Szz
end do
end do
hess=hesel+hesrep
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
~~~~