X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fold_F%2FMREMD.F.safe;fp=source%2Funres%2Fsrc_MD%2Fold_F%2FMREMD.F.safe;h=0000000000000000000000000000000000000000;hb=0a11a2c4ccee14ed99ae44f2565b270ba8d4bbb6;hp=110dea3f18d31ce41032f4a61de5501a4713c184;hpb=5eb407964903815242c59de10960f42761139e10;p=unres.git diff --git a/source/unres/src_MD/old_F/MREMD.F.safe b/source/unres/src_MD/old_F/MREMD.F.safe deleted file mode 100644 index 110dea3..0000000 --- a/source/unres/src_MD/old_F/MREMD.F.safe +++ /dev/null @@ -1,1756 +0,0 @@ - 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' - 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) - 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./ - -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 - 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(i),i=0,nodes-1) - write(iout,*) (i2set(i),i=0,nodes-1) - 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) 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 - -c write (iout,'(a6,100i4)') "ifirst", -c & (ifirst(i),i=1,remd_m(1)) -c do il=1,nodes -c write (iout,'(a6,i4,a1,100i4)') "nupa",il,":", -c & (nupa(i,il),i=1,nupa(0,il)) -c write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":", -c & (ndowna(i,il),i=1,ndowna(0,il)) -c 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) then - iset=i2set(me) - 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 - call rescale_weights(t_bath) - endif - -#ifdef MPI - t_MDsetup = MPI_Wtime() -#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 - call setup_fricmat - if (lang.eq.2) then - call sd_verlet_p_setup - else - call sd_verlet_ciccotti_setup - endif -#ifndef LANG0 - 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 -#endif - flag_stoch(0)=.true. - do i=1,maxflag_stoch - flag_stoch(i)=.false. - enddo - 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 - call flush(iout) -#ifdef MPI - t_langsetup=MPI_Wtime() - 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 - call setup_fricmat - if (lang.eq.2) then - call sd_verlet_p_setup - else - call sd_verlet_ciccotti_setup - endif -#ifndef LANG0 - 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 -#endif - flag_stoch(0)=.true. - do i=1,maxflag_stoch - flag_stoch(i)=.false. - enddo - 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/(dimen*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 - call brown_step(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) 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) 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) 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. - time00=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-time00 - 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 - - call flush(iout) - endif - if(synflag.and..not.end_of_run) then - time02=MPI_WTIME() - synflag=.false. - -cd 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 - 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 - 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 - 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 - 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 -cd write(iout,*) 'REMD exchange temp,ene' -c do i=1,nodes -co write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i) -c write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene) -c enddo - endif -c------------------------------------- - IF(.not.usampl) THEN - write (iout,*) "Enter exchnge, remd_m",remd_m(1), - & " nodes",nodes - call flush(iout) - do irr=1,remd_m(1) - i=ifirst(iran_num(1,remd_m(1))) - call flush(iout) - - do ii=1,nodes-1 - - if(i.gt.0.and.nupa(0,i).gt.0) then - iex=nupa(iran_num(1,int(nupa(0,i))),i) - 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) -cd write (iout,*) "rescaling weights with temperature", -cd & remd_t_bath(i) -cd call flush(iout) - call rescale_weights(remd_t_bath(i)) - -cd write (iout,*) "0,iex",remd_t_bath(i) -cd call enerprint(remd_ene(0,iex)) - - 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) - -cd write (iout,*) "0,i",remd_t_bath(i) -cd call enerprint(remd_ene(0,i)) - - call sum_energy(remd_ene(0,i),.false.) -cd write (iout,*) "ene_i_i",remd_ene(0,i) -cd call flush(iout) -cd write (iout,*) "rescaling weights with temperature", -cd & remd_t_bath(iex) - if (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 - call rescale_weights(remd_t_bath(iex)) - -cd write (iout,*) "0,i",remd_t_bath(iex) -cd call enerprint(remd_ene(0,i)) - - call sum_energy(remd_ene(0,i),.false.) -cd write (iout,*) "ene_i_iex",remd_ene(0,i) -cd call flush(iout) - ene_i_iex=remd_ene(0,i) - -cd write (iout,*) "0,iex",remd_t_bath(iex) -cd 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 - write (iout,*) "ERROR: inconsistent energies:",iex, - & ene_iex_iex,remd_ene(0,iex) - endif -cd write (iout,*) "ene_iex_iex",remd_ene(0,iex) -cd 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 -cd 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 -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)) - 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) -cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx -cd 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 - -cd write(iout,*) 'exchange',i,iex -cd write (iout,'(a8,100i4)') "@ ifirst", -cd & (ifirst(k),k=1,remd_m(1)) -cd do il=1,nodes -cd write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":", -cd & (nupa(k,il),k=1,nupa(0,il)) -cd write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":", -cd & (ndowna(k,il),k=1,ndowna(0,il)) -cd enddo - 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) - ELSE - 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 - - 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_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 - 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 - - - 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 - - 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) - & 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 - 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 - -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) - 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() -#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 ' - endif - return - end - -c----------------------------------------------------------------------- - - subroutine write1rst_oldtxt - 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 - - - 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 - open(irest1,file=mremd_rst_name,status='unknown') - write (irest1,*) (i2rep(i),i=0,nodes-1) - write (irest1,*) (ifirst(i),i=1,remd_m(1)) - do il=1,nodes - write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il)) - write (irest1,*) ndowna(0,il), - & (ndowna(i,il),i=1,ndowna(0,il)) - enddo - - do il=1,nodes - write(irest1,*) (t_restart1(j,il),j=1,4) - enddo - - do il=0,nodes-1 - do i=1,2*nres - write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3) - enddo - enddo - do il=0,nodes-1 - do i=1,2*nres - write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3) - enddo - enddo - close(irest1) - endif - - return - end - - 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) - - - 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 -c open(irest1,file=mremd_rst_name,status='unknown') - call xdrfopen(ixdrf,mremd_rst_name, "w", iret) -c write (irest1,*) (i2rep(i),i=0,nodes-1) - do i=0,nodes-1 - call xdrfint(ixdrf, i2rep(i), iret) - enddo -c write (irest1,*) (ifirst(i),i=1,remd_m(1)) - do i=1,remd_m(1) - call xdrfint(ixdrf, ifirst(i), iret) - enddo - do il=1,nodes -c write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il)) - do i=0,nupa(0,il) - call xdrfint(ixdrf, nupa(i,il), iret) - enddo - -c write (irest1,*) ndowna(0,il), -c & (ndowna(i,il),i=1,ndowna(0,il)) - do i=0,ndowna(0,il) - call xdrfint(ixdrf, ndowna(i,il), iret) - enddo - enddo - - do il=1,nodes -c write(irest1,*) (t_restart1(j,il),j=1,4) - do j=1,4 - call xdrffloat(ixdrf, t_restart1(j,il), iret) - enddo - enddo - - do il=0,nodes-1 - do i=1,2*nres -c write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3) - 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 -c write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3) - do j=1,3 - call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret) - enddo - enddo - enddo - - - if(usampl) then -c write (irest1,*) nset - call xdrfint(ixdrf, nset, iret) -c write (irest1,*) (mset(i),i=1,nset) - do i=1,nset - call xdrfint(ixdrf,mset(i), iret) - enddo -c write (irest1,*) (i2set(i),i=0,nodes-1) - do i=0,nodes-1 - call xdrfint(ixdrf,i2set(i), iret) - enddo -c write (irest1,*) "i_index" - do il=1,nset - do il1=1,mset(il) - do i=1,nrep -c write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i)) - do j=1,remd_m(i) - itmp=i_index(i,j,il,il1) - call xdrfint(ixdrf,itmp, iret) - enddo - enddo - enddo - enddo - - endif - -c close(irest1) - call xdrfclose(ixdrf, iret) - 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) - - call mpi_bcast(ii_write,1,mpi_integer, - & king,CG_COMM,ierr) - -c debugging - print *,'traj1file',me,ii_write,ntwx_cache -c end debugging - - if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret) - 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 - 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 - - 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 - enddo - if(me.eq.king) call xdrfclose(ixdrf, iret) - 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) - 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 - call xdrfopen(ixdrf,mremd_rst_name, "r", iret) - -c read (irest2,*) (i2rep(i),i=0,nodes-1) - do i=0,nodes-1 - call xdrfint(ixdrf, i2rep(i), iret) - enddo -c read (irest2,*) (ifirst(i),i=1,remd_m(1)) - do i=1,remd_m(1) - call xdrfint(ixdrf, ifirst(i), iret) - enddo - do il=1,nodes -c read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il)) - call xdrfint(ixdrf, nupa(0,il), iret) - do i=1,nupa(0,il) - call xdrfint(ixdrf, nupa(i,il), iret) - enddo - -c read (irest2,*) ndowna(0,il), -c & (ndowna(i,il),i=1,ndowna(0,il)) - 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 -c read(irest2,*) (t_restart1(j,il),j=1,4) - do j=1,4 - call xdrffloat(ixdrf, t_restart1(j,il), iret) - enddo - 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 -c read(irest2,'(3e15.5)') -c & (d_restart1(j,i+2*nres*il),j=1,3) - do j=1,3 - call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret) - 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 - call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret) - 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 - 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 - - 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) - - 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 -