correction in new ions and clusterfix
[unres4.git] / source / unres / energy.F90
index 376d318..63a76fb 100644 (file)
            zj=c(3,nres+j)
               call to_box(xj,yj,zj)
               call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
-              write (iout,*) "KWA2", itypi,itypj
+!              write (iout,*) "KWA2", itypi,itypj
               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 &
         iabs(itype(jjj,1)).eq.1) then
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
-!d          write (iout,*) "eij",eij
+!          write (iout,*) "eij",eij,iii,jjj
          endif
         else if (ii.gt.nres .and. jj.gt.nres) then
 !c Restraints from contact prediction
       itypj=iabs(itype(j,1))
 !      dscj_inv=dsc_inv(itypj)
       dscj_inv=vbld_inv(nres+j)
-      xj=c(1,nres+j)-xi
-      yj=c(2,nres+j)-yi
-      zj=c(3,nres+j)-zi
+      xj=c(1,nres+j)
+      yj=c(2,nres+j)
+      zj=c(3,nres+j)
           call to_box(xj,yj,zj)
+      xj=boxshift(xj-xi,boxxsize)
+      yj=boxshift(yj-yi,boxysize)
+      zj=boxshift(zj-zi,boxzsize)
       dxj=dc_norm(1,nres+j)
       dyj=dc_norm(2,nres+j)
       dzj=dc_norm(3,nres+j)
       eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2) &
         +akct*deltad*deltat12 &
         +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi+ebr
-!      write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
-!     &  " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
-!     &  " deltat12",deltat12," eij",eij 
+!      write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth, &
+!       " akct",akct," deltad",deltad," deltat",deltat1,deltat2, &
+!       " deltat12",deltat12," eij",eij 
       ed=2*akcm*deltad+akct*deltat12
       pom1=akct*deltad
       pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
           call to_box(xj,yj,zj)
+      xj=boxshift(xj-xi,boxxsize)
+      yj=boxshift(yj-yi,boxysize)
+      zj=boxshift(zj-zi,boxzsize)
+
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             fac_augm=rrij**expon
             e_augm=augm(itypi,itypj)*fac_augm
@@ -22271,6 +22278,10 @@ chip1=chip(itypi)
       real(kind=8) ::  facd4, adler, Fgb, facd3
       integer troll,jj,istate
       real (kind=8) :: dcosom1(3),dcosom2(3)
+      real(kind=8) ::locbox(3)
+      locbox(1)=boxxsize
+          locbox(2)=boxysize
+      locbox(3)=boxzsize
 
       evdw=0.0D0
       if (nres_molec(5).eq.0) return
@@ -22313,6 +22324,8 @@ chip1=chip(itypi)
          zj=c(3,j)
  
       call to_box(xj,yj,zj)
+!      write(iout,*) "xi,yi,zi,xj,yj,zj", xi,yi,zi,xj,yj,zj
+
 !      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
 !      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
 !       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
@@ -22321,6 +22334,7 @@ chip1=chip(itypi)
       xj=boxshift(xj-xi,boxxsize)
       yj=boxshift(yj-yi,boxysize)
       zj=boxshift(zj-zi,boxzsize)
+!      write(iout,*) "xj,yj,zj", xj,yj,zj,boxxsize
 
 !          dxj = dc_norm( 1, nres+j )
 !          dyj = dc_norm( 2, nres+j )
@@ -22364,11 +22378,13 @@ chip1=chip(itypi)
       ctail(k,1)=c(k,i+nres)
       ctail(k,2)=c(k,j)
        END DO
+      call to_box(ctail(1,1),ctail(2,1),ctail(3,1))
+      call to_box(ctail(1,2),ctail(2,2),ctail(3,2))
 !c! tail distances will be themselves usefull elswhere
 !c1 (in Gcav, for example)
-       Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
-       Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
-       Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
+       do k=1,3
+       Rtail_distance(k) = boxshift(ctail(k,2) - ctail(k,1),locbox(k))
+       enddo 
        Rtail = dsqrt( &
         (Rtail_distance(1)*Rtail_distance(1)) &
       + (Rtail_distance(2)*Rtail_distance(2)) &
@@ -22383,10 +22399,15 @@ chip1=chip(itypi)
 ! see unres publications for very informative images
       chead(k,1) = c(k, i+nres) + d1 * dc_norm(k, i+nres)
       chead(k,2) = c(k, j)
+      enddo
+      call to_box(chead(1,1),chead(2,1),chead(3,1))
+      call to_box(chead(1,2),chead(2,2),chead(3,2))
+
 ! distance 
 !        Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
-!        Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
-      Rhead_distance(k) = chead(k,2) - chead(k,1)
+!         Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+      do k=1,3
+      Rhead_distance(k) = boxshift(chead(k,2) - chead(k,1),locbox(k))
        END DO
 ! pitagoras (root of sum of squares)
        Rhead = dsqrt( &
@@ -22408,6 +22429,7 @@ chip1=chip(itypi)
        dPOLdOM1 = 0.0d0
        dPOLdOM2 = 0.0d0
         Fcav = 0.0d0
+        Fisocav=0.0d0
         dFdR = 0.0d0
         dCAVdOM1  = 0.0d0
         dCAVdOM2  = 0.0d0
@@ -22635,6 +22657,10 @@ chip1=chip(itypi)
          yj=c(2,j)
          zj=c(3,j)
         call to_box(xj,yj,zj)
+      xj=boxshift(xj-xi,boxxsize)
+      yj=boxshift(yj-yi,boxysize)
+      zj=boxshift(zj-zi,boxzsize)
+
         dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 
         dxj = 0.0d0! dc_norm( 1, nres+j )
@@ -22679,11 +22705,16 @@ chip1=chip(itypi)
       ctail(k,1)=(c(k,i)+c(k,i+1))/2.0
       ctail(k,2)=c(k,j)
        END DO
+      call to_box(ctail(1,1),ctail(2,1),ctail(3,1))
+      call to_box(ctail(1,2),ctail(2,2),ctail(3,2))
+!c! tail distances will be themselves usefull elswhere
+!c1 (in Gcav, for example)
+       do k=1,3
+       Rtail_distance(k) = boxshift(ctail(k,2) - ctail(k,1),locbox(k))
+       enddo
+
 !c! tail distances will be themselves usefull elswhere
 !c1 (in Gcav, for example)
-       Rtail_distance(1) = ctail( 1, 2 ) - ctail( 1,1 )
-       Rtail_distance(2) = ctail( 2, 2 ) - ctail( 2,1 )
-       Rtail_distance(3) = ctail( 3, 2 ) - ctail( 3,1 )
        Rtail = dsqrt( &
         (Rtail_distance(1)*Rtail_distance(1)) &
       + (Rtail_distance(2)*Rtail_distance(2)) &
@@ -22700,11 +22731,20 @@ chip1=chip(itypi)
 ! see unres publications for very informative images
       chead(k,1) = (c(k, i)+c(k,i+1))/2.0 + d1 * dc_norm(k, i)
       chead(k,2) = c(k, j)
+       ENDDO
 ! distance 
 !        Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
 !        Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
-      Rhead_distance(k) = chead(k,2) - chead(k,1)
+      call to_box(chead(1,1),chead(2,1),chead(3,1))
+      call to_box(chead(1,2),chead(2,2),chead(3,2))
+
+! distance 
+!        Rsc_distance(k) = dabs(c(k, i+nres) - c(k, j+nres))
+!         Rsc(k) = Rsc_distance(k) * Rsc_distance(k)
+      do k=1,3
+      Rhead_distance(k) = boxshift(chead(k,2) - chead(k,1),locbox(k))
        END DO
+
 ! pitagoras (root of sum of squares)
        Rhead = dsqrt( &
         (Rhead_distance(1)*Rhead_distance(1)) &
@@ -23453,6 +23493,7 @@ chip1=chip(itypi)
              yj=c(2,j)
              zj=c(3,j)
       call to_box(xj,yj,zj)
+!      write(iout,*) "xi,yi,zi,xj,yj,zj", xi,yi,zi,xj,yj,zj
 !      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
 !      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
 !       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
@@ -23461,7 +23502,7 @@ chip1=chip(itypi)
       xj=boxshift(xj-xi,boxxsize)
       yj=boxshift(yj-yi,boxysize)
       zj=boxshift(zj-zi,boxzsize)
-
+!       write(iout,*) 'after shift', xj,yj,zj
              dist_init=xj**2+yj**2+zj**2
 
              itypi=itype(i,2)
@@ -25240,28 +25281,7 @@ chip1=chip(itypi)
       zi=c(3,nres+i)
         call to_box(xi,yi,zi)
         call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
-       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
+!       endif
 !       print *, sslipi,ssgradlipi
       dxi=dc_norm(1,nres+i)
       dyi=dc_norm(2,nres+i)
@@ -25318,7 +25338,7 @@ chip1=chip(itypi)
          zj=c(3,j+nres)
      call to_box(xj,yj,zj)
      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
-      write(iout,*) "KRUWA", i,j
+!      write(iout,*) "KRUWA", i,j
       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 &