c ....................................................................
c
c Force Field Program UFF
c
c A.K.Rappe,C.J.Casewitt,K.S.Colwell,W.A.Goddard and W.M.Skiff:
c UFF, a Full Periodic Table Force Field for Molecular
c Mechanic and Molecular Dynamic Simulations,
c J.Am.Chem.Soc, 114, 10024 (1992)
c
c Author : K.May and R.Ahlrichs
c Email : Klaus.May@chemie.uni-karlsruhe.de
c Date : 01.06.1999
c Last Change : 04.02.2002
c
c ....................................................................
c
c Calculates the energy, gradient (analytically) and the
c hessian (analytically) of a molecule with anlytical potentials.
c It optimizes the structure with the newton method and a
c linesearch algorithm.
c
c All quantitites are in atomic units.
c ....................................................................
Program UFF
implicit none
external
, fnatoms,rdparm,nol2,iargc
integer
, icoord,igrad,igradx,iforce,ihess,iforcer,ihessr,
, icontrol,itopo,itmp,ifx,ifapprox,nbdim,
, fnatoms,natoms,natomsfx,nhess,nhessfx,
, i,jj,k,l,id,jd,ja,je,jn,idebug,iconv,inotconv,
, nbond,nangle,ntorsion,ninversion,nnonbond,
, nbondv,nanglev,ntorsionv,ninversionv,nnonbondv,
, nbondmx,nanglemx,ntorsionmx,ninversionmx,nnonbondmx,
, imodus,nper,nafix,irc,ierr,ios,ndiis,dbgflg,iauffd,
, iargc,nargs
character*80, allocatable ::
, argstr(:)
character*8, allocatable ::
, mtype(:),ufftype(:)
integer, allocatable ::
, bond(:,:),nonbond(:,:),angle(:,:),torsion(:,:),
, inversion(:,:),
, natt(:),gpnr(:,:),
, nxtnei12(:),atoms12(:,:),
, nxtnei13(:),atoms13(:,:),
, nxtnei14(:),atoms14(:,:),
, ipq(:)
real*8, allocatable ::
, wmass(:),eig(:),eigvec(:,:),
, xyz(:,:),xyzold(:,:),valbond(:,:),bobond(:),
, valangle(:,:),
, valtorsion(:,:),
, valinversion(:,:),
, valnonbond(:,:),
, deges(:,:),degesold(:,:),d2eges(:),d2enum(:),fmat(:),
, q(:),disp(:),emat(:,:),gold(:,:),qold(:,:),bmat(:,:)
character
, non*2,cdbgflg*80
parameter
, (nbdim=6,idebug=0,nper=110,non=' ')
character
, fmd*20,fdummy*20,
, fgrad_def*20,fhess_def*20,ftopo_def*20,
, fcoord*20,fgrad*20,fgradx*20,
, fhess*20,fforce*20,ffapprox*20,
, fhesso*40,fforceo*40,
, fcontrol*20,ftopo*20,fstop*20,
, fconv*20,fnotconv*20,zeile*80,ftmp*20,
, date*8,daytime*10,persys(nper)*2
parameter
, ( fmd='mdmaster',
, ifapprox=4 ,ffapprox='fapprox',
, icoord=7 ,fcoord='coord',
, icontrol=8 ,fcontrol='control',
, itopo=9 ,ftopo_def='ufftopology',
, iconv=10 ,fconv='uffconverged',
, inotconv=11,fnotconv='not.uffconverged',
, fstop='stop',
, igrad=18 ,fgrad_def='uffgradient',
, igradx=19 ,fgradx='uffgradx',
, ihess=20 ,fhess_def='uffhessian0-0',
, iforce=30 ,fforce='uffforceapprox',
, itmp=99 ,ftmp='ufftmp')
parameter
, (nbondv=50,nanglev=100,ntorsionv=300,ninversionv=30,
, nnonbondv=150)
real*8
, eges,egesold,deabs,dipol(3),epsdiag,
, qtot,qtot_def,gnorm,gnormold,dfac,rot(3,3),traeg(3),
, pi,nol2,epssearch,epssteep,
, gconv,econv,gconv_def,econv_def,dfac_def,
, epssteepsmall_def,epssteepbig_def,
, dqmax,dqmax_def,
, alpha,alpha_def,beta,beta_def,gamma,gamma_def,
, epssearchsmall_def,epssearchbig_def,
, xm,ym,zm,atrad(nper),size,dhls,
, dhls_def,ahls,ahls_def
integer
, maxcycle,maxcycle_def,cycle,cycle0,
, modus,modus_def,topmod,inxt12,inxt13,iterm,
, iterm_def,idummy,nqeq_def,nqeq,
, mxls,mxls_def
logical
, lb,la,lt,li,lnv,lne,
, harmonic,crystal,lennard,transform,lnumhess,lmd,
, transform_def,lnumhess_def,lmd_def,
, harmonic_def,crystal_def,lennard_def,
, ex,ex1,lconv,loh,lnext,lsearch,lsteep,excoord,lbfgs,
, interactiv
parameter
, (iterm_def=111111,nqeq_def=0,
, harmonic_def=.true.,crystal_def=.false.,
, maxcycle_def=1,transform_def=.false.,
, lnumhess_def=.false.,lmd_def=.false.,
, gconv_def=1.d-5,econv_def=1.d-08,
, alpha_def=1.d0,beta_def=0.d0,gamma_def=0.d0,
, dfac_def=1.10d0,qtot_def=0.d0,
, epssearchsmall_def=gconv_def,
, epssearchbig_def=1.d+4,
, epssteepsmall_def=1.d+2,
, epssteepbig_def=gconv_def*0.1,
, dqmax_def=0.3,
, lennard_def=.true.,modus_def=1,
, epsdiag=1.d-14,mxls_def=25,dhls_def=0.1d0,
, ahls_def=0.0d0)
c ... symmetry stuff
character sflies*4
integer ndi14,ntrans,ngen,nhessred,natomsred,motz
real*8 thrsym
parameter (ndi14=120)
integer mi(ndi14),mj(ndi14),invt(ndi14)
real*8 trans(9,ndi14),gen(9,3)
integer,allocatable :: ict(:,:),ictnr(:,:),jkseq(:,:),lieff(:),
, nsymeq(:)
real*8,allocatable :: detsym(:,:),symd2eges(:,:,:),xyzeq(:,:,:),
, sumcoo(:,:),chg(:),chgeq(:),
, detgrd(:,:)
character*8, allocatable :: mtypeq(:)
c ... project trans. and rot out of hessian stuff
real*8, allocatable :: trd2eges(:)
c ... Periodic System of elements
data persys
,/'h_',non,non,non,non,non,non,non,non,non,
, non,non,non,non,non,non,non,'he',
, 'li','be',non,non,non,non,non,non,non,non,non,non,
, 'b_','c_','n_','o_','f_','ne',
, 'na','mg',non,non,non,non,non,non,non,non,non,non,
, 'al','si','p_','s_','cl','ar',
, 'k_','ca','sc','ti','v_','cr','mn','fe','co','ni','cu','zn',
, 'ga','ge','as','se','br','kr',
, 'rb','sr','y_','zr','nb','mo','tc','ru','rh','pd','ag','cd',
, 'in','sn','sb','te','i_','xe',
, 'cs','ba','la','hf','ta','w_','re','os','ir','pt','au','hg',
, 'tl','pb','bi','po','at','rn',
, 'fr','ra'/
c ... atrad contains covalent radii (in pm)
c ... interatomic distances (metals)
c
data atrad
,/ 37.0d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,
, 0.d0,0.d0,0.d0,0.d0,0.d0,19.0d0,
, 157.0d0,112.0d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,
, 0.d0,89.0d0,77.0d0,74.0d0,74.0d0,72.0d0,70.0d0,
, 191.0d0,160.0d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,
, 0.d0,143.0d0,117.0d0,110.0d0,104.0d0,99.0d0,95.0d0,
, 235.0d0,197.0d0,164.0d0,147.0d0,135.0d0,129.0d0,137.0d0,
, 126.0d0,125.0d0,125.0d0,128.0d0,137.0d0,
, 153.0d0,139.0d0,121.0d0,117.0d0,114.0d0,111.0d0,
, 250.0d0,215.0d0,182.0d0,160.0d0,147.0d0,140.0d0,135.0d0,
, 134.0d0,134.0d0,137.0d0,144.0d0,152.0d0,
, 167.0d0,158.0d0,161.0d0,137.0d0,133.0d0,129.0d0,
, 272.0d0,224.0d0,187.0d0,159.0d0,147.0d0,141.0d0,137.0d0,
, 135.0d0,136.0d0,139.0d0,144.0d0,155.0d0,
, 171.0d0,175.0d0,182.0d0,164.0d0,145.0d0,140.0d0,
, 282.0d0,235.0d0/
c ......................................................................
c ... Force Field parameter ...
c ... from Rappe and Goddard, JACS, Vol. 114, 1024 (1992) ...
c ......................................................................
integer
, nelementsmx,nelements,rdparm,inparm
character
, finparm*50
parameter
, (nelementsmx=250,inparm=88,finparm='parms.in')
character
, atomtypes(nelementsmx)*5
real*8
, radii(nelementsmx),z(nelementsmx),chi(nelementsmx),
, th0(nelementsmx),x(nelementsmx),d(nelementsmx),
, zeta(nelementsmx),
, J(nelementsmx),Rqeq(nelementsmx),
, Qmin(nelementsmx),Qmax(nelementsmx),
, vsp3(nelementsmx),vsp2(nelementsmx),dconst,epsilon
data dconst/1.11552D+00/,epsilon/1.0D+00/
c ......................................................................
call timex('init','UFF (total)')
c ... check if the program has been called with the argument interactiv
nargs=iargc()
interactiv=.false.
if (nargs.lt.1) then
c ... argstr is part of the next call statement, so it needs a valid
c ... address, even if it will not be used at all
allocate (argstr(1),stat=ierr)
else
allocate (argstr(nargs),stat=ierr)
if (ierr.ne.0)
, STOP ' Error while trying to allocate argstr!'
c ... get the program arguments
call argus(argstr,nargs)
c ... UNIX : analyze procedure arguments
c ... accepted arguments : interactiv
do i=1,nargs
if(index(argstr(i),'interactiv').gt.0) then
write(*,'(/,a,/)') ' UFF will be used interactively'
interactiv=.true.
end if
end do
c free memory for argument string
deallocate (argstr, stat=ierr)
if (ierr.ne.0)
, STOP ' Error in deallocation of argstr!'
end if ! if nargs < 1
call timex('init','Pre.Proc.')
call releas ('uff')
write (*,'(15x,a)') ' '
write (*,'(15x,a)') ' u f f'
write (*,'(15x,a)') ' '
write (*,'(15x,a)') ' by K. May, R. Ahlrichs'
write (*,'(15x,a)') ' '
write (*,'(15x,a)') ' quantum chemistry group '
write (*,'(15x,a)') ' university karlsruhe'
write (*,'(15x,a)') ' germany'
pi=4.d0*atan(1.d0)
c ... read the parameters from file finparm
nelements=rdparm(inparm,finparm,
, atomtypes,radii,th0,x,d,zeta,z,vsp3,vsp2,
, chi,J,rqeq,qmin,qmax,nelementsmx)
c ... number of atoms
natoms=fnatoms(icoord,fcoord)
nhess=3*natoms
ndiis=5
sflies='c1 '
call rdinput
, (icontrol,fcontrol,natoms,
, fgrad,fgrad_def,fhess,fhess_def,ftopo,ftopo_def,
, econv_def,gconv_def,dfac_def,
, harmonic_def,crystal_def,maxcycle_def,lennard_def,
, modus_def,iterm_def,transform_def,lnumhess_def,lmd_def,
, epssearchsmall_def,epssearchbig_def,
, epssteepsmall_def,epssteepbig_def,
, alpha_def,beta_def,gamma_def,dqmax_def,
, mxls_def,dhls_def,ahls_def,qtot,qtot_def,
, econv,gconv,dfac,harmonic,crystal,lennard,
, iterm,modus,transform,lnumhess,lmd,
, nqeq,nqeq_def,maxcycle,sflies,epssearch,epssteep,
, alpha,beta,gamma,mxls,dhls,ahls,dqmax,
, lb,la,lt,li,lnv,lne,excoord)
c ... get debug output flag dbgflg
dbgflg=0
open(icontrol,file=fcontrol,iostat=ios)
if (ios.ne.0) STOP ' Unable to open control file'
do
read (icontrol,'(a)',iostat=ios) zeile
if (ios.ne.0) exit
iauffd=index(zeile,'$uffdebug ')
if (iauffd.gt.0) then
cdbgflg=zeile(iauffd+10:)
read (cdbgflg,*) dbgflg
write (*,'(/,a,i2)') ' Debug print level is set to ',
, dbgflg
exit
end if
end do
close (icontrol)
c ... print the used terms
write (*,'(a)')
, ' The following terms are included in the Force Field !'
if (lb) write (*,'(a)')
, ' Bond terms are included !'
if (la) write (*,'(a)')
, ' Angle terms are included !'
if (lt) write (*,'(a)')
, ' Torsion terms are included !'
if (li) write (*,'(a)')
, ' Inversion terms are included !'
if (lnv) write (*,'(a)')
, ' Nonbonded terms (vdW) are included !'
if (lne) write (*,'(a)')
, ' Nonbonded terms (el) are included !'
c ................................................................
c ..... ....
c ..... Allocate required memory ....
c ..... ....
c ................................................................
c ... get amount of required disc space
nbondmx=nbondv*natoms
nanglemx=nanglev*natoms
ntorsionmx=ntorsionv*natoms
ninversionmx=ninversionv*natoms
nnonbondmx=natoms*natoms
size=0.d0
allocate (nsymeq(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for nsymeq !'
size=size+4.d0*natoms
allocate (sumcoo(3,natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for sumcoo !'
size=size+8.d0*3*natoms
allocate (xyz(3,natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for xyz !'
size=size+8.d0*3*natoms
allocate (xyzold(3,natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for xyzold !'
size=size+8.d0*3*natoms
allocate (xyzeq(3,natoms,natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for xyzeq !'
size=size+8.d0*3*natoms*natoms
allocate (ict(natoms,ndi14),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for ict !'
size=size+4.d0*natoms*ndi14
allocate (ictnr(natoms,ndi14),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for ictnr !'
size=size+4.d0*natoms*ndi14
allocate (jkseq(2,natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for jkseq !'
size=size+4.d0*2.d0*natoms
allocate (wmass(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for wmass !'
size=size+8.d0*natoms
allocate (mtype(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for mtype !'
size=size+8.d0*natoms
allocate (mtypeq(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for mtypeq !'
size=size+8.d0*natoms
allocate (ufftype(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for ufftype !'
size=size+8.d0*natoms
allocate (q(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for q !'
size=size+8.d0*natoms
allocate (chg(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for chg !'
size=size+8.d0*natoms
allocate (chgeq(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for chgeq !'
size=size+8.d0*natoms
allocate (natt(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for natt !'
size=size+4.d0*natoms
allocate (gpnr(natoms,3),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for gpnr !'
size=size+8.d0*3.d0*natoms
allocate (nxtnei12(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for nxtnei12 !'
size=size+4.d0*natoms
allocate (atoms12(natoms,natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for atoms12 !'
size=size+4.d0*natoms*natoms
allocate (nxtnei13(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for nxtnei13 !'
size=size+4.d0*natoms
allocate (atoms13(natoms,natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for atoms13 !'
size=size+4.d0*natoms*natoms
allocate (nxtnei14(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for nxtnei14 !'
size=size+4.d0*natoms
allocate (atoms14(natoms,natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for atoms14 !'
size=size+4.d0*natoms*natoms
allocate (disp(nhess),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for disp !'
size=size+8.d0*nhess
allocate (deges(3,natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for deges !'
size=size+8.d0*3.d0*natoms
allocate (degesold(3,natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for degesold !'
size=size+8.d0*3.d0*natoms
allocate (d2eges(nhess*(nhess+1)/2),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for d2eges !'
size=size+8.d0*nhess*(nhess+1)/2
allocate (d2enum(nhess*(nhess+1)/2),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for d2enum !'
size=size+8.d0*nhess*(nhess+1)/2
allocate (fmat(nhess*(nhess+1)/2),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for fmat !'
size=size+8.d0*nhess*(nhess+1)/2
allocate (bmat(nhess,nbdim),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for bmat !'
size=size+8.d0*6.d0*nhess
allocate (symd2eges(3,3,natoms*(natoms+1)/2),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for symd2eges !'
size=size+8.d0*9.d0*natoms*(natoms+1)/2
allocate (detsym(9,natoms*(natoms+1)/2),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for detsym !'
size=size+8.d0*9.d0*natoms*(natoms+1)/2
allocate (detgrd(3,natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for detgrd!'
size=size+8.d0*3.d0*natoms
allocate (ipq(3*natoms+1),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for ipq !'
size=size+4.d0*(3.d0*natoms+1)
allocate (emat(nhess,ndiis),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for emat !'
size=size+8.d0*nhess*ndiis
allocate (qold(nhess,ndiis),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for qold !'
size=size+8.d0*nhess*ndiis
allocate (gold(nhess,ndiis),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for gold !'
size=size+8.d0*nhess*ndiis
allocate (bond(nbondmx,2),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for bond !'
size=size+4.d0*2.d0*nbondmx
allocate (valbond(nbondmx,4),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for valbond !'
size=size+8.d0*4.d0*nbondmx
allocate (bobond(nbondmx),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for bobond !'
size=size+8.d0*nbondmx
allocate (angle(nanglemx,4),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for angle !'
size=size+4.d0*4.d0*nanglemx
allocate (valangle(nanglemx,7),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for valangle !'
size=size+8.d0*7.d0*nanglemx
allocate (torsion(ntorsionmx,6),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for torsion !'
size=size+4.d0*6.d0*ntorsionmx
allocate (valtorsion(ntorsionmx,6),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for valtorsion !'
size=size+8.d0*6.d0*ntorsionmx
allocate (inversion(ninversionmx,7),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for inversion !'
size=size+4.d0*7.d0*ninversionmx
allocate (valinversion(ninversionmx,15),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for valinversion !'
size=size+8.d0*15.d0*ninversionmx
allocate (nonbond(nnonbondmx,2),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for nonbond !'
size=size+4.d0*2.d0*nnonbondmx
allocate (valnonbond(nnonbondmx,6),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for valnonbond !'
size=size+8.d0*6.d0*nnonbondmx
allocate (eig(nhess),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for eig !'
size=size+8.d0*nhess
allocate (eigvec(nhess,nhess),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for eigvec !'
size=size+8.d0*nhess*nhess
allocate (lieff(natoms),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for lieff !'
size=size+4.d0*natoms
allocate (trd2eges(nhess*(nhess+1)/2),stat=ierr)
if (ierr.ne.0) STOP ':Allocation failure for proj. !'
size=size+8.d0*nhess*(nhess+1)/2
c ... disp space for automatic fields
c ... svd,v and w in uffrelax
size=size+2.d0*8.d0*nhess*nhess+8.d0*nhess
c ... automatic fields C,Jij in qeq
size=size+8.d0*natoms*natoms+8.d0*natoms*(natoms+1)/2
c ... automatic fields degesp,degesm and d2egespm in numhess
size=size+8.d0*2.d0*3.d0*natoms+8.d0*nhess*(nhess+1)/2
c ... automatic fields qmp1,gmp1,qnew,rs
size=size+8.d0*nhess+8.d0*nhess+8.d0*nhess+8.d0*(ndiis+1)+
, 8.d0*(ndiis+1)*(ndiis+1)
c ... static fields for ff parmaters
size=size+13.d0*8.d0*nelementsmx+5.d0*nelementsmx
c ... get amount of required disc space in MB
write (*,'(/,a,f8.2,a,/)')
, ' Required memory :',size/(1024.d0*1024.d0),' MB'
c ... Initializiation
fmat=0.d0
bmat=0.d0
qold=0.d0
gold=0.d0
emat=0.d0
nsymeq=0
sumcoo=0.d0
lieff=0
eig=0.d0
eigvec=0.d0
trd2eges=0.d0
detsym=0.d0
detgrd=0.d0
symd2eges=0.d0
xyz=0.d0
xyzold=0.d0
xyzeq=0.d0
deges=0.d0
degesold=0.d0
d2eges=0.d0
d2enum=0.d0
q=0.d0
chg=0.d0
chgeq=0.d0
mtype=' '
mtypeq=' '
ufftype=' '
natt=0
gpnr=0
nxtnei12=0
nxtnei13=0
nxtnei14=0
atoms12=0
atoms13=0
atoms14=0
bond=0
valbond=0.d0
bobond=0.d0
angle=0
valangle=0.d0
torsion=0
valtorsion=0.d0
inversion=0
valinversion=0.d0
nonbond=0
valnonbond=0.d0
ict=0
ictnr=0
jkseq=0
disp=0.d0
gnorm=+9.d9
ipq=0
c ................................................................
c ..... ....
c ..... get force field coordinates ....
c ..... ....
c ................................................................
c ... look if MD-Run
inquire (file=fmd,exist=ex)
if (ex) lmd=.true.
c ... read cartesian coordinates from file fcoord
call rdcoord
, (icoord,fcoord,
, xyz,mtype,natoms)
c ... arrange xyz according fixed or not
c ... (if there aren't fixed atoms then natomsfx=natoms)
call getflist
, (xyz,mtype,natoms,natomsfx,nafix)
nhessfx=3*natomsfx
c ... get nuclear charges
call getchg
, (chg,mtype,natoms)
c ... symmetry of molecule --> ict,trans
thrsym=1.d-6
call symmetry
, (sflies,thrsym,xyz,natoms,
, ict,ictnr,trans,ntrans,mi,mj,invt,ndi14,
, gen,ngen,jkseq,lieff,nhessred,natomsred)
c ... get group and period of element -> gpnr(.,1)=group
c ... gpnr(.,2)=period
c ... gpnr(.,3)=index in persys
call getgp
, (gpnr,mtype,natoms,persys)
c ... if file ftopo exists => quit
topmod=0
inquire (file=ftopo,exist=ex)
if (ex) then
open(itopo,file=ftopo)
do
read (itopo,'(a)',iostat=ios) zeile
if (ios.ne.0) exit
inxt12=index(zeile,'nxtnei12')
inxt13=index(zeile,'nxtnei13')
if (inxt12.gt.0) then
topmod=2
read (itopo,*) idummy,idummy
read (itopo,*) idummy
if (idummy.eq.-9) then
lnext=.true.
else
lnext=.false.
end if
else if (inxt13.gt.0) then
topmod=3
end if
end do
close(itopo)
if (topmod.eq.0) then
write (*,'(a)') ' : topology-file is bad !'
call flush(6)
STOP
end if
else
topmod=1
end if
modus=modus*topmod
c ... get 1-2 neigbourhood nxtnei12,atoms12,bond
if (abs(modus).eq.1) then
call getnxt12
, (xyz,mtype,natoms,nxtnei12,atoms12,
, gpnr,atrad,dfac,
, ictnr,ntrans,
, bond,valbond,nbond,nbondmx,interactiv)
else
call rdnxt12
, (itopo,ftopo,lnext,
, xyz,natoms,nxtnei12,atoms12,
, ictnr,ntrans,
, bond,valbond,nbond,nbondmx)
end if
c ... partial automatic transformation from mtype to ufftype
call autouff
, (xyz,mtype,ufftype,natoms,nxtnei12,gpnr)
c ... get force field atom type natt
call atomtype
, (ufftype,natoms,atomtypes,nelements,
, natt)
c ... bond analyse
call boana
, (mtype,ufftype,natoms,nxtnei12)
call timex('measure','Pre.Proc.')
call timex('init','Topology')
if (abs(modus).eq.1 .or. abs(modus).eq.2) then
c ... get topology of molecule
call cptopo
, (xyz,natoms,ufftype,ictnr,ntrans,ndi14,gpnr,natt,
, nxtnei12,atoms12,nxtnei13,atoms13,
, nxtnei14,atoms14,nbond,nbondmx,bond,bobond,
, nangle,nanglemx,angle,valangle,ntorsion,ntorsionmx,
, torsion,valtorsion,ninversion,ninversionmx,inversion,
, valinversion,nnonbond,nnonbondmx,nonbond,valnonbond,
, pi,nelements,th0)
c ... get partial charges of the atoms => q
call qeq
, (xyz,gpnr,mtype,natoms,natt,
, chi,J,Rqeq,Qmin,Qmax,nelements,
, q,qtot)
else
call intopo
, (itopo,ftopo,pi,modus,
, sflies,q,natoms,
, nxtnei12,atoms12,nxtnei13,atoms13,nxtnei14,atoms14,
, nbond,bond,valbond,bobond,nbondmx,
, nangle,angle,valangle,nanglemx,
, valtorsion,ntorsion,torsion,ntorsionmx,
, ninversion,inversion,valinversion,ninversionmx,
, nnonbond,nonbond,valnonbond,nnonbondmx)
end if
c ... Print the topology of the molecule
write (*,'(a,i5)') ' Molecule: natoms=',natoms
do i=1,natoms
write (*,'(a,1x,a,3(1x,f20.14),x,4(i2))')
, mtype(i),ufftype(i),xyz(1,i),xyz(2,i),xyz(3,i)
end do
write (*,'(/,a,a,i4)') ' Number of fixed atoms is ',
, 'nafix=',nafix
call outtopo
, (itopo,ftopo,pi,
, sflies,q,natoms,
, nxtnei12,atoms12,nxtnei13,atoms13,nxtnei14,atoms14,
, nbond,bond,valbond,bobond,nbondmx,
, nangle,angle,valangle,nanglemx,
, valtorsion,ntorsion,torsion,ntorsionmx,
, ninversion,inversion,valinversion,ninversionmx,
, nnonbond,nonbond,valnonbond,nnonbondmx)
call timex('measure','Topology')
c ... modus<0 => No Optimization, only print the topology
if (modus.lt.0) G O T O 222
c ................................................................
c ..... get the values of the coordinates ....
c ..... get the energy-expression parameters ....
c ................................................................
call timex('init','FF val/prms')
call getval
, (xyz,natoms,
, bond,angle,torsion,inversion,nonbond,
, valbond,valangle,valtorsion,valinversion,valnonbond,
, nbond,nangle,ntorsion,ninversion,nnonbond,
, nbondmx,nanglemx,ntorsionmx,ninversionmx,nnonbondmx)
call geteprm
, (bond,valbond,bobond,nbond,nbondmx,
, angle,valangle,nangle,nanglemx,
, torsion,valtorsion,ntorsion,ntorsionmx,
, inversion,valinversion,ninversion,ninversionmx,
, nonbond,valnonbond,nnonbond,nnonbondmx,
, lb,la,lt,li,lnv,
, nxtnei12,natt,natoms,
, radii,th0,z,d,x,chi,vsp2,vsp3,
, atomtypes,zeta,nelements,
, dconst,pi,harmonic,crystal)
call timex('measure','FF val/prms')
c ................................................................
c ..... ....
c ..... GEOMETRYOPTIMIZATION ....
c ..... ....
c ..... energy, gradients and hessian ....
c ..... ....
c ................................................................
write (*,*)
write (*,*) '.................................................'
write (*,*) ' Optimization'
write (*,*) '.................................................'
write (*,*)
c ... create control-file for relax
call creatcon
, (icontrol,fcontrol,fcoord,fhess,fgrad,
, mtype,natoms,excoord)
c ... read number of performed optimization cycles
call rdcy
, (igrad,fgrad,
, cycle0)
call setipq(ipq,nhess)
call defmss(wmass,mtype,natoms,0)
ex=.false.
lconv=.false.
c ...................................................................
c ................. Optimization Loop ...........................
c ...................................................................
do cycle=1,maxcycle
c ... get partial charges of the atoms => q
if (nqeq.gt.0) then
if (mod(cycle-1,nqeq).eq.0) then
write (*,'(a)')
, ' QEq: Calculation of partial charges !'
call timex('init','QEq')
call qeq
, (xyz,gpnr,mtype,natoms,natt,
, chi,J,Rqeq,Qmin,Qmax,nelements,
, q,qtot)
call timex('measure','QEq')
end if
end if
c ... get dipole moment from partial charges
cc call getdip
cc , (xyz,q,natoms,
cc , dipol)
c ... initialize gradient and hessian
deges=0.d0
d2eges=0.d0
c ... get total energy, cartesian gradients and cartesian hessian
call timex('init','cpegh')
call cpegh
, (eges,xyz,deges,d2eges,q,ipq,natoms,nhess,
, harmonic,crystal,lennard,
, lb,la,lt,li,lnv,lne,
, nbond,nbondmx,bond,nangle,nanglemx,angle,
, ntorsion,ntorsionmx,torsion,
, ninversion,ninversionmx,inversion,
, nnonbond,nnonbondmx,nonbond,
, pi,dconst,epsilon,
, valbond,bobond,valangle,valtorsion,
, valinversion,valnonbond,.true.,dbgflg)
call timex('measure','cpegh')
c ... delete fixed atoms in hessian
if (nafix.gt.0) then
ifx=0
do i=1,nhessfx
do jj=1,i
ifx=ifx+1
trd2eges(ifx)=d2eges(ifx)
end do
end do
else
trd2eges=d2eges
end if
c ... L2 norm GNORM of deges
gnormold=gnorm
gnorm=nol2(deges,3,natomsfx)
c ... Print the energy of the molecule
deabs=eges-egesold
write (*,*)
write (*,*)
, cycle+cycle0,' -cycle .................................'
cc write (*,'(2(a,e12.6),/,a,3(1x,d12.6))')'UFF energy = ',eges,
cc , ' |dE/dxyz| = ',gnorm,
cc , 'Dipole moment =',(dipol(i),i=1,3)
write (*,'(2(a,d12.6))')'UFF energy = ',eges,
, ' |dE/dxyz| = ',gnorm
write (*,'(a,d12.6,/)') 'UFF Delta E=',deabs
c ... Print the gradient on file fgrad
call outgrad
, (igrad,fgrad,igradx,fgradx,
, natoms,xyz,deges,mtype,
, eges,gnorm,cycle+cycle0)
c ... single point calculation
if (maxcycle.eq.1) exit
c .....................................................................
c ..... Relax geometry ....
c .....................................................................
if (gnorm.gt.epssteep) then
lsteep=.true.
loh=.false.
else
lsteep=.false.
if (gnorm.lt.1.d-3) then
loh=mod(cycle-1,10).eq.0
else
loh=.true.
end if
end if
if (loh) then
call timex('init','dsyprj')
c ... project translational & rotational space out of hessian
write (*,'(a)')
, ' Project trans. and rot. out of hessian '
trd2eges=d2eges
call dsyprj(nhessfx,nbdim,bmat,nhessfx,trd2eges)
call timex('measure','dsyprj')
end if
call timex('init','uffrelax')
lsearch=gnorm.gt.epssearch
if (lsearch .and. gnorm.gt.gnormold) then
if (ahls.lt.0.d0) then
write (*,'(a)') ' Set ahls to zero !'
ahls=0.d0
else
write (*,'(a)') ' Set ahls to -1.d0 !'
ahls=-1.d0
end if
end if
call uffrelax
, (xyz,deges,xyzold,degesold,
, fmat,ifapprox,ffapprox,
, ndiis,cycle,emat,qold,gold,
, lsearch,lsteep,mxls,ahls,dhls,dqmax,
, sflies,ntrans,ict,trans,
, trd2eges,disp,ipq,nhess,natoms,nhessfx,q,
, harmonic,crystal,lennard,
, lb,la,lt,li,lnv,lne,
, nbond,nbondmx,bond,nangle,nanglemx,angle,
, ntorsion,ntorsionmx,torsion,
, ninversion,ninversionmx,inversion,
, nnonbond,nnonbondmx,nonbond,
, pi,dconst,epsilon,
, valbond,bobond,valangle,valtorsion,
, valinversion,valnonbond,lbfgs)
call getval
, (xyz,natoms,
, bond,angle,torsion,inversion,nonbond,
, valbond,valangle,valtorsion,valinversion,valnonbond,
, nbond,nangle,ntorsion,ninversion,nnonbond,
, nbondmx,nanglemx,ntorsionmx,ninversionmx,nnonbondmx)
call timex('measure','uffrelax')
c ... test for convergence
if (gnorm.le.gconv .and. abs(deabs).le.econv) then
inquire (file=fnotconv,exist=ex1)
if (ex1) call system ('rm '//fnotconv)
open(iconv,file=fconv)
write (iconv,'(a)') 'your optimization has converged !'
write (iconv,'(a,i6,a)')
, 'convergence criteria are fulfilled in ',cycle0+cycle,
, '-cycles'
write (iconv,'(a,e16.7,a,e12.7)') 'dE/au=',deabs,
, ' |dE/dxyz|=',gnorm
close (iconv)
lconv=.true.
write (*,'(/,a,/)')
, ' Optimization sucessfully ................'
exit
end if
egesold=eges
c ... if file fstop exists => quit
inquire (file=fstop,exist=ex)
if (ex) exit
call flush(6)
end do
c ...................................................................
c .............. End of Optimization Loop .........................
c ...................................................................
if (ex) call system('rm '//fstop)
c ... write approx. inverse of cartesian hessian fmat
if (lbfgs .and. .not.lconv .and. maxcycle.gt.1) then
write (*,'(a,a)') ' Write fmat on file ',ffapprox
call outhess
, (ifapprox,ffapprox,idummy,fdummy,
, nhess,fmat,fmat,ipq,.false.)
end if
if (lmd .or. nafix.gt.0) then
write (*,'(a)')
, ' Because of presence of fixed atoms or MD run'
write (*,'(a)') ' => No postprocessing'
c ... dump current cart. coordinates on file
call outcoord
, (icoord,fcoord,itmp,ftmp,
, xyz,mtype,natoms)
c ... calculation of the hessian for the last geometry
else
c ... symmetrize xyz
if (sflies.ne.'c1 ') then
call timex('init','crtsym')
thrsym=1.d-3
motz=1
write (*,'(a,a)') ' Symmetrize coordinates to ',sflies
call crtsym(xyz,xyzeq,sumcoo,trans,nsymeq,ntrans,
, natoms,natoms,
, thrsym,chg,chgeq,mtype,mtypeq,ierr,motz)
if (ierr.ne.0) then
write(6,'(/,2x,a,/)')
, ' can''t symmetrize cartesian coordinates !'
end if
call timex('measure','crtsym')
end if
c ... transform coordinates into principal axis system
if (transform) then
write (*,*)'transforming optimized geometry into ',
, 'prinicpal axis system'
call inert (xyz,wmass,rot,traeg,natoms,mtype,.false.)
end if
c ... dump current cart. coordinates on file
call outcoord
, (icoord,fcoord,itmp,ftmp,
, xyz,mtype,natoms)
c ... calculation of the hessian for the last geometry
loh=.true.
deges=0.d0
d2eges=0.d0
call timex('init','cpegh')
call cpegh
, (eges,xyz,deges,d2eges,q,ipq,natoms,nhess,
, harmonic,crystal,lennard,
, lb,la,lt,li,lnv,lne,
, nbond,nbondmx,bond,nangle,nanglemx,angle,
, ntorsion,ntorsionmx,torsion,
, ninversion,ninversionmx,inversion,
, nnonbond,nnonbondmx,nonbond,
, pi,dconst,epsilon,
, valbond,bobond,valangle,valtorsion,
, valinversion,valnonbond,loh,dbgflg)
call timex('measure','cpegh')
c ... check hessian by comparing with numerical one
call timex('init','numhess')
if (lconv .or. maxcycle.eq.1) then
if (lnumhess) then
call timex('init','numhess')
call numhess
, (eges,xyz,deges,d2eges,d2enum,q,ipq,natoms,nhess,
, harmonic,crystal,lennard,
, lb,la,lt,li,lnv,lne,
, nbond,nbondmx,bond,nangle,nanglemx,angle,
, ntorsion,ntorsionmx,torsion,
, ninversion,ninversionmx,inversion,
, nnonbond,nnonbondmx,nonbond,
, pi,dconst,epsilon,
, valbond,bobond,valangle,valtorsion,
, valinversion,valnonbond)
call timex('measure','numhess')
write (*,'(a,a)') ' Use numerical hessian instead of',
, ' analytical one !'
d2eges=d2enum
fforceo='uffnumforceapprox'
fhesso='uffnumhessian'
else
fforceo=fforce
fhesso=fhess
end if
c ... symmetrize hessian
if (sflies.ne.'c1 ') then
detsym=0.d0
write (*,'(a,a)') ' Symmetrize hessian to ',sflies
call timex('init','symhess')
call symhess
, (symd2eges,d2eges,
, detsym,ntrans,natoms,nhess,ndi14,ict,trans,ipq)
call timex('measure','symhess')
detgrd=0.d0
write (*,'(a,a)') ' Symmetrize gradient to ',sflies
call timex('init','grdsym')
call grdsym(deges,detgrd,natoms,ntrans,natoms,ict,trans)
call timex('measure','grdsym')
end if
if (beta.ne.0.d0) then
c ... calculate eigenvalues of d2eg (d2eges will be destroyed)
call timex('init','eighess')
call eighess
, (d2eges,'d2eges',eig,eigvec,
, nhess,epsdiag)
call timex('measure','eighess')
c ... modify eigenvalues and calculate new hessian
d2eges=0.d0
call modeig
, (alpha,beta,gamma,ipq,
, eig,eigvec,d2eges,nhess)
end if
c ... project translational & rotational space out of hessian with
c ... TURBOMOLE routine dsyprj
call timex('init','dsyprj')
write (*,'(a)')
, ' Project trans. and rot. out of hessian '
trd2eges=d2eges
call dsyprj(nhessfx,nbdim,bmat,nhessfx,trd2eges)
call timex('measure','dsyprj')
c ... Print the hessian on file fhess
call outhess
, (ihess,fhesso,iforce,fforceo,
, nhess,d2eges,trd2eges,ipq,.true.)
c ... calculate eigenvalues of trd2eg (trd2eges will be destroyed)
call timex('init','eighess')
call eighess
, (trd2eges,'Hesse matric without Trans. and Rot.',eig,
, eigvec,
, nhess,epsdiag)
call timex('measure','eighess')
end if
c ... endif lconv
if (.not. lconv.and.maxcycle.ne.1) then
open (inotconv,file=fnotconv)
write (inotconv,'(a)')'your optimization has NOT converged !'
write (inotconv,'(a,i6,a)')
, 'convergence criteria are NOT fulfilled in ',
, cycle0+cycle,'-cycles'
write (inotconv,'(a,e16.7,a,e12.7)') 'dE/au=',deabs,
, ' |dE/dxyz|=',gnorm
close (inotconv)
end if
c ... end if lmd or nafix>0
end if
222 CONTINUE
call timex('init','Postproc.')
c ... print topology on file ftopo
call outtopo
, (itopo,ftopo,pi,
, sflies,q,natoms,
, nxtnei12,atoms12,nxtnei13,atoms13,nxtnei14,atoms14,
, nbond,bond,valbond,bobond,nbondmx,
, nangle,angle,valangle,nanglemx,
, valtorsion,ntorsion,torsion,ntorsionmx,
, ninversion,inversion,valinversion,ninversionmx,
, nnonbond,nonbond,valnonbond,nnonbondmx)
c ... free allocated memory
deallocate (nsymeq,stat=ierr)
deallocate (sumcoo,stat=ierr)
deallocate (lieff,stat=ierr)
deallocate (xyz,xyzold,stat=ierr)
deallocate (xyzeq,stat=ierr)
deallocate (deges,degesold,stat=ierr)
deallocate (d2eges,d2enum,symd2eges,fmat,bmat,stat=ierr)
deallocate (ipq,stat=ierr)
deallocate (mtype,stat=ierr)
deallocate (mtypeq,stat=ierr)
deallocate (ufftype,stat=ierr)
deallocate (q,stat=ierr)
deallocate (chg,stat=ierr)
deallocate (chgeq,stat=ierr)
deallocate (natt,stat=ierr)
deallocate (gpnr,stat=ierr)
deallocate (nxtnei12,stat=ierr)
deallocate (atoms12,stat=ierr)
deallocate (nxtnei13,stat=ierr)
deallocate (atoms13,stat=ierr)
deallocate (nxtnei14,stat=ierr)
deallocate (atoms14,stat=ierr)
deallocate (wmass,stat=ierr)
deallocate (ict,stat=ierr)
deallocate (ictnr,stat=ierr)
deallocate (jkseq,stat=ierr)
deallocate (bond,stat=ierr)
deallocate (valbond,stat=ierr)
deallocate (bobond,stat=ierr)
deallocate (angle,stat=ierr)
deallocate (valangle,stat=ierr)
deallocate (torsion,stat=ierr)
deallocate (valtorsion,stat=ierr)
deallocate (inversion,stat=ierr)
deallocate (valinversion,stat=ierr)
deallocate (nonbond,stat=ierr)
deallocate (valnonbond,stat=ierr)
deallocate (disp,stat=ierr)
deallocate (detsym,stat=ierr)
deallocate (detgrd,stat=ierr)
deallocate (eig,stat=ierr)
deallocate (eigvec,stat=ierr)
deallocate (trd2eges,stat=ierr)
deallocate (emat,qold,gold,stat=ierr)
call timex('measure','Postproc.')
call timex('measure','UFF (total)')
call timex('display_all','UFF')
write (*,'(a)') 'UFF ended normally'
End
KlsDMay@aol.com
Zurück