module MREMDyn !----------------------------------------------------------------------------- use io_units use names use MPI_data use md_data use remd_data use geometry_data use energy_data ! use control_data, only:maxprocs use MDyn implicit none !----------------------------------------------------------------------------- ! commom.remd ! common /remdrestart/ integer(kind=2),dimension(:),allocatable :: i2set !(0:maxprocs) integer(kind=2),dimension(:),allocatable :: ifirst !(maxprocs) integer(kind=2),dimension(:,:),allocatable :: nupa,& ndowna !(0:maxprocs/4,0:maxprocs) real(kind=4),dimension(:,:),allocatable :: t_restart1 !(5,maxprocs) integer,dimension(:),allocatable :: iset_restart1 !(maxprocs) ! common /traj1cache/ real(kind=4),dimension(:),allocatable :: totT_cache,EK_cache,& potE_cache,t_bath_cache,Uconst_cache !(max_cache_traj) real(kind=4),dimension(:,:),allocatable :: qfrag_cache !(50,max_cache_traj) real(kind=4),dimension(:,:),allocatable :: qpair_cache !(100,max_cache_traj) real(kind=4),dimension(:,:),allocatable :: ugamma_cache,& utheta_cache,uscdiff_cache !(maxfrag_back,max_cache_traj) real(kind=4),dimension(:,:,:),allocatable :: c_cache !(3,maxres2+2,max_cache_traj) integer :: ntwx_cache,ii_write !,max_cache_traj_use integer,dimension(:),allocatable :: iset_cache !(max_cache_traj) !----------------------------------------------------------------------------- ! common /przechowalnia/ real(kind=4),dimension(:,:),allocatable :: d_restart1 !(3,2*nres*maxprocs) real(kind=4),dimension(:,:),allocatable :: d_restart2 !(3,2*nres*maxprocs) real(kind=4),dimension(:,:),allocatable :: p_c !(3,(nres2+2)*maxprocs) !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- contains !----------------------------------------------------------------------------- ! !----------------------------------------------------------------------------- ! MREMD.F !----------------------------------------------------------------------------- subroutine MREMD use comm_gucio use control, only:tcpu,ovrtim use io_base, only:ilen use control_data use geometry_data use random, only: iran_num,ran_number use compare, only:hairpin,secondary2 use io, only:cartout,statout ! 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 real(kind=8),dimension(3) :: L,vcm real(kind=8) :: energia(0:n_ene) real(kind=8) :: remd_t_bath(Nprocs) !(maxprocs) integer :: iremd_iset(Nprocs) !(maxprocs) integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200) ! (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200) real(kind=8) :: remd_ene(0:n_ene+4,Nprocs) !(0:n_ene+4,maxprocs) integer :: iremd_acc(Nprocs),iremd_tot(Nprocs) !(maxprocs) integer :: iremd_acc_usa(Nprocs),iremd_tot_usa(Nprocs) !(maxprocs) integer :: rstcount !el ilen, !el external ilen character(len=50) :: tytul !el common /gucio/ cm integer :: itime !old integer nup(0:maxprocs),ndown(0:maxprocs) integer :: rep2i(0:Nprocs),ireqi(Nprocs) !(maxprocs) integer :: icache_all(Nprocs) !(maxprocs) integer :: status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,Nprocs)! (MPI_STATUS_SIZE,maxprocs) logical :: synflag, end_of_run, file_exist = .false.!, ovrtim real(kind=8) :: delta,time00,time01,time001,time02,time03,time04,& time05,time06,time07,time08,tt0,scalfac,ene_iex_iex,& ene_i_i,ene_iex_i,ene_i_iex,xxx,tmp,econstr_temp_i,& econstr_temp_iex integer :: k,il,il1,i,j,nharp,ii,ierr,itime_master,irr,iex,& i_set_temp,itmp,i_temp,i_mult,i_iset,i_mset,i_dir,i_temp1,& i_mult1,i_iset1,i_mset1,ierror integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3) !deb imin_itime_old=0 integer :: nres2 !el nres2=2*nres time001=0.0d0 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" !d print *,'MREMD',nodes !d print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep) !de 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 ! print *,'i2rep',me,i2rep(me) ! print *,'rep2i',(rep2i(i),i=0,nrep) !old if (i2rep(me).eq.nrep) then !old nup(0)=0 !old else !old nup(0)=remd_m(i2rep(me)+1) !old k=rep2i(int(i2rep(me)))+1 !old do i=1,nup(0) !old nup(i)=k !old k=k+1 !old enddo !old endif !d print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0)) !old if (i2rep(me).eq.1) then !old ndown(0)=0 !old else !old ndown(0)=remd_m(i2rep(me)-1) !old k=rep2i(i2rep(me)-2)+1 !old do i=1,ndown(0) !old ndown(i)=k !old k=k+1 !old enddo !old endif !d print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0)) !el common /przechowalnia/ if(.not.allocated(d_restart1)) allocate(d_restart1(3,nres2*nodes)) if(.not.allocated(d_restart2)) allocate(d_restart2(3,nres2*nodes)) if(.not.allocated(p_c)) allocate(p_c(3,(nres2+2)*nodes)) !el------------- write (*,*) "Processor",me," rest",rest,& "restart1fie",restart1file if(rest.and.restart1file) then if (me.eq.king) & inquire(file=mremd_rst_name,exist=file_exist) !d write (*,*) me," Before broadcast: file_exist",file_exist call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,& IERR) !d 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 ! ! 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 !d 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 ! stdfp=dsqrt(2*Rb*t_bath/d_time) do i=1,ntyp stdfsc(i)=dsqrt(2*Rb*t_bath/d_time) enddo ! 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 !------copy MD-------------- ! The driver for molecular dynamics subroutines !------------------------------------------------ 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 ! Determine the inverse of the inertia matrix. call setup_MD_matrices ! 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 ! 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) !d 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 ! Time-reversible RESPA algorithm ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992) call RESPA_step(itime) else ! 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 ! REMD - exchange ! 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) !deb if (out1file.or.traj1file) then !deb call mpi_gather(itime,1,mpi_integer, !deb & icache_all,1,mpi_integer,king, !deb & 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. !time call flush(iout) endif endif ! 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) !d write(iout,*) 'REMD synchro with',i !d 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 !deb call mpi_gather(itime,1,mpi_integer, !deb & itime_all,1,mpi_integer,king, !deb & CG_COMM,ierr) !deb write(iout,'(a19,8000i8)') ' REMD synchro itime', !deb & (itime_all(i),i=1,nodes) if(traj1file) then !deb imin_itime=itime_all(1) !deb do i=2,nodes !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i) !deb enddo !deb ii_write=(imin_itime-imin_itime_old)/ntwx !deb imin_itime_old=int(imin_itime/ntwx)*ntwx !deb 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) ! 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 !time 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 ! 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 ! call mpi_gather(t_bath,1,mpi_double_precision, ! & remd_t_bath,1,mpi_double_precision,king, ! & 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 !d debugging !deb call mpi_gather(ntwx_cache,1,mpi_integer, !deb & icache_all,1,mpi_integer,king, !deb & CG_COMM,ierr) !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file', !deb & (icache_all(i),i=1,nodes) !d 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 #ifdef DEBUG if(lmuca) then !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh' do i=1,nodes write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),& elowi(i),ehighi(i) enddo else write(iout,*) 'REMD exchange temp,ene' do i=1,nodes write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i) write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene) enddo endif #endif !------------------------------------- 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 #ifdef DEBUG write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i)) #endif if(i.gt.0.and.nupa(0,i).gt.0) then iex=i ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then ! write (iout,*) ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD" ! call flush(iout) ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr) ! endif ! do while (iex.eq.i) ! write (iout,*) "upper",nupa(int(nupa(0,i)),i) iex=nupa(iran_num(1,int(nupa(0,i))),i) ! enddo ! 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 ! Swap temperatures between conformations i and iex with recalculating the free energies ! following temperature changes. ene_iex_iex=remd_ene(0,iex) ene_i_i=remd_ene(0,i) ! write (iout,*) "i",i," ene_i_i",ene_i_i, ! & " iex",iex," ene_iex_iex",ene_iex_iex ! write (iout,*) "rescaling weights with temperature", ! & remd_t_bath(i) ! call flush(iout) call rescale_weights(remd_t_bath(i)) ! write (iout,*) "0,iex",remd_t_bath(i) ! call enerprint(remd_ene(0,iex)) call sum_energy(remd_ene(0,iex),.false.) ene_iex_i=remd_ene(0,iex) ! write (iout,*) "ene_iex_i",remd_ene(0,iex) ! write (iout,*) "0,i",remd_t_bath(i) ! call enerprint(remd_ene(0,i)) call sum_energy(remd_ene(0,i),.false.) ! write (iout,*) "ene_i_i",remd_ene(0,i) ! call flush(iout) ! write (iout,*) "rescaling weights with temperature", ! & 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)) ! write (iout,*) "0,i",remd_t_bath(iex) ! call enerprint(remd_ene(0,i)) call sum_energy(remd_ene(0,i),.false.) ! write (iout,*) "ene_i_iex",remd_ene(0,i) ! call flush(iout) ene_i_iex=remd_ene(0,i) ! write (iout,*) "0,iex",remd_t_bath(iex) ! 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 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex) ! write (iout,*) "i",i," iex",iex ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i, ! & " ene_i_iex",ene_i_iex, ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex ! 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 ! write(iout,*) 'delta',delta ! delta=(remd_t_bath(i)-remd_t_bath(iex))* ! & (remd_ene(i)-remd_ene(iex))/Rb/ ! & (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) ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx ! 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 ! write(iout,*) 'exchange',i,iex ! write (iout,'(a8,100i4)') "@ ifirst", ! & (ifirst(k),k=1,remd_m(1)) ! do il=1,nodes ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":", ! & (nupa(k,il),k=1,nupa(0,il)) ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":", ! & (ndowna(k,il),k=1,ndowna(0,il)) ! 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 !d write (iout,*) "exchange completed" !d call flush(iout) ELSE do ii=1,nodes !d 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) !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset i_dir=iran_num(1,3) !d 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 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1 call flush(iout) ! Swap temperatures between conformations i and iex with recalculating the free energies ! following temperature changes. ene_iex_iex=remd_ene(0,iex) ene_i_i=remd_ene(0,i) !o write (iout,*) "rescaling weights with temperature", !o & 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) !d write (iout,*) "ene_iex_i",remd_ene(0,iex) ! call sum_energy(remd_ene(0,i),.false.) !d write (iout,*) "ene_i_i",remd_ene(0,i) ! write (iout,*) "rescaling weights with temperature", ! & 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)) call sum_energy(remd_ene(0,i),.false.) !d write (iout,*) "ene_i_iex",remd_ene(0,i) ene_i_iex=remd_ene(0,i) ! 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 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex) ! write (iout,*) "i",i," iex",iex !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i, !d & " ene_i_iex",ene_i_iex, !d & " 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 !d write(iout,*) 'delta',delta ! delta=(remd_t_bath(i)-remd_t_bath(iex))* ! & (remd_ene(i)-remd_ene(iex))/Rb/ ! & (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) !d 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 !d do il=1,nset !d do il1=1,mset(il) !d do i=1,nrep !d do j=1,remd_m(i) !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1) !d enddo !d enddo !d enddo !d enddo 444 continue enddo ENDIF !------------------------------------- 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) !d write (iout,'(a6,100i4)') "ifirst", !d & (ifirst(i),i=1,remd_m(1)) !d do il=1,nodes !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":", !d & (nupa(i,il),i=1,nupa(0,il)) !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":", !d & (ndowna(i,il),i=1,ndowna(0,il)) !d enddo endif time06=MPI_WTIME() !d write (iout,*) "Before scatter" !d call flush(iout) #ifdef DEBUG if (me.eq.king) then write (iout,*) "t_bath before scatter",remd_t_bath call flush(iout) endif #endif call mpi_scatter(remd_t_bath,1,mpi_double_precision,& t_bath,1,mpi_double_precision,king,& CG_COMM,ierr) !d write (iout,*) "After scatter" !d 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) !o write (iout,*) "Processor",me, !o & " 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 !de write(iout,*) 'REMD after',me,t_bath time08=MPI_WTIME() if (me.eq.king .or. .not. out1file) then write(iout,*) 'REMD exchange time 8-0=',time08-time00 write(iout,*) 'REMD exchange time 8-7=',time08-time07 write(iout,*) 'REMD exchange time 7-6=',time07-time06 write(iout,*) 'REMD exchange time 6-5=',time06-time05 write(iout,*) 'REMD exchange time 5-4=',time05-time04 write(iout,*) 'REMD exchange time 4-3=',time04-time03 write(iout,*) 'REMD exchange time 3-2=',time03-time02 write(iout,*) 'REMD exchange time 2-1=',time02-time01 write(iout,*) 'REMD exchange time 1-0=',time01-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 !d debugging !deb call mpi_gather(ntwx_cache,1,mpi_integer, !deb & icache_all,1,mpi_integer,king, !deb & CG_COMM,ierr) !deb write(iout,'(a40,8000i8)') !deb & ' ntwx_cache after traj1file at the end', !deb & (icache_all(i),i=1,nodes) !d 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 !el common /przechowalnia/ ! deallocate(d_restart1) ! deallocate(d_restart2) ! deallocate(p_c) !el-------------- return end subroutine MREMD !----------------------------------------------------------------------------- subroutine write1rst(i_index) use control_data ! 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' !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),& !el d_restart2(3,2*nres*maxprocs) real(kind=4) :: r_d(3,2*nres) real(kind=4) :: t5_restart1(5) integer :: iret,itmp integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200) !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200) !el common /przechowalnia/ d_restart1,d_restart2 integer :: i,j,il,il1,ierr,ixdrf 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 write1rst !----------------------------------------------------------------------------- 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(kind=4) :: t5_restart1(5) integer :: iret,itmp real(kind=4) :: xcoord(3,2*nres+2),prec real(kind=4) :: r_qfrag(50),r_qpair(100) real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100) real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs) real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),& p_uscdiff(100*Nprocs) !(100*maxprocs) !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs) real(kind=4) :: r_c(3,2*nres+2) !el common /przechowalnia/ p_c integer :: i,j,il,ierr,ii,ixdrf call mpi_bcast(ii_write,1,mpi_integer,& king,CG_COMM,ierr) ! debugging print *,'traj1file',me,ii_write,ntwx_cache ! 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 ! write (iout,*) "before gather write1traj: from node",ii ! call flush(iout) ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii) ! call flush(iout) 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) ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4) call flush(iout) call mpi_gather(t5_restart1,5,mpi_real,& t_restart1,5,mpi_real,king,CG_COMM,ierr) ! do il=1,nodes ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il) ! enddo 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 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 #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) ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il) 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 write1traj !----------------------------------------------------------------------------- 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' !el real(kind=4) :: d_restart1(3,2*nres*maxprocs) real(kind=4) :: r_d(3,2*nres),t5_restart1(5) integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200) !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200) !el common /przechowalnia/ d_restart1 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf 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 ! read(irest2,'(3e15.5)') ! & (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 ! read(irest2,'(3e15.5)') ! & (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 !----------------------------------------------------------------------------- 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' !el real(kind=4) :: d_restart1(3,2*nres*maxprocs) real(kind=4) :: r_d(3,2*nres),t5_restart1(5) !el common /przechowalnia/ d_restart1 integer :: i,j,il,ierr 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 subroutine read1restart_old !---------------------------------------------------------------- subroutine alloc_MREMD_arrays ! if(.not.allocated(mset)) allocate(mset(max0(nset,1))) if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1)) !(ntyp1)) ! commom.remd ! common /remdcommon/ in io: read_REMDpar ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs) ! integer,dimension(:),allocatable :: remd_m !(maxprocs) ! common /remdrestart/ if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes)) allocate(i2set(0:2*nodes)) !(0:maxprocs) allocate(ifirst(0:nodes)) !(maxprocs) allocate(nupa(0:nodes,0:2*nodes)) allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs) allocate(t_restart1(5,nodes)) !(5,maxprocs) allocate(iset_restart1(nodes)) !(maxprocs) ! common /traj1cache/ allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj)) allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj)) allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj) allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj) allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj) allocate(ugamma_cache(nfrag_back,max_cache_traj)) allocate(utheta_cache(nfrag_back,max_cache_traj)) allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj) allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj) allocate(iset_cache(max_cache_traj)) !(max_cache_traj) return end subroutine alloc_MREMD_arrays !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module MREMDyn