respa cleaning (my bug) and shielding MPI and previous bug cleaning
[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F
index f557ebc..934bfb9 100644 (file)
@@ -10,6 +10,8 @@ cMS$ATTRIBUTES C ::  proc_proc
 #ifdef MPI
       include "mpif.h"
       double precision weights_(n_ene)
+      integer IERR
+      integer status(MPI_STATUS_SIZE)
 #endif
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
@@ -25,6 +27,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.CONTROL'
       include 'COMMON.TIME1'
       include 'COMMON.SPLITELE'
+      include 'COMMON.SHIELD'
 #ifdef MPI      
 c      print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
 c     & " nfgtasks",nfgtasks
@@ -149,6 +152,63 @@ C      write (iout,*) "shield_mode",shield_mode
        call set_shield_fac
       else if  (shield_mode.eq.2) then
        call set_shield_fac2
+      if (nfgtasks.gt.1) then
+C#define DEBUG
+#ifdef DEBUG
+       write(iout,*) "befor reduce fac_shield reduce"
+       do i=1,nres
+        write(2,*) "fac",itype(i),fac_shield(i),grad_shield(1,i)
+        write(2,*) "list", shield_list(1,i),ishield_list(i),
+     &  grad_shield_side(1,1,i),grad_shield_loc(1,1,i)
+       enddo
+#endif
+       call MPI_Allgatherv(fac_shield(ivec_start),ivec_count(fg_rank1),
+     &  MPI_DOUBLE_PRECISION,fac_shield(1),ivec_count(0),ivec_displ(0),
+     &  MPI_DOUBLE_PRECISION,FG_COMM,IERR)
+       call MPI_Allgatherv(shield_list(1,ivec_start),
+     &  ivec_count(fg_rank1),
+     &  MPI_I50,shield_list(1,1),ivec_count(0),
+     &  ivec_displ(0),
+     &  MPI_I50,FG_COMM,IERR)
+       call MPI_Allgatherv(ishield_list(ivec_start),
+     &  ivec_count(fg_rank1),
+     &  MPI_INTEGER,ishield_list(1),ivec_count(0),
+     &  ivec_displ(0),
+     &  MPI_INTEGER,FG_COMM,IERR)
+       call MPI_Allgatherv(grad_shield(1,ivec_start),
+     &  ivec_count(fg_rank1),
+     &  MPI_UYZ,grad_shield(1,1),ivec_count(0),
+     &  ivec_displ(0),
+     &  MPI_UYZ,FG_COMM,IERR)
+       call MPI_Allgatherv(grad_shield_side(1,1,ivec_start),
+     &  ivec_count(fg_rank1),
+     &  MPI_SHI,grad_shield_side(1,1,1),ivec_count(0),
+     &  ivec_displ(0),
+     &  MPI_SHI,FG_COMM,IERR)
+       call MPI_Allgatherv(grad_shield_loc(1,1,ivec_start),
+     &  ivec_count(fg_rank1),
+     &  MPI_SHI,grad_shield_loc(1,1,1),ivec_count(0),
+     &  ivec_displ(0),
+     &  MPI_SHI,FG_COMM,IERR)
+#ifdef DEBUG
+       write(iout,*) "after reduce fac_shield reduce"
+       do i=1,nres
+        write(2,*) "fac",itype(i),fac_shield(i),grad_shield(1,i)
+        write(2,*) "list", shield_list(1,i),ishield_list(i),
+     &  grad_shield_side(1,1,i),grad_shield_loc(1,1,i)
+       enddo
+#endif
+C#undef DEBUG
+      endif
+#ifdef DEBUG
+      do i=1,nres
+      write(iout,*) fac_shield(i),ishield_list(i),i,grad_shield(1,i)
+        do j=1,ishield_list(i)
+         write(iout,*) "grad", grad_shield_side(1,j,i),
+     &   grad_shield_loc(1,j,i)
+        enddo
+      enddo
+#endif
       endif
 c      print *,"Processor",myrank," left VEC_AND_DERIV"
       if (ipot.lt.6) then
@@ -300,6 +360,8 @@ C based on partition function
 C      print *,"przed lipidami"
       if (wliptran.gt.0) then
         call Eliptransfer(eliptran)
+      else
+       eliptran=0.0d0
       endif
 C      print *,"za lipidami"
       if (AFMlog.gt.0) then
@@ -413,11 +475,11 @@ cMS$ATTRIBUTES C ::  proc_proc
         time00=MPI_Wtime()
         call MPI_Reduce(enebuff(0),energia(0),n_ene+1,
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
-#ifdef DEBUG
+C#ifdef DEBUG
         write (iout,*) "energies after REDUCE"
         call enerprint(energia)
         call flush(iout)
-#endif
+C#endif
         time_Reduce=time_Reduce+MPI_Wtime()-time00
       endif
       if (fg_rank.eq.0) then
@@ -4696,11 +4758,13 @@ C        else
 C        fac_shield(i)=0.4
 C        fac_shield(j)=0.6
         endif
+C         if (j.eq.78)
+C     &   write(iout,*) i,j,fac_shield(i),fac_shield(j)
         eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
      &  *fac_shield(i)*fac_shield(j)
         eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
      &  *fac_shield(i)*fac_shield(j)
-C#ifdef NEWCORR
+#ifdef NEWCORR
 C Derivatives in theta
         gloc(nphi+i,icg)=gloc(nphi+i,icg)
      &  +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3
@@ -4708,7 +4772,7 @@ C Derivatives in theta
         gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)
      &  +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3
      &   *fac_shield(i)*fac_shield(j)
-C#endif
+#endif
 
 C        if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
 C Derivatives in shield mode
@@ -11527,19 +11591,29 @@ C first for shielding is setting of function of side-chains
       include 'COMMON.IOUNITS'
       include 'COMMON.SHIELD'
       include 'COMMON.INTERACT'
+      include 'COMMON.LOCAL'
+
 C this is the squar root 77 devided by 81 the epislion in lipid (in protein)
       double precision div77_81/0.974996043d0/,
      &div4_81/0.2222222222d0/,sh_frac_dist_grad(3)
-
+  
 C the vector between center of side_chain and peptide group
        double precision pep_side(3),long,side_calf(3),
      &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
      &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
+C      write(2,*) "ivec",ivec_start,ivec_end
+      do i=1,nres
+        fac_shield(i)=0.0d0
+        do j=1,3
+        grad_shield(j,i)=0.0d0
+        enddo
+      enddo
 C the line belowe needs to be changed for FGPROC>1
-      do i=iatscp_s,iatscp_e
+      do i=ivec_start,ivec_end
 C      do i=1,nres-1
-      if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
+C      if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
       ishield_list(i)=0
+      if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
 Cif there two consequtive dummy atoms there is no peptide group between them
 C the line below has to be changed for FGPROC>1
       VolumeTotal=0.0
@@ -11572,6 +11646,7 @@ C now sscale fraction
 C       print *,buff_shield,"buff"
 C now sscale
         if (sh_frac_dist.le.0.0) cycle
+C        print *,ishield_list(i),i
 C If we reach here it means that this side chain reaches the shielding sphere
 C Lets add him to the list for gradient       
         ishield_list(i)=ishield_list(i)+1
@@ -11686,7 +11761,7 @@ C grad_shield_side is Cbeta sidechain gradient
       VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
       enddo
       fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
-C      write(2,*) "TOTAL VOLUME",i,VolumeTotal,fac_shield(i)
+C      write(2,*) "TOTAL VOLUME",i,itype(i),fac_shield(i)
       enddo
       return
       end
@@ -11861,6 +11936,7 @@ C first sugare-phosphate group for NARES this would be peptide group
 C for UNRES
       do i=1,nres
 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)