update
[unres.git] / source / unres / src_MD-M / MREMD.F
index 1df3f0a..afecaa5 100644 (file)
@@ -29,7 +29,7 @@
       integer iremd_iset(maxprocs)
       integer*2 i_index
      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
-      double precision remd_ene(0:n_ene+4,maxprocs)
+      double precision remd_ene(0:n_ene+8,maxprocs)
       integer iremd_acc(maxprocs),iremd_tot(maxprocs)
       integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
       integer ilen,rstcount
@@ -272,6 +272,7 @@ cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
           iset=i2set(me)
           if(me.eq.king.or..not.out1file) 
      &     write(iout,*) me,"iset=",iset,"t_bath=",t_bath
+          call flush(iout)
        endif        
 c
        stdfp=dsqrt(2*Rb*t_bath/d_time)
@@ -525,7 +526,9 @@ c Variable time step algorithm.
               ugamma_cache(i,ntwx_cache)=ugamma(i)
               uscdiff_cache(i,ntwx_cache)=uscdiff(i)
             enddo
-
+C            print *,'przed returnbox'
+            call returnbox
+C            call enerprint(remd_ene(0,i))
             do i=1,nres*2
              do j=1,3
               c_cache(j,i,ntwx_cache)=c(j,i)
@@ -697,19 +700,78 @@ c     &             CG_COMM,ierr)
                i_set_temp=iset
                iset=iset+1
                call EconstrQ
-               potEcomp(n_ene+3)=Uconst
+               if (loc_qlike) then
+                 call Econstr_back_qlike
+               else
+                 call Econstr_back
+               endif
+               potEcomp(n_ene+3)=Uconst+Uconst_back
                iset=i_set_temp
              endif
              if (iset.gt.1) then
                i_set_temp=iset
                iset=iset-1
                call EconstrQ
-               potEcomp(n_ene+4)=Uconst 
+               if (loc_qlike) then
+                 call Econstr_back_qlike
+               else
+                 call Econstr_back
+               endif
+               potEcomp(n_ene+4)=Uconst+Uconst_back
                iset=i_set_temp
              endif
+             if (adaptive) then
+c 9/11/17 Adaptive US
+               itt=i2rep(me)
+#ifdef DEBUG
+               write (iout,*) "me ",me," itt",itt
+#endif
+               if (itt.lt.nrep) then
+                 t_bath_temp=t_bath
+                 t_bath=remd_t(itt+1)
+                 call EconstrQ
+                 potEcomp(n_ene+5)=Uconst
+#ifdef DEBUG
+                 write (iout,*) "t_bath",t_bath_temp,t_bath,
+     &                   " Uconst",Uconst
+#endif
+                 if (iset.lt.nset) then
+                   i_set_temp=iset
+                   iset=iset+1
+                   call EconstrQ
+                   potEcomp(n_ene+7)=Uconst
+#ifdef DEBUG
+                   write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
+#endif
+                   iset=i_set_temp
+                 endif
+                 t_bath=t_bath_temp
+               endif
+               if (itt.gt.1) then
+                 t_bath_temp=t_bath
+                 t_bath=remd_t(itt-1)
+                 call EconstrQ
+                 potEcomp(n_ene+6)=Uconst
+#ifdef DEBUG
+                 write (iout,*) "t_bath",t_bath_temp,t_bath,
+     &                   " Uconst",Uconst
+#endif
+                 if (iset.gt.1) then
+                   i_set_temp=iset
+                   iset=iset-1
+                   call EconstrQ
+                   potEcomp(n_ene+8)=Uconst
+#ifdef DEBUG
+                   write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
+#endif
+                   iset=i_set_temp
+                 endif
+                 t_bath=t_bath_temp
+               endif
+             endif
            endif
-           call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
-     &             remd_ene(0,1),n_ene+5,mpi_double_precision,king,
+           call mpi_gather(potEcomp(0),n_ene+9,mpi_double_precision,
+     &             remd_ene(0,1),n_ene+9,mpi_double_precision,king,
      &             CG_COMM,ierr)
            if(lmuca) then 
             call mpi_gather(elow,1,mpi_double_precision,
@@ -755,6 +817,7 @@ cd end
                remd_t_bath(i)=remd_ene(n_ene+1,i)
                iremd_iset(i)=remd_ene(n_ene+2,i)
             enddo
+            if (mremd_dec) then
             if(lmuca) then
 co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
              do i=1,nodes
@@ -768,20 +831,21 @@ co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
                 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
               enddo
             endif
+            endif
 c-------------------------------------           
            IF(.not.usampl) THEN
+            if (mremd_dec) then
             write (iout,*) "Enter exchnge, remd_m",remd_m(1),
      &        " nodes",nodes
-            call flush(iout)
             write (iout,*) "remd_m(1)",remd_m(1)
+            endif
             do irr=1,remd_m(1)
                i=ifirst(iran_num(1,remd_m(1)))
-             write (iout,*) "i",i
-             call flush(iout)
 
              do ii=1,nodes-1
 
-              write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
+              if (mremd_dec) 
+     &          write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
              if(i.gt.0.and.nupa(0,i).gt.0) then
               iex=i
 c              if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
@@ -795,6 +859,11 @@ c                write (iout,*) "upper",nupa(int(nupa(0,i)),i)
                 iex=nupa(iran_num(1,int(nupa(0,i))),i)
 c              enddo
 c              write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
+#ifdef DEBUG
+c 8/21/17 AL : debug
+              write (iout,*) "i",i,"iex",iex," temperatures",
+     &            remd_t_bath(i),remd_t_bath(iex)
+#endif
               if (lmuca) then
                call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
               else
@@ -846,12 +915,14 @@ c               call enerprint(remd_ene(0,iex))
                 write (iout,*) "ERROR: inconsistent energies:",iex,
      &            ene_iex_iex,remd_ene(0,iex)
                endif
-c               write (iout,*) "ene_iex_iex",remd_ene(0,iex)
-c               write (iout,*) "i",i," iex",iex
-c               write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
-c     &           " ene_i_iex",ene_i_iex,
-c     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
-c               call flush(iout)
+#ifdef DEBUG
+               write (iout,*) "ene_iex_iex",remd_ene(0,iex)
+               write (iout,*) "i",i," iex",iex
+               write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
+     &           " ene_i_iex",ene_i_iex,
+     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
+               call flush(iout)
+#endif
                delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
                delta=-delta
@@ -925,18 +996,18 @@ c              call flush(iout)
                 itmp=i2rep(i-1)
                 i2rep(i-1)=i2rep(iex-1)
                 i2rep(iex-1)=itmp
-
-c                write(iout,*) 'exchange',i,iex
-c                write (iout,'(a8,100i4)') "@ ifirst",
-c     &                    (ifirst(k),k=1,remd_m(1))
-c                do il=1,nodes
-c                 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
-c     &                    (nupa(k,il),k=1,nupa(0,il))
-c                 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
-c     &                    (ndowna(k,il),k=1,ndowna(0,il))
-c                enddo
-c                call flush(iout) 
-
+#ifdef DEBUG
+                write(iout,*) 'exchange',i,iex
+                write (iout,'(a8,100i4)') "@ ifirst",
+     &                    (ifirst(k),k=1,remd_m(1))
+                do il=1,nodes
+                 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
+     &                    (nupa(k,il),k=1,nupa(0,il))
+                 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
+     &                    (ndowna(k,il),k=1,ndowna(0,il))
+                enddo
+                call flush(iout) 
+#endif
               else
                remd_ene(0,iex)=ene_iex_iex
                remd_ene(0,i)=ene_i_i
@@ -959,7 +1030,7 @@ cd            write(iout,*) "########",ii
 
 cd            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
 
-            i_dir=iran_num(1,3)
+             i_dir=iran_num(1,3)
 cd            write(iout,*) "i_dir=",i_dir
 
             if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then            
@@ -969,7 +1040,15 @@ cd            write(iout,*) "i_dir=",i_dir
                i_iset1=i_iset
                i_mset1=iran_num(1,mset(i_iset1))
                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
-
+c 9/1/17 AL: Correction; otherwise the restraint energies are mis-assigned 
+c       on failed replica exchange attempt
+               econstr_temp_i=remd_ene(20,i)
+               econstr_temp_iex=remd_ene(20,iex)
+c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
+               if (adaptive) then
+                 remd_ene(20,i)=remd_ene(n_ene+5,i)
+                 remd_ene(20,iex)=remd_ene(n_ene+6,iex)
+               endif
             elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
 
                i_temp1=i_temp
@@ -991,9 +1070,13 @@ cd            write(iout,*) "i_dir=",i_dir
                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
                econstr_temp_i=remd_ene(20,i)
                econstr_temp_iex=remd_ene(20,iex)
-               remd_ene(20,i)=remd_ene(n_ene+3,i)
-               remd_ene(20,iex)=remd_ene(n_ene+4,iex)
-
+               if (adaptive) then
+                 remd_ene(20,i)=remd_ene(n_ene+7,i)
+                 remd_ene(20,iex)=remd_ene(n_ene+8,iex)
+               else
+                 remd_ene(20,i)=remd_ene(n_ene+3,i)
+                 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
+               endif
             else
                goto 444 
             endif
@@ -1003,6 +1086,10 @@ cd            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
 
 c Swap temperatures between conformations i and iex with recalculating the free energies
 c following temperature changes.
+#ifdef DEBUG
+              write (iout,*) "i_dir",i_dir," i",i," iex",iex,
+     &           " t_bath",remd_t_bath(i),remd_t_bath(iex)
+#endif
               ene_iex_iex=remd_ene(0,iex)
               ene_i_i=remd_ene(0,i)
 co              write (iout,*) "rescaling weights with temperature",
@@ -1170,7 +1257,7 @@ co     &    " rescaling weights with temperature",t_bath
          stdfp=dsqrt(2*Rb*t_bath/d_time)
          do i=1,ntyp
            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
-         enddo 
+         enddo
 
 cde         write(iout,*) 'REMD after',me,t_bath
            time08=MPI_WTIME()
@@ -1422,7 +1509,7 @@ c-----------------------------------------------------------------------
      &           king,CG_COMM,ierr)
 
 c debugging
-      print *,'traj1file',me,ii_write,ntwx_cache
+c      print *,'traj1file',me,ii_write,ntwx_cache
 c end debugging
 
 #ifdef AIX
@@ -1500,8 +1587,13 @@ c end debugging
           call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
           call xdrfint_(ixdrf, nss, iret) 
           do j=1,nss
-           call xdrfint_(ixdrf, ihpb(j), iret)
-           call xdrfint_(ixdrf, jhpb(j), iret)
+           if (dyn_ss) then
+            call xdrfint(ixdrf, idssb(j)+nres, iret)
+            call xdrfint(ixdrf, jdssb(j)+nres, iret)
+           else
+            call xdrfint_(ixdrf, ihpb(j), iret)
+            call xdrfint_(ixdrf, jhpb(j), iret)
+           endif
           enddo
           call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
           call xdrfint_(ixdrf, iset_restart1(il), iret)
@@ -1538,8 +1630,13 @@ c end debugging
           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
           call xdrfint(ixdrf, nss, iret) 
           do j=1,nss
-           call xdrfint(ixdrf, ihpb(j), iret)
-           call xdrfint(ixdrf, jhpb(j), iret)
+           if (dyn_ss) then
+            call xdrfint(ixdrf, idssb(j)+nres, iret)
+            call xdrfint(ixdrf, jdssb(j)+nres, iret)
+           else
+            call xdrfint(ixdrf, ihpb(j), iret)
+            call xdrfint(ixdrf, jhpb(j), iret)
+           endif
           enddo
           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
           call xdrfint(ixdrf, iset_restart1(il), iret)
@@ -1782,7 +1879,7 @@ c     &                (d_restart1(j,i+2*nres*il),j=1,3)
               enddo
              endif
 #endif
-c Corrected AL 8/19/2014: each processor needs whole iset array not only its
+Corrected AL 8/19/2014: each processor needs whole iset array not only its
 c own element
 c              call mpi_scatter(i2set,1,mpi_integer,
 c     &           iset,1,mpi_integer,king,
@@ -1866,4 +1963,4 @@ c     &           CG_COMM,ierr)
         if(me.eq.king) close(irest2)
         return
         end
-          
+c------------------------------------------