--- /dev/null
+ 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'
+ include 'COMMON.HAIRPIN'
+ 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)
+ 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
+
+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
+
+ 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) 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()-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
+ 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(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 hairpin(.true.,nharp,iharp)
+ call secondary2(.true.)
+ 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 .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.
+ 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) then
+c Update the time safety limiy
+ if (time001-time00.gt.safety) then
+ safety=time001-time00+600
+ 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.
+
+ 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
+ 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
+c-------------------------------------
+ IF(.not.usampl) THEN
+ write (iout,*) "Enter exchnge, remd_m",remd_m(1),
+ & " nodes",nodes
+ call flush(iout)
+ write (iout,*) "remd_m(1)",remd_m(1)
+ do irr=1,remd_m(1)
+ i=ifirst(iran_num(1,remd_m(1)))
+ write (iout,*) "i",i
+ call flush(iout)
+
+ do ii=1,nodes-1
+
+ write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
+ if(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 (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))
+
+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 (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
+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)
+ 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()-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 '
+ 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
+ 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
+C do j=1,nss
+ if (dyn_ss) then
+ call xdrfint(ixdrf, idssb(j)+nres, iret)
+ call xdrfint(ixdrf, jdssb(j)+nres, iret)
+ else
+ call xdrfint_(ixdrf, ihpb(j), iret)
+ call xdrfint_(ixdrf, jhpb(j), iret)
+ endif
+ enddo
+ 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
+ if (dyn_ss) then
+ call xdrfint(ixdrf, idssb(j)+nres, iret)
+ call xdrfint(ixdrf, jdssb(j)+nres, iret)
+ else
+ call xdrfint(ixdrf, ihpb(j), iret)
+ call xdrfint(ixdrf, jhpb(j), iret)
+ endif
+ enddo
+ call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
+ call xdrfint(ixdrf, iset_restart1(il), iret)
+ 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
+
projects / unres.git / blobdiff