X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2FMREMD.F;h=afecaa5ff6aa6b6aa0a5e60254f3bdb25d3d7b55;hb=34d3ad3987785642be58fb2f26557d3314215577;hp=1df3f0a9810911dae70dc9f6d5d61a0b78cf14f2;hpb=f690e8b70bab14132839afebf080d4a28363b226;p=unres.git diff --git a/source/unres/src_MD-M/MREMD.F b/source/unres/src_MD-M/MREMD.F index 1df3f0a..afecaa5 100644 --- a/source/unres/src_MD-M/MREMD.F +++ b/source/unres/src_MD-M/MREMD.F @@ -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------------------------------------------