correction in tube for parallel code
[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F
index 6078e06..38efafa 100644 (file)
@@ -671,7 +671,26 @@ c      enddo
 
 
         enddo
-      enddo 
+      enddo
+      j=1
+      i=0
+      print *,"KUPA2",gradbufc(j,i),wsc*gvdwc(j,i),
+     &                wscp*gvdwc_scp(j,i),gvdwc_scpp(j,i),
+     &                welec*gelc_long(j,i),wvdwpp*gvdwpp(j,i),
+     &                wel_loc*gel_loc_long(j,i),
+     &                wcorr*gradcorr_long(j,i),
+     &                wcorr5*gradcorr5_long(j,i),
+     &                wcorr6*gradcorr6_long(j,i),
+     &                wturn6*gcorr6_turn_long(j,i),
+     &                wstrain*ghpbc(j,i)
+     &                ,wliptran*gliptranc(j,i)
+     &                ,gradafm(j,i)
+     &                 ,welec*gshieldc(j,i)
+     &                 ,wcorr*gshieldc_ec(j,i)
+     &                 ,wturn3*gshieldc_t3(j,i)
+     &                 ,wturn4*gshieldc_t4(j,i)
+     &                 ,wel_loc*gshieldc_ll(j,i)
+     &                ,wtube*gg_tube(j,i) 
 #else
       do i=0,nct
         do j=1,3
@@ -726,7 +745,7 @@ c      call flush(iout)
 #ifdef TIMING
 c      time_allreduce=time_allreduce+MPI_Wtime()-time00
 #endif
-      do i=nnt,nres
+      do i=0,nres
         do k=1,3
           gradbufc(k,i)=0.0d0
         enddo
@@ -902,7 +921,42 @@ C          print *,gradafm(1,13),"AFM"
 
 
         enddo
-      enddo 
+      enddo
+C       i=0
+C       j=1
+C       print *,"KUPA",    gradbufc(j,i),welec*gelc(j,i),
+C     &                wel_loc*gel_loc(j,i),
+C     &                0.5d0*wscp*gvdwc_scpp(j,i),
+C     &                welec*gelc_long(j,i),wvdwpp*gvdwpp(j,i),
+C     &                wel_loc*gel_loc_long(j,i),
+C     &                wcorr*gradcorr_long(j,i),
+C     &                wcorr5*gradcorr5_long(j,i),
+C     &                wcorr6*gradcorr6_long(j,i),
+C     &                wturn6*gcorr6_turn_long(j,i),
+C     &                wbond*gradb(j,i),
+C     &                wcorr*gradcorr(j,i),
+C     &                wturn3*gcorr3_turn(j,i),
+C     &                wturn4*gcorr4_turn(j,i),
+C     &                wcorr5*gradcorr5(j,i),
+C     &                wcorr6*gradcorr6(j,i),
+C     &                wturn6*gcorr6_turn(j,i),
+C     &                wsccor*gsccorc(j,i)
+C     &               ,wscloc*gscloc(j,i)
+C     &               ,wliptran*gliptranc(j,i)
+C     &                ,gradafm(j,i)
+C     &                 +welec*gshieldc(j,i)
+C     &                 +welec*gshieldc_loc(j,i)
+C     &                 +wcorr*gshieldc_ec(j,i)
+C     &                 +wcorr*gshieldc_loc_ec(j,i)
+C     &                 +wturn3*gshieldc_t3(j,i)
+C     &                 +wturn3*gshieldc_loc_t3(j,i)
+C     &                 +wturn4*gshieldc_t4(j,i)
+C     &                 ,wturn4*gshieldc_loc_t4(j,i)
+C     &                 ,wel_loc*gshieldc_ll(j,i)
+C     &                 ,wel_loc*gshieldc_loc_ll(j,i)
+C     &                ,wtube*gg_tube(j,i)
+
+C      print *,gg_tube(1,0),"TU3" 
 #ifdef DEBUG
       write (iout,*) "gloc before adding corr"
       do i=1,4*nres
@@ -954,7 +1008,7 @@ c#undef DEBUG
         call MPI_Barrier(FG_COMM,IERR)
         time_barrier_g=time_barrier_g+MPI_Wtime()-time00
         time00=MPI_Wtime()
-        call MPI_Reduce(gradbufc(1,1),gradc(1,1,icg),3*nres,
+        call MPI_Reduce(gradbufc(1,0),gradc(1,0,icg),3*nres+3,
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         call MPI_Reduce(gradbufx(1,1),gradx(1,1,icg),3*nres,
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
@@ -3636,6 +3690,7 @@ 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
@@ -12095,20 +12150,38 @@ C simple Kihara potential
       include 'COMMON.SBRIDGE'
       double precision tub_r,vectube(3),enetube(maxres*2)
       Etube=0.0d0
-      do i=1,2*nres
+      do i=itube_start,itube_end
         enetube(i)=0.0d0
+        enetube(i+nres)=0.0d0
       enddo
 C first we calculate the distance from tube center
 C first sugare-phosphate group for NARES this would be peptide group 
 C for UNRES
-      do i=1,nres
+       do i=itube_start,itube_end
 C lets ommit dummy atoms for now
        if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
 C now calculate distance from center of tube and direction vectors
-      vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
-          if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
-      vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
-          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+      xmin=boxxsize
+      ymin=boxysize
+        do j=-1,1
+         vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+         vectube(1)=vectube(1)+boxxsize*j
+         vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+         vectube(2)=vectube(2)+boxysize*j
+       
+         xminact=abs(vectube(1)-tubecenter(1))
+         yminact=abs(vectube(2)-tubecenter(2))
+           if (xmin.gt.xminact) then
+            xmin=xminact
+            xtemp=vectube(1)
+           endif
+           if (ymin.gt.yminact) then
+             ymin=yminact
+             ytemp=vectube(2)
+            endif
+         enddo
+      vectube(1)=xtemp
+      vectube(2)=ytemp
       vectube(1)=vectube(1)-tubecenter(1)
       vectube(2)=vectube(2)-tubecenter(2)
 
@@ -12128,12 +12201,12 @@ C calculte rdiffrence between r and r0
 C and its 6 power
       rdiff6=rdiff**6.0d0
 C for vectorization reasons we will sumup at the end to avoid depenence of previous
-       enetube(i)=pep_aa_tube/rdiff6**2.0d0-pep_bb_tube/rdiff6
+       enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6
 C       write(iout,*) "TU13",i,rdiff6,enetube(i)
 C       print *,rdiff,rdiff6,pep_aa_tube
 C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
 C now we calculate gradient
-       fac=(-12.0d0*pep_aa_tube/rdiff6+
+       fac=(-12.0d0*pep_aa_tube/rdiff6-
      &       6.0d0*pep_bb_tube)/rdiff6/rdiff
 C       write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
 C     &rdiff,fac
@@ -12145,7 +12218,10 @@ C now direction of gg_tube vector
         enddo
         enddo
 C basically thats all code now we split for side-chains (REMEMBER to sum up at the END)
-        do i=1,nres
+C        print *,gg_tube(1,0),"TU"
+
+
+       do i=itube_start,itube_end
 C Lets not jump over memory as we use many times iti
          iti=itype(i)
 C lets ommit dummy atoms for now
@@ -12153,13 +12229,29 @@ C lets ommit dummy atoms for now
 C in UNRES uncomment the line below as GLY has no side-chain...
 C      .or.(iti.eq.10)
      &   ) cycle
-          vectube(1)=c(1,i+nres)
-          vectube(1)=mod(vectube(1),boxxsize)
-          if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
-          vectube(2)=c(2,i+nres)
-          vectube(2)=mod(vectube(2),boxysize)
-          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
-
+      xmin=boxxsize
+      ymin=boxysize
+        do j=-1,1
+         vectube(1)=mod((c(1,i+nres)),boxxsize)
+         vectube(1)=vectube(1)+boxxsize*j
+         vectube(2)=mod((c(2,i+nres)),boxysize)
+         vectube(2)=vectube(2)+boxysize*j
+
+         xminact=abs(vectube(1)-tubecenter(1))
+         yminact=abs(vectube(2)-tubecenter(2))
+           if (xmin.gt.xminact) then
+            xmin=xminact
+            xtemp=vectube(1)
+           endif
+           if (ymin.gt.yminact) then
+             ymin=yminact
+             ytemp=vectube(2)
+            endif
+         enddo
+      vectube(1)=xtemp
+      vectube(2)=ytemp
+C          write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2),
+C     &     tubecenter(2)
       vectube(1)=vectube(1)-tubecenter(1)
       vectube(2)=vectube(2)-tubecenter(2)
 
@@ -12171,6 +12263,7 @@ C now calculte the distance
 C now normalize vector
       vectube(1)=vectube(1)/tub_r
       vectube(2)=vectube(2)/tub_r
+
 C calculte rdiffrence between r and r0
       rdiff=tub_r-tubeR0
 C and its 6 power
@@ -12178,10 +12271,10 @@ C and its 6 power
 C for vectorization reasons we will sumup at the end to avoid depenence of previous
        sc_aa_tube=sc_aa_tube_par(iti)
        sc_bb_tube=sc_bb_tube_par(iti)
-       enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0-sc_bb_tube/rdiff6
+       enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6
 C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
 C now we calculate gradient
-       fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff+
+       fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff-
      &       6.0d0*sc_bb_tube/rdiff6/rdiff
 C now direction of gg_tube vector
          do j=1,3
@@ -12189,8 +12282,8 @@ C now direction of gg_tube vector
           gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
          enddo
         enddo
-        do i=1,2*nres
-          Etube=Etube+enetube(i)
+        do i=itube_start,itube_end
+          Etube=Etube+enetube(i)+enetube(i+nres)
         enddo
 C        print *,"ETUBE", etube
         return
@@ -12230,13 +12323,14 @@ C simple Kihara potential
       include 'COMMON.SBRIDGE'
       double precision tub_r,vectube(3),enetube(maxres*2)
       Etube=0.0d0
-      do i=1,2*nres
+      do i=itube_start,itube_end
         enetube(i)=0.0d0
+        enetube(i+nres)=0.0d0
       enddo
 C first we calculate the distance from tube center
 C first sugare-phosphate group for NARES this would be peptide group 
 C for UNRES
-      do i=1,nres
+       do i=itube_start,itube_end
 C lets ommit dummy atoms for now
        
        if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle
@@ -12264,12 +12358,12 @@ C calculte rdiffrence between r and r0
 C and its 6 power
       rdiff6=rdiff**6.0d0
 C for vectorization reasons we will sumup at the end to avoid depenence of previous
-       enetube(i)=pep_aa_tube/rdiff6**2.0d0-pep_bb_tube/rdiff6
+       enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6
 C       write(iout,*) "TU13",i,rdiff6,enetube(i)
 C       print *,rdiff,rdiff6,pep_aa_tube
 C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
 C now we calculate gradient
-       fac=(-12.0d0*pep_aa_tube/rdiff6+
+       fac=(-12.0d0*pep_aa_tube/rdiff6-
      &       6.0d0*pep_bb_tube)/rdiff6/rdiff
 C       write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
 C     &rdiff,fac
@@ -12281,7 +12375,8 @@ C now direction of gg_tube vector
         enddo
         enddo
 C basically thats all code now we split for side-chains (REMEMBER to sum up at the END)
-        do i=1,nres
+C        print *,gg_tube(1,0),"TU"
+        do i=itube_start,itube_end
 C Lets not jump over memory as we use many times iti
          iti=itype(i)
 C lets ommit dummy atoms for now
@@ -12359,11 +12454,11 @@ C and its 6 power
 C for vectorization reasons we will sumup at the end to avoid depenence of previous
        sc_aa_tube=sc_aa_tube_par(iti)
        sc_bb_tube=sc_bb_tube_par(iti)
-       enetube(i+nres)=(sc_aa_tube/rdiff6**2.0d0-sc_bb_tube/rdiff6)
+       enetube(i+nres)=(sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6)
      &                 *sstube+enetube(i+nres)
 C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
 C now we calculate gradient
-       fac=(-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff+
+       fac=(-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff-
      &       6.0d0*sc_bb_tube/rdiff6/rdiff)*sstube
 C now direction of gg_tube vector
          do j=1,3
@@ -12376,8 +12471,8 @@ C now direction of gg_tube vector
      &+ssgradtube*enetube(i+nres)/sstube
 
         enddo
-        do i=1,2*nres
-          Etube=Etube+enetube(i)
+        do i=itube_start,itube_end
+          Etube=Etube+enetube(i)+enetube(i+nres)
         enddo
 C        print *,"ETUBE", etube
         return