X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fsrc_MD-M%2FssMD.F;h=2e0ace739ee7f581ebb4ec496ce552ff79910d8f;hb=d96a47c690dabd9c300b08916cd2811a37c3e5eb;hp=15800ae427e473b78ecf1775da3081a838d84eb9;hpb=e3a9e82f4769b262c53a0b3d0aa297880600d179;p=unres.git diff --git a/source/unres/src_MD-M/ssMD.F b/source/unres/src_MD-M/ssMD.F index 15800ae..2e0ace7 100644 --- a/source/unres/src_MD-M/ssMD.F +++ b/source/unres/src_MD-M/ssMD.F @@ -150,11 +150,119 @@ 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 +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)-c(1,nres+i) - yj=c(2,nres+j)-c(2,nres+i) - zj=c(3,nres+j)-c(3,nres+i) + 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- + & ((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) dxj=dc_norm(1,nres+j) dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) @@ -172,6 +280,8 @@ c-------END TESTING CODE 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 @@ -187,9 +297,9 @@ 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 + ljA=ljA*aa + ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0) ssXs=d0cm deltat1=1.0d0-om1 @@ -223,7 +333,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/aa if (ssm.lt.ljm .and. & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then nicheck=1000 @@ -248,17 +358,18 @@ 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 + e2=fac*bb eij=eps1*eps2rt*eps3rt*(e1+e2) eps2der=eij*eps3rt eps3der=eij*eps2rt - eij=eij*eps2rt*eps3rt + eij=eij*eps2rt*eps3rt*sss 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 @@ -267,8 +378,9 @@ 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 @@ -309,13 +421,15 @@ 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(itypi,itypj)/aa(itypi,itypj) - d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB + ljm=-0.25D0*ljB*bb/aa + d_ljm(1)=-0.5D0*bb/aa*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) @@ -331,6 +445,8 @@ 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) @@ -433,6 +549,8 @@ 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 @@ -442,10 +560,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) + gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(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) + gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(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 @@ -456,8 +574,8 @@ cgrad enddo cgrad enddo do l=1,3 - gvdwc(l,i)=gvdwc(l,i)-gg(l) - gvdwc(l,j)=gvdwc(l,j)+gg(l) + gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k) + gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k) enddo return @@ -624,9 +742,9 @@ cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss) 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 @@ -639,9 +757,9 @@ cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss) 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 @@ -655,38 +773,7 @@ cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss) 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----------------------------------------------------------------------------- @@ -1949,3 +2036,151 @@ c$$$ return c$$$ end c$$$ c$$$C----------------------------------------------------------------------------- +c$$$C----------------------------------------------------------------------------- + subroutine triple_ssbond_ene(resi,resj,resk,eij) + include 'DIMENSIONS' + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' +#ifndef CLUST +#ifndef WHAM + include 'COMMON.MD' +#endif +#endif + +c External functions + double precision h_base + external h_base + +c Input arguments + integer resi,resj,resk + +c Output arguments + double precision eij,eij1,eij2,eij3 + +c Local variables + logical havebond +c integer itypi,itypj,k,l + double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi + double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij + double precision xik,yik,zik,xjk,yjk,zjk + double precision sig0ij,ljd,sig,fac,e1,e2 + double precision dcosom1(3),dcosom2(3),ed + double precision pom1,pom2 + double precision ljA,ljB,ljXs + double precision d_ljB(1:3) + 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) + if (dtriss.eq.0) return + i=resi + j=resj + k=resk +C write(iout,*) resi,resj,resk + itypi=itype(i) + dxi=dc_norm(1,nres+i) + 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) + + itypj=itype(j) + xj=c(1,nres+j) + yj=c(2,nres+j) + zj=c(3,nres+j) + + dxj=dc_norm(1,nres+j) + dyj=dc_norm(2,nres+j) + dzj=dc_norm(3,nres+j) + dscj_inv=vbld_inv(j+nres) + itypk=itype(k) + xk=c(1,nres+k) + yk=c(2,nres+k) + zk=c(3,nres+k) + + dxk=dc_norm(1,nres+k) + dyk=dc_norm(2,nres+k) + dzk=dc_norm(3,nres+k) + dscj_inv=vbld_inv(k+nres) + xij=xj-xi + xik=xk-xi + xjk=xk-xj + yij=yj-yi + yik=yk-yi + yjk=yk-yj + zij=zj-zi + zik=zk-zi + zjk=zk-zj + rrij=(xij*xij+yij*yij+zij*zij) + rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse + rrik=(xik*xik+yik*yik+zik*zik) + rik=dsqrt(rrik) + rrjk=(xjk*xjk+yjk*yjk+zjk*zjk) + rjk=dsqrt(rrjk) +C there are three combination of distances for each trisulfide bonds +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) +C second case jth atom is center + eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss) +C the third case kth atom is the center + eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss) +C eij2=0.0 +C eij3=0.0 +C eij1=0.0 + eij=eij1+eij2+eij3 +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)) + gg(1)=xij*fac/rij + gg(2)=yij*fac/rij + gg(3)=zij*fac/rij + do m=1,3 + gvdwx(m,i)=gvdwx(m,i)-gg(m) + gvdwx(m,j)=gvdwx(m,j)+gg(m) + enddo + do l=1,3 + gvdwc(l,i)=gvdwc(l,i)-gg(l) + 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)) + gg(1)=xik*fac/rik + gg(2)=yik*fac/rik + gg(3)=zik*fac/rik + do m=1,3 + gvdwx(m,i)=gvdwx(m,i)-gg(m) + gvdwx(m,k)=gvdwx(m,k)+gg(m) + enddo + do l=1,3 + gvdwc(l,i)=gvdwc(l,i)-gg(l) + 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)) + gg(1)=xjk*fac/rjk + gg(2)=yjk*fac/rjk + gg(3)=zjk*fac/rjk + do m=1,3 + gvdwx(m,j)=gvdwx(m,j)-gg(m) + gvdwx(m,k)=gvdwx(m,k)+gg(m) + enddo + do l=1,3 + gvdwc(l,j)=gvdwc(l,j)-gg(l) + gvdwc(l,k)=gvdwc(l,k)+gg(l) + enddo + return + end