X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-NEWSC%2FMREMD.F;fp=source%2Funres%2Fsrc_MD-NEWSC%2FMREMD.F;h=576e43dda288d438ed3d57ce0ba1837bad607749;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_MD-NEWSC/MREMD.F b/source/unres/src_MD-NEWSC/MREMD.F new file mode 100644 index 0000000..576e43d --- /dev/null +++ b/source/unres/src_MD-NEWSC/MREMD.F @@ -0,0 +1,2102 @@ +#ifdef MPI + subroutine MREMD + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' +#ifndef LANG0 + include 'COMMON.LANGEVIN' +#else + include 'COMMON.LANGEVIN.lang0' +#endif + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + include 'COMMON.TIME1' + include 'COMMON.REMD' + include 'COMMON.SETUP' + include 'COMMON.MUCA' + integer ERRCODE + double precision cm(3),L(3),vcm(3) + double precision energia(0:n_ene) + double precision remd_t_bath(maxprocs) + 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),t_bath_old + integer iremd_acc(maxprocs),iremd_tot(maxprocs) + integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs) + integer ilen,rstcount + external ilen + character*50 tytul + common /gucio/ cm + integer itime +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 + real ene_tol /1.0e-5/ + +cdeb imin_itime_old=0 + ntwx_cache=0 + time00=MPI_WTIME() + time01=time00 + if(me.eq.king.or..not.out1file) then + write (iout,*) 'MREMD',nodes,'time before',time00-walltime + write (iout,*) "NREP=",nrep + endif + + synflag=.false. + if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then + call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst") + endif + mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst" + +cd print *,'MREMD',nodes +cd print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep) +cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath + + if(hremd.gt.0) then + nset=hremd + do i=1,nset + mset(i)=1 + enddo + endif + + k=0 + rep2i(k)=-1 + do il=1,max0(nset,1) + do il1=1,max0(mset(il),1) + do i=1,nrep + iremd_acc(i)=0 + iremd_acc_usa(i)=0 + iremd_tot(i)=0 + do j=1,remd_m(i) + i2rep(k)=i + i2set(k)=il + rep2i(i)=k + k=k+1 + i_index(i,j,il,il1)=k + enddo + enddo + enddo + enddo + + if(me.eq.king.or..not.out1file) then + 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 + do j=1,remd_m(i) + write(iout,*) i,j,il,il1,i_index(i,j,il,il1) + enddo + enddo + enddo + enddo + endif + +c print *,'i2rep',me,i2rep(me) +c print *,'rep2i',(rep2i(i),i=0,nrep) + +cold if (i2rep(me).eq.nrep) then +cold nup(0)=0 +cold else +cold nup(0)=remd_m(i2rep(me)+1) +cold k=rep2i(int(i2rep(me)))+1 +cold do i=1,nup(0) +cold nup(i)=k +cold k=k+1 +cold enddo +cold endif + +cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0)) + +cold if (i2rep(me).eq.1) then +cold ndown(0)=0 +cold else +cold ndown(0)=remd_m(i2rep(me)-1) +cold k=rep2i(i2rep(me)-2)+1 +cold do i=1,ndown(0) +cold ndown(i)=k +cold k=k+1 +cold enddo +cold endif + +cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0)) + + + write (*,*) "Processor",me," rest",rest," + & restart1fie",restart1file + if(rest.and.restart1file) then + if (me.eq.king) + & inquire(file=mremd_rst_name,exist=file_exist) +cd write (*,*) me," Before broadcast: file_exist",file_exist + call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM, + & IERR) +cd write (*,*) me," After broadcast: file_exist",file_exist + if(file_exist) then + if(me.eq.king.or..not.out1file) + & write (iout,*) 'Master is reading restart1file' + call read1restart(i_index) + else + if(me.eq.king.or..not.out1file) + & write (iout,*) 'WARNING : no restart1file' + endif + + if(me.eq.king.or..not.out1file) then + write(iout,*) "i2set",(i2set(i),i=0,nodes-1) + write(iout,*) "i_index" + do il=1,nset + do il1=1,mset(il) + do i=1,nrep + do j=1,remd_m(i) + write(iout,*) i,j,il,il1,i_index(i,j,il,il1) + enddo + enddo + enddo + enddo + endif + endif + + if(me.eq.king) then + if (rest.and..not.restart1file) + & inquire(file=mremd_rst_name,exist=file_exist) + if(.not.file_exist.and.rest.and..not.restart1file) + & write(iout,*) 'WARNING : no restart file',mremd_rst_name + IF (rest.and.file_exist.and..not.restart1file) THEN + write (iout,*) 'Master is reading restart file', + & mremd_rst_name + open(irest2,file=mremd_rst_name,status='unknown') + read (irest2,*) + read (irest2,*) (i2rep(i),i=0,nodes-1) + read (irest2,*) + read (irest2,*) (ifirst(i),i=1,remd_m(1)) + do il=1,nodes + read (irest2,*) + read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il)) + read (irest2,*) + read (irest2,*) ndowna(0,il), + & (ndowna(i,il),i=1,ndowna(0,il)) + enddo + if(usampl.or.hremd.gt.0) then + read (irest2,*) + read (irest2,*) nset + read (irest2,*) + read (irest2,*) (mset(i),i=1,nset) + read (irest2,*) + read (irest2,*) (i2set(i),i=0,nodes-1) + read (irest2,*) + do il=1,nset + do il1=1,mset(il) + do i=1,nrep + read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i)) + enddo + enddo + enddo + + write(iout,*) "i2set",(i2set(i),i=0,nodes-1) + write(iout,*) "i_index" + do il=1,nset + do il1=1,mset(il) + do i=1,nrep + do j=1,remd_m(i) + write(iout,*) i,j,il,il1,i_index(i,j,il,il1) + enddo + enddo + enddo + enddo + endif + + close(irest2) + + write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1) + write (iout,'(a6,1000i5)') "ifirst", + & (ifirst(i),i=1,remd_m(1)) + do il=1,nodes + write (iout,'(a6,i4,a1,100i4)') "nupa",il,":", + & (nupa(i,il),i=1,nupa(0,il)) + write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":", + & (ndowna(i,il),i=1,ndowna(0,il)) + enddo + ELSE IF (.not.(rest.and.file_exist)) THEN + do il=1,remd_m(1) + ifirst(il)=il + enddo + + do il=1,nodes + if (i2rep(il-1).eq.nrep) then + nupa(0,il)=0 + else + nupa(0,il)=remd_m(i2rep(il-1)+1) + k=rep2i(int(i2rep(il-1)))+1 + do i=1,nupa(0,il) + nupa(i,il)=k+1 + k=k+1 + enddo + endif + if (i2rep(il-1).eq.1) then + ndowna(0,il)=0 + else + ndowna(0,il)=remd_m(i2rep(il-1)-1) + k=rep2i(i2rep(il-1)-2)+1 + do i=1,ndowna(0,il) + ndowna(i,il)=k+1 + k=k+1 + enddo + endif + enddo + + write (iout,'(a6,100i4)') "ifirst", + & (ifirst(i),i=1,remd_m(1)) + do il=1,nodes + write (iout,'(a6,i4,a1,100i4)') "nupa",il,":", + & (nupa(i,il),i=1,nupa(0,il)) + write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":", + & (ndowna(i,il),i=1,ndowna(0,il)) + enddo + + ENDIF + endif +c +c t_bath=retmin+(retmax-retmin)*me/(nodes-1) + if(.not.(rest.and.file_exist.and.restart1file)) then + if (me .eq. king) then + t_bath=retmin + else + t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep)) + endif +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.or.hremd.gt.0) then + iset=i2set(me) + if (hremd.gt.0) call set_hweights(iset) + if(me.eq.king.or..not.out1file) + & write(iout,*) me,"iset=",iset,"t_bath=",t_bath + endif +c + stdfp=dsqrt(2*Rb*t_bath/d_time) + do i=1,ntyp + stdfsc(i)=dsqrt(2*Rb*t_bath/d_time) + enddo + +c print *,'irep',me,t_bath + if (.not.rest) then + if (me.eq.king .or. .not. out1file) + & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath + call rescale_weights(t_bath) + endif + + +c------copy MD-------------- +c The driver for molecular dynamics subroutines +c------------------------------------------------ + t_MDsetup=0.0d0 + t_langsetup=0.0d0 + t_MD=0.0d0 + t_enegrad=0.0d0 + t_sdsetup=0.0d0 + if(me.eq.king.or..not.out1file) + & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started" +#ifdef MPI + tt0 = MPI_Wtime() +#else + tt0 = tcpu() +#endif +c Determine the inverse of the inertia matrix. + call setup_MD_matrices +c Initialize MD + call init_MD + if (rest) then + if (me.eq.king .or. .not. out1file) + & write (iout,'(a60,f10.5)') "REMD restart 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 + if (lang.gt.0 .and. .not.surfarea) then + do i=nnt,nct-1 + stdforcp(i)=stdfp*dsqrt(gamp) + enddo + do i=nnt,nct + stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i))) + enddo + elseif (lang.gt.0 .and. surfarea ) then + call setup_fricmat + endif + call rescale_weights(t_bath) + endif + +#ifdef MPI + t_MDsetup = MPI_Wtime()-tt0 +#else + t_MDsetup = tcpu()-tt0 +#endif + rstcount=0 +c Entering the MD loop +#ifdef MPI + tt0 = MPI_Wtime() +#else + tt0 = tcpu() +#endif + if (lang.eq.2 .or. lang.eq.3) then +#ifndef LANG0 + call setup_fricmat + if (lang.eq.2) then + call sd_verlet_p_setup + else + call sd_verlet_ciccotti_setup + endif + do i=1,dimen + do j=1,dimen + pfric0_mat(i,j,0)=pfric_mat(i,j) + afric0_mat(i,j,0)=afric_mat(i,j) + vfric0_mat(i,j,0)=vfric_mat(i,j) + prand0_mat(i,j,0)=prand_mat(i,j) + vrand0_mat1(i,j,0)=vrand_mat1(i,j) + vrand0_mat2(i,j,0)=vrand_mat2(i,j) + enddo + enddo + flag_stoch(0)=.true. + do i=1,maxflag_stoch + flag_stoch(i)=.false. + enddo +#else + write (iout,*) + & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0" +#ifdef MPI + call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) +#endif + stop +#endif + else if (lang.eq.1 .or. lang.eq.4) then + call setup_fricmat + endif + time00=MPI_WTIME() + if (me.eq.king .or. .not. out1file) + & write(iout,*) 'Setup time',time00-walltime +ctime call flush(iout) +#ifdef MPI + t_langsetup=MPI_Wtime()-tt0 + tt0=MPI_Wtime() +#else + t_langsetup=tcpu()-tt0 + tt0=tcpu() +#endif + itime=0 + end_of_run=.false. + do while(.not.end_of_run) + itime=itime+1 + if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true. + if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true. + rstcount=rstcount+1 + if (lang.gt.0 .and. surfarea .and. + & mod(itime,reset_fricmat).eq.0) then + if (lang.eq.2 .or. lang.eq.3) then +#ifndef LANG0 + call setup_fricmat + if (lang.eq.2) then + call sd_verlet_p_setup + else + call sd_verlet_ciccotti_setup + endif + do i=1,dimen + do j=1,dimen + pfric0_mat(i,j,0)=pfric_mat(i,j) + afric0_mat(i,j,0)=afric_mat(i,j) + vfric0_mat(i,j,0)=vfric_mat(i,j) + prand0_mat(i,j,0)=prand_mat(i,j) + vrand0_mat1(i,j,0)=vrand_mat1(i,j) + vrand0_mat2(i,j,0)=vrand_mat2(i,j) + enddo + enddo + flag_stoch(0)=.true. + do i=1,maxflag_stoch + flag_stoch(i)=.false. + enddo +#endif + else if (lang.eq.1 .or. lang.eq.4) then + call setup_fricmat + endif + write (iout,'(a,i10)') + & "Friction matrix reset based on surface area, itime",itime + endif + if (reset_vel .and. tbf .and. lang.eq.0 + & .and. mod(itime,count_reset_vel).eq.0) then + call random_vel + if (me.eq.king .or. .not. out1file) + & write(iout,'(a,f20.2)') + & "Velocities reset to random values, time",totT + do i=0,2*nres + do j=1,3 + d_t_old(j,i)=d_t(j,i) + enddo + enddo + endif + if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then + call inertia_tensor + call vcm_vel(vcm) + do j=1,3 + d_t(j,0)=d_t(j,0)-vcm(j) + enddo + call kinetic(EK) + kinetic_T=2.0d0/(dimen3*Rb)*EK + scalfac=dsqrt(T_bath/kinetic_T) +cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT + do i=0,2*nres + do j=1,3 + d_t_old(j,i)=scalfac*d_t(j,i) + enddo + enddo + endif + if (lang.ne.4) then + if (RESPA) then +c Time-reversible RESPA algorithm +c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992) + call RESPA_step(itime) + else +c Variable time step algorithm. + call velverlet_step(itime) + endif + else +#ifdef BROWN + call brown_step(itime) +#else + print *,"Brown dynamics not here!" +#ifdef MPI + call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) +#endif + stop +#endif + endif + if(hmc.gt.0 .and. mod(itime,hmc).eq.0) then + call statout(itime) + call hmc_test(itime) + endif + if(ntwe.ne.0) then + if (mod(itime,ntwe).eq.0) call statout(itime) + endif + if (mod(itime,ntwx).eq.0.and..not.traj1file) then + write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath + if(mdpdb) then + call pdbout(potE,tytul,ipdb) + else + call cartout(totT) + endif + endif + if (mod(itime,ntwx).eq.0.and.traj1file) then + if(ntwx_cache.lt.max_cache_traj_use) then + ntwx_cache=ntwx_cache+1 + else + if (max_cache_traj_use.ne.1) + & print *,itime,"processor ",me," over cache ",ntwx_cache + do i=1,ntwx_cache-1 + + totT_cache(i)=totT_cache(i+1) + EK_cache(i)=EK_cache(i+1) + potE_cache(i)=potE_cache(i+1) + t_bath_cache(i)=t_bath_cache(i+1) + Uconst_cache(i)=Uconst_cache(i+1) + iset_cache(i)=iset_cache(i+1) + + do ii=1,nfrag + qfrag_cache(ii,i)=qfrag_cache(ii,i+1) + enddo + do ii=1,npair + qpair_cache(ii,i)=qpair_cache(ii,i+1) + enddo + do ii=1,nfrag_back + utheta_cache(ii,i)=utheta_cache(ii,i+1) + ugamma_cache(ii,i)=ugamma_cache(ii,i+1) + uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1) + enddo + + + do ii=1,nres*2 + do j=1,3 + c_cache(j,ii,i)=c_cache(j,ii,i+1) + enddo + enddo + enddo + endif + + totT_cache(ntwx_cache)=totT + EK_cache(ntwx_cache)=EK + potE_cache(ntwx_cache)=potE + t_bath_cache(ntwx_cache)=t_bath + Uconst_cache(ntwx_cache)=Uconst + iset_cache(ntwx_cache)=iset + + do i=1,nfrag + qfrag_cache(i,ntwx_cache)=qfrag(i) + enddo + do i=1,npair + qpair_cache(i,ntwx_cache)=qpair(i) + enddo + do i=1,nfrag_back + utheta_cache(i,ntwx_cache)=utheta(i) + ugamma_cache(i,ntwx_cache)=ugamma(i) + uscdiff_cache(i,ntwx_cache)=uscdiff(i) + enddo + + do i=1,nres*2 + do j=1,3 + c_cache(j,i,ntwx_cache)=c(j,i) + enddo + enddo + + endif + if ((rstcount.eq.1000.or.itime.eq.n_timestep) + & .and..not.restart1file) then + + if(me.eq.king) then + open(irest1,file=mremd_rst_name,status='unknown') + write (irest1,*) "i2rep" + write (irest1,*) (i2rep(i),i=0,nodes-1) + write (irest1,*) "ifirst" + write (irest1,*) (ifirst(i),i=1,remd_m(1)) + do il=1,nodes + write (irest1,*) "nupa",il + write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il)) + write (irest1,*) "ndowna",il + write (irest1,*) ndowna(0,il), + & (ndowna(i,il),i=1,ndowna(0,il)) + enddo + if(usampl.or.hremd.gt.0) then + write (irest1,*) "nset" + write (irest1,*) nset + write (irest1,*) "mset" + write (irest1,*) (mset(i),i=1,nset) + write (irest1,*) "i2set" + write (irest1,*) (i2set(i),i=0,nodes-1) + write (irest1,*) "i_index" + do il=1,nset + do il1=1,mset(il) + do i=1,nrep + write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i)) + enddo + enddo + enddo + + endif + close(irest1) + endif + open(irest2,file=rest2name,status='unknown') + write(irest2,*) totT,EK,potE,totE,t_bath + do i=1,2*nres + write (irest2,'(3e15.5)') (d_t(j,i),j=1,3) + enddo + do i=1,2*nres + write (irest2,'(3e15.5)') (dc(j,i),j=1,3) + enddo + if(usampl.or.hremd.gt.0) then + write (irest2,*) iset + endif + close(irest2) + rstcount=0 + endif + +c REMD - exchange +c forced synchronization + if (mod(itime,i_sync_step).eq.0 .and. me.ne.king + & .and. .not. mremdsync) then + synflag=.false. + call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr) + if (synflag) then + call mpi_recv(itime_master, 1, MPI_INTEGER, + & 0,101,CG_COMM, status, ierr) + call mpi_barrier(CG_COMM, ierr) +cdeb if (out1file.or.traj1file) then +cdeb call mpi_gather(itime,1,mpi_integer, +cdeb & icache_all,1,mpi_integer,king, +cdeb & CG_COMM,ierr) + if(traj1file) + & call mpi_gather(ntwx_cache,1,mpi_integer, + & icache_all,1,mpi_integer,king, + & CG_COMM,ierr) + if (.not.out1file) + & write(iout,*) 'REMD synchro at',itime_master,itime + if (itime_master.ge.n_timestep .or. ovrtim()) + & end_of_run=.true. +ctime call flush(iout) + endif + endif + +c REMD - exchange + if ((mod(itime,nstex).eq.0.and.me.eq.king + & .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) +cd write(iout,*) 'REMD synchro with',i +cd call flush(iout) + enddo + 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-time01_ + if (out1file.or.traj1file) then +cdeb call mpi_gather(itime,1,mpi_integer, +cdeb & itime_all,1,mpi_integer,king, +cdeb & CG_COMM,ierr) +cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime', +cdeb & (itime_all(i),i=1,nodes) + if(traj1file) then +cdeb imin_itime=itime_all(1) +cdeb do i=2,nodes +cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i) +cdeb enddo +cdeb ii_write=(imin_itime-imin_itime_old)/ntwx +cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx +cdeb write(iout,*) imin_itime,imin_itime_old,ii_write + call mpi_gather(ntwx_cache,1,mpi_integer, + & icache_all,1,mpi_integer,king, + & CG_COMM,ierr) +c write(iout,'(a19,8000i8)') ' ntwx_cache', +c & (icache_all(i),i=1,nodes) + ii_write=icache_all(1) + do i=2,nodes + if(icache_all(i).lt.ii_write) ii_write=icache_all(i) + enddo +c write(iout,*) "MIN ii_write=",ii_write + endif + endif +ctime call flush(iout) + endif + if(mremdsync .and. mod(itime,nstex).eq.0) then + synflag=.true. + if (me.eq.king .or. .not. out1file) + & write(iout,*) 'REMD synchro at',itime + + if(traj1file) then + call mpi_gather(ntwx_cache,1,mpi_integer, + & icache_all,1,mpi_integer,king, + & CG_COMM,ierr) + if (me.eq.king) then + write(iout,'(a19,8000i8)') ' ntwx_cache', + & (icache_all(i),i=1,nodes) + ii_write=icache_all(1) + do i=2,nodes + if(icache_all(i).lt.ii_write) ii_write=icache_all(i) + enddo + write(iout,*) "MIN ii_write=",ii_write + endif + endif +ctest call flush(iout) + endif + if (synflag) then +c Update the time safety limiy + if (time001-time00.gt.safety) then + safety=time001-time00+600 + if (me.eq.king .or. .not. out1file) + & write (iout,*) "****** SAFETY increased to",safety," s" + endif + if (ovrtim()) end_of_run=.true. + endif + if(synflag.and..not.end_of_run) then + time02=MPI_WTIME() + synflag=.false. + +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 + t_bath_old=t_bath + if (usampl) 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 + iset=i_set_temp + endif + if (iset.gt.1) then + i_set_temp=iset + iset=iset-1 + call EconstrQ + potEcomp(n_ene+4)=Uconst + iset=i_set_temp + endif + endif + if(hremd.gt.0) potEcomp(n_ene+2)=iset + call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision, + & remd_ene(0,1),n_ene+5,mpi_double_precision,king, + & CG_COMM,ierr) + if(lmuca) then + call mpi_gather(elow,1,mpi_double_precision, + & elowi,1,mpi_double_precision,king, + & CG_COMM,ierr) + call mpi_gather(ehigh,1,mpi_double_precision, + & ehighi,1,mpi_double_precision,king, + & CG_COMM,ierr) + endif + + time03=MPI_WTIME() + if (me.eq.king .or. .not. out1file) then + write(iout,*) 'REMD gather times=',time03-time01 + & ,time03-time02 + endif + + if (restart1file) call write1rst(i_index) + + time04=MPI_WTIME() + if (me.eq.king .or. .not. out1file) then + write(iout,*) 'REMD writing rst time=',time04-time03 + endif + + if (traj1file) call write1traj +cd debugging +cdeb call mpi_gather(ntwx_cache,1,mpi_integer, +cdeb & icache_all,1,mpi_integer,king, +cdeb & CG_COMM,ierr) +cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file', +cdeb & (icache_all(i),i=1,nodes) +cd end + + + time05=MPI_WTIME() + if (me.eq.king .or. .not. out1file) then + write(iout,*) 'REMD writing traj time=',time05-time04 +ctime call flush(iout) + endif + + + if (me.eq.king) then + do i=1,nodes + remd_t_bath(i)=remd_ene(n_ene+1,i) + iremd_iset(i)=remd_ene(n_ene+2,i) + enddo +#ifdef DEBUG + if(lmuca) then +co write(iout,*) 'REMD exchange temp,ene,elow,ehigh' + do i=1,nodes + write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i), + & elowi(i),ehighi(i) + enddo + else + write(iout,*) 'REMD exchange temp,ene' + do i=1,nodes + write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i) + write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene) + enddo + endif +#endif +c------------------------------------- + IF(.not.usampl.and.hremd.eq.0) THEN +#ifdef DEBUG + write (iout,*) "Enter exchnge, remd_m",remd_m(1), + & " nodes",nodes +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 +#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 write (iout,*) +c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD" +c call flush(iout) +c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr) +c endif +c do while (iex.eq.i) +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 + if (lmuca) then + call muca_delta(remd_t_bath,remd_ene,i,iex,delta) + else +c Swap temperatures between conformations i and iex with recalculating the free energies +c following temperature changes. + ene_iex_iex=remd_ene(0,iex) + ene_i_i=remd_ene(0,i) +c write (iout,*) "i",i," ene_i_i",ene_i_i, +c & " iex",iex," ene_iex_iex",ene_iex_iex +c write (iout,*) "rescaling weights with temperature", +c & remd_t_bath(i) +c call flush(iout) + call rescale_weights(remd_t_bath(i)) + +c write (iout,*) "0,iex",remd_t_bath(i) +c call enerprint(remd_ene(0,iex)) + + call sum_energy(remd_ene(0,iex),.false.) + ene_iex_i=remd_ene(0,iex) +c write (iout,*) "ene_iex_i",remd_ene(0,iex) + +c write (iout,*) "0,i",remd_t_bath(i) +c call enerprint(remd_ene(0,i)) + + call sum_energy(remd_ene(0,i),.false.) +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 (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 + call rescale_weights(remd_t_bath(iex)) + +c write (iout,*) "0,i",remd_t_bath(iex) +c call enerprint(remd_ene(0,i)) + + call sum_energy(remd_ene(0,i),.false.) +c write (iout,*) "ene_i_iex",remd_ene(0,i) +c call flush(iout) + ene_i_iex=remd_ene(0,i) + +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 (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 +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) + 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 +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)) + endif + if (delta .gt. 50.0d0) then + delta=0.0d0 + else +#ifdef OSF + if(isnan(delta))then + delta=0.0d0 + else if (delta.lt.-50.0d0) then + delta=dexp(50.0d0) + else + delta=dexp(-delta) + endif +#else + delta=dexp(-delta) +#endif + endif + iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1 + xxx=ran_number(0.0d0,1.0d0) +c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx +c call flush(iout) + if (delta .gt. xxx) then + tmp=remd_t_bath(i) + remd_t_bath(i)=remd_t_bath(iex) + remd_t_bath(iex)=tmp + remd_ene(0,i)=ene_i_iex + remd_ene(0,iex)=ene_iex_i + if(lmuca) then + tmp=elowi(i) + elowi(i)=elowi(iex) + elowi(iex)=tmp + tmp=ehighi(i) + ehighi(i)=ehighi(iex) + ehighi(iex)=tmp + endif + + + do k=0,nodes + itmp=nupa(k,i) + nupa(k,i)=nupa(k,iex) + nupa(k,iex)=itmp + itmp=ndowna(k,i) + ndowna(k,i)=ndowna(k,iex) + ndowna(k,iex)=itmp + enddo + do il=1,nodes + if (ifirst(il).eq.i) ifirst(il)=iex + do k=1,nupa(0,il) + if (nupa(k,il).eq.i) then + nupa(k,il)=iex + elseif (nupa(k,il).eq.iex) then + nupa(k,il)=i + endif + enddo + do k=1,ndowna(0,il) + if (ndowna(k,il).eq.i) then + ndowna(k,il)=iex + elseif (ndowna(k,il).eq.iex) then + ndowna(k,il)=i + endif + enddo + enddo + + iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1 + 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) + + else + remd_ene(0,iex)=ene_iex_iex + remd_ene(0,i)=ene_i_i + i=iex + endif + endif + enddo + enddo +cd write (iout,*) "exchange completed" +cd call flush(iout) + ELSEIF (usampl) THEN + do ii=1,nodes +cd write(iout,*) "########",ii + + i_temp=iran_num(1,nrep) + i_mult=iran_num(1,remd_m(i_temp)) + i_iset=iran_num(1,nset) + 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 + + if(t_exchange_only)then + i_dir=1 + else + i_dir=iran_num(1,3) + endif +cd write(iout,*) "i_dir=",i_dir + + if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then + + i_temp1=i_temp+1 + i_mult1=iran_num(1,remd_m(i_temp1)) + i_iset1=i_iset + i_mset1=iran_num(1,mset(i_iset1)) + iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1) + + elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then + + i_temp1=i_temp + i_mult1=iran_num(1,remd_m(i_temp1)) + 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) + + elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then + + i_temp1=i_temp+1 + i_mult1=iran_num(1,remd_m(i_temp1)) + 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) + + else + goto 444 + endif + +cd 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. + ene_iex_iex=remd_ene(0,iex) + ene_i_i=remd_ene(0,i) +co write (iout,*) "rescaling weights with temperature", +co & remd_t_bath(i) + call rescale_weights(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) +c call sum_energy(remd_ene(0,i),.false.) +cd 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 & ene_i_i,remd_ene(0,i) +c endif + 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) +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 & ene_iex_iex,remd_ene(0,iex) +c endif +cd 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 + 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 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)) + if (delta .gt. 50.0d0) then + delta=0.0d0 + else + delta=dexp(-delta) + endif + if (i_dir.eq.1.or.i_dir.eq.3) + & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1 + if (i_dir.eq.2.or.i_dir.eq.3) + & 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 + if (delta .gt. xxx) then + tmp=remd_t_bath(i) + remd_t_bath(i)=remd_t_bath(iex) + remd_t_bath(iex)=tmp + + itmp=iremd_iset(i) + iremd_iset(i)=iremd_iset(iex) + iremd_iset(iex)=itmp + + remd_ene(0,i)=ene_i_iex + remd_ene(0,iex)=ene_iex_i + + if (i_dir.eq.1.or.i_dir.eq.3) + & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1 + + itmp=i2rep(i-1) + i2rep(i-1)=i2rep(iex-1) + i2rep(iex-1)=itmp + + if (i_dir.eq.2.or.i_dir.eq.3) + & iremd_acc_usa(int(i2set(i-1)))= + & iremd_acc_usa(int(i2set(i-1)))+1 + + itmp=i2set(i-1) + i2set(i-1)=i2set(iex-1) + i2set(iex-1)=itmp + + itmp=i_index(i_temp,i_mult,i_iset,i_mset) + i_index(i_temp,i_mult,i_iset,i_mset)= + & i_index(i_temp1,i_mult1,i_iset1,i_mset1) + i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp + + 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 + endif + +cd do il=1,nset +cd do il1=1,mset(il) +cd do i=1,nrep +cd do j=1,remd_m(i) +cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1) +cd enddo +cd enddo +cd enddo +cd enddo + + 444 continue + + enddo + + ELSEIF (hremd.gt.0) THEN + do ii=1,nodes +cd write(iout,*) "########",ii + + i_temp=iran_num(1,nrep) + i_mult=iran_num(1,remd_m(i_temp)) + i_iset=iran_num(1,nset) + i_mset=1 + i=i_index(i_temp,i_mult,i_iset,i_mset) + +cd 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) + endif + +cd write(iout,*) "i_dir=",i_dir + + if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then + + i_temp1=i_temp+1 + i_mult1=iran_num(1,remd_m(i_temp1)) + i_iset1=i_iset + i_mset1=1 + iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1) + + elseif(i_dir.eq.2)then + + i_temp1=i_temp + i_mult1=iran_num(1,remd_m(i_temp1)) + i_iset1=iran_num(1,hremd) + do while(i_iset1.eq.i_iset) + i_iset1=iran_num(1,hremd) + enddo + i_mset1=1 + iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1) + + elseif(remd_m(i_temp+1).gt.0)then + + i_temp1=i_temp+1 + i_mult1=iran_num(1,remd_m(i_temp1)) + i_iset1=iran_num(1,hremd) + do while(i_iset1.eq.i_iset) + i_iset1=iran_num(1,hremd) + enddo + i_mset1=1 + iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1) + + else + goto 445 + endif + +cd 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. + ene_iex_iex=remd_ene(0,iex) + ene_i_i=remd_ene(0,i) + + call set_hweights(i_iset) + call rescale_weights(remd_t_bath(i)) + call sum_energy(remd_ene(0,iex),.false.) + ene_iex_i=remd_ene(0,iex) + + call set_hweights(i_iset1) + call rescale_weights(remd_t_bath(iex)) + call sum_energy(remd_ene(0,i),.false.) + ene_i_iex=remd_ene(0,i) + +cd write(iout,*) ene_iex_iex,ene_i_i,ene_iex_i,ene_i_iex + + 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 + + if (delta .gt. 50.0d0) then + delta=0.0d0 + else + delta=dexp(-delta) + endif + + if (i_dir.eq.1.or.i_dir.eq.3) + & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1 + if (i_dir.eq.2.or.i_dir.eq.3) + & 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 + if (delta .gt. xxx) then + +cd write (iout,*) "exchange" + tmp=remd_t_bath(i) + remd_t_bath(i)=remd_t_bath(iex) + remd_t_bath(iex)=tmp + + itmp=iremd_iset(i) + iremd_iset(i)=iremd_iset(iex) + iremd_iset(iex)=itmp + + remd_ene(0,i)=ene_i_iex + remd_ene(0,iex)=ene_iex_i + + if (i_dir.eq.1.or.i_dir.eq.3) + & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1 + + itmp=i2rep(i-1) + i2rep(i-1)=i2rep(iex-1) + i2rep(iex-1)=itmp + + if (i_dir.eq.2.or.i_dir.eq.3) + & iremd_acc_usa(int(i2set(i-1)))= + & iremd_acc_usa(int(i2set(i-1)))+1 + + itmp=i2set(i-1) + i2set(i-1)=i2set(iex-1) + i2set(iex-1)=itmp + + itmp=i_index(i_temp,i_mult,i_iset,i_mset) + i_index(i_temp,i_mult,i_iset,i_mset)= + & i_index(i_temp1,i_mult1,i_iset1,i_mset1) + i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp + +cd do il=1,nset +cd do il1=1,mset(il) +cd do i=1,nrep +cd do j=1,remd_m(i) +cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1) +cd enddo +cd enddo +cd enddo +cd enddo + + else + remd_ene(0,iex)=ene_iex_iex + remd_ene(0,i)=ene_i_i + endif + + + + 445 continue + + enddo + + ENDIF + +c------------------------------------- + write (iout,*) "NREP",nrep + do i=1,nrep + if(iremd_tot(i).ne.0) + & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) + & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i) + enddo + + if(usampl) then + do i=1,nset + if(iremd_tot_usa(i).ne.0) + & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i, + & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i) + enddo + endif + + if(hremd.gt.0) then + do i=1,nset + if(iremd_tot_usa(i).ne.0) + & write(iout,'(a10,i4,f12.5,i8)') 'ACC_hremd',i, + & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i) + enddo + endif + + +ctime call flush(iout) + +cd write (iout,'(a6,100i4)') "ifirst", +cd & (ifirst(i),i=1,remd_m(1)) +cd do il=1,nodes +cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":", +cd & (nupa(i,il),i=1,nupa(0,il)) +cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":", +cd & (ndowna(i,il),i=1,ndowna(0,il)) +cd enddo + endif + + time06=MPI_WTIME() +cd write (iout,*) "Before scatter" +cd call flush(iout) + call mpi_scatter(remd_t_bath,1,mpi_double_precision, + & t_bath,1,mpi_double_precision,king, + & CG_COMM,ierr) +cd write (iout,*) "After scatter" +cd call flush(iout) + if(usampl.or.hremd.gt.0) + & call mpi_scatter(iremd_iset,1,mpi_integer, + & iset,1,mpi_integer,king, + & CG_COMM,ierr) + + time07=MPI_WTIME() + if (me.eq.king .or. .not. out1file) then + write(iout,*) 'REMD scatter time=',time07-time06 + endif + + if(lmuca) then + call mpi_scatter(elowi,1,mpi_double_precision, + & elow,1,mpi_double_precision,king, + & CG_COMM,ierr) + call mpi_scatter(ehighi,1,mpi_double_precision, + & ehigh,1,mpi_double_precision,king, + & CG_COMM,ierr) + endif + + if(hremd.gt.0) call set_hweights(iset) + call rescale_weights(t_bath) +co write (iout,*) "Processor",me, +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 + if (lang.gt.0) then + do i=nnt,nct-1 + stdforcp(i)=stdforcp(i)*sqrt(t_bath/t_bath_old) + enddo + do i=nnt,nct + stdforcsc(i)=stdforcsc(i)*sqrt(t_bath/t_bath_old) + 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-time02 +ctime call flush(iout) + endif + endif + enddo + + if (restart1file) then + if (me.eq.king .or. .not. out1file) + & write(iout,*) 'writing restart at the end of run' + call write1rst(i_index) + endif + + if (traj1file) call write1traj +cd debugging +cdeb call mpi_gather(ntwx_cache,1,mpi_integer, +cdeb & icache_all,1,mpi_integer,king, +cdeb & CG_COMM,ierr) +cdeb write(iout,'(a40,8000i8)') +cdeb & ' ntwx_cache after traj1file at the end', +cdeb & (icache_all(i),i=1,nodes) +cd end + + +#ifdef MPI + t_MD=MPI_Wtime()-tt0 +#else + t_MD=tcpu()-tt0 +#endif + if (me.eq.king .or. .not. out1file) then + write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') + & ' Timing ', + & 'MD calculations setup:',t_MDsetup, + & 'Energy & gradient evaluation:',t_enegrad, + & 'Stochastic MD setup:',t_langsetup, + & 'Stochastic MD step setup:',t_sdsetup, + & 'MD steps:',t_MD + write (iout,'(/28(1h=),a25,27(1h=))') + & ' End of MD calculation ' + if(hmc.gt.0) write (iout,*) 'HMC acceptance ratio', + & n_timestep*1.0d0/hmc/hmc_acc + endif + return + end + +c----------------------------------------------------------------------- + subroutine write1rst(i_index) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.MD' + include 'COMMON.IOUNITS' + include 'COMMON.REMD' + include 'COMMON.SETUP' + include 'COMMON.CHAIN' + include 'COMMON.SBRIDGE' + include 'COMMON.INTERACT' + + real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres), + & d_restart2(3,2*maxres*maxprocs) + real t5_restart1(5) + integer iret,itmp + integer*2 i_index + & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200) + common /przechowalnia/ d_restart1,d_restart2 + + t5_restart1(1)=totT + t5_restart1(2)=EK + t5_restart1(3)=potE + t5_restart1(4)=t_bath + t5_restart1(5)=Uconst + + call mpi_gather(t5_restart1,5,mpi_real, + & t_restart1,5,mpi_real,king,CG_COMM,ierr) + + + do i=1,2*nres + do j=1,3 + r_d(j,i)=d_t(j,i) + enddo + enddo + call mpi_gather(r_d,3*2*nres,mpi_real, + & d_restart1,3*2*nres,mpi_real,king, + & CG_COMM,ierr) + + + do i=1,2*nres + do j=1,3 + r_d(j,i)=dc(j,i) + enddo + enddo + call mpi_gather(r_d,3*2*nres,mpi_real, + & d_restart2,3*2*nres,mpi_real,king, + & CG_COMM,ierr) + + if(me.eq.king) then +#ifdef AIX + call xdrfopen_(ixdrf,mremd_rst_name, "w", iret) + do i=0,nodes-1 + call xdrfint_(ixdrf, i2rep(i), iret) + enddo + do i=1,remd_m(1) + call xdrfint_(ixdrf, ifirst(i), iret) + enddo + do il=1,nodes + do i=0,nupa(0,il) + call xdrfint_(ixdrf, nupa(i,il), iret) + enddo + + do i=0,ndowna(0,il) + call xdrfint_(ixdrf, ndowna(i,il), iret) + enddo + enddo + + do il=1,nodes + do j=1,4 + call xdrffloat_(ixdrf, t_restart1(j,il), iret) + enddo + enddo + + do il=0,nodes-1 + do i=1,2*nres + do j=1,3 + call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret) + enddo + enddo + enddo + do il=0,nodes-1 + do i=1,2*nres + do j=1,3 + call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret) + enddo + enddo + enddo + + if(usampl) then + call xdrfint_(ixdrf, nset, iret) + do i=1,nset + call xdrfint_(ixdrf,mset(i), iret) + enddo + do i=0,nodes-1 + call xdrfint_(ixdrf,i2set(i), iret) + enddo + do il=1,nset + do il1=1,mset(il) + do i=1,nrep + do j=1,remd_m(i) + itmp=i_index(i,j,il,il1) + call xdrfint_(ixdrf,itmp, iret) + enddo + enddo + enddo + enddo + + endif + call xdrfclose_(ixdrf, iret) +#else + call xdrfopen(ixdrf,mremd_rst_name, "w", iret) + do i=0,nodes-1 + call xdrfint(ixdrf, i2rep(i), iret) + enddo + do i=1,remd_m(1) + call xdrfint(ixdrf, ifirst(i), iret) + enddo + do il=1,nodes + do i=0,nupa(0,il) + call xdrfint(ixdrf, nupa(i,il), iret) + enddo + + do i=0,ndowna(0,il) + call xdrfint(ixdrf, ndowna(i,il), iret) + enddo + enddo + + do il=1,nodes + do j=1,4 + call xdrffloat(ixdrf, t_restart1(j,il), iret) + enddo + enddo + + do il=0,nodes-1 + do i=1,2*nres + do j=1,3 + call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret) + enddo + enddo + enddo + do il=0,nodes-1 + do i=1,2*nres + do j=1,3 + call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret) + enddo + enddo + enddo + + + if(usampl) then + call xdrfint(ixdrf, nset, iret) + do i=1,nset + call xdrfint(ixdrf,mset(i), iret) + enddo + do i=0,nodes-1 + call xdrfint(ixdrf,i2set(i), iret) + enddo + do il=1,nset + do il1=1,mset(il) + do i=1,nrep + do j=1,remd_m(i) + itmp=i_index(i,j,il,il1) + call xdrfint(ixdrf,itmp, iret) + enddo + enddo + enddo + enddo + + endif + call xdrfclose(ixdrf, iret) +#endif + endif + return + end + + + subroutine write1traj + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.MD' + include 'COMMON.IOUNITS' + include 'COMMON.REMD' + include 'COMMON.SETUP' + include 'COMMON.CHAIN' + include 'COMMON.SBRIDGE' + include 'COMMON.INTERACT' + + real t5_restart1(5) + integer iret,itmp + real xcoord(3,maxres2+2),prec + real r_qfrag(50),r_qpair(100) + real r_utheta(50),r_ugamma(100),r_uscdiff(100) + real p_qfrag(50*maxprocs),p_qpair(100*maxprocs) + real p_utheta(50*maxprocs),p_ugamma(100*maxprocs), + & p_uscdiff(100*maxprocs) + real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2) + common /przechowalnia/ p_c + + call mpi_bcast(ii_write,1,mpi_integer, + & king,CG_COMM,ierr) + +c debugging + print *,'traj1file',me,ii_write,ntwx_cache +c end debugging + +#ifdef AIX + if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret) +#else + if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret) +#endif + do ii=1,ii_write + t5_restart1(1)=totT_cache(ii) + t5_restart1(2)=EK_cache(ii) + t5_restart1(3)=potE_cache(ii) + t5_restart1(4)=t_bath_cache(ii) + t5_restart1(5)=Uconst_cache(ii) + call mpi_gather(t5_restart1,5,mpi_real, + & t_restart1,5,mpi_real,king,CG_COMM,ierr) + + call mpi_gather(iset_cache(ii),1,mpi_integer, + & iset_restart1,1,mpi_integer,king,CG_COMM,ierr) + + do i=1,nfrag + r_qfrag(i)=qfrag_cache(i,ii) + enddo + do i=1,npair + r_qpair(i)=qpair_cache(i,ii) + enddo + do i=1,nfrag_back + r_utheta(i)=utheta_cache(i,ii) + r_ugamma(i)=ugamma_cache(i,ii) + r_uscdiff(i)=uscdiff_cache(i,ii) + enddo + + call mpi_gather(r_qfrag,nfrag,mpi_real, + & p_qfrag,nfrag,mpi_real,king, + & CG_COMM,ierr) + call mpi_gather(r_qpair,npair,mpi_real, + & p_qpair,npair,mpi_real,king, + & CG_COMM,ierr) + call mpi_gather(r_utheta,nfrag_back,mpi_real, + & p_utheta,nfrag_back,mpi_real,king, + & CG_COMM,ierr) + call mpi_gather(r_ugamma,nfrag_back,mpi_real, + & p_ugamma,nfrag_back,mpi_real,king, + & CG_COMM,ierr) + call mpi_gather(r_uscdiff,nfrag_back,mpi_real, + & p_uscdiff,nfrag_back,mpi_real,king, + & CG_COMM,ierr) + +#ifdef DEBUG + write (iout,*) "p_qfrag" + do i=1,nodes + write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag) + enddo + write (iout,*) "p_qpair" + do i=1,nodes + write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair) + enddo +ctime call flush(iout) +#endif + do i=1,nres*2 + do j=1,3 + r_c(j,i)=c_cache(j,i,ii) + enddo + enddo + + call mpi_gather(r_c,3*2*nres,mpi_real, + & p_c,3*2*nres,mpi_real,king, + & CG_COMM,ierr) + + if(me.eq.king) then +#ifdef AIX + do il=1,nodes + call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret) + call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret) + call xdrffloat_(ixdrf, real(t_restart1(5,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) + enddo + call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret) + call xdrfint_(ixdrf, iset_restart1(il), iret) + do i=1,nfrag + call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret) + enddo + do i=1,npair + call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret) + enddo + do i=1,nfrag_back + call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret) + call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret) + call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret) + enddo + prec=10000.0 + do i=1,nres + do j=1,3 + xcoord(j,i)=p_c(j,i+(il-1)*nres*2) + enddo + enddo + do i=nnt,nct + do j=1,3 + xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2) + enddo + enddo + itmp=nres+nct-nnt+1 + call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret) + enddo +#else + do il=1,nodes + call xdrffloat(ixdrf, real(t_restart1(1,il)), iret) + call xdrffloat(ixdrf, real(t_restart1(3,il)), iret) + call xdrffloat(ixdrf, real(t_restart1(5,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) + enddo + call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret) + call xdrfint(ixdrf, iset_restart1(il), iret) + do i=1,nfrag + call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret) + enddo + do i=1,npair + call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret) + enddo + do i=1,nfrag_back + call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret) + call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret) + call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret) + enddo + prec=10000.0 + do i=1,nres + do j=1,3 + xcoord(j,i)=p_c(j,i+(il-1)*nres*2) + enddo + enddo + do i=nnt,nct + do j=1,3 + xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2) + enddo + enddo + itmp=nres+nct-nnt+1 + call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret) + enddo +#endif + endif + enddo +#ifdef AIX + if(me.eq.king) call xdrfclose_(ixdrf, iret) +#else + if(me.eq.king) call xdrfclose(ixdrf, iret) +#endif + do i=1,ntwx_cache-ii_write + + totT_cache(i)=totT_cache(ii_write+i) + EK_cache(i)=EK_cache(ii_write+i) + potE_cache(i)=potE_cache(ii_write+i) + t_bath_cache(i)=t_bath_cache(ii_write+i) + Uconst_cache(i)=Uconst_cache(ii_write+i) + iset_cache(i)=iset_cache(ii_write+i) + + do ii=1,nfrag + qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i) + enddo + do ii=1,npair + qpair_cache(ii,i)=qpair_cache(ii,ii_write+i) + enddo + do ii=1,nfrag_back + utheta_cache(ii,i)=utheta_cache(ii,ii_write+i) + ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i) + uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i) + enddo + + do ii=1,nres*2 + do j=1,3 + c_cache(j,ii,i)=c_cache(j,ii,ii_write+i) + enddo + enddo + enddo + ntwx_cache=ntwx_cache-ii_write + return + end + + + subroutine read1restart(i_index) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.MD' + include 'COMMON.IOUNITS' + include 'COMMON.REMD' + include 'COMMON.SETUP' + include 'COMMON.CHAIN' + include 'COMMON.SBRIDGE' + include 'COMMON.INTERACT' + real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres), + & t5_restart1(5) + integer*2 i_index + & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200) + common /przechowalnia/ d_restart1 + write (*,*) "Processor",me," called read1restart" + + if(me.eq.king)then + open(irest2,file=mremd_rst_name,status='unknown') + read(irest2,*,err=334) i + write(iout,*) "Reading old rst in ASCI format" + close(irest2) + call read1restart_old + return + 334 continue +#ifdef AIX + call xdrfopen_(ixdrf,mremd_rst_name, "r", iret) + + do i=0,nodes-1 + call xdrfint_(ixdrf, i2rep(i), iret) + enddo + do i=1,remd_m(1) + call xdrfint_(ixdrf, ifirst(i), iret) + enddo + do il=1,nodes + call xdrfint_(ixdrf, nupa(0,il), iret) + do i=1,nupa(0,il) + call xdrfint_(ixdrf, nupa(i,il), iret) + enddo + + call xdrfint_(ixdrf, ndowna(0,il), iret) + do i=1,ndowna(0,il) + call xdrfint_(ixdrf, ndowna(i,il), iret) + enddo + enddo + do il=1,nodes + do j=1,4 + call xdrffloat_(ixdrf, t_restart1(j,il), iret) + enddo + enddo +#else + call xdrfopen(ixdrf,mremd_rst_name, "r", iret) + + do i=0,nodes-1 + call xdrfint(ixdrf, i2rep(i), iret) + enddo + do i=1,remd_m(1) + call xdrfint(ixdrf, ifirst(i), iret) + enddo + do il=1,nodes + call xdrfint(ixdrf, nupa(0,il), iret) + do i=1,nupa(0,il) + call xdrfint(ixdrf, nupa(i,il), iret) + enddo + + call xdrfint(ixdrf, ndowna(0,il), iret) + do i=1,ndowna(0,il) + call xdrfint(ixdrf, ndowna(i,il), iret) + enddo + enddo + do il=1,nodes + do j=1,4 + call xdrffloat(ixdrf, t_restart1(j,il), iret) + enddo + enddo +#endif + endif + call mpi_scatter(t_restart1,5,mpi_real, + & t5_restart1,5,mpi_real,king,CG_COMM,ierr) + totT=t5_restart1(1) + EK=t5_restart1(2) + potE=t5_restart1(3) + t_bath=t5_restart1(4) + + if(me.eq.king)then + do il=0,nodes-1 + do i=1,2*nres +c read(irest2,'(3e15.5)') +c & (d_restart1(j,i+2*nres*il),j=1,3) + do j=1,3 +#ifdef AIX + call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret) +#else + call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret) +#endif + enddo + enddo + enddo + 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 j=1,3 + d_t(j,i)=r_d(j,i) + enddo + enddo + if(me.eq.king)then + do il=0,nodes-1 + do i=1,2*nres +c read(irest2,'(3e15.5)') +c & (d_restart1(j,i+2*nres*il),j=1,3) + do j=1,3 +#ifdef AIX + call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret) +#else + call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret) +#endif + enddo + enddo + enddo + 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 j=1,3 + dc(j,i)=r_d(j,i) + enddo + enddo + + + if(usampl) then +#ifdef AIX + if(me.eq.king)then + call xdrfint_(ixdrf, nset, iret) + do i=1,nset + call xdrfint_(ixdrf,mset(i), iret) + enddo + do i=0,nodes-1 + call xdrfint_(ixdrf,i2set(i), iret) + enddo + do il=1,nset + do il1=1,mset(il) + do i=1,nrep + do j=1,remd_m(i) + call xdrfint_(ixdrf,itmp, iret) + i_index(i,j,il,il1)=itmp + enddo + enddo + enddo + enddo + endif +#else + if(me.eq.king)then + call xdrfint(ixdrf, nset, iret) + do i=1,nset + call xdrfint(ixdrf,mset(i), iret) + enddo + do i=0,nodes-1 + call xdrfint(ixdrf,i2set(i), iret) + enddo + do il=1,nset + do il1=1,mset(il) + do i=1,nrep + do j=1,remd_m(i) + call xdrfint(ixdrf,itmp, iret) + i_index(i,j,il,il1)=itmp + enddo + enddo + enddo + enddo + endif +#endif + call mpi_scatter(i2set,1,mpi_integer, + & iset,1,mpi_integer,king, + & CG_COMM,ierr) + + endif + + + if(me.eq.king) close(irest2) + return + end + + subroutine read1restart_old + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.MD' + include 'COMMON.IOUNITS' + include 'COMMON.REMD' + include 'COMMON.SETUP' + include 'COMMON.CHAIN' + include 'COMMON.SBRIDGE' + include 'COMMON.INTERACT' + real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres), + & t5_restart1(5) + common /przechowalnia/ d_restart1 + if(me.eq.king)then + open(irest2,file=mremd_rst_name,status='unknown') + read (irest2,*) (i2rep(i),i=0,nodes-1) + read (irest2,*) (ifirst(i),i=1,remd_m(1)) + do il=1,nodes + read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il)) + read (irest2,*) ndowna(0,il), + & (ndowna(i,il),i=1,ndowna(0,il)) + enddo + do il=1,nodes + read(irest2,*) (t_restart1(j,il),j=1,4) + enddo + endif + call mpi_scatter(t_restart1,5,mpi_real, + & t5_restart1,5,mpi_real,king,CG_COMM,ierr) + totT=t5_restart1(1) + EK=t5_restart1(2) + potE=t5_restart1(3) + t_bath=t5_restart1(4) + + if(me.eq.king)then + do il=0,nodes-1 + do i=1,2*nres + read(irest2,'(3e15.5)') + & (d_restart1(j,i+2*nres*il),j=1,3) + enddo + enddo + 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 j=1,3 + d_t(j,i)=r_d(j,i) + enddo + enddo + if(me.eq.king)then + do il=0,nodes-1 + do i=1,2*nres + read(irest2,'(3e15.5)') + & (d_restart1(j,i+2*nres*il),j=1,3) + enddo + enddo + 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 j=1,3 + dc(j,i)=r_d(j,i) + enddo + enddo + if(me.eq.king) close(irest2) + return + end +c------------------------------------------------------------------- + subroutine set_hweights(iiset) + implicit real*8 (a-h,o-z) + integer i + include 'DIMENSIONS' + include 'COMMON.FFIELD' + include 'COMMON.REMD' + + do i=1,n_ene + weights(i)=hweights(iiset,i) + enddo + + wsc =weights(1) + wscp =weights(2) + welec =weights(3) + wcorr =weights(4) + wcorr5 =weights(5) + wcorr6 =weights(6) + wel_loc=weights(7) + wturn3 =weights(8) + wturn4 =weights(9) + wturn6 =weights(10) + wang =weights(11) + wscloc =weights(12) + wtor =weights(13) + wtor_d =weights(14) + wstrain=weights(15) + wvdwpp =weights(16) + wbond =weights(17) + scal14 =weights(18) + wsccor =weights(21) + + return + end +#endif