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
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
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
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)
enddo
#ifndef CLUST
#ifndef WHAM
- if (.not.found.and.fg_rank.eq.0)
- & write(iout,'(a15,f12.2,f8.1,2i5)')
- & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
+c if (.not.found.and.fg_rank.eq.0)
+c & write(iout,'(a15,f12.2,f8.1,2i5)')
+c & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
#endif
#endif
enddo
enddo
#ifndef CLUST
#ifndef WHAM
- if (.not.found.and.fg_rank.eq.0)
- & write(iout,'(a15,f12.2,f8.1,2i5)')
- & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
+c if (.not.found.and.fg_rank.eq.0)
+c & write(iout,'(a15,f12.2,f8.1,2i5)')
+c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
#endif
#endif
enddo
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-----------------------------------------------------------------------------
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
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
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
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
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
projects / unres.git / blobdiff