c c c ################################################### c ## COPYRIGHT (C) 1990 by Jay William Ponder ## c ## All Rights Reserved ## c ################################################### c c ################################################################## c ## ## c ## subroutine gradient -- find energy & gradient components ## c ## ## c ################################################################## c c c "gradient" calls subroutines to calculate the potential energy c and first derivatives with respect to Cartesian coordinates c c subroutine gradient (energy,derivs) implicit none include 'sizes.i' include 'atoms.i' include 'bound.i' include 'deriv.i' include 'energi.i' include 'inter.i' include 'potent.i' include 'rigid.i' include 'vdwpot.i' include 'virial.i' integer i,j real*8 energy,cutoff real*8 derivs(3,maxatm) integer nat0ai parameter (nat0ai=500) logical laic real*8 fcar,xeq,x1,x2,y1,y2,z1,z2,sum integer i1,i2,jp common /aic/fcar(3*nat0ai,3*nat0ai),xeq(3*nat0ai),laic c c spectral gradient: logical lexsp real*8 gradsp,xold,yold,zold,ensp,spalpha,sumspold,rshold common /spegrad/gradsp(3*maxatm),xold(maxatm),yold(maxatm), 1zold(maxatm),rshold(maxatm),ensp,lexsp,spalpha,sumspold real*8 rsh(3*maxatm) c c c zero out each of the potential energy components c ensp = 0.0d0 eb = 0.0d0 ea = 0.0d0 eba = 0.0d0 eub = 0.0d0 eaa = 0.0d0 eopb = 0.0d0 eopd = 0.0d0 eid = 0.0d0 eit = 0.0d0 et = 0.0d0 ept = 0.0d0 ebt = 0.0d0 ett = 0.0d0 ev = 0.0d0 ec = 0.0d0 ecd = 0.0d0 ed = 0.0d0 em = 0.0d0 ep = 0.0d0 er = 0.0d0 es = 0.0d0 elf = 0.0d0 eg = 0.0d0 ex = 0.0d0 c c zero out each of the first derivative components c do i = 1, n do j = 1, 3 deb(j,i) = 0.0d0 dea(j,i) = 0.0d0 deba(j,i) = 0.0d0 deub(j,i) = 0.0d0 deaa(j,i) = 0.0d0 deopb(j,i) = 0.0d0 deopd(j,i) = 0.0d0 deid(j,i) = 0.0d0 deit(j,i) = 0.0d0 det(j,i) = 0.0d0 dept(j,i) = 0.0d0 debt(j,i) = 0.0d0 dett(j,i) = 0.0d0 dev(j,i) = 0.0d0 dec(j,i) = 0.0d0 decd(j,i) = 0.0d0 ded(j,i) = 0.0d0 dem(j,i) = 0.0d0 dep(j,i) = 0.0d0 der(j,i) = 0.0d0 des(j,i) = 0.0d0 delf(j,i) = 0.0d0 deg(j,i) = 0.0d0 dex(j,i) = 0.0d0 end do end do c c zero out the virial and the intermolecular energy c do i = 1, 3 do j = 1, 3 vir(j,i) = 0.0d0 end do end do einter = 0.0d0 c c if harmonic ab initio FF: if(laic)then do j=1,n i1=3*(j-1) x1=x(j)-xeq(i1+1) y1=y(j)-xeq(i1+2) z1=z(j)-xeq(i1+3) ex=ex +fcar(i1+1,i1+1)*x1*x1 1 +fcar(i1+2,i1+2)*y1*y1 2 +fcar(i1+3,i1+3)*z1*z1 3 +2.0d0*( fcar(i1+1,i1+2)*x1*y1 4 +fcar(i1+1,i1+3)*x1*z1 5 +fcar(i1+2,i1+3)*y1*z1) dex(1,j)=dex(1,j)+fcar(i1+1,i1+1)*x1 1 +fcar(i1+1,i1+2)*y1 2 +fcar(i1+1,i1+3)*z1 dex(2,j)=dex(2,j)+fcar(i1+2,i1+1)*x1 1 +fcar(i1+2,i1+2)*y1 2 +fcar(i1+2,i1+3)*z1 dex(3,j)=dex(3,j)+fcar(i1+3,i1+1)*x1 1 +fcar(i1+3,i1+2)*y1 2 +fcar(i1+3,i1+3)*z1 sum=0.0d0 do jp=j+1,n i2=3*(jp-1) x2=x(jp)-xeq(i2+1) y2=y(jp)-xeq(i2+2) z2=z(jp)-xeq(i2+3) sum=sum+x1*(fcar(i1+1,i2+1)*x2 1 +fcar(i1+1,i2+2)*y2 2 +fcar(i1+1,i2+3)*z2) 3 +y1*(fcar(i1+2,i2+1)*x2 4 +fcar(i1+2,i2+2)*y2 5 +fcar(i1+2,i2+3)*z2) 6 +z1*(fcar(i1+3,i2+1)*x2 7 +fcar(i1+3,i2+2)*y2 8 +fcar(i1+3,i2+3)*z2) dex(1,j)=dex(1,j)+fcar(i1+1,i2+1)*x2 1 +fcar(i1+1,i2+2)*y2 2 +fcar(i1+1,i2+3)*z2 dex(2,j)=dex(2,j)+fcar(i1+2,i2+1)*x2 1 +fcar(i1+2,i2+2)*y2 2 +fcar(i1+2,i2+3)*z2 dex(3,j)=dex(3,j)+fcar(i1+3,i2+1)*x2 1 +fcar(i1+3,i2+2)*y2 2 +fcar(i1+3,i2+3)*z2 enddo do jp=1,j-1 i2=3*(jp-1) x2=x(jp)-xeq(i2+1) y2=y(jp)-xeq(i2+2) z2=z(jp)-xeq(i2+3) dex(1,j)=dex(1,j)+fcar(i1+1,i2+1)*x2 1 +fcar(i1+1,i2+2)*y2 2 +fcar(i1+1,i2+3)*z2 dex(2,j)=dex(2,j)+fcar(i1+2,i2+1)*x2 1 +fcar(i1+2,i2+2)*y2 2 +fcar(i1+2,i2+3)*z2 dex(3,j)=dex(3,j)+fcar(i1+3,i2+1)*x2 1 +fcar(i1+3,i2+2)*y2 2 +fcar(i1+3,i2+3)*z2 enddo ex=ex+2.0d0*sum enddo ex=ex/2.0d0 endif if(laic) goto 9999 c c maintain any periodic boundary conditions c if (use_bounds .and. .not.use_rigid) call bounds c c remove any previous use of the replicates method c cutoff = 0.0d0 if (use_image) call replica (cutoff) c c alter bond and torsion constants for pisystem c if (use_orbit) call piscf c c call the local geometry energy and gradient routines c c gradient from spectra comparison if(lexsp)then ensp=spalpha*sumspold endif if (use_bond) call ebond1 if (use_angle) call eangle1 if (use_strbnd) call estrbnd1 if (use_urey) call eurey1 if (use_angang) call eangang1 if (use_opbend) call eopbend1 if (use_opdist) call eopdist1 if (use_improp) call eimprop1 if (use_imptor) call eimptor1 if (use_tors) call etors1 if (use_pitors) call epitors1 if (use_strtor) call estrtor1 if (use_tortor) call etortor1 c c call the van der Waals energy and gradient routines c if (use_vdw) then if (vdwtyp .eq. 'LENNARD-JONES') call elj1 if (vdwtyp .eq. 'BUCKINGHAM') call ebuck1 if (vdwtyp .eq. 'MM3-HBOND') call emm3hb1 if (vdwtyp .eq. 'BUFFERED-14-7') call ehal1 if (vdwtyp .eq. 'GAUSSIAN') call egauss1 end if c c call the electrostatic energy and gradient routines c if (use_charge) call echarge1 if (use_chgdpl) call echgdpl1 if (use_dipole) call edipole1 if (use_mpole .or. use_polar) call empole1 if (use_rxnfld) call erxnfld1 c c call any miscellaneous energy and gradient routines c if (use_solv) call esolv1 if (use_metal) call emetal1 if (use_geom) call egeom1 if (use_extra) call extra1 c c sum up to get the total energy and first derivatives c 9999 energy = eb + ea + eba + eub + eaa + eopb + eopd + eid + eit & + et + ept + ebt + ett + ev + ec + ecd + ed + em & + ep + er + es + elf + eg + ex + ensp do i = 1, n do j = 1, 3 desum(j,i) = deb(j,i) + dea(j,i) + deba(j,i) & + deub(j,i) + deaa(j,i) + deopb(j,i) & + deopd(j,i) + deid(j,i) + deit(j,i) & + det(j,i) + dept(j,i) + debt(j,i) & + dett(j,i) + dev(j,i) + dec(j,i) & + decd(j,i) + ded(j,i) + dem(j,i) & + dep(j,i) + der(j,i) + des(j,i)+ & gradsp(j+3*(i-1))*spalpha & + delf(j,i) + deg(j,i) + dex(j,i) derivs(j,i) = desum(j,i) end do end do return end