common /sschecks/ checkstop,transgrad
integer icheck,nicheck,jcheck,njcheck
- double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi
+ double precision echeck(-1:1),deps,ssx0,ljx0
c-------END TESTING CODE
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
dsci_inv=vbld_inv(i+nres)
- xi=c(1,nres+i)
- yi=c(2,nres+i)
- zi=c(3,nres+i)
- xi=dmod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=dmod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=dmod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
-C define scaling factor for lipids
-
-C if (positi.le.0) positi=positi+boxzsize
-C print *,i
-C first for peptide groups
-c for each residue check if it is in lipid or lipid water border area
- if ((zi.gt.bordlipbot)
- &.and.(zi.lt.bordliptop)) then
-C the energy transfer exist
- if (zi.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((positi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipi=1.0d0
- ssgradlipi=0.0
- endif
- else
- sslipi=0.0d0
- ssgradlipi=0.0
- endif
+
itypj=itype(j)
- xj=c(1,nres+j)
- yj=c(2,nres+j)
- zj=c(3,nres+j)
- xj=dmod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=dmod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=dmod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- if ((zj.gt.bordlipbot)
- &.and.(zj.lt.bordliptop)) then
-C the energy transfer exist
- if (zj.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((positi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipj=1.0d0
- ssgradlipj=0.0
- endif
- else
- sslipj=0.0d0
- ssgradlipj=0.0
- endif
- aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
- bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
-
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
-
-C xj=c(1,nres+j)-c(1,nres+i)
-C yj=c(2,nres+j)-c(2,nres+i)
-C zj=c(3,nres+j)-c(3,nres+i)
+ xj=c(1,nres+j)-c(1,nres+i)
+ yj=c(2,nres+j)-c(2,nres+i)
+ zj=c(3,nres+j)-c(3,nres+i)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
- sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
- sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
c The following are set in sc_angular
c erij(1)=xj*rij
c erij(2)=yj*rij
ljXs=sig-sig0ij
ljA=eps1*eps2rt**2*eps3rt**2
- ljB=ljA*bb
- ljA=ljA*aa
- ljxm=ljXs+(-2.0D0*aa/bb)**(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/aa
+ 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
- e2=fac*bb
+ 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
- eij=eij*eps2rt*eps3rt*sss
+ eij=eij*eps2rt*eps3rt
sigder=-sig/sigsq
e1=e1*eps1*eps2rt**2*eps3rt**2
ed=-expon*(e1+eij)/ljd
sigder=ed*sigder
- ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
eom12=eij*eps1_om12+eps2der*eps2rt_om12
havebond=.true.
ssd=rij-ssXs
eij=ssA*ssd*ssd+ssB*ssd+ssC
- eij=eij*sss
+
ed=2*akcm*ssd+akct*deltat12
- ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
pom1=akct*ssd
pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
eom1=-2*akth*deltat1-pom1-om2*pom2
fac1=deltasq_inv*fac*(xm-rij)
fac2=deltasq_inv*fac*(rij-ssxm)
ed=delta_inv*(Ht*hd2-ssm*hd1)
- eij=eij*sss
- ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
else
havebond=.false.
- ljm=-0.25D0*ljB*bb/aa
- d_ljm(1)=-0.5D0*bb/aa*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)
fac1=deltasq_inv*fac*(ljxm-rij)
fac2=deltasq_inv*fac*(rij-xm)
ed=delta_inv*(ljm*hd2-Ht*hd1)
- eij=eij*sss
- ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
checkstop=.false.
endif
c-------END TESTING CODE
- gg_lipi(3)=ssgradlipi*eij
- gg_lipj(3)=ssgradlipj*eij
do k=1,3
dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
enddo
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)
& +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
& +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
- gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)
& +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
& +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
enddo
cgrad enddo
do l=1,3
- gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
- gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
enddo
return
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
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
& +ft(1)*welec*ees
& +ft(1)*wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
-#ifndef CLUST
-#ifndef WHAM
+C#ifndef CLUST
+C#ifndef WHAM
C include 'COMMON.MD'
-#endif
-#endif
+C#endif
+C#endif
c External functions
double precision h_base
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
dsci_inv=vbld_inv(i+nres)
- xi=c(1,nres+i)
- yi=c(2,nres+i)
- zi=c(3,nres+i)
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
- if ((zi.gt.bordlipbot)
- &.and.(zi.lt.bordliptop)) then
-C the energy transfer exist
- if (zi.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipi=1.0d0
- ssgradlipi=0.0
- endif
- else
- sslipi=0.0d0
- ssgradlipi=0.0
- endif
+
itypj=itype(j)
- xj=c(1,nres+j)
- yj=c(2,nres+j)
- zj=c(3,nres+j)
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- if ((zj.gt.bordlipbot)
- &.and.(zj.lt.bordliptop)) then
-C the energy transfer exist
- if (zj.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((zj-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zj.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipj=1.0d0
- ssgradlipj=0.0
- endif
- else
- sslipj=0.0d0
- ssgradlipj=0.0
- endif
- aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
- bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
- xj=xj-xi
- yj=yj-yi
- zj=zj-zi
+ xj=c(1,nres+j)-c(1,nres+i)
+ yj=c(2,nres+j)-c(2,nres+i)
+ zj=c(3,nres+j)-c(3,nres+i)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
ljXs=sig-sig0ij
ljA=eps1*eps2rt**2*eps3rt**2
- ljB=ljA*bb
- ljA=ljA*aa
- ljxm=ljXs+(-2.0D0*aa/bb)**(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/aa
+ 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
- e2=fac*bb
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
eij=eps1*eps2rt*eps3rt*(e1+e2)
-C write(iout,*) eij,'TU?1'
eps2der=eij*eps3rt
eps3der=eij*eps2rt
eij=eij*eps2rt*eps3rt
havebond=.true.
ssd=rij-ssXs
eij=ssA*ssd*ssd+ssB*ssd+ssC
-C write(iout,*) 'TU?2',ssc,ssd
+
ed=2*akcm*ssd+akct*deltat12
pom1=akct*ssd
pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
h1=h_base(f1,hd1)
h2=h_base(f2,hd2)
eij=ssm*h1+Ht*h2
-C write(iout,*) eij,'TU?3'
delta_inv=1.0d0/(xm-ssxm)
deltasq_inv=delta_inv*delta_inv
fac=ssm*hd1-Ht*hd2
eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
else
havebond=.false.
- ljm=-0.25D0*ljB*bb/aa
- d_ljm(1)=-0.5D0*bb/aa*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)
h1=h_base(f1,hd1)
h2=h_base(f2,hd2)
eij=Ht*h1+ljm*h2
-C write(iout,*) 'TU?4',ssA
delta_inv=1.0d0/(ljxm-xm)
deltasq_inv=delta_inv*delta_inv
fac=Ht*hd1-ljm*hd2
c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
endif
- write(iout,*) 'havebond',havebond
+
if (havebond) then
#ifndef CLUST
#ifndef WHAM
& allihpb(maxdim),alljhpb(maxdim),
& newnss,newihpb(maxdim),newjhpb(maxdim)
logical found
- integer i_newnss(1024),displ(0:1024)
+ integer i_newnss(256),displ(0:256)
integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
allnss=0
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-----------------------------------------------------------------------------
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
& +ft(1)*welec*ees
& +ft(1)*wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(1)*welec*ees
& +ft(1)*wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(1)*welec*ees
& +ft(1)*wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
projects / unres.git / commitdiff