X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fsrc_MD-M%2FMREMD.F;h=7a70f1db18105bd964917152be784b46a22d41f9;hb=26591bffd0d706caeed1e3b49f76215dc2411363;hp=026c1b6f389f585a2036ddef2b397c03dbe62ca4;hpb=c6335222f02abf32e0cc697ec8a4666c754655b9;p=unres.git diff --git a/source/unres/src_MD-M/MREMD.F b/source/unres/src_MD-M/MREMD.F index 026c1b6..7a70f1d 100644 --- a/source/unres/src_MD-M/MREMD.F +++ b/source/unres/src_MD-M/MREMD.F @@ -1,3 +1,4 @@ +#define DEBUG subroutine MREMD implicit real*8 (a-h,o-z) include 'DIMENSIONS' @@ -29,7 +30,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+4,maxprocs),t_bath_old,e_tmp integer iremd_acc(maxprocs),iremd_tot(maxprocs) integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs) integer ilen,rstcount @@ -42,6 +43,7 @@ cold integer nup(0:maxprocs),ndown(0:maxprocs) integer icache_all(maxprocs) integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs) logical synflag,end_of_run,file_exist /.false./,ovrtim + real ene_tol /1.0e-5/ cdeb imin_itime_old=0 ntwx_cache=0 @@ -58,9 +60,19 @@ cdeb imin_itime_old=0 endif mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst" -cd print *,'MREMD',nodes +cd print *,'MREMD',nodes,homol_nset cd print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep) cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath + 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) @@ -81,8 +93,9 @@ cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath enddo if(me.eq.king.or..not.out1file) then - write(iout,*) (i2rep(i),i=0,nodes-1) - write(iout,*) (i2set(i),i=0,nodes-1) + write(iout,*) "i2rep",(i2rep(i),i=0,nodes-1) + write(iout,*) "i2set",(i2set(i),i=0,nodes-1) + write(iout,*) "i,j,il,il1,i_index(i,j,il,il1)" do il=1,nset do il1=1,mset(il) do i=1,nrep @@ -177,7 +190,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) then + if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then read (irest2,*) read (irest2,*) nset read (irest2,*) @@ -268,8 +281,13 @@ 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) then + if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then iset=i2set(me) +c broadcast iset to slaves + if (nfgtasks.gt.1) then + call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR) + call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR) + endif if(me.eq.king.or..not.out1file) & write(iout,*) me,"iset=",iset,"t_bath=",t_bath endif @@ -525,7 +543,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) @@ -549,7 +569,7 @@ c Variable time step algorithm. write (irest1,*) ndowna(0,il), & (ndowna(i,il),i=1,ndowna(0,il)) enddo - if(usampl) then + if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then write (irest1,*) "nset" write (irest1,*) nset write (irest1,*) "mset" @@ -576,7 +596,7 @@ c Variable time step algorithm. do i=1,2*nres write (irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo - if(usampl) then + if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then write (irest2,*) iset endif close(irest2) @@ -614,6 +634,7 @@ c REMD - exchange & .or.end_of_run.and.me.eq.king ) & .and. .not. mremdsync ) then synflag=.true. + time01_=MPI_WTIME() do i=1,nodes-1 call mpi_isend(itime,1,MPI_INTEGER,i,101, & CG_COMM, ireqi(i), ierr) @@ -623,7 +644,7 @@ cd call flush(iout) call mpi_waitall(nodes-1,ireqi,statusi,ierr) call mpi_barrier(CG_COMM, ierr) time01=MPI_WTIME() - write(iout,*) 'REMD synchro at',itime,'time=',time01-time00 + write(iout,*) 'REMD synchro at',itime,'time=',time01-time01_ if (out1file.or.traj1file) then cdeb call mpi_gather(itime,1,mpi_integer, cdeb & itime_all,1,mpi_integer,king, @@ -679,7 +700,8 @@ c Update the time safety limiy safety=time001-time00+600 write (iout,*) "****** SAFETY increased to",safety," s" endif - if (ovrtim()) end_of_run=.true. + if (ovrtim() .and. me.eq.king) end_of_run=.true. + call MPI_Bcast(end_of_run,1,MPI_LOGICAL,king,CG_COMM,IERR) endif if(synflag.and..not.end_of_run) then time02=MPI_WTIME() @@ -691,21 +713,68 @@ c call mpi_gather(t_bath,1,mpi_double_precision, c & remd_t_bath,1,mpi_double_precision,king, c & CG_COMM,ierr) potEcomp(n_ene+1)=t_bath - if (usampl) then + t_bath_old=t_bath + if (usampl.or.homol_nset.gt.1) then potEcomp(n_ene+2)=iset 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 +c broadcast iset to slaves and reduce energy + if (nfgtasks.gt.1) then + call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR) + call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR) + call e_modeller(e_tmp) +c write(iout,*) "iset+1 before reduce",e_tmp + call MPI_Barrier(FG_COMM,IERR) + call MPI_Reduce(e_tmp,potEcomp(n_ene+3),1, + & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) + else + call e_modeller(potEcomp(n_ene+3)) + endif +c write(iout,*) "iset+1",potEcomp(n_ene+3) + else + call EconstrQ + potEcomp(n_ene+3)=Uconst + endif iset=i_set_temp +c broadcast iset to slaves + if (nfgtasks.gt.1) then + call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR) + call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR) + endif + else + potEcomp(n_ene+3)=0.0 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 +c broadcast iset to slaves and reduce energy + if (nfgtasks.gt.1) then + call MPI_Bcast(12,1,MPI_INTEGER,king,FG_COMM,IERROR) + call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR) + call e_modeller(e_tmp) +c write(iout,*) "iset-1 before reduce",e_tmp + call MPI_Barrier(FG_COMM,IERR) + call MPI_Reduce(e_tmp,potEcomp(n_ene+4),1, + & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) + else + call e_modeller(potEcomp(n_ene+4)) + endif +c write(iout,*) "iset-1",potEcomp(n_ene+4) + else + call EconstrQ + potEcomp(n_ene+4)=Uconst + endif iset=i_set_temp +c broadcast iset to slaves + if (nfgtasks.gt.1) then + call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR) + call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR) + endif + else + potEcomp(n_ene+4)=0.0 endif endif call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision, @@ -751,10 +820,18 @@ cd end if (me.eq.king) then + if(homol_nset.gt.1) 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) + if(homol_nset.gt.1) + & 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 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh' do i=1,nodes @@ -768,20 +845,27 @@ 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(.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 - call flush(iout) +ctime 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))) +#ifdef DEBUG write (iout,*) "i",i - call flush(iout) +#endif +ctime call flush(iout) 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 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then @@ -824,7 +908,7 @@ c write (iout,*) "ene_i_i",remd_ene(0,i) c call flush(iout) c write (iout,*) "rescaling weights with temperature", c & remd_t_bath(iex) - if (real(ene_i_i).ne.real(remd_ene(0,i))) then + if (abs(ene_i_i-remd_ene(0,i)).gt.ene_tol) then write (iout,*) "ERROR: inconsistent energies:",i, & ene_i_i,remd_ene(0,i) endif @@ -842,7 +926,7 @@ c write (iout,*) "0,iex",remd_t_bath(iex) c call enerprint(remd_ene(0,iex)) call sum_energy(remd_ene(0,iex),.false.) - if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then + if (abs(ene_iex_iex-remd_ene(0,iex)).gt.ene_tol) then write (iout,*) "ERROR: inconsistent energies:",iex, & ene_iex_iex,remd_ene(0,iex) endif @@ -947,9 +1031,9 @@ c call flush(iout) enddo cd write (iout,*) "exchange completed" cd call flush(iout) - ELSE + 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)) @@ -957,10 +1041,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) -cd write(iout,*) "i_dir=",i_dir + endif +c write(iout,*) "i_dir=",i_dir if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then @@ -977,10 +1065,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 @@ -989,17 +1078,17 @@ 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 - call flush(iout) +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 c following temperature changes. @@ -1011,33 +1100,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)) @@ -1052,7 +1147,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) @@ -1088,8 +1183,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 @@ -1117,7 +1212,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, @@ -1145,10 +1240,17 @@ cd call flush(iout) & CG_COMM,ierr) cd write (iout,*) "After scatter" cd call flush(iout) - if(usampl) - & call mpi_scatter(iremd_iset,1,mpi_integer, + if(usampl.or.hremd.gt.0.or.homol_nset.gt.1) then + call mpi_scatter(iremd_iset,1,mpi_integer, & iset,1,mpi_integer,king, & CG_COMM,ierr) +c 8/31/2015 Correction by AL: send new iset to slaves + if (nfgtasks.gt.1) then + call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR) + call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR) + endif + + endif time07=MPI_WTIME() if (me.eq.king .or. .not. out1file) then @@ -1175,8 +1277,8 @@ co & " rescaling weights with temperature",t_bath cde 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 - call flush(iout) + write(iout,*) 'REMD exchange time=',time08-time02 +ctime call flush(iout) endif endif enddo @@ -1229,6 +1331,7 @@ c----------------------------------------------------------------------- include 'COMMON.CHAIN' include 'COMMON.SBRIDGE' include 'COMMON.INTERACT' + include 'COMMON.CONTROL' real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres), & d_restart2(3,2*maxres*maxprocs) @@ -1307,7 +1410,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) @@ -1368,7 +1471,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) @@ -1501,8 +1604,8 @@ c end debugging call xdrfint_(ixdrf, nss, iret) do j=1,nss if (dyn_ss) then - call xdrfint(ixdrf, idssb(j)+nres, iret) - call xdrfint(ixdrf, jdssb(j)+nres, iret) + 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) @@ -1629,6 +1732,7 @@ c end debugging include 'COMMON.CHAIN' include 'COMMON.SBRIDGE' include 'COMMON.INTERACT' + include 'COMMON.CONTROL' real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres), & t5_restart1(5) integer*2 i_index @@ -1716,6 +1820,14 @@ c & (d_restart1(j,i+2*nres*il),j=1,3) #endif enddo enddo +#ifdef DEBUG + write (iout,*) "Conformation read",il + do i=1,nres + write (iout,'(i5,3f10.5,5x,3f10.5)') + & i,(d_restart1(j,i+2*nres*il),j=1,3), + & (d_restart1(j,nres+i+2*nres*il),j=1,3) + enddo +#endif enddo endif call mpi_scatter(d_restart1,3*2*nres,mpi_real, @@ -1750,7 +1862,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) @@ -1792,10 +1904,19 @@ c & (d_restart1(j,i+2*nres*il),j=1,3) enddo endif #endif - call mpi_scatter(i2set,1,mpi_integer, - & iset,1,mpi_integer,king, - & CG_COMM,ierr) - +c 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, +c & CG_COMM,ierr) + call mpi_bcast(i2set(0),nodes,mpi_integer,king, + & CG_COMM,ierr) + iset=i2set(me) +c broadcast iset to slaves + if (nfgtasks.gt.1) then + call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR) + call MPI_Bcast(iset,1,MPI_INTEGER,king,FG_COMM,IERROR) + endif endif @@ -1871,4 +1992,4 @@ c & (d_restart1(j,i+2*nres*il),j=1,3) if(me.eq.king) close(irest2) return end - +c------------------------------------------