X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2FssMD.F;h=387225737e260f124164e858f7d0b09fae3b7a26;hb=7d64cc3ff0edffb6aa37e309e4375f58bd5875a2;hp=f59466c9099cc5d4e29cad7402a494c927491c44;hpb=0e3666de2f89a5360d30e6ae5e18f8d7856daacd;p=unres.git diff --git a/source/unres/src_MD-M/ssMD.F b/source/unres/src_MD-M/ssMD.F index f59466c..3872257 100644 --- a/source/unres/src_MD-M/ssMD.F +++ b/source/unres/src_MD-M/ssMD.F @@ -187,9 +187,10 @@ c om12=dxi*dxj+dyi*dyj+dzi*dzj ljXs=sig-sig0ij ljA=eps1*eps2rt**2*eps3rt**2 - ljB=ljA*bb(itypi,itypj) - ljA=ljA*aa(itypi,itypj) - ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) + ljB=ljA*bb_aq(itypi,itypj) + ljA=ljA*aa_aq(itypi,itypj) + ljxm=ljXs+(-2.0D0*aa_aq(itypi,itypj)/ + & bb_aq(itypi,itypj))**(1.0D0/6.0D0) ssXs=d0cm deltat1=1.0d0-om1 @@ -223,7 +224,7 @@ c-------TESTING CODE c Stop and plot energy and derivative as a function of distance if (checkstop) then ssm=ssC-0.25D0*ssB*ssB/ssA - ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) + ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj) if (ssm.lt.ljm .and. & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then nicheck=1000 @@ -248,8 +249,8 @@ c-------END TESTING CODE havebond=.false. ljd=rij-ljXs fac=(1.0D0/ljd)**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) + e1=fac*fac*aa_aq(itypi,itypj) + e2=fac*bb_aq(itypi,itypj) eij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=eij*eps3rt eps3der=eij*eps2rt @@ -314,8 +315,8 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3) else havebond=.false. - ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) - d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB + ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj) + d_ljm(1)=-0.5D0*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)*ljB d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt) d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt- + alf12/eps3rt) @@ -655,38 +656,7 @@ c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i) return end -c---------------------------------------------------------------------------- - -#ifdef WHAM - subroutine read_ssHist - implicit none - -c Includes - include 'DIMENSIONS' - include "DIMENSIONS.FREE" - include 'COMMON.FREE' - -c Local variables - integer i,j - character*80 controlcard - - do i=1,dyn_nssHist - call card_concat(controlcard,.true.) - read(controlcard,*) - & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0)) - enddo - - return - end -#endif -c---------------------------------------------------------------------------- - - -C----------------------------------------------------------------------------- -C----------------------------------------------------------------------------- -C----------------------------------------------------------------------------- -C----------------------------------------------------------------------------- C----------------------------------------------------------------------------- C----------------------------------------------------------------------------- C----------------------------------------------------------------------------- @@ -1990,7 +1960,8 @@ c integer itypi,itypj,k,l double precision ssA,ssB,ssC,ssXs double precision ssxm,ljxm,ssm,ljm double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3) - + eij=0.0 + if (dtriss.eq.0) return i=resi j=resj k=resk @@ -2042,11 +2013,23 @@ C The first case the ith atom is the center C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first C distance y is second distance the a,b,c,d are parameters derived for C this problem d parameter was set as a penalty currenlty set to 1. - eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss) + if ((iabs(j-i).le.2).or.(iabs(i-k).le.2)) then + eij1=0.0d0 + else + eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**6+ctriss) + endif C second case jth atom is center - eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss) + if ((iabs(j-i).le.2).or.(iabs(j-k).le.2)) then + eij2=0.0d0 + else + eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**6+ctriss) + endif C the third case kth atom is the center - eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss) + if ((iabs(i-k).le.2).or.(iabs(j-k).le.2)) then + eij3=0.0d0 + else + eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**6+ctriss) + endif C eij2=0.0 C eij3=0.0 C eij1=0.0 @@ -2054,8 +2037,8 @@ C eij1=0.0 C write(iout,*)i,j,k,eij C The energy penalty calculated now time for the gradient part C derivative over rij - fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik)) - &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk)) + fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) + &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5) gg(1)=xij*fac/rij gg(2)=yij*fac/rij gg(3)=zij*fac/rij @@ -2068,8 +2051,9 @@ C derivative over rij gvdwc(l,j)=gvdwc(l,j)+gg(l) enddo C now derivative over rik - fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik)) - &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk)) + fac=-eij1**2/dtriss* + &(-2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) + &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5) gg(1)=xik*fac/rik gg(2)=yik*fac/rik gg(3)=zik*fac/rik @@ -2082,8 +2066,9 @@ C now derivative over rik gvdwc(l,k)=gvdwc(l,k)+gg(l) enddo C now derivative over rjk - fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))- - &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk)) + fac=-eij2**2/dtriss* + &(-2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)- + &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5) gg(1)=xjk*fac/rjk gg(2)=yjk*fac/rjk gg(3)=zjk*fac/rjk