small chanegs in nanotube + working wham for lipid

[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 1390fe3..e3a17dd 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -4647,8 +4647,8 @@ c     &   a33*gmuji2(4)
 #endif
 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 
-          if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
-     &            'eelloc',i,j,eel_loc_ij
+          if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2f7.3)')
+     &            'eelloc',i,j,eel_loc_ij,a22*muij(1),a23*muij(2)
 c           if (eel_loc_ij.ne.0)
 c     &      write (iout,'(a4,2i4,8f9.5)')'chuj',
 c     &     i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
@@ -12611,17 +12611,17 @@ C now calculate distance from center of tube and direction vectors
       zmin=boxzsize
 
         do j=-1,1
-         vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+         vectube(1)=dmod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
          vectube(1)=vectube(1)+boxxsize*j
-         vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+         vectube(2)=dmod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
          vectube(2)=vectube(2)+boxysize*j
-         vectube(3)=mod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
+         vectube(3)=dmod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
          vectube(3)=vectube(3)+boxzsize*j
 
 
-         xminact=abs(vectube(1)-tubecenter(1))
-         yminact=abs(vectube(2)-tubecenter(2))
-         zminact=abs(vectube(3)-tubecenter(3))
+         xminact=dabs(vectube(1)-tubecenter(1))
+         yminact=dabs(vectube(2)-tubecenter(2))
+         zminact=dabs(vectube(3)-tubecenter(3))
 
            if (xmin.gt.xminact) then
             xmin=xminact
@@ -12674,13 +12674,13 @@ C go to 667
          enecavtube(i)=0.0
          faccav=0.0
          else
-         denominator=(1.0+dcavtubpep*rdiff6*rdiff6)
+         denominator=(1.0d0+dcavtubpep*rdiff6*rdiff6)
          enecavtube(i)=
-     &   (bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)+ccavtubpep)
+     &   (bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)+ccavtubpep)
      &   /denominator
          enecavtube(i)=0.0
-         faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/sqrt(rdiff))
-     &   *denominator-(bcavtubpep*rdiff+acavtubpep*sqrt(rdiff)
+         faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/dsqrt(rdiff))
+     &   *denominator-(bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)
      &   +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0)
      &   /denominator**2.0d0
 C         faccav=0.0
@@ -12701,7 +12701,7 @@ C now direction of gg_tube vector
         enddo
 
        do i=itube_start,itube_end
-        enecavtube(i)=0.0 
+        enecavtube(i)=0.0d0
 C Lets not jump over memory as we use many times iti
          iti=itype(i)
 C lets ommit dummy atoms for now
@@ -12713,17 +12713,17 @@ C      .or.(iti.eq.10)
       ymin=boxysize
       zmin=boxzsize
         do j=-1,1
-         vectube(1)=mod((c(1,i+nres)),boxxsize)
+         vectube(1)=dmod((c(1,i+nres)),boxxsize)
          vectube(1)=vectube(1)+boxxsize*j
-         vectube(2)=mod((c(2,i+nres)),boxysize)
+         vectube(2)=dmod((c(2,i+nres)),boxysize)
          vectube(2)=vectube(2)+boxysize*j
-         vectube(3)=mod((c(3,i+nres)),boxzsize)
+         vectube(3)=dmod((c(3,i+nres)),boxzsize)
          vectube(3)=vectube(3)+boxzsize*j
 
 
-         xminact=abs(vectube(1)-tubecenter(1))
-         yminact=abs(vectube(2)-tubecenter(2))
-         zminact=abs(vectube(3)-tubecenter(3))
+         xminact=dabs(vectube(1)-tubecenter(1))
+         yminact=dabs(vectube(2)-tubecenter(2))
+         zminact=dabs(vectube(3)-tubecenter(3))
 
            if (xmin.gt.xminact) then
             xmin=xminact
@@ -12772,16 +12772,16 @@ C now direction of gg_tube vector
 C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12)
          if (acavtub(iti).eq.0.0d0) then
 C go to 667
-         enecavtube(i+nres)=0.0
-         faccav=0.0
+         enecavtube(i+nres)=0.0d0
+         faccav=0.0d0
          else
-         denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6)
+         denominator=(1.0d0+dcavtub(iti)*rdiff6*rdiff6)
          enecavtube(i+nres)=
-     &   (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+     &   (bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)+ccavtub(iti))
      &   /denominator
 C         enecavtube(i)=0.0
-         faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/sqrt(rdiff))
-     &   *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)
+         faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/dsqrt(rdiff))
+     &   *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)
      &   +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)*12.0d0)
      &   /denominator**2.0d0
 C         faccav=0.0