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
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)
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)
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,
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
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
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
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
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
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
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
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
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",
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()
& 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
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)
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)
enddo
endif
#endif
- call mpi_scatter(i2set,1,mpi_integer,
- & iset,1,mpi_integer,king,
- & CG_COMM,ierr)
+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)
endif
if(me.eq.king) close(irest2)
return
end
-
+c------------------------------------------