MD changes micorcanononical not working
[unres4.git] / source / unres / MREMD.f90
index f33432f..350ab52 100644 (file)
@@ -7,7 +7,7 @@
       use remd_data
       use geometry_data
       use energy_data
-      use control_data, only:maxprocs
+!      use control_data, only:maxprocs
       use MDyn
 
       implicit none
       integer(kind=2),dimension(:),allocatable :: ifirst !(maxprocs)
       integer(kind=2),dimension(:,:),allocatable :: nupa,&
        ndowna !(0:maxprocs/4,0:maxprocs)
-      real(kind=8),dimension(:,:),allocatable :: t_restart1 !(5,maxprocs)
+      real(kind=4),dimension(:,:),allocatable :: t_restart1 !(5,maxprocs)
       integer,dimension(:),allocatable :: iset_restart1 !(maxprocs)
 !      common /traj1cache/
-      real(kind=8),dimension(:),allocatable :: totT_cache,EK_cache,&
+      real(kind=4),dimension(:),allocatable :: totT_cache,EK_cache,&
        potE_cache,t_bath_cache,Uconst_cache !(max_cache_traj)
-      real(kind=8),dimension(:,:),allocatable :: qfrag_cache !(50,max_cache_traj)
-      real(kind=8),dimension(:,:),allocatable :: qpair_cache !(100,max_cache_traj)
-      real(kind=8),dimension(:,:),allocatable :: ugamma_cache,&
+      real(kind=4),dimension(:,:),allocatable :: qfrag_cache !(50,max_cache_traj)
+      real(kind=4),dimension(:,:),allocatable :: qpair_cache !(100,max_cache_traj)
+      real(kind=4),dimension(:,:),allocatable :: ugamma_cache,&
        utheta_cache,uscdiff_cache !(maxfrag_back,max_cache_traj)
-      real(kind=8),dimension(:,:,:),allocatable :: c_cache !(3,maxres2+2,max_cache_traj)
+      real(kind=4),dimension(:,:,:),allocatable :: c_cache !(3,maxres2+2,max_cache_traj)
       integer :: ntwx_cache,ii_write   !,max_cache_traj_use
       integer,dimension(:),allocatable :: iset_cache !(max_cache_traj)
 !-----------------------------------------------------------------------------
       integer :: ERRCODE
       real(kind=8),dimension(3) :: L,vcm
       real(kind=8) :: energia(0:n_ene)
-      real(kind=8) :: remd_t_bath(maxprocs)
-      integer :: iremd_iset(maxprocs)
-      integer(kind=2) :: i_index(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
-      real(kind=8) :: remd_ene(0:n_ene+4,maxprocs)
-      integer :: iremd_acc(maxprocs),iremd_tot(maxprocs)
-      integer :: iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
+      real(kind=8) :: remd_t_bath(Nprocs) !(maxprocs)
+      integer :: iremd_iset(Nprocs) !(maxprocs)
+      integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
+! (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
+      real(kind=8) :: remd_ene(0:n_ene+4,Nprocs) !(0:n_ene+4,maxprocs)
+      integer :: iremd_acc(Nprocs),iremd_tot(Nprocs) !(maxprocs)
+      integer :: iremd_acc_usa(Nprocs),iremd_tot_usa(Nprocs) !(maxprocs)
       integer :: rstcount      !el ilen,
 !el      external ilen
       character(len=50) :: tytul
 !el      common /gucio/ cm
       integer :: itime
 !old      integer nup(0:maxprocs),ndown(0:maxprocs)
-      integer :: rep2i(0:maxprocs),ireqi(maxprocs)
-      integer :: icache_all(maxprocs)
-      integer :: status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
+      integer :: rep2i(0:Nprocs),ireqi(Nprocs) !(maxprocs)
+      integer :: icache_all(Nprocs) !(maxprocs)
+      integer :: status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,Nprocs)! (MPI_STATUS_SIZE,maxprocs)
       logical :: synflag, end_of_run, file_exist = .false.!, ovrtim
 
       real(kind=8) :: delta,time00,time01,time001,time02,time03,time04,&
       nres2=2*nres
       time001=0.0d0
 
-write(iout,*) "jestesmy na poczatku MREMD"
       ntwx_cache=0
       time00=MPI_WTIME()
       time01=time00
@@ -121,7 +121,6 @@ write(iout,*) "jestesmy na poczatku MREMD"
        write (iout,*) "NREP=",nrep
       endif
 
-write(iout,*) "jestesmy na poczatku MREMD"
       synflag=.false.
       if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
         call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
@@ -194,8 +193,8 @@ write(iout,*) "jestesmy na poczatku MREMD"
 !d       print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
 
 !el  common /przechowalnia/
-       if(.not.allocated(d_restart1)) allocate(d_restart1(3,nres2*nodes))
-       if(.not.allocated(d_restart2)) allocate(d_restart2(3,nres2*nodes))
+       if(.not.allocated(d_restart1)) allocate(d_restart1(3,(nres2+1)*nodes))
+       if(.not.allocated(d_restart2)) allocate(d_restart2(3,(nres2+1)*nodes))
        if(.not.allocated(p_c)) allocate(p_c(3,(nres2+2)*nodes))
 !el-------------
 
@@ -348,11 +347,15 @@ write(iout,*) "jestesmy na poczatku MREMD"
           if(me.eq.king.or..not.out1file) &
            write(iout,*) me,"iset=",iset,"t_bath=",t_bath
        endif        
-!
-       stdfp=dsqrt(2*Rb*t_bath/d_time)
+!    
+       do i=1,5
+       stdfp(i)=dsqrt(2*Rb*t_bath/d_time)
+       enddo
+       do j=1,5
        do i=1,ntyp
-          stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
+          stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
        enddo 
+       enddo
 
 !      print *,'irep',me,t_bath
        if (.not.rest) then  
@@ -384,10 +387,12 @@ write(iout,*) "jestesmy na poczatku MREMD"
       if (rest) then  
        if (me.eq.king .or. .not. out1file) &
         write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
-       stdfp=dsqrt(2*Rb*t_bath/d_time)
+       do j=1,5
+       stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
        do i=1,ntyp
-          stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
+          stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
        enddo 
+       enddo
        call rescale_weights(t_bath)
       endif
 
@@ -831,6 +836,7 @@ write(iout,*) "jestesmy na poczatku MREMD"
                remd_t_bath(i)=remd_ene(n_ene+1,i)
                iremd_iset(i)=remd_ene(n_ene+2,i)
             enddo
+#ifdef DEBUG
             if(lmuca) then
 !o             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
              do i=1,nodes
@@ -844,6 +850,7 @@ write(iout,*) "jestesmy na poczatku MREMD"
                 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
               enddo
             endif
+#endif
 !-------------------------------------           
            IF(.not.usampl) THEN
             write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
@@ -857,7 +864,9 @@ write(iout,*) "jestesmy na poczatku MREMD"
 
              do ii=1,nodes-1
 
+#ifdef DEBUG
               write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
+#endif
              if(i.gt.0.and.nupa(0,i).gt.0) then
               iex=i
 !              if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
@@ -1216,6 +1225,12 @@ write(iout,*) "jestesmy na poczatku MREMD"
          time06=MPI_WTIME()
 !d         write (iout,*) "Before scatter"
 !d         call flush(iout)
+#ifdef DEBUG
+         if (me.eq.king) then
+           write (iout,*) "t_bath before scatter",remd_t_bath
+           call flush(iout)
+         endif
+#endif
          call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
                  t_bath,1,mpi_double_precision,king,&
                  CG_COMM,ierr) 
@@ -1242,16 +1257,25 @@ write(iout,*) "jestesmy na poczatku MREMD"
          call rescale_weights(t_bath)
 !o         write (iout,*) "Processor",me,
 !o     &    " rescaling weights with temperature",t_bath
-
-         stdfp=dsqrt(2*Rb*t_bath/d_time)
+         do j=1,5
+         stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
          do i=1,ntyp
-           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
+           stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
          enddo 
+         enddo
 
 !de         write(iout,*) 'REMD after',me,t_bath
            time08=MPI_WTIME()
            if (me.eq.king .or. .not. out1file) then
-            write(iout,*) 'REMD exchange time=',time08-time00
+            write(iout,*) 'REMD exchange time 8-0=',time08-time00
+            write(iout,*) 'REMD exchange time 8-7=',time08-time07
+            write(iout,*) 'REMD exchange time 7-6=',time07-time06
+            write(iout,*) 'REMD exchange time 6-5=',time06-time05
+            write(iout,*) 'REMD exchange time 5-4=',time05-time04
+            write(iout,*) 'REMD exchange time 4-3=',time04-time03
+            write(iout,*) 'REMD exchange time 3-2=',time03-time02
+            write(iout,*) 'REMD exchange time 2-1=',time02-time01
+            write(iout,*) 'REMD exchange time 1-0=',time01-time00
             call flush(iout)
            endif
         endif
@@ -1291,9 +1315,9 @@ write(iout,*) "jestesmy na poczatku MREMD"
        '  End of MD calculation  '
       endif
 !el  common /przechowalnia/
-      deallocate(d_restart1)
-      deallocate(d_restart2)
-      deallocate(p_c)
+!      deallocate(d_restart1)
+!      deallocate(d_restart2)
+!      deallocate(p_c)
 !el--------------
       return
       end subroutine MREMD
@@ -1314,10 +1338,11 @@ write(iout,*) "jestesmy na poczatku MREMD"
                
 !el      real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
 !el           d_restart2(3,2*nres*maxprocs)
-      real(kind=4) :: r_d(3,2*nres)
+      real(kind=4) :: r_d(3,0:2*nres)
       real(kind=4) :: t5_restart1(5)
       integer :: iret,itmp
-      integer(kind=2) :: i_index(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
+      integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
+      !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
 !el       common /przechowalnia/ d_restart1,d_restart2
       integer :: i,j,il,il1,ierr,ixdrf
 
@@ -1331,23 +1356,26 @@ write(iout,*) "jestesmy na poczatku MREMD"
             t_restart1,5,mpi_real,king,CG_COMM,ierr)
 
 
-       do i=1,2*nres
+       do i=0,2*nres
          do j=1,3
            r_d(j,i)=d_t(j,i)
          enddo
        enddo
-       call mpi_gather(r_d,3*2*nres,mpi_real,&
-                 d_restart1,3*2*nres,mpi_real,king,&
+       call mpi_gather(r_d,3*2*nres+3,mpi_real,&
+                 d_restart1,3*2*nres+3,mpi_real,king,&
                  CG_COMM,ierr)
+       do j=1,3
+       dc(j,0)=c(j,1)
+       enddo
 
 
-       do i=1,2*nres
+       do i=0,2*nres
          do j=1,3
            r_d(j,i)=dc(j,i)
          enddo
        enddo
-       call mpi_gather(r_d,3*2*nres,mpi_real,&
-                 d_restart2,3*2*nres,mpi_real,king,&
+       call mpi_gather(r_d,3*2*nres+3,mpi_real,&
+                 d_restart2,3*2*nres+3,mpi_real,king,&
                  CG_COMM,ierr)
 
        if(me.eq.king) then
@@ -1376,16 +1404,16 @@ write(iout,*) "jestesmy na poczatku MREMD"
          enddo
 
          do il=0,nodes-1
-           do i=1,2*nres
+           do i=0,2*nres
             do j=1,3
-             call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
             enddo
            enddo
          enddo
          do il=0,nodes-1
-           do i=1,2*nres
+           do i=0,2*nres
             do j=1,3
-             call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
+             call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
             enddo
            enddo
          enddo
@@ -1436,16 +1464,16 @@ write(iout,*) "jestesmy na poczatku MREMD"
          enddo
 
          do il=0,nodes-1
-           do i=1,2*nres
+           do i=0,2*nres
             do j=1,3
-             call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
             enddo
            enddo
          enddo
          do il=0,nodes-1
-           do i=1,2*nres
+           do i=0,2*nres
             do j=1,3
-             call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
+             call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
             enddo
            enddo
          enddo
@@ -1495,9 +1523,9 @@ write(iout,*) "jestesmy na poczatku MREMD"
       real(kind=4) :: xcoord(3,2*nres+2),prec
       real(kind=4) :: r_qfrag(50),r_qpair(100)
       real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
-      real(kind=4) :: p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
-      real(kind=4) :: p_utheta(50*maxprocs),p_ugamma(100*maxprocs),&
-           p_uscdiff(100*maxprocs)
+      real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
+      real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
+           p_uscdiff(100*Nprocs) !(100*maxprocs)
 !el      real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
       real(kind=4) :: r_c(3,2*nres+2)
 !el      common /przechowalnia/ p_c
@@ -1517,13 +1545,22 @@ write(iout,*) "jestesmy na poczatku MREMD"
       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
 #endif
       do ii=1,ii_write
+!       write (iout,*) "before gather write1traj: from node",ii
+!       call flush(iout)
+!       write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
+!       call flush(iout)
        t5_restart1(1)=totT_cache(ii)
        t5_restart1(2)=EK_cache(ii)
        t5_restart1(3)=potE_cache(ii)
        t5_restart1(4)=t_bath_cache(ii)
        t5_restart1(5)=Uconst_cache(ii)
+!       write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
+       call flush(iout)
        call mpi_gather(t5_restart1,5,mpi_real,&
             t_restart1,5,mpi_real,king,CG_COMM,ierr)
+!       do il=1,nodes
+!       write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
+!       enddo
 
        call mpi_gather(iset_cache(ii),1,mpi_integer,&
             iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
@@ -1627,6 +1664,7 @@ write(iout,*) "jestesmy na poczatku MREMD"
           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
+!          write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
           call xdrfint(ixdrf, nss, iret) 
           do j=1,nss
            if (dyn_ss) then
@@ -1716,8 +1754,9 @@ write(iout,*) "jestesmy na poczatku MREMD"
 !      include 'COMMON.SBRIDGE'
 !      include 'COMMON.INTERACT'
 !el      real(kind=4) :: d_restart1(3,2*nres*maxprocs)
-      real(kind=4) :: r_d(3,2*nres),t5_restart1(5)
-      integer(kind=2) :: i_index(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
+      real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
+      integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
+      !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
 !el      common /przechowalnia/ d_restart1
       integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
 
@@ -1792,45 +1831,45 @@ write(iout,*) "jestesmy na poczatku MREMD"
 
          if(me.eq.king)then
               do il=0,nodes-1
-               do i=1,2*nres
+               do i=0,2*nres
 !                read(irest2,'(3e15.5)') 
 !     &                (d_restart1(j,i+2*nres*il),j=1,3)
             do j=1,3
 #ifdef AIX
-             call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
 #else
-             call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
 #endif
             enddo
                enddo
               enddo
          endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
-                 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
+         call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
+                 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
 
-         do i=1,2*nres
+         do i=0,2*nres
            do j=1,3
             d_t(j,i)=r_d(j,i)
            enddo
          enddo
          if(me.eq.king)then 
               do il=0,nodes-1
-               do i=1,2*nres
+               do i=0,2*nres
 !                read(irest2,'(3e15.5)') 
 !     &                (d_restart1(j,i+2*nres*il),j=1,3)
             do j=1,3
 #ifdef AIX
-             call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
 #else
-             call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
+             call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
 #endif
             enddo
                enddo
               enddo
          endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
-                 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
-         do i=1,2*nres
+         call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
+                 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
+         do i=0,2*nres
            do j=1,3
             dc(j,i)=r_d(j,i)
            enddo
@@ -1902,7 +1941,7 @@ write(iout,*) "jestesmy na poczatku MREMD"
 !      include 'COMMON.SBRIDGE'
 !      include 'COMMON.INTERACT'
 !el      real(kind=4) :: d_restart1(3,2*nres*maxprocs)
-      real(kind=4) :: r_d(3,2*nres),t5_restart1(5)
+      real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
 !el      common /przechowalnia/ d_restart1
 
       integer :: i,j,il,ierr
@@ -1929,16 +1968,16 @@ write(iout,*) "jestesmy na poczatku MREMD"
 
          if(me.eq.king)then
               do il=0,nodes-1
-               do i=1,2*nres
+               do i=0,2*nres
                 read(irest2,'(3e15.5)')  &
-                      (d_restart1(j,i+2*nres*il),j=1,3)
+                      (d_restart1(j,i+(2*nres+1)*il),j=1,3)
                enddo
               enddo
          endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
-                 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
+         call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
+                 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
 
-         do i=1,2*nres
+         do i=0,2*nres
            do j=1,3
             d_t(j,i)=r_d(j,i)
            enddo
@@ -1947,13 +1986,13 @@ write(iout,*) "jestesmy na poczatku MREMD"
               do il=0,nodes-1
                do i=1,2*nres
                 read(irest2,'(3e15.5)') &
-                      (d_restart1(j,i+2*nres*il),j=1,3)
+                      (d_restart1(j,i+(2*nres+1)*il),j=1,3)
                enddo
               enddo
          endif
-         call mpi_scatter(d_restart1,3*2*nres,mpi_real,&
-                 r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
-         do i=1,2*nres
+         call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
+                 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
+         do i=0,2*nres
            do j=1,3
             dc(j,i)=r_d(j,i)
            enddo
@@ -1965,7 +2004,7 @@ write(iout,*) "jestesmy na poczatku MREMD"
       subroutine alloc_MREMD_arrays
 
 !      if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
-      if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1)) !(ntyp1))
+      if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1,5)) !(ntyp1))
 ! commom.remd
 !      common /remdcommon/ in io: read_REMDpar
 !      real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)