From: Adam Sieradzan Date: Wed, 13 Jul 2016 08:41:31 +0000 (+0200) Subject: poprawka w make (usuniecie nss) X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=661698da20bbce5a332222e67a7bb4c67df3ce6c poprawka w make (usuniecie nss) --- diff --git a/source/unres/src_MD-M/ssMD.F b/source/unres/src_MD-M/ssMD.F index aa938b5..3872257 100644 --- a/source/unres/src_MD-M/ssMD.F +++ b/source/unres/src_MD-M/ssMD.F @@ -138,7 +138,7 @@ c-------TESTING CODE 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 @@ -150,119 +150,11 @@ 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) @@ -280,8 +172,6 @@ C zj=c(3,nres+j)-c(3,nres+i) 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 @@ -297,9 +187,10 @@ c om12=dxi*dxj+dyi*dyj+dzi*dzj 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 @@ -333,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/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 @@ -358,18 +249,17 @@ c-------END TESTING CODE 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 @@ -378,9 +268,8 @@ c-------END TESTING CODE 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 @@ -421,15 +310,13 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE 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) @@ -445,8 +332,6 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE 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) @@ -549,8 +434,6 @@ c-------TESTING CODE 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 @@ -560,10 +443,10 @@ c-------END TESTING CODE 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 @@ -574,8 +457,8 @@ cgrad 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 @@ -2077,6 +1960,7 @@ 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 @@ -2129,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 @@ -2141,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 @@ -2155,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 @@ -2169,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 diff --git a/source/wham/src-M/make_ensemble1.F b/source/wham/src-M/make_ensemble1.F index 959c61a..80538f6 100644 --- a/source/wham/src-M/make_ensemble1.F +++ b/source/wham/src-M/make_ensemble1.F @@ -169,7 +169,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +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 @@ -179,7 +179,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 @@ -191,7 +191,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 @@ -201,7 +201,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 diff --git a/source/wham/src-M/ssMD.F b/source/wham/src-M/ssMD.F index f298c98..3811a69 100644 --- a/source/wham/src-M/ssMD.F +++ b/source/wham/src-M/ssMD.F @@ -96,11 +96,11 @@ c Includes 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 @@ -150,76 +150,11 @@ 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=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) @@ -252,9 +187,10 @@ c om12=dxi*dxj+dyi*dyj+dzi*dzj 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 @@ -288,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/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 @@ -313,10 +249,9 @@ c-------END TESTING CODE 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 @@ -333,7 +268,7 @@ C write(iout,*) eij,'TU?1' 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 @@ -369,7 +304,6 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE 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 @@ -381,8 +315,8 @@ C write(iout,*) eij,'TU?3' 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) @@ -392,7 +326,6 @@ C write(iout,*) eij,'TU?3' 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 @@ -464,7 +397,7 @@ c$$$ if (ed.gt.0.0d0) havebond=.true. c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE endif - write(iout,*) 'havebond',havebond + if (havebond) then #ifndef CLUST #ifndef WHAM @@ -600,7 +533,7 @@ c Local variables & 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 @@ -723,41 +656,10 @@ 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----------------------------------------------------------------------------- c$$$c----------------------------------------------------------------------------- c$$$ @@ -2058,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 @@ -2110,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 @@ -2122,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 @@ -2136,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 @@ -2150,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 diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F index f0180ba..a3d1658 100644 --- a/source/wham/src-M/wham_calc1.F +++ b/source/wham/src-M/wham_calc1.F @@ -338,7 +338,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +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 @@ -348,7 +348,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 @@ -360,7 +360,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 @@ -370,7 +370,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 @@ -699,7 +699,7 @@ C edihcnstr=0.0d0 & +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 @@ -709,7 +709,7 @@ C edihcnstr=0.0d0 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 @@ -721,7 +721,7 @@ C edihcnstr=0.0d0 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 @@ -731,7 +731,7 @@ C edihcnstr=0.0d0 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 @@ -1081,7 +1081,7 @@ c write (iout,*) "ftbis",ftbis & +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 @@ -1112,7 +1112,7 @@ C & +ftprim(6)*evdw_t 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 @@ -1137,7 +1137,7 @@ C & +ftprim(6)*evdw_t 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 @@ -1164,7 +1164,7 @@ C & +ftprim(6)*evdw_t 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