bug fix in respa and readpdb

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 3815c13..0a638ef 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
@@ -662,6 +662,29 @@ c     if (icall.eq.0) lprn=.false.
           if (yi.lt.0) yi=yi+boxysize
           zi=mod(zi,boxzsize)
           if (zi.lt.0) zi=zi+boxzsize
+       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-
+     &        ((zi-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-zi)/lipbufthick)
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipi=1.0d0
+         ssgradlipi=0.0
+        endif
+       else
+         sslipi=0.0d0
+         ssgradlipi=0.0
+       endif
+
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -707,12 +730,12 @@ C the energy transfer exist
         if (zj.lt.buflipbot) then
 C what fraction I am in
          fracinbuf=1.0d0-
-     &        ((positi-bordlipbot)/lipbufthick)
+     &        ((zj-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)
+        elseif (zj.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
          sslipj=sscalelip(fracinbuf)
          ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
         else
@@ -723,10 +746,11 @@ C lipbufthick is thickenes of lipid buffore
          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
+     &  +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
+     &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
 
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
@@ -788,6 +812,7 @@ cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 c---------------------------------------------------------------
               rij_shift=1.0D0/rij_shift 
               fac=rij_shift**expon
+              faclip=fac
               e1=fac*fac*aa
               e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
@@ -822,8 +847,12 @@ C Calculate the radial part of the gradient
               gg(1)=xj*fac
               gg(2)=yj*fac
               gg(3)=zj*fac
-              gg_lipi(3)=ssgradlipi*evdwij
-              gg_lipj(3)=ssgradlipj*evdwij
+            gg_lipi(3)=eps1*(eps2rt*eps2rt)
+     &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+     & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+     &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+            gg_lipj(3)=ssgradlipj*gg_lipi(3)
+            gg_lipi(3)=gg_lipi(3)*ssgradlipi
 C Calculate angular part of the gradient.
               call sc_grad_scale(1.0d0-sss)
             endif
@@ -853,6 +882,7 @@ C
       include 'COMMON.CALC'
       include 'COMMON.CONTROL'
       logical lprn
+      integer xshift,yshift,zshift
       evdw=0.0D0
 ccccc      energy_dec=.false.
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -873,6 +903,29 @@ c     if (icall.eq.0) lprn=.false.
           if (yi.lt.0) yi=yi+boxysize
           zi=mod(zi,boxzsize)
           if (zi.lt.0) zi=zi+boxzsize
+       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-
+     &        ((zi-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-zi)/lipbufthick)
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipi=1.0d0
+         ssgradlipi=0.0
+        endif
+       else
+         sslipi=0.0d0
+         ssgradlipi=0.0
+       endif
+
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -918,12 +971,12 @@ C the energy transfer exist
         if (zj.lt.buflipbot) then
 C what fraction I am in
          fracinbuf=1.0d0-
-     &        ((positi-bordlipbot)/lipbufthick)
+     &        ((zj-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)
+        elseif (zj.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
          sslipj=sscalelip(fracinbuf)
          ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
         else
@@ -935,9 +988,9 @@ C lipbufthick is thickenes of lipid buffore
          ssgradlipj=0.0
        endif
       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
-     &  +aa_aq(itypi,itypj)*(2.0d0-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
+     &  +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
@@ -998,6 +1051,7 @@ cd     &          rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 c---------------------------------------------------------------
               rij_shift=1.0D0/rij_shift 
               fac=rij_shift**expon
+              faclip=fac
               e1=fac*fac*aa
               e2=fac*bb
               evdwij=eps1*eps2rt*eps3rt*(e1+e2)
@@ -1032,15 +1086,19 @@ C Calculate the radial part of the gradient
               gg(1)=xj*fac
               gg(2)=yj*fac
               gg(3)=zj*fac
-              gg_lipi(3)=ssgradlipi*evdwij
-              gg_lipj(3)=ssgradlipj*evdwij
+            gg_lipi(3)=eps1*(eps2rt*eps2rt)
+     &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+     & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+     &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+            gg_lipj(3)=ssgradlipj*gg_lipi(3)
+            gg_lipi(3)=gg_lipi(3)*ssgradlipi
 C Calculate angular part of the gradient.
               call sc_grad_scale(sss)
             endif
           enddo      ! j
         enddo        ! iint
       enddo          ! i
-c      write (iout,*) "Number of loop steps in EGB:",ind
+C      write (iout,*) "Number of loop steps in EGB:",ind
 cccc      energy_dec=.false.
       return
       end
@@ -1317,6 +1375,8 @@ c      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
       enddo
       do k=1,3
         gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
+        gg_lipi(k)=gg_lipi(k)*scalfac
+        gg_lipj(k)=gg_lipj(k)*scalfac
       enddo 
 c      write (iout,*) "gg",(gg(k),k=1,3)
       do k=1,3
@@ -1589,6 +1649,7 @@ C 13-go grudnia roku pamietnego...
       double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
      &                   0.0d0,1.0d0,0.0d0,
      &                   0.0d0,0.0d0,1.0d0/
+       integer xhift,yshift,zshift
 c          time00=MPI_Wtime()
 cd      write (iout,*) "eelecij",i,j
 C      print *,"WCHODZE2"
@@ -1652,8 +1713,8 @@ C      print *,"WCHODZE2"
           rij=dsqrt(rij)
           rmij=1.0D0/rij
 c For extracting the short-range part of Evdwpp
-          sss=sscale(rij/rpp(iteli,itelj))
-          sssgrad=sscagrad(rij/rpp(iteli,itelj))
+          sss=sscale(rij)
+          sssgrad=sscagrad(rij)
           r3ij=rrmij*rmij
           r6ij=r3ij*r3ij  
           cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
@@ -1778,9 +1839,9 @@ cgrad            do l=1,3
 cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
 cgrad            enddo
 cgrad          enddo
-          ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
-          ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
-          ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
+          ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj
+          ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj
+          ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj
 c          do k=1,3
 c            ghalf=0.5D0*ggg(k)
 c            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
@@ -1835,9 +1896,9 @@ c 9/28/08 AL Gradient compotents will be summed only at the end
 C          ggg(1)=facvdw*xj
 C          ggg(2)=facvdw*yj
 C          ggg(3)=facvdw*zj
-          ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
-          ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
-          ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
+          ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj
+          ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj
+          ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj
           do k=1,3
             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
@@ -2499,8 +2560,8 @@ c     &   ' ielend',ielend_vdw(i)
           rij=xj*xj+yj*yj+zj*zj
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
-          sss=sscale(rij/rpp(iteli,itelj))
-            sssgrad=sscagrad(rij/rpp(iteli,itelj))
+          sss=sscale(rij)
+            sssgrad=sscagrad(rij)
           if (sss.gt.0.0d0) then
             rmij=1.0D0/rij
             r3ij=rrmij*rmij
@@ -2518,9 +2579,9 @@ C
 C Calculate contributions to the Cartesian gradient.
 C
             facvdw=-6*rrmij*(ev1+evdwij)*sss
-          ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
-          ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
-          ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
+          ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+          ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+          ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
 C            ggg(1)=facvdw*xj
 C            ggg(2)=facvdw*yj
 C            ggg(3)=facvdw*zj
@@ -2552,6 +2613,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
       dimension ggg(3)
+      integer xshift,yshift,zshift
       evdw2=0.0D0
       evdw2_14=0.0d0
 CD        print '(a)','Enter ESCP KURWA'
@@ -2581,6 +2643,13 @@ C Uncomment following three lines for Ca-p interactions
           xj=c(1,j)
           yj=c(2,j)
           zj=c(3,j)
+         xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zi,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
       yj_safe=yj
@@ -2615,8 +2684,8 @@ C Uncomment following three lines for Ca-p interactions
 
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 
-          sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
-          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
+          sss=sscale(1.0d0/(dsqrt(rrij)))
+          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
           if (sss.lt.1.0d0) then
             fac=rrij**expon2
             e1=fac*fac*aad(itypj,iteli)
@@ -2635,7 +2704,7 @@ C Calculate contributions to the gradient in the virtual-bond and SC vectors.
 C
              
             fac=-(evdwij+e1)*rrij*(1.0d0-sss)
-            fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
+            fac=fac-(evdwij)*sssgrad*dsqrt(rrij)
             ggg(1)=xj*fac
             ggg(2)=yj*fac
             ggg(3)=zj*fac
@@ -2691,6 +2760,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
       dimension ggg(3)
+      integer xshift,yshift,zshift
       evdw2=0.0D0
       evdw2_14=0.0d0
 cd    print '(a)','Enter ESCP'
@@ -2721,6 +2791,12 @@ C Uncomment following three lines for Ca-p interactions
           xj=c(1,j)
           yj=c(2,j)
           zj=c(3,j)
+         xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zi,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
       yj_safe=yj
@@ -2753,8 +2829,8 @@ C Uncomment following three lines for Ca-p interactions
           zj=zj_safe-zi
        endif
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
-          sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
-          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
+          sss=sscale(1.0d0/(dsqrt(rrij)))
+          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
           if (sss.gt.0.0d0) then
 
             fac=rrij**expon2
@@ -2773,7 +2849,7 @@ C
 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
 C
             fac=-(evdwij+e1)*rrij*sss
-            fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
+            fac=fac+(evdwij)*sssgrad*dsqrt(rrij)
             ggg(1)=xj*fac
             ggg(2)=yj*fac
             ggg(3)=zj*fac

diff --git a/source/unres/src_MD-M/energy_split-sep.F b/source/unres/src_MD-M/energy_split-sep.F
index e15ab61..b3f924b 100644 (file)
--- a/source/unres/src_MD-M/energy_split-sep.F
+++ b/source/unres/src_MD-M/energy_split-sep.F
@@ -246,7 +246,12 @@ c     &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6
       if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then
          call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
       endif
-C 
+C
+C      if (wliptran.gt.0) then
+C        call Eliptransfer(eliptran)
+C      else
+C       eliptran=0.0d0
+C      endif 
 C If performing constraint dynamics, call the constraint energy
 C  after the equilibration time
       if(usampl.and.totT.gt.eq_time) then
@@ -287,7 +292,7 @@ C
       energia(20)=Uconst+Uconst_back
       call sum_energy(energia,.true.)
 C      call enerprint
-      write (iout,*) "Exit ETOTAL_LONG"
+C      write (iout,*) "Exit ETOTAL_LONG"
       call flush(iout)
       return
       end
@@ -575,6 +580,12 @@ C
       else
         esccor=0.0d0
       endif
+      if (wliptran.gt.0) then
+        call Eliptransfer(eliptran)
+      else
+       eliptran=0.0d0
+      endif
+C       print *,eliptran,wliptran
 C
 C Put energy components into an array
 C
@@ -611,10 +622,11 @@ C
       energia(17)=estr
       energia(19)=edihcnstr
       energia(21)=esccor
-      write (iout,*) "ETOTAL_SHORT before SUM_ENERGY"
+      energia(22)=eliptran
+C      write (iout,*) "ETOTAL_SHORT before SUM_ENERGY"
       call flush(iout)
       call sum_energy(energia,.true.)
-      write (iout,*) "Exit ETOTAL_SHORT"
+C      write (iout,*) "Exit ETOTAL_SHORT"
       call flush(iout)
       return
       end

diff --git a/source/unres/src_MD-M/readpdb.F b/source/unres/src_MD-M/readpdb.F
index 7aa8fd4..48bb75a 100644 (file)
--- a/source/unres/src_MD-M/readpdb.F
+++ b/source/unres/src_MD-M/readpdb.F
@@ -310,7 +310,9 @@ c        enddo
 c enddiagnostic       
 C makes copy of chains
         write (iout,*) "symetr", symetr
-       
+      do j=1,3
+      dc(j,0)=c(j,1)
+      enddo 
       if (symetr.gt.1) then
        call permut(symetr)
        nperm=1