common /sschecks/ checkstop,transgrad
integer icheck,nicheck,jcheck,njcheck
- double precision echeck(-1:1),deps,ssx0,ljx0
+ double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi
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)-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=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)
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(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
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
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
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(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)
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)
+ 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
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
& allihpb(maxdim),alljhpb(maxdim),
& newnss,newihpb(maxdim),newjhpb(maxdim)
logical found
- integer i_newnss(max_fg_procs),displ(max_fg_procs)
+ integer i_newnss(max_fg_procs),displ(0:max_fg_procs)
integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
allnss=0
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-----------------------------------------------------------------------------
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
projects / unres.git / blobdiff