X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2FMREMD.F;h=f22e2f6acc0f08f35e39e28652a5e9266d8a9712;hb=57038e4bdff4cc9534106b25bfbd4b9a844d47fd;hp=38db8a836ca5b6ce497a1978a98518c978a83043;hpb=29dca221ada246dbf1a26c9e0e4fe12ee1e6b978;p=unres.git diff --git a/source/unres/src-HCD-5D/MREMD.F b/source/unres/src-HCD-5D/MREMD.F index 38db8a8..f22e2f6 100644 --- a/source/unres/src-HCD-5D/MREMD.F +++ b/source/unres/src-HCD-5D/MREMD.F @@ -57,11 +57,12 @@ cold integer nup(0:maxprocs),ndown(0:maxprocs) integer rep2i(0:maxprocs),ireqi(maxprocs) integer icache_all(maxprocs) integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs) - logical synflag,end_of_run,file_exist /.false./,ovrtim + logical synflag,end_of_run,file_exist /.false./,ovrtim,first_pass double precision t_bath_temp,delta,ene_iex_iex,ene_i_i,ene_iex_i, & ene_i_iex,xxx,tmp,econstr_temp_iex,econstr_temp_i integer iran_num double precision ran_number + integer i_econstr/20/ cdeb imin_itime_old=0 ntwx_cache=0 @@ -144,8 +145,8 @@ cold endif cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0)) - write (*,*) "Processor",me," rest",rest," - & restart1fie",restart1file +c write (*,*) "Processor",me," rest",rest," +c & restart1fie",restart1file if(rest.and.restart1file) then if (me.eq.king) & inquire(file=mremd_rst_name,exist=file_exist) @@ -186,8 +187,9 @@ cd write (*,*) me," After broadcast: file_exist",file_exist do i=nnt,nct if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i))) & *dsqrt(gamsc(iabs(itype(i)))) - enddo + enddo endif + endif if(me.eq.king) then @@ -250,20 +252,20 @@ cd write (*,*) me," After broadcast: file_exist",file_exist write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":", & (ndowna(i,il),i=1,ndowna(0,il)) enddo + stdfp=dsqrt(2*Rb*t_bath/d_time) do i=1,ntyp stdfsc(i)=dsqrt(2*Rb*t_bath/d_time) enddo if (lang.gt.0 .and. .not.surfarea) then - do i=nnt,nct-1 - stdforcp(i)=stdfp*dsqrt(gamp) - enddo + 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 - + enddo + endif ELSE IF (.not.(rest.and.file_exist)) THEN do il=1,remd_m(1) ifirst(il)=il @@ -319,7 +321,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) +c call flush(iout) endif c stdfp=dsqrt(2*Rb*t_bath/d_time) @@ -422,6 +424,8 @@ c Entering the MD loop #endif itime=0 end_of_run=.false. + first_pass=.not.rest +c write (iout,*) "first_pass",first_pass do while(.not.end_of_run) itime=itime+1 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true. @@ -721,7 +725,7 @@ ctime call flush(iout) write(iout,*) "MIN ii_write=",ii_write endif endif - call flush(iout) +c call flush(iout) endif if (synflag) then c Update the time safety limiy @@ -855,7 +859,7 @@ cd end time05=MPI_WTIME() if (me.eq.king .or. .not. out1file) then write(iout,*) 'REMD writing traj time=',time05-time04 - call flush(iout) +c call flush(iout) endif @@ -940,7 +944,10 @@ 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 +c write (iout,*) "first_pass",first_pass + if (.not.first_pass.and. + & real(ene_i_i).ne.real(remd_ene(0,i))) + & then write (iout,*) "ERROR: inconsistent energies:",i, & ene_i_i,remd_ene(0,i) endif @@ -958,7 +965,8 @@ 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 (.not.first_pass.and. + & real(ene_iex_iex).ne.real(remd_ene(0,iex))) then write (iout,*) "ERROR: inconsistent energies:",iex, & ene_iex_iex,remd_ene(0,iex) endif @@ -1063,6 +1071,7 @@ c call flush(iout) endif enddo enddo + first_pass=.false. cd write (iout,*) "exchange completed" cd call flush(iout) ELSE @@ -1089,12 +1098,12 @@ cd write(iout,*) "i_dir=",i_dir 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) + econstr_temp_i=remd_ene(i_econstr,i) + econstr_temp_iex=remd_ene(i_econstr,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) + remd_ene(i_econstr,i)=remd_ene(n_ene+5,i) + remd_ene(i_econstr,iex)=remd_ene(n_ene+6,iex) endif elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then @@ -1103,10 +1112,10 @@ c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials) 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 @@ -1115,14 +1124,14 @@ c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials) 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) + econstr_temp_i=remd_ene(i_econstr,i) + econstr_temp_iex=remd_ene(i_econstr,iex) if (adaptive) then - remd_ene(20,i)=remd_ene(n_ene+7,i) - remd_ene(20,iex)=remd_ene(n_ene+8,iex) + remd_ene(i_econstr,i)=remd_ene(n_ene+7,i) + remd_ene(i_econstr,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) + remd_ene(i_econstr,i)=remd_ene(n_ene+3,i) + remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex) endif else goto 444 @@ -1222,8 +1231,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 @@ -1259,7 +1268,7 @@ c------------------------------------- enddo endif - call flush(iout) +c call flush(iout) cd write (iout,'(a6,100i4)') "ifirst", cd & (ifirst(i),i=1,remd_m(1)) @@ -1305,7 +1314,6 @@ co & " rescaling weights with temperature",t_bath do i=1,ntyp stdfsc(i)=dsqrt(2*Rb*t_bath/d_time) 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 @@ -1317,11 +1325,11 @@ c if the friction coefficients do not depend on surface area & *dsqrt(gamsc(iabs(itype(i)))) enddo endif -cde write(iout,*) 'REMD after',me,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) +c call flush(iout) endif endif enddo @@ -1382,7 +1390,7 @@ c----------------------------------------------------------------------- include 'COMMON.SBRIDGE' include 'COMMON.INTERACT' - real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres), + real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1), & d_restart2(3,2*maxres*maxprocs) real t5_restart1(5) integer iret,itmp @@ -1402,7 +1410,7 @@ c----------------------------------------------------------------------- & t_restart1,5,mpi_real,king,CG_COMM,ierr) - do i=1,2*nres + do i=0,2*nres-1 do j=1,3 r_d(j,i)=d_t(j,i) enddo @@ -1412,7 +1420,7 @@ c----------------------------------------------------------------------- & CG_COMM,ierr) - do i=1,2*nres + do i=0,2*nres-1 do j=1,3 r_d(j,i)=dc(j,i) enddo @@ -1793,14 +1801,14 @@ c end debugging include 'COMMON.CHAIN' include 'COMMON.SBRIDGE' include 'COMMON.INTERACT' - real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres), + real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1), & t5_restart1(5) integer*2 i_index & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200) common /przechowalnia/ d_restart1 integer i,j,il,il1,ixdrf,iret,itmp integer ierr - write (*,*) "Processor",me," called read1restart" +c write (*,*) "Processor",me," called read1restart" if(me.eq.king)then open(irest2,file=mremd_rst_name,status='unknown') @@ -1887,7 +1895,7 @@ c & (d_restart1(j,i+2*nres*il),j=1,3) call mpi_scatter(d_restart1,3*2*nres,mpi_real, & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr) - do i=1,2*nres + do i=0,2*nres-1 do j=1,3 d_t(j,i)=r_d(j,i) enddo @@ -1909,7 +1917,7 @@ c & (d_restart1(j,i+2*nres*il),j=1,3) endif call mpi_scatter(d_restart1,3*2*nres,mpi_real, & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr) - do i=1,2*nres + do i=0,2*nres-1 do j=1,3 dc(j,i)=r_d(j,i) enddo @@ -1990,7 +1998,7 @@ c & CG_COMM,ierr) include 'COMMON.CHAIN' include 'COMMON.SBRIDGE' include 'COMMON.INTERACT' - real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres), + real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1), & t5_restart1(5) common /przechowalnia/ d_restart1 integer i,j,il,itmp @@ -2026,7 +2034,7 @@ c & CG_COMM,ierr) call mpi_scatter(d_restart1,3*2*nres,mpi_real, & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr) - do i=1,2*nres + do i=0,2*nres-1 do j=1,3 d_t(j,i)=r_d(j,i) enddo @@ -2041,7 +2049,7 @@ c & CG_COMM,ierr) endif call mpi_scatter(d_restart1,3*2*nres,mpi_real, & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr) - do i=1,2*nres + do i=0,2*nres-1 do j=1,3 dc(j,i)=r_d(j,i) enddo