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)
& 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)
+ eij=0.0
+ 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.
+ 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
+ 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
+ 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
+ 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)+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
+ 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)+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
+ 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)+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
+ 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