+c#define DEBUG
subroutine MREMD
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
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
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
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)
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
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,*)
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
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)
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"
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)
& .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)
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,
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()
synflag=.false.
- write(iout,*) 'REMD before',me,t_bath
+c write(iout,*) 'REMD before',me,t_bath
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,
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
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
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
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
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))
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
- i_dir=iran_num(1,3)
-cd write(iout,*) "i_dir=",i_dir
+ if(t_exchange_only)then
+ i_dir=1
+ else
+ i_dir=iran_num(1,3)
+ endif
+c write(iout,*) "i_dir=",i_dir
if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
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
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.
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))
& 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)
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
& ,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,
& 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
stdfp=dsqrt(2*Rb*t_bath/d_time)
do i=1,ntyp
stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
- enddo
+ enddo
+
+c Compute the standard deviations of stochastic forces for Langevin dynamics
+c if the friction coefficients do not depend on surface area
+ if (lang.gt.0 .and. .not.surfarea) then
+ do i=nnt,nct-1
+ stdforcp(i)=stdfp*dsqrt(gamp)
+ enddo
+ do i=nnt,nct
+ if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
+ & *dsqrt(gamsc(iabs(itype(i))))
+ enddo
+ endif
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
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)
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)
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)
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)
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
#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,
enddo
- if(usampl) then
+ if(usampl.or.homol_nset.gt.1) then
#ifdef AIX
if(me.eq.king)then
call xdrfint_(ixdrf, nset, iret)
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
if(me.eq.king) close(irest2)
return
end
-
+c------------------------------------------
projects / unres.git / blobdiff