#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(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 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.and.hremd.eq.0) THEN write (iout,*) "Enter exchnge, remd_m",remd_m(1), & " nodes",nodes ctime 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 ctime 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 (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 ' 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