Fixed eello5, eello6, eturn6, and shortrange RESPA

[unres.git] / source / unres / src_MD-M / energy_p_new-sep_barrier.F

diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
index c0b4c84..1c6470d 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
@@ -885,14 +885,16 @@ C Calculate SC interaction energy.
 C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
+c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+c            call flush(iout)
             ind=ind+1
             itypj=itype(j)
             if (itypj.eq.ntyp1) cycle
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
-c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
+c            write (iout,*) "j",j,itypi,itypj,dsc_inv(itypj),dscj_inv,
 c     &       1.0d0/vbld(j+nres)
-c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+c            call flush(iout)
             sig0ij=sigma(itypi,itypj)
             chi1=chi(itypi,itypj)
             chi2=chi(itypj,itypi)
@@ -911,6 +913,8 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
           yj=mod(yj,boxysize)
           if (yj.lt.0) yj=yj+boxysize
           zj=mod(zj,boxzsize)
+c          write (iout,*) "After box"
+c          call flush(iout)
           if (zj.lt.0) zj=zj+boxzsize
        if ((zj.gt.bordlipbot)
      &.and.(zj.lt.bordliptop)) then
@@ -934,6 +938,8 @@ C lipbufthick is thickenes of lipid buffore
          sslipj=0.0d0
          ssgradlipj=0.0
        endif
+c          write (iout,*) "After lipid"
+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
@@ -981,6 +987,8 @@ C lipbufthick is thickenes of lipid buffore
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
               call sc_angular
+c          write (iout,*) "After sc_angular"
+c          call flush(iout)
               sigsq=1.0D0/sigsq
               sig=sig0ij*dsqrt(sigsq)
               rij_shift=1.0D0/rij-sig+sig0ij
@@ -989,9 +997,9 @@ c              rij_shift=1.2*sig0ij
 C I hate to put IF's in the loops, but here don't have another choice!!!!
               if (rij_shift.le.0.0D0) then
                 evdw=1.0D20
-cd                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
-cd     &          restyp(itypi),i,restyp(itypj),j,
-cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
+c                write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+c     &          restyp(itypi),i,restyp(itypj),j,
+c     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
                 return
               endif
               sigder=-sig*sigsq
@@ -1005,6 +1013,7 @@ c---------------------------------------------------------------
               eps3der=evdwij*eps2rt
 c              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
+c              call flush(iout)
               evdwij=evdwij*eps2rt*eps3rt
               evdw=evdw+evdwij*sss
               if (lprn) then
@@ -1016,10 +1025,13 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
      &          eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
      &          om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
      &          evdwij
+                call flush(iout)
               endif
 
-              if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
-     &                        'evdw',i,j,evdwij
+              if (energy_dec) then
+                write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij
+                call flush(iout)
+              endif
 
 C Calculate gradient components.
               e1=e1*eps1*eps2rt**2*eps3rt**2
@@ -1034,8 +1046,12 @@ C Calculate the radial part of the gradient
               gg(3)=zj*fac
               gg_lipi(3)=ssgradlipi*evdwij
               gg_lipj(3)=ssgradlipj*evdwij
+c              write (iout,*) "Calling sc_grad_scale"
+c              call flush(iout)
 C Calculate angular part of the gradient.
               call sc_grad_scale(sss)
+c              write (iout,*) "After sc_grad_scale"
+c              call flush(iout)
             endif
           enddo      ! j
         enddo        ! iint
@@ -1298,6 +1314,8 @@ C----------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       double precision dcosom1(3),dcosom2(3)
       double precision scalfac
+c      write (iout,*) "sc_grad_scale",i,j,k,l
+c      call flush(iout)
       eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
       eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
       eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
@@ -1335,8 +1353,8 @@ C
 C Calculate the components of the gradient in DC and X
 C
       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
       end