X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.f90;h=197ca8cd7d8c5d137438e484a20ef17677024d2a;hb=726ec392dd7f57d3faac20616e7701620d0454f2;hp=1d581367284325d47a81aa33e1010cbe392dcc70;hpb=c43a5a1d5cc1d3f94c53de28dd1fd93bea770790;p=unres4.git

diff --git a/source/unres/energy.f90 b/source/unres/energy.f90
index 1d58136..197ca8c 100644
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -2691,6 +2691,12 @@
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=dmod(xmedi,boxxsize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=dmod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=dmod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
         num_conti=0
         call eelecij(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
@@ -2709,6 +2715,12 @@
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=dmod(xmedi,boxxsize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=dmod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=dmod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
         num_conti=num_cont_hb(i)
         call eelecij(i,i+3,ees,evdw1,eel_loc)
         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) &
@@ -2729,6 +2741,13 @@
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=dmod(xmedi,boxxsize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=dmod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=dmod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
+
 !        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
         do j=ielstart(i),ielend(i)
@@ -2775,6 +2794,9 @@
       real(kind=8),dimension(2,2) :: acipa !el,a_temp
 !el      real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
       real(kind=8),dimension(4) :: muij
+      real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+                    dist_temp, dist_init
+      integer xshift,yshift,zshift
 !el      integer :: num_conti,j1,j2
 !el      real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,&
 !el        dz_normi,xmedi,ymedi,zmedi
@@ -2796,7 +2818,7 @@
                                              0.0d0,0.0d0,1.0d0/),shape(unmat)) 
 !      integer :: maxconts=nres/4
 !el local variables
-      integer :: k,i,j,iteli,itelj,kkk,l,kkll,m
+      integer :: k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap
       real(kind=8) :: ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
       real(kind=8) :: ees,evdw1,eel_loc,aaa,bbb,ael3i
       real(kind=8) :: dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,&
@@ -2827,12 +2849,59 @@
           dx_normj=dc_norm(1,j)
           dy_normj=dc_norm(2,j)
           dz_normj=dc_norm(3,j)
-          xj=c(1,j)+0.5D0*dxj-xmedi
-          yj=c(2,j)+0.5D0*dyj-ymedi
-          zj=c(3,j)+0.5D0*dzj-zmedi
+!          xj=c(1,j)+0.5D0*dxj-xmedi
+!          yj=c(2,j)+0.5D0*dyj-ymedi
+!          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
+      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      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
+!C          print *,i,j
+          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
+
           rij=xj*xj+yj*yj+zj*zj
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
+            sss_ele_cut=sscale_ele(rij)
+            sss_ele_grad=sscagrad_ele(rij)
+!            print *,sss_ele_cut,sss_ele_grad,&
+!            (rij),r_cut_ele,rlamb_ele
+            if (sss_ele_cut.le.0.0) go to 128
+
           rmij=1.0D0/rij
           r3ij=rrmij*rmij
           r6ij=r3ij*r3ij  
@@ -2852,8 +2921,8 @@
           eesij=el1+el2
 ! 12/26/95 - for the evaluation of multi-body H-bonding interactions
           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
-          ees=ees+eesij
-          evdw1=evdw1+evdwij
+          ees=ees+eesij*sss_ele_cut
+          evdw1=evdw1+evdwij*sss_ele_cut
 !d          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
 !d     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
 !d     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
@@ -2870,8 +2939,8 @@
 ! Calculate contributions to the Cartesian gradient.
 !
 #ifdef SPLITELE
-          facvdw=-6*rrmij*(ev1+evdwij)
-          facel=-3*rrmij*(el1+eesij)
+          facvdw=-6*rrmij*(ev1+evdwij)*sss_ele_cut
+          facel=-3*rrmij*(el1+eesij)*sss_ele_cut
           fac1=fac
           erij(1)=xj*rmij
           erij(2)=yj*rmij
@@ -2900,9 +2969,9 @@
 !grad              gelc(l,k)=gelc(l,k)+ggg(l)
 !grad            enddo
 !grad          enddo
-          ggg(1)=facvdw*xj
-          ggg(2)=facvdw*yj
-          ggg(3)=facvdw*zj
+          ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj
+          ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj
+          ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj
 !          do k=1,3
 !            ghalf=0.5D0*ggg(k)
 !            gvdwpp(k,i)=gvdwpp(k,i)+ghalf
@@ -2922,8 +2991,8 @@
 !grad            enddo
 !grad          enddo
 #else
-          facvdw=ev1+evdwij 
-          facel=el1+eesij  
+          facvdw=(ev1+evdwij)*sss_ele_cut
+          facel=(el1+eesij)*sss_ele_cut
           fac1=fac
           fac=-3*rrmij*(facvdw+facvdw+facel)
           erij(1)=xj*rmij
@@ -2932,9 +3001,9 @@
 !
 ! Radial derivatives. First process both termini of the fragment (i,j)
 ! 
-          ggg(1)=fac*xj
-          ggg(2)=fac*yj
-          ggg(3)=fac*zj
+          ggg(1)=fac*xj+sss_ele_grad*rmij*(eesij+evdwij)*xj
+          ggg(2)=fac*yj+sss_ele_grad*rmij*(eesij+evdwij)*yj
+          ggg(3)=fac*zj+sss_ele_grad*rmij*(eesij+evdwij)*zj
 !          do k=1,3
 !            ghalf=0.5D0*ggg(k)
 !            gelc(k,i)=gelc(k,i)+ghalf
@@ -3003,6 +3072,7 @@
             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
           enddo
+ 128      continue
           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 &
               .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 &
               .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
@@ -10999,6 +11069,20 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       endif
       return
       end function sscale
+      real(kind=8) function sscale_grad(r)
+!      include "COMMON.SPLITELE"
+      real(kind=8) :: r,gamm
+      if(r.lt.r_cut-rlamb) then
+        sscale_grad=0.0d0
+      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+        gamm=(r-(r_cut-rlamb))/rlamb
+        sscale_grad=gamm*(6*gamm-6.0d0)/rlamb
+      else
+        sscale_grad=0d0
+      endif
+      return
+      end function sscale_grad
+
 !!!!!!!!!! PBCSCALE
       real(kind=8) function sscale_ele(r)
 !      include "COMMON.SPLITELE"
@@ -11648,7 +11732,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !el local variables
       integer :: iint,itypi,itypi1,itypj,subchap
       real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig,sig0ij,rij_shift
-      real(kind=8) :: sss,e1,e2,evdw
+      real(kind=8) :: sss,e1,e2,evdw,sss_grad
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
                     dist_temp, dist_init
 
@@ -11752,6 +11836,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
             sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
             sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+            sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
             if (sss_ele_cut.le.0.0) cycle
             if (sss.lt.1.0d0) then
 
@@ -11805,8 +11890,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
               fac=-expon*(e1+evdwij)*rij_shift
               sigder=fac*sigder
               fac=rij*fac
-              fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
-            /sigma(itypi,itypj)*rij
+              fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+            /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij  &
+            /sigmaii(itypi,itypj))
 !              fac=0.0d0
 ! Calculate the radial part of the gradient
               gg(1)=xj*fac
@@ -11845,7 +11931,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !el local variables
       integer :: iint,itypi,itypi1,itypj,subchap
       real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig0ij,sig
-      real(kind=8) :: sss,e1,e2,evdw,rij_shift
+      real(kind=8) :: sss,e1,e2,evdw,rij_shift,sss_grad
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
                     dist_temp, dist_init
       evdw=0.0D0
@@ -11955,6 +12041,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
             sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+            sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
             sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
             sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
             if (sss_ele_cut.le.0.0) cycle
@@ -11989,7 +12076,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !              write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 !     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
               evdwij=evdwij*eps2rt*eps3rt
-              evdw=evdw+evdwij*sss
+              evdw=evdw+evdwij*sss*sss_ele_cut
               if (lprn) then
               sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
               epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@ -12011,8 +12098,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
               fac=-expon*(e1+evdwij)*rij_shift
               sigder=fac*sigder
               fac=rij*fac
-              fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
-            /sigma(itypi,itypj)*rij
+              fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+            /sigma(itypi,itypj)*rij+sss_grad/sss*rij  &
+            /sigmaii(itypi,itypj))
 
 !              fac=0.0d0
 ! Calculate the radial part of the gradient