po przyspieszeniu przed checkgradem i sympcheckiem

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 252cff0..c8acad2 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -152,9 +152,9 @@ c      print *,"Processor",myrank," left VEC_AND_DERIV"
             eello_turn4=0.0d0
          endif
       else
-c        write (iout,*) "Soft-spheer ELEC potential"
-        call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
-     &   eello_turn4)
+        write (iout,*) "Soft-spheer ELEC potential"
+c        call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3,
+c     &   eello_turn4)
       endif
 c      print *,"Processor",myrank," computed UELEC"
 C
@@ -1425,9 +1425,9 @@ c     if (icall.eq.0) lprn=.false.
       ind=0
 C the loop over all 27 posible neigbours (for xshift=0,yshift=0,zshift=0
 C we have the original box)
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
+C      do xshift=-1,1
+C      do yshift=-1,1
+C      do zshift=-1,1
       do i=iatsc_s,iatsc_e
         itypi=iabs(itype(i))
         if (itypi.eq.ntyp1) cycle
@@ -1466,9 +1466,9 @@ c        endif
           if (yi.lt.0) yi=yi+boxysize
           zi=mod(zi,boxzsize)
           if (zi.lt.0) zi=zi+boxzsize
-          xi=xi+xshift*boxxsize
-          yi=yi+yshift*boxysize
-          zi=zi+zshift*boxzsize
+C          xi=xi+xshift*boxxsize
+C          yi=yi+yshift*boxysize
+C          zi=zi+zshift*boxzsize
 
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
@@ -1544,12 +1544,43 @@ c        endif
           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
+      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
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
-            xj=xj-xi
-            yj=yj-yi
-            zj=zj-zi
+C            xj=xj-xi
+C            yj=yj-yi
+C            zj=zj-zi
 c            write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
 c            write (iout,*) "j",j," dc_norm",
 c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
@@ -1623,9 +1654,9 @@ C Calculate angular part of the gradient.
           enddo      ! j
         enddo        ! iint
       enddo          ! i
-      enddo          ! zshift
-      enddo          ! yshift
-      enddo          ! xshift
+C      enddo          ! zshift
+C      enddo          ! yshift
+C      enddo          ! xshift
 c      write (iout,*) "Number of loop steps in EGB:",ind
 cccc      energy_dec=.false.
       return
@@ -2989,9 +3020,9 @@ c        endif
         num_cont_hb(i)=num_conti
       enddo   ! i
 C Loop over all neighbouring boxes
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
+C      do xshift=-1,1
+C      do yshift=-1,1
+C      do zshift=-1,1
 c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
@@ -3015,9 +3046,9 @@ c
           if (ymedi.lt.0) ymedi=ymedi+boxysize
           zmedi=mod(zmedi,boxzsize)
           if (zmedi.lt.0) zmedi=zmedi+boxzsize
-          xmedi=xmedi+xshift*boxxsize
-          ymedi=ymedi+yshift*boxysize
-          zmedi=zmedi+zshift*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
@@ -3057,9 +3088,9 @@ c          write (iout,*) i,j,itype(i),itype(j)
         enddo ! j
         num_cont_hb(i)=num_conti
       enddo   ! i
-      enddo   ! zshift
-      enddo   ! yshift
-      enddo   ! xshift
+C     enddo   ! zshift
+C      enddo   ! yshift
+C      enddo   ! xshift
 
 c      write (iout,*) "Number of loop steps in EELEC:",ind
 cd      do i=1,nres
@@ -3138,7 +3169,37 @@ C          zj=c(3,j)+0.5D0*dzj-zmedi
           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
+      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-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
+            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
@@ -3165,9 +3226,9 @@ c     &       (zj.lt.((-0.5d0)*boxzsize))) then
 c        go to 176
 c        endif
 C        endif !endPBC condintion
-        xj=xj-xmedi
-        yj=yj-ymedi
-        zj=zj-zmedi
+C        xj=xj-xmedi
+C        yj=yj-ymedi
+C        zj=zj-zmedi
           rij=xj*xj+yj*yj+zj*zj
 
             sss=sscale(sqrt(rij))
@@ -4110,9 +4171,9 @@ C
       r0_scp=4.5d0
 cd    print '(a)','Enter ESCP'
 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
+C      do xshift=-1,1
+C      do yshift=-1,1
+C      do zshift=-1,1
       do i=iatscp_s,iatscp_e
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         iteli=itel(i)
@@ -4150,9 +4211,9 @@ c        endif
           if (yi.lt.0) yi=yi+boxysize
           zi=mod(zi,boxzsize)
           if (zi.lt.0) zi=zi+boxzsize
-          xi=xi+xshift*boxxsize
-          yi=yi+yshift*boxysize
-          zi=zi+zshift*boxzsize
+C          xi=xi+xshift*boxxsize
+C          yi=yi+yshift*boxysize
+C          zi=zi+zshift*boxzsize
         do iint=1,nscp_gr(i)
 
         do j=iscpstart(i,iint),iscpend(i,iint)
@@ -4195,10 +4256,41 @@ c        go to 176
           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
+      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
 c c       endif
-          xj=xj-xi
-          yj=yj-yi
-          zj=zj-zi
+C          xj=xj-xi
+C          yj=yj-yi
+C          zj=zj-zi
           rij=xj*xj+yj*yj+zj*zj
 
           r0ij=r0_scp
@@ -4251,9 +4343,9 @@ cgrad          enddo
 
         enddo ! iint
       enddo ! i
-      enddo !zshift
-      enddo !yshift
-      enddo !xshift
+C      enddo !zshift
+C      enddo !yshift
+C      enddo !xshift
       return
       end
 C-----------------------------------------------------------------------------
@@ -4278,11 +4370,12 @@ C
       dimension ggg(3)
       evdw2=0.0D0
       evdw2_14=0.0d0
+c        print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
 cd    print '(a)','Enter ESCP'
 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
-      do xshift=-1,1
-      do yshift=-1,1
-      do zshift=-1,1
+C      do xshift=-1,1
+C      do yshift=-1,1
+C      do zshift=-1,1
       do i=iatscp_s,iatscp_e
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         iteli=itel(i)
@@ -4295,9 +4388,10 @@ cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
           if (yi.lt.0) yi=yi+boxysize
           zi=mod(zi,boxzsize)
           if (zi.lt.0) zi=zi+boxzsize
-          xi=xi+xshift*boxxsize
-          yi=yi+yshift*boxysize
-          zi=zi+zshift*boxzsize
+c          xi=xi+xshift*boxxsize
+c          yi=yi+yshift*boxysize
+c          zi=zi+zshift*boxzsize
+c        print *,xi,yi,zi,'polozenie i'
 C Return atom into box, boxxsize is size of box in x dimension
 c  134   continue
 c        if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
@@ -4368,13 +4462,46 @@ c        if ((zj.gt.((0.5d0)*boxzsize)).or.
 c     &       (zj.lt.((-0.5d0)*boxzsize))) then
 c        go to 176
 c        endif
-          xj=xj-xi
-          yj=yj-yi
-          zj=zj-zi
+CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      subchap=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-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
+c          print *,xj,yj,zj,'polozenie j'
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+c          print *,rrij
           sss=sscale(1.0d0/(dsqrt(rrij)))
+c          print *,r_cut,1.0d0/dsqrt(rrij),sss,'tu patrz'
+c          if (sss.eq.0) print *,'czasem jest OK'
+          if (sss.le.0.0d0) cycle
           sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
-          if (sss.gt.0.0d0) then
           fac=rrij**expon2
           e1=fac*fac*aad(itypj,iteli)
           e2=fac*bad(itypj,iteli)
@@ -4427,14 +4554,14 @@ cgrad          enddo
             gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
           enddo
-        endif !endif for sscale cutoff
+c        endif !endif for sscale cutoff
         enddo ! j
 
         enddo ! iint
       enddo ! i
-      enddo !zshift
-      enddo !yshift
-      enddo !xshift
+c      enddo !zshift
+c      enddo !yshift
+c      enddo !xshift
       do i=1,nct
         do j=1,3
           gvdwc_scp(j,i)=expon*gvdwc_scp(j,i)