projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
2D replica exchange with homology constraints
[unres.git] / source / unres / src_MD / MREMD.F
diff --git a/source/unres/src_MD/MREMD.F b/source/unres/src_MD/MREMD.F
index be6af9c..9301122 100644 (file)
--- a/source/unres/src_MD/MREMD.F
+++ b/source/unres/src_MD/MREMD.F
@@ -70,6 +70,16 @@ cde      write (iout,*) "Start MREMD: me",me," t_bath",t_bath
          enddo
       endif
 
+      if(homol_nset.gt.1) then
+         i_econstr=24
+         nset=homol_nset
+         do i=1,nset
+          mset(i)=1
+         enddo
+      endif
+
+      if(usampl) i_econstr=20
+
       k=0
       rep2i(k)=-1
       do il=1,max0(nset,1)
@@ -187,7 +197,7 @@ cd           write (*,*) me," After broadcast: file_exist",file_exist
               read (irest2,*) ndowna(0,il),
      &                    (ndowna(i,il),i=1,ndowna(0,il))
              enddo
-             if(usampl.or.hremd.gt.0) then
+             if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
               read (irest2,*)
               read (irest2,*) nset
               read (irest2,*) 
@@ -278,7 +288,7 @@ cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
          if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
 
        endif
-       if(usampl.or.hremd.gt.0) then
+       if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
           iset=i2set(me)
           if (hremd.gt.0) call set_hweights(iset)
           if(me.eq.king.or..not.out1file) 
@@ -572,7 +582,7 @@ c Variable time step algorithm.
               write (irest1,*) ndowna(0,il),
      &                   (ndowna(i,il),i=1,ndowna(0,il))
              enddo
-             if(usampl.or.hremd.gt.0) then
+             if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
               write (irest1,*) "nset"
               write (irest1,*) nset
               write (irest1,*) "mset"
@@ -599,7 +609,7 @@ c Variable time step algorithm.
            do i=1,2*nres
             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
            enddo
-           if(usampl.or.hremd.gt.0) then
+           if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then
              write (irest2,*) iset
            endif
           close(irest2)
@@ -717,20 +727,31 @@ c     &             remd_t_bath,1,mpi_double_precision,king,
 c     &             CG_COMM,ierr)
            potEcomp(n_ene+1)=t_bath
            t_bath_old=t_bath
-           if (usampl) then
+           if (usampl.or.homol_nset.gt.1) then
              potEcomp(n_ene+2)=iset
+c             write(iout,*) potEcomp(n_ene+1),potEcomp(n_ene+2),iset,nset
              if (iset.lt.nset) then
                i_set_temp=iset
                iset=iset+1
-               call EconstrQ
-               potEcomp(n_ene+3)=Uconst
+               if (homol_nset.gt.1) then
+                call e_modeller(potEcomp(n_ene+3))
+c                write(iout,*) "iset+1",potEcomp(n_ene+3)
+               else
+                call EconstrQ
+                potEcomp(n_ene+3)=Uconst
+               endif
                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 (homol_nset.gt.1) then
+                call e_modeller(potEcomp(n_ene+4))
+c                write(iout,*) "iset-1",potEcomp(n_ene+4)
+               else
+                call EconstrQ
+                potEcomp(n_ene+4)=Uconst
+               endif
                iset=i_set_temp
              endif
            endif
@@ -778,9 +799,15 @@ ctime            call flush(iout)
 
 
           if (me.eq.king) then
+            write(iout,*) 
+     &     'energy_c temperature iset energy_c(iset+1) energy_c(iset-1)'
             do i=1,nodes
                remd_t_bath(i)=remd_ene(n_ene+1,i)
                iremd_iset(i)=remd_ene(n_ene+2,i)
+               write(iout,'(i4,f10.3,f6.0,i3,2f10.3)') 
+     &                i,remd_ene(i_econstr,i),
+     &                remd_ene(n_ene+1,i),iremd_iset(i),
+     &                remd_ene(n_ene+3,i),remd_ene(n_ene+4,i)
             enddo
 #ifdef DEBUG
             if(lmuca) then
@@ -798,7 +825,7 @@ co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
             endif
 #endif
 c-------------------------------------           
-           IF(.not.usampl.and.hremd.eq.0) THEN
+           IF(.not.usampl.and.hremd.eq.0.and.homol_nset.le.1) THEN
 #ifdef DEBUG
             write (iout,*) "Enter exchnge, remd_m",remd_m(1),
      &        " nodes",nodes
@@ -982,9 +1009,9 @@ c                call flush(iout)
            enddo
 cd           write (iout,*) "exchange completed"
 cd           call flush(iout) 
-        ELSEIF (usampl) THEN
+        ELSEIF (usampl.or.homol_nset.gt.1) THEN
           do ii=1,nodes  
-cd            write(iout,*) "########",ii
+c            write(iout,*) "########",ii
 
             i_temp=iran_num(1,nrep)
             i_mult=iran_num(1,remd_m(i_temp))
@@ -992,14 +1019,14 @@ cd            write(iout,*) "########",ii
             i_mset=iran_num(1,mset(i_iset))
             i=i_index(i_temp,i_mult,i_iset,i_mset)
 
-cd            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
+c            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
 
             if(t_exchange_only)then
              i_dir=1
             else
              i_dir=iran_num(1,3)
             endif
-cd            write(iout,*) "i_dir=",i_dir
+c            write(iout,*) "i_dir=",i_dir
 
             if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then            
                
@@ -1016,10 +1043,11 @@ cd            write(iout,*) "i_dir=",i_dir
                i_iset1=i_iset+1
                i_mset1=iran_num(1,mset(i_iset1))
                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)
+                
+               econstr_temp_i=remd_ene(i_econstr,i)
+               econstr_temp_iex=remd_ene(i_econstr,iex)
+               remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
+               remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
 
             elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
 
@@ -1028,16 +1056,16 @@ cd            write(iout,*) "i_dir=",i_dir
                i_iset1=i_iset+1
                i_mset1=iran_num(1,mset(i_iset1))
                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)
+               econstr_temp_i=remd_ene(i_econstr,i)
+               econstr_temp_iex=remd_ene(i_econstr,iex)
+               remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
+               remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
 
             else
                goto 444 
             endif
  
-cd            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
+c            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
 ctime            call flush(iout)
 
 c Swap temperatures between conformations i and iex with recalculating the free energies
@@ -1050,33 +1078,39 @@ co     &          remd_t_bath(i)
               
               call sum_energy(remd_ene(0,iex),.false.)
               ene_iex_i=remd_ene(0,iex)
-cd              write (iout,*) "ene_iex_i",remd_ene(0,iex)
+cdebug
+c ERROR only makes sense for dir =1
+c              write (iout,*) "ene_iex_i",remd_ene(0,iex)
 c              call sum_energy(remd_ene(0,i),.false.)
-cd              write (iout,*) "ene_i_i",remd_ene(0,i)
+c              write (iout,*) "ene_i_i",remd_ene(0,i)
 c              write (iout,*) "rescaling weights with temperature",
 c     &          remd_t_bath(iex)
 c              if (real(ene_i_i).ne.real(remd_ene(0,i))) then
-c                write (iout,*) "ERROR: inconsistent energies:",i,
+c                write (iout,*) "ERROR: inconsistent energies i:",i,
 c     &            ene_i_i,remd_ene(0,i)
 c              endif
+cdebug_end
               call rescale_weights(remd_t_bath(iex))
               call sum_energy(remd_ene(0,i),.false.)
 cd              write (iout,*) "ene_i_iex",remd_ene(0,i)
               ene_i_iex=remd_ene(0,i)
+cdebug
+c ERROR only makes sense for dir =1
 c              call sum_energy(remd_ene(0,iex),.false.)
 c              if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
-c                write (iout,*) "ERROR: inconsistent energies:",iex,
+c                write (iout,*) "ERROR: inconsistent energies iex:",iex,
 c     &            ene_iex_iex,remd_ene(0,iex)
 c              endif
-cd              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
+c              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
 c              write (iout,*) "i",i," iex",iex
-cd              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
-cd     &           " ene_i_iex",ene_i_iex,
-cd     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_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
+cdebug_end
               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
-cd              write(iout,*) 'delta',delta
+c              write(iout,*) 'delta',delta
 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
 c     &              (remd_t_bath(i)*remd_t_bath(iex))
@@ -1091,7 +1125,7 @@ c     &              (remd_t_bath(i)*remd_t_bath(iex))
      &          iremd_tot_usa(int(i2set(i-1)))=
      &                 iremd_tot_usa(int(i2set(i-1)))+1
               xxx=ran_number(0.0d0,1.0d0)
-cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
+c              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
               if (delta .gt. xxx) then
                 tmp=remd_t_bath(i)       
                 remd_t_bath(i)=remd_t_bath(iex)
@@ -1127,8 +1161,8 @@ cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
               else
                remd_ene(0,iex)=ene_iex_iex
                remd_ene(0,i)=ene_i_i
-               remd_ene(20,iex)=econstr_temp_iex
-               remd_ene(20,i)=econstr_temp_i
+               remd_ene(i_econstr,iex)=econstr_temp_iex
+               remd_ene(i_econstr,i)=econstr_temp_i
               endif
 
 cd      do il=1,nset
@@ -1301,7 +1335,7 @@ c-------------------------------------
      &           ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
              enddo
 
-             if(usampl) then
+             if(usampl.or.homol_nset.gt.1) then
               do i=1,nset
                if(iremd_tot_usa(i).ne.0)
      &           write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
@@ -1338,7 +1372,7 @@ cd         call flush(iout)
      &           CG_COMM,ierr) 
 cd         write (iout,*) "After scatter"
 cd         call flush(iout)
-         if(usampl.or.hremd.gt.0)
+         if(usampl.or.hremd.gt.0.or.homol_nset.gt.1)
      &    call mpi_scatter(iremd_iset,1,mpi_integer,
      &           iset,1,mpi_integer,king,
      &           CG_COMM,ierr) 
@@ -1511,7 +1545,7 @@ c-----------------------------------------------------------------------
            enddo
          enddo
 
-         if(usampl) then
+         if(usampl.or.homol_nset.gt.1) then
            call xdrfint_(ixdrf, nset, iret)
            do i=1,nset
              call xdrfint_(ixdrf,mset(i), iret)
@@ -1572,7 +1606,7 @@ c-----------------------------------------------------------------------
          enddo
 
 
-             if(usampl) then
+             if(usampl.or.homol_nset.gt.1) then
               call xdrfint(ixdrf, nset, iret)
               do i=1,nset
                 call xdrfint(ixdrf,mset(i), iret)
@@ -1954,7 +1988,7 @@ c     &                (d_restart1(j,i+2*nres*il),j=1,3)
          enddo
        
 
-           if(usampl) then
+           if(usampl.or.homol_nset.gt.1) then
 #ifdef AIX
              if(me.eq.king)then
               call xdrfint_(ixdrf, nset, iret)
UNRESPACK 3.x