mobcal2.f

c     *****************************************************
c
      subroutine dljpot(x,y,z,pot,dpotx,dpoty,dpotz,dmax)
c
c     Subroutine to calculate L-J + ion-dipole potential.
c
      implicit double precision (a-h,m-z)
      parameter (len=1000)
      common/printswitch/ip,it,iu1,iu2,iu3,iv,im2,im4,igs
      common/constants/mu,ro,eo,pi,cang,ro2,dipol,emax,m1,m2,
     ?xe,xeo,xk,xn,mconst,correct,romax,inatom,icoord,iic
      common/charge/pcharge(len)
      common/coordinates/fx(len),fy(len),fz(len),
     ?ox(len),oy(len),oz(len)
      common/ljparameters/eolj(len),rolj(len),eox4(len),
     ?ro6lj(len),ro12lj(len),dro6(len),dro12(len)
      common/hsparameters/rhs(len),rhs2(len)
      common/trajectory/sw1,sw2,dtsf1,dtsf2,cmin,ifail,ifailc,inwr
      common/angles/theta,phi,gamma
      common/xrandom/i1,i2,i3,i4,i5,i6
      REAL pottry(2,3)
      REAL pot_mol(3)
      REAL dpotxtry(2,3)
      REAL dpotx_mol(3)
      REAL dpotytry(2,3)
      REAL dpoty_mol(3)
      REAL dpotztry(2,3)
      REAL dpotz_mol(3)
c
c     nitrogen : five charge model (Allen Tildesley, 14 page)
c     Every data from B3LYP//aug-cc-pVDZ
      dmax=2.d0*romax
      bond=1.0976d-10
      Ptfn=0.d0
      xkT=500.d0*xk
      pc=-0.4825d0
c     -1.447 D A 
      pc_center=-(pc)
c     B3LYP/aug-cc-pVDZ
c      dipolzz=2.263d-30/(2.d0*4.d0*pi*xeo)*1.0
      dipolzz=1.710d-30/(2.d0*4.d0*pi*xeo)
c      dipolzz=1.7543d-30/(2.d0*4.d0*pi*xeo)
      dipolzz=dipolzz*xe*xe 
      dipolxx=1.710d-30/(2.d0*4.d0*pi*xeo)
c      dipolxx=1.500d-30/(2.d0*4.d0*pi*xeo)*1.0
c      dipolxx=1.7543d-30/(2.d0*4.d0*pi*xeo)
      dipolxx=dipolxx*xe*xe 
      pot_min=1.0d8
c
      do 1300 isamp=1,3
      do 1200 ibatom=1,2
      rx=0.d0
      ry=0.d0
      rz=0.d0
      e00=0.d0
      de00x=0.d0
      de00y=0.d0
      de00z=0.d0
      sum1=0.d0
      sum2=0.d0
      sum3=0.d0
      sum4=0.d0
      sum5=0.d0
      sum6=0.d0
      qpol=0.d0
      dqpolx=0.d0
      dqpoly=0.d0
      dqpolz=0.d0
c
      do 1100 iatom=1,inatom
      if(isamp.eq.1) then
      xc=(bond/2.d0)*(2.d0*ibatom-3.d0)
      yc=0.d0
      zc=0.d0
      dpolx=dipolzz
      dpoly=dipolxx
      dpolz=dipolxx
      endif
      if(isamp.eq.2) then
      xc=0.d0
      yc=(bond/2.d0)*(2.d0*ibatom-3.d0)
      zc=0.d0
      dpolx=dipolxx
      dpoly=dipolzz
      dpolz=dipolxx
      endif 
      if(isamp.eq.3) then
      xc=0.d0
      yc=0.d0
      zc=(bond/2.d0)*(2.d0*ibatom-3.d0)
      dpolx=dipolxx
      dpoly=dipolxx
      dpolz=dipolzz
      endif
c
      xx_center=x-fx(iatom)
      xx=xx_center+xc
      xx_center2=xx_center*xx_center
      xx2=xx*xx
c
      yy_center=y-fy(iatom)
      yy=yy_center+yc
      yy_center2=yy_center*yy_center
      yy2=yy*yy
c
      zz_center=z-fz(iatom)
      zz=zz_center+zc
      zz_center2=zz_center*zz_center
      zz2=zz*zz
c
      rxyz_center2=xx_center2+yy_center2+zz_center2
      rxyz2=xx2+yy2+zz2
c
      rxyz_center=dsqrt(rxyz_center2)
      rxyz=dsqrt(rxyz2)
      if(rxyz.lt.dmax) dmax=rxyz
      rxyz3=rxyz2*rxyz
      rxyz5=rxyz3*rxyz2
      rxyz6=rxyz5*rxyz
      rxyz8=rxyz5*rxyz3
      rxyz12=rxyz6*rxyz6
      rxyz14=rxyz12*rxyz2
      rxyz_center3=rxyz_center2*rxyz_center
      rxyz_center5=rxyz_center3*rxyz_center2
c     LJ potential 
      e00=e00+(eox4(iatom)*((ro12lj(iatom)/rxyz12)-
     ?(ro6lj(iatom)/rxyz6)))
c     LJ derivative
      de00=eox4(iatom)*((dro6(iatom)/rxyz8)-
     ?(dro12(iatom)/rxyz14))
      de00x=de00x+(de00*xx)
      de00y=de00y+(de00*yy)
      de00z=de00z+(de00*zz)
c     ion-induced dipole potential
      if(pcharge(iatom).eq.0.d0) goto 1100
      rxyz3i=pcharge(iatom)/rxyz_center3
      rxyz5i=-3.d0*pcharge(iatom)/rxyz_center5
      rx=rx+(xx_center*rxyz3i)
      ry=ry+(yy_center*rxyz3i)
      rz=rz+(zz_center*rxyz3i)
c     ion-induced dipole derivative
      sum1=sum1+(rxyz3i+(xx_center2*rxyz5i))
      sum2=sum2+(xx_center*yy_center*rxyz5i)
      sum3=sum3+(xx_center*zz_center*rxyz5i)
      sum4=sum4+(rxyz3i+(yy_center2*rxyz5i))
      sum5=sum5+(yy_center*zz_center*rxyz5i)
      sum6=sum6+(rxyz3i+(zz_center2*rxyz5i))
c     ion-partial charge coulomb potential(quadrupole)       
      const_k=pcharge(iatom)*(xe*xe)/(4.d0*pi*xeo)
      qpol=qpol+(pc_center*const_k/rxyz_center)
      qpol=qpol+(pc*const_k/rxyz)
c     ion-partial charge coulomb derivative(quadrupole)       
      dqpolx=dqpolx-((pc_center*const_k/rxyz_center3)*(xx_center))
      dqpoly=dqpoly-((pc_center*const_k/rxyz_center3)*(yy_center))
      dqpolz=dqpolz-((pc_center*const_k/rxyz_center3)*(zz_center))
      dqpolx=dqpolx-((pc*const_k/rxyz3)*(xx))
      dqpoly=dqpoly-((pc*const_k/rxyz3)*(yy))
      dqpolz=dqpolz-((pc*const_k/rxyz3)*(zz))
c
 1100 continue
c
      pottry(ibatom,isamp)=e00-0.5*(((dpolx*rx*rx)+(dpoly*ry*ry)
     ?+(dpolz*rz*rz)))+qpol
      dpotxtry(ibatom,isamp)=de00x-0.5*((dpolx*2.d0*rx*sum1)
     ?+(dpoly*2.d0*ry*sum2)+(dpolz*2.d0*rz*sum3))+dqpolx
      dpotytry(ibatom,isamp)=de00y-0.5*((dpolx*2.d0*rx*sum2)
     ?+(dpoly*2.d0*ry*sum4)+(dpolz*2.d0*rz*sum5))+dqpoly
      dpotztry(ibatom,isamp)=de00z-0.5*((dpolx*2.d0*rx*sum3)
     ?+(dpoly*2.d0*ry*sum5)+(dpolz*2.d0*rz*sum6))+dqpolz
c 
c
 1200 continue
      pot_mol(isamp)=pottry(1,isamp)+pottry(2,isamp)
      tpot=pot_mol(isamp)
      if(pot_min.ge.pot_mol(isamp)) pot_min=pot_mol(isamp)
      dpotx_mol(isamp)=dpotxtry(1,isamp)+dpotxtry(2,isamp)
      dpoty_mol(isamp)=dpotytry(1,isamp)+dpotytry(2,isamp)
      dpotz_mol(isamp)=dpotztry(1,isamp)+dpotztry(2,isamp)
 1300 continue
c
      do 1410 isamp=1,3
      temp_pot=pot_mol(isamp)-pot_min
      Ptfn=Ptfn+dexp(-temp_pot/xkT)
 1410 continue
c
      pot=0.d0
      dpotx=0.d0
      dpoty=0.d0
      dpotz=0.d0
      do 1400 isamp=1,3
      temp_pot=pot_mol(isamp)-pot_min
      weight=dexp(-temp_pot/xkT)/Ptfn
      pot=pot+weight*pot_mol(isamp)
      dpotx=dpotx+weight*dpotx_mol(isamp)
      dpoty=dpoty+weight*dpoty_mol(isamp)
      dpotz=dpotz+weight*dpotz_mol(isamp)
 1400 continue
c
      return
      end