changes

[unres.git] / source / unres / src-HCD-5D / ssMD.F

diff --git a/source/unres/src-HCD-5D/ssMD.F b/source/unres/src-HCD-5D/ssMD.F
index 26807a0..67e9f10 100644 (file)
--- a/source/unres/src-HCD-5D/ssMD.F
+++ b/source/unres/src-HCD-5D/ssMD.F
@@ -100,6 +100,7 @@ C-----------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
       include 'COMMON.NAMES'
+      include 'COMMON.SPLITELE'
 #ifndef CLUST
 #ifndef WHAM
       include 'COMMON.MD'
@@ -153,6 +154,8 @@ c-------END TESTING CODE
       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)') 
@@ -177,6 +180,8 @@ c-------END TESTING CODE
       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
@@ -184,12 +189,18 @@ C        print *,i
 C first for peptide groups
 c for each residue check if it is in lipid or lipid water border area
       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)
       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
@@ -197,6 +208,8 @@ c for each residue check if it is in lipid or lipid water border area
       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)
@@ -214,8 +227,8 @@ c for each residue check if it is in lipid or lipid water border area
 
       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
@@ -223,7 +236,11 @@ c      erij(3)=zj*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)
@@ -265,6 +282,8 @@ c-------END TESTING CODE
 
 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
@@ -447,6 +466,8 @@ c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
 
       endif
 
+c      write (iout,*) "havebond",havebond
+c      call flush(iout)
       if (havebond) then
 #ifndef CLUST
 #ifndef WHAM
@@ -509,8 +530,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)+gg_lipi(l)
+        gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
       enddo
 
       return