C-----------------------------------------------------------------------------
subroutine dyn_ssbond_ene(resi,resj,eij)
-c implicit none
-
-c Includes
+ implicit none
include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+ include 'COMMON.SETUP'
+ integer ierr
+#endif
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.NAMES'
+ include 'COMMON.SPLITELE'
#ifndef CLUST
#ifndef WHAM
include 'COMMON.MD'
double precision omega,delta_inv,deltasq_inv,fac1,fac2
c-------FIRST METHOD
double precision xm,d_xm(1:3)
+ double precision sslipi,sslipj,ssgradlipi,ssgradlipj
+ integer ici,icj,itypi,itypj
+ double precision boxshift,sscale,sscagrad
c-------END FIRST METHOD
c-------SECOND METHOD
c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
i=resi
j=resj
-
+ ici=icys(i)
+ icj=icys(j)
+c write (iout,*) "dyn_ssbond",resi,resj,ici,icj
+c call flush(iout)
+ if (ici.eq.0 .or. icj.eq.0) then
+#ifdef MPI
+ write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)')
+ & "Processor",me," attempt to create",
+ & " a disulfide link between non-cysteine residues ",restyp(i),i,
+ & restyp(j),j
+ call MPI_Abort(MPI_COMM_WORLD,ierr)
+#else
+ write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)')
+ & "Processor",me," attempt to create",
+ & " a disulfide link between non-cysteine residues ",restyp(i),i,
+ & restyp(j),j
+ stop
+#endif
+ endif
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)
- 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
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+c write (iout,*) "After to_box i",xi,yi,zi
+c call flush(iout)
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
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+c write (iout,*) "After lipid_layer"
+c call flush(iout)
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
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+c write (iout,*) "After to_box j",xj,yj,zj
+c call flush(iout)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+c write (iout,*) "After lipid_layer"
+c call flush(iout)
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=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+c write (iout,*) "After boxshift"
+c call flush(iout)
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))
+ sss=sscale((1.0d0/rij)/sigma(itypi,itypj),r_cut_int)
+ sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj),r_cut_int)
c The following are set in sc_angular
c erij(1)=xj*rij
c erij(2)=yj*rij
c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
c om12=dxi*dxj+dyi*dyj+dzi*dzj
+c write (iout,*) "Calling sc_angular"
+c call flush(iout)
call sc_angular
+c write (iout,*) "After sc_angular"
+c call flush(iout)
rij=1.0D0/rij ! Reset this so it makes sense
sig0ij=sigma(itypi,itypj)
c-------TESTING CODE
c Stop and plot energy and derivative as a function of distance
+c write (iout,*) "checkstop",checkstop
+c call flush(iout)
if (checkstop) then
ssm=ssC-0.25D0*ssB*ssB/ssA
ljm=-0.25D0*ljB*bb/aa
endif
+c write (iout,*) "havebond",havebond
+c call flush(iout)
if (havebond) then
#ifndef CLUST
#ifndef WHAM
c endif
#endif
#endif
- dyn_ssbond_ij(i,j)=eij
- else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
- dyn_ssbond_ij(i,j)=1.0d300
+ dyn_ssbond_ij(ici,icj)=eij
+ else if (.not.havebond .and. dyn_ssbond_ij(ici,icj).lt.1.0d300)
+ &then
+ dyn_ssbond_ij(ici,icj)=1.0d300
#ifndef CLUST
#ifndef WHAM
c write(iout,'(a15,f12.2,f8.1,2i5)')
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)+gg_lipi(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
enddo
return
c Local variables
double precision emin
integer i,j,imin
- integer diff,allflag(maxdim),allnss,
- & allihpb(maxdim),alljhpb(maxdim),
- & newnss,newihpb(maxdim),newjhpb(maxdim)
+ integer diff,allflag(maxss),allnss,
+ & allihpb(maxss),alljhpb(maxss),
+ & newnss,newihpb(maxss),newjhpb(maxss)
logical found
integer i_newnss(max_fg_procs),displ(0:max_fg_procs)
- integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
+ integer g_newihpb(maxss),g_newjhpb(maxss),g_newnss
allnss=0
- do i=1,nres-1
- do j=i+1,nres
+ do i=1,ns-1
+ do j=i+1,ns
if (dyn_ssbond_ij(i,j).lt.1.0d300) then
allnss=allnss+1
allflag(allnss)=0
c$$$
c$$$C-----------------------------------------------------------------------------
c$$$C-----------------------------------------------------------------------------
- subroutine triple_ssbond_ene(resi,resj,resk,eij)
+ subroutine triple_ssbond_ene(resi,resj,resk,eij)
include 'DIMENSIONS'
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
projects / unres.git / blobdiff