X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-HCD%2Fenergy_p_new.F;fp=source%2Fwham%2Fsrc-HCD%2Fenergy_p_new.F;h=ce7a6a7ad1cdb402f08724cc05283149fecbd3c7;hb=57038e4bdff4cc9534106b25bfbd4b9a844d47fd;hp=610515637abdc807973f7f28c7ca524cc9d960dc;hpb=32caa3b64eb94b90fa9fd402b77263ea89efffa1;p=unres.git

diff --git a/source/wham/src-HCD/energy_p_new.F b/source/wham/src-HCD/energy_p_new.F
index 6105156..ce7a6a7 100644
--- a/source/wham/src-HCD/energy_p_new.F
+++ b/source/wham/src-HCD/energy_p_new.F
@@ -159,6 +159,7 @@ c         write (iout,*) "Calling multibody_hbond"
          call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
       endif
 #endif
+c      write (iout,*) "nsaxs",nsaxs
 c      write (iout,*) "From Esaxs: Esaxs_constr",Esaxs_constr
       if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
         call e_saxs(Esaxs_constr)
@@ -192,8 +193,12 @@ c      write(iout,*)'edfan is finished!', wdfa_nei,edfanei
       edfabet=0.0d0
       if (wdfa_beta.gt.0) call edfab(edfabet)
 c      write(iout,*)'edfab is finished!', wdfa_beta,edfabet
+#else
+      edfadis=0.0d0
+      edfator=0.0d0
+      edfanei=0.0d0
+      edfabet=0.0d0
 #endif
-
 c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
 #ifdef SPLITELE
       if (shield_mode.gt.0) then
@@ -515,6 +520,9 @@ C     Bartek
       edfator = energia(29)
       edfanei = energia(30)
       edfabet = energia(31)
+      Eafmforc=0.0d0
+      etube=0.0d0
+      Uconst=0.0d0
 #ifdef SPLITELE
       write(iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,wvdwpp,
      &  estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*fact(1),
@@ -680,6 +688,7 @@ cROZNICA
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C Change 12/1/95
         num_conti=0
 C
@@ -694,6 +703,10 @@ cd   &                  'iend=',iend(i,iint)
             xj=c(1,nres+j)-xi
             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)
 C Change 12/1/95 to calculate four-body interactions
             rij=xj*xj+yj*yj+zj*zj
             rrij=1.0D0/rij
@@ -865,6 +878,7 @@ c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C
 C Calculate SC interaction energy.
 C
@@ -875,6 +889,10 @@ C
             xj=c(1,nres+j)-xi
             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
@@ -982,6 +1000,7 @@ c     endif
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -1014,9 +1033,13 @@ c           chip12=0.0D0
 c           alf1=0.0D0
 c           alf2=0.0D0
 c           alf12=0.0D0
-            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)
@@ -1128,35 +1151,8 @@ c      if (icall.gt.0) lprn=.true.
         yi=c(2,nres+i)
         zi=c(3,nres+i)
 C returning the ith atom to box
-          xi=mod(xi,boxxsize)
-          if (xi.lt.0) xi=xi+boxxsize
-          yi=mod(yi,boxysize)
-          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
-
+        call to_box(xi,yi,zi)
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -1214,80 +1210,15 @@ c           alf12=0.0D0
             yj=c(2,nres+j)
             zj=c(3,nres+j)
 C returning jth atom to box
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-       if ((zj.gt.bordlipbot)
-     &.and.(zj.lt.bordliptop)) then
-C the energy transfer exist
-        if (zj.lt.buflipbot) then
-C what fraction I am in
-         fracinbuf=1.0d0-
-     &        ((zj-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
-         sslipj=sscalelip(fracinbuf)
-         ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
-        elseif (zj.gt.bufliptop) then
-         fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
-         sslipj=sscalelip(fracinbuf)
-         ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
-        else
-         sslipj=1.0d0
-         ssgradlipj=0.0
-        endif
-       else
-         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
-      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
-     &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
-C       if (aa.ne.aa_aq(itypi,itypj)) then
-       
-C      write(iout,*) "tu,", i,j,aa_aq(itypi,itypj)-aa,
-C     & bb_aq(itypi,itypj)-bb,
-C     & sslipi,sslipj
-C         endif
-
-C        write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj)
-C checking the distance
-      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      subchap=0
-C finding the closest
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            subchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (subchap.eq.1) then
-          xj=xj_temp-xi
-          yj=yj_temp-yi
-          zj=zj_temp-zi
-       else
-          xj=xj_safe-xi
-          yj=yj_safe-yi
-          zj=zj_safe-zi
-       endif
-
+            call to_box(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
+            bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &        +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+            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)
@@ -1413,6 +1344,8 @@ c      if (icall.gt.0) lprn=.true.
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -1447,9 +1380,21 @@ c           chip12=0.0D0
 c           alf1=0.0D0
 c           alf2=0.0D0
 c           alf12=0.0D0
-            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)
+            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
+            bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &        +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+C      if (aa.ne.aa_aq(itypi,itypj)) write(63,'2e10.5')
+C     &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
+C      write(iout,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj
+            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)
@@ -2271,12 +2216,7 @@ c        end if
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+        call to_box(xmedi,ymedi,zmedi)
         num_conti=0
         call eelecij(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
@@ -2306,37 +2246,7 @@ c     &    .or. itype(i-1).eq.ntyp1
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-C Return atom into box, boxxsize is size of box in x dimension
-c  194   continue
-c        if (xmedi.gt.((0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
-c        if (xmedi.lt.((-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
-C Condition for being inside the proper box
-c        if ((xmedi.gt.((0.5d0)*boxxsize)).or.
-c     &       (xmedi.lt.((-0.5d0)*boxxsize))) then
-c        go to 194
-c        endif
-c  195   continue
-c        if (ymedi.gt.((0.5d0)*boxysize)) ymedi=ymedi-boxysize
-c        if (ymedi.lt.((-0.5d0)*boxysize)) ymedi=ymedi+boxysize
-C Condition for being inside the proper box
-c        if ((ymedi.gt.((0.5d0)*boxysize)).or.
-c     &       (ymedi.lt.((-0.5d0)*boxysize))) then
-c        go to 195
-c        endif
-c  196   continue
-c        if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
-c        if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
-C Condition for being inside the proper box
-c        if ((zmedi.gt.((0.5d0)*boxzsize)).or.
-c     &       (zmedi.lt.((-0.5d0)*boxzsize))) then
-c        go to 196
-c        endif
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+        call to_box(xmedi,ymedi,zmedi)
 #ifdef FOURBODY
         num_conti=num_cont_hb(i)
 #endif
@@ -2376,43 +2286,7 @@ c     &  .or. itype(i-1).eq.ntyp1
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-          xmedi=mod(xmedi,boxxsize)
-          if (xmedi.lt.0) xmedi=xmedi+boxxsize
-          ymedi=mod(ymedi,boxysize)
-          if (ymedi.lt.0) ymedi=ymedi+boxysize
-          zmedi=mod(zmedi,boxzsize)
-          if (zmedi.lt.0) zmedi=zmedi+boxzsize
-C          xmedi=xmedi+xshift*boxxsize
-C          ymedi=ymedi+yshift*boxysize
-C          zmedi=zmedi+zshift*boxzsize
-
-C Return tom into box, boxxsize is size of box in x dimension
-c  164   continue
-c        if (xmedi.gt.((xshift+0.5d0)*boxxsize)) xmedi=xmedi-boxxsize
-c        if (xmedi.lt.((xshift-0.5d0)*boxxsize)) xmedi=xmedi+boxxsize
-C Condition for being inside the proper box
-c        if ((xmedi.gt.((xshift+0.5d0)*boxxsize)).or.
-c     &       (xmedi.lt.((xshift-0.5d0)*boxxsize))) then
-c        go to 164
-c        endif
-c  165   continue
-c        if (ymedi.gt.((yshift+0.5d0)*boxysize)) ymedi=ymedi-boxysize
-c        if (ymedi.lt.((yshift-0.5d0)*boxysize)) ymedi=ymedi+boxysize
-C Condition for being inside the proper box
-c        if ((ymedi.gt.((yshift+0.5d0)*boxysize)).or.
-c     &       (ymedi.lt.((yshift-0.5d0)*boxysize))) then
-c        go to 165
-c        endif
-c  166   continue
-c        if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
-c        if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
-cC Condition for being inside the proper box
-c        if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or.
-c     &       (zmedi.lt.((zshift-0.5d0)*boxzsize))) then
-c        go to 166
-c        endif
-
-c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
+        call to_box(xmedi,ymedi,zmedi)
 #ifdef FOURBODY
         num_conti=num_cont_hb(i)
 #endif
@@ -2518,73 +2392,10 @@ C          zj=c(3,j)+0.5D0*dzj-zmedi
           xj=c(1,j)+0.5D0*dxj
           yj=c(2,j)+0.5D0*dyj
           zj=c(3,j)+0.5D0*dzj
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,boxzsize)
-          if (zj.lt.0) zj=zj+boxzsize
-          if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
-      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-      xj_safe=xj
-      yj_safe=yj
-      zj_safe=zj
-      isubchap=0
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            isubchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (isubchap.eq.1) then
-          xj=xj_temp-xmedi
-          yj=yj_temp-ymedi
-          zj=zj_temp-zmedi
-       else
-          xj=xj_safe-xmedi
-          yj=yj_safe-ymedi
-          zj=zj_safe-zmedi
-       endif
-C        if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC
-c  174   continue
-c        if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c        if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c        if ((xj.gt.((0.5d0)*boxxsize)).or.
-c     &       (xj.lt.((-0.5d0)*boxxsize))) then
-c        go to 174
-c        endif
-c  175   continue
-c        if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c        if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-C Condition for being inside the proper box
-c        if ((yj.gt.((0.5d0)*boxysize)).or.
-c     &       (yj.lt.((-0.5d0)*boxysize))) then
-c        go to 175
-c        endif
-c  176   continue
-c        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c        if ((zj.gt.((0.5d0)*boxzsize)).or.
-c     &       (zj.lt.((-0.5d0)*boxzsize))) then
-c        go to 176
-c        endif
-C        endif !endPBC condintion
-C        xj=xj-xmedi
-C        yj=yj-ymedi
-C        zj=zj-zmedi
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xmedi,boxxsize)
+          yj=boxshift(yj-ymedi,boxysize)
+          zj=boxshift(zj-zmedi,boxzsize)
           rij=xj*xj+yj*yj+zj*zj
 
           sss=sscale(sqrt(rij))
@@ -4075,13 +3886,7 @@ c     &   " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
         zi=0.5D0*(c(3,i)+c(3,i+1))
-C Returning the ith atom to box
-          xi=mod(xi,boxxsize)
-          if (xi.lt.0) xi=xi+boxxsize
-          yi=mod(yi,boxysize)
-          if (yi.lt.0) yi=yi+boxysize
-          zi=mod(zi,boxzsize)
-          if (zi.lt.0) zi=zi+boxzsize
+        call to_box(xi,yi,zi)
         do iint=1,nscp_gr(i)
 
         do j=iscpstart(i,iint),iscpend(i,iint)
@@ -4096,44 +3901,10 @@ C Uncomment following three lines for Ca-p interactions
           yj=c(2,j)
           zj=c(3,j)
 C returning the jth atom to box
-          xj=mod(xj,boxxsize)
-          if (xj.lt.0) xj=xj+boxxsize
-          yj=mod(yj,boxysize)
-          if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zj,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
-      zj_safe=zj
-      subchap=0
-C Finding the closest jth atom
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
-          xj=xj_safe+xshift*boxxsize
-          yj=yj_safe+yshift*boxysize
-          zj=zj_safe+zshift*boxzsize
-          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          if(dist_temp.lt.dist_init) then
-            dist_init=dist_temp
-            xj_temp=xj
-            yj_temp=yj
-            zj_temp=zj
-            subchap=1
-          endif
-       enddo
-       enddo
-       enddo
-       if (subchap.eq.1) then
-          xj=xj_temp-xi
-          yj=yj_temp-yi
-          zj=zj_temp-zi
-       else
-          xj=xj_safe-xi
-          yj=yj_safe-yi
-          zj=zj_safe-zi
-       endif
+          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)
 C sss is scaling function for smoothing the cutoff gradient otherwise
 C the gradient would not be continuouse