X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Fcluster%2Fwham%2Fsrc-HCD%2Fenergy_p_new.F;h=119bad6fe1c7f0036202fc7da276f0a8a070c8fa;hb=57038e4bdff4cc9534106b25bfbd4b9a844d47fd;hp=5d07d5db79749520c5530904889c4b22b6c7dd5c;hpb=32caa3b64eb94b90fa9fd402b77263ea89efffa1;p=unres.git

diff --git a/source/cluster/wham/src-HCD/energy_p_new.F b/source/cluster/wham/src-HCD/energy_p_new.F
index 5d07d5d..119bad6 100644
--- a/source/cluster/wham/src-HCD/energy_p_new.F
+++ b/source/cluster/wham/src-HCD/energy_p_new.F
@@ -180,14 +180,23 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
 c      write(iout,*) "TEST_ENE1 ehomology_constr=",ehomology_constr
 #ifdef DFA
 C     BARTEK for dfa test!
+      edfadis=0.0d0
       if (wdfa_dist.gt.0) call edfad(edfadis)
 c      write(iout,*)'edfad is finished!', wdfa_dist,edfadis
+      edfator=0.0d0
       if (wdfa_tor.gt.0) call edfat(edfator)
 c      write(iout,*)'edfat is finished!', wdfa_tor,edfator
+      edfanei=0.0d0
       if (wdfa_nei.gt.0) call edfan(edfanei)
 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
 
 #ifdef SPLITELE
@@ -511,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),
@@ -674,6 +686,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
@@ -688,6 +701,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
@@ -857,6 +874,7 @@ c      enddo
         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
@@ -867,6 +885,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
@@ -972,6 +994,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)
@@ -1004,9 +1027,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)
@@ -1110,35 +1137,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)
@@ -1196,80 +1196,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)
@@ -1391,6 +1326,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)
@@ -1425,9 +1362,15 @@ 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
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -2256,12 +2199,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)
@@ -2291,38 +2229,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
@@ -2362,43 +2269,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
@@ -2503,75 +2374,11 @@ 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))
           if (sss.eq.0.0d0) return
           sssgrad=sscagrad(sqrt(rij))
@@ -4055,12 +3862,7 @@ c     &   " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
         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)
@@ -4075,44 +3877,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
@@ -10217,7 +9985,6 @@ c         min_odl=minval(distancek)
      &              min_odl=distancek(kk)
            enddo
          endif
-
 c        write (iout,* )"min_odl",min_odl
 #ifdef DEBUG
          write (iout,*) "ij dij",i,j,dij