2 !-----------------------------------------------------------------------------
10 ! use control_data, only:maxprocs
14 !-----------------------------------------------------------------------------
16 ! common /remdrestart/
17 integer(kind=2),dimension(:),allocatable :: i2set !(0:maxprocs)
18 integer(kind=2),dimension(:),allocatable :: ifirst !(maxprocs)
19 integer(kind=2),dimension(:,:),allocatable :: nupa,&
20 ndowna !(0:maxprocs/4,0:maxprocs)
21 real(kind=4),dimension(:,:),allocatable :: t_restart1 !(5,maxprocs)
22 integer,dimension(:),allocatable :: iset_restart1 !(maxprocs)
24 real(kind=4),dimension(:),allocatable :: totT_cache,EK_cache,&
25 potE_cache,t_bath_cache,Uconst_cache !(max_cache_traj)
26 real(kind=4),dimension(:,:),allocatable :: qfrag_cache !(50,max_cache_traj)
27 real(kind=4),dimension(:,:),allocatable :: qpair_cache !(100,max_cache_traj)
28 real(kind=4),dimension(:,:),allocatable :: ugamma_cache,&
29 utheta_cache,uscdiff_cache !(maxfrag_back,max_cache_traj)
30 real(kind=4),dimension(:,:,:),allocatable :: c_cache !(3,maxres2+2,max_cache_traj)
31 integer :: ntwx_cache,ii_write !,max_cache_traj_use
32 integer,dimension(:),allocatable :: iset_cache !(max_cache_traj)
33 !-----------------------------------------------------------------------------
34 ! common /przechowalnia/
35 real(kind=4),dimension(:,:),allocatable :: d_restart1 !(3,2*nres*maxprocs)
36 real(kind=4),dimension(:,:),allocatable :: d_restart2 !(3,2*nres*maxprocs)
37 real(kind=4),dimension(:,:),allocatable :: p_c !(3,(nres2+2)*maxprocs)
38 !-----------------------------------------------------------------------------
41 !-----------------------------------------------------------------------------
43 !-----------------------------------------------------------------------------
45 !-----------------------------------------------------------------------------
47 !-----------------------------------------------------------------------------
52 use control, only:tcpu,ovrtim
53 use io_base, only:ilen
56 use random, only: iran_num,ran_number
57 use compare, only:hairpin,secondary2
58 use io, only:cartout,statout
59 ! implicit real*8 (a-h,o-z)
60 ! include 'DIMENSIONS'
62 ! include 'COMMON.CONTROL'
63 ! include 'COMMON.VAR'
66 ! include 'COMMON.LANGEVIN'
68 ! include 'COMMON.LANGEVIN.lang0'
70 ! include 'COMMON.CHAIN'
71 ! include 'COMMON.DERIV'
72 ! include 'COMMON.GEO'
73 ! include 'COMMON.LOCAL'
74 ! include 'COMMON.INTERACT'
75 ! include 'COMMON.IOUNITS'
76 ! include 'COMMON.NAMES'
77 ! include 'COMMON.TIME1'
78 ! include 'COMMON.REMD'
79 ! include 'COMMON.SETUP'
80 ! include 'COMMON.MUCA'
81 ! include 'COMMON.HAIRPIN'
83 real(kind=8),dimension(3) :: L,vcm
84 real(kind=8) :: energia(0:n_ene)
86 real(kind=8) :: remd_t_bath(Nprocs) !(maxprocs)
87 integer :: iremd_iset(Nprocs) !(maxprocs)
88 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
89 ! (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
90 real(kind=8) :: remd_ene(0:n_ene+4,Nprocs) !(0:n_ene+4,maxprocs)
91 integer :: iremd_acc(Nprocs),iremd_tot(Nprocs) !(maxprocs)
92 integer :: iremd_acc_usa(Nprocs),iremd_tot_usa(Nprocs) !(maxprocs)
93 integer :: rstcount !el ilen,
95 character(len=50) :: tytul
98 !old integer nup(0:maxprocs),ndown(0:maxprocs)
99 integer :: rep2i(0:Nprocs),ireqi(Nprocs) !(maxprocs)
100 integer :: icache_all(Nprocs) !(maxprocs)
101 integer :: status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,Nprocs)! (MPI_STATUS_SIZE,maxprocs)
102 logical :: synflag, end_of_run, file_exist = .false.!, ovrtim
104 real(kind=8) :: delta,time00,time01,time001,time02,time03,time04,&
105 time05,time06,time07,time08,tt0,scalfac,ene_iex_iex,&
106 ene_i_i,ene_iex_i,ene_i_iex,xxx,tmp,econstr_temp_i,&
108 integer :: k,il,il1,i,j,nharp,ii,ierr,itime_master,irr,iex,&
109 i_set_temp,itmp,i_temp,i_mult,i_iset,i_mset,i_dir,i_temp1,&
110 i_mult1,i_iset1,i_mset1,ierror
111 integer,dimension(4,nres/3) :: iharp !(4,nres/3)(4,maxres/3)
112 !deb imin_itime_old=0
114 WRITE(iout,*) "JUST AFTER CALL"
115 ! if (.not.allocated(remd_ene)) allocate(remd_ene(0:n_ene+4,Nprocs))
122 if(me.eq.king.or..not.out1file) then
123 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
124 write (iout,*) "NREP=",nrep
128 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
129 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
131 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
133 print *,'MREMD',nodes
134 print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
135 write (iout,*) "Start MREMD: me",me," t_bath",t_bath
136 print *,"NSET=",nset, "MSET=", mset
140 do il1=1,max0(mset(il),1)
150 i_index(i,j,il,il1)=k
156 if(me.eq.king.or..not.out1file) then
157 write(iout,*) (i2rep(i),i=0,nodes-1)
158 write(iout,*) (i2set(i),i=0,nodes-1)
163 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
170 ! print *,'i2rep',me,i2rep(me)
171 ! print *,'rep2i',(rep2i(i),i=0,nrep)
173 !old if (i2rep(me).eq.nrep) then
176 !old nup(0)=remd_m(i2rep(me)+1)
177 !old k=rep2i(int(i2rep(me)))+1
184 !d print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
186 !old if (i2rep(me).eq.1) then
189 !old ndown(0)=remd_m(i2rep(me)-1)
190 !old k=rep2i(i2rep(me)-2)+1
197 !d print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
199 !el common /przechowalnia/
200 if(.not.allocated(d_restart1)) allocate(d_restart1(3,0:(nres2+1)*nodes))
201 if(.not.allocated(d_restart2)) allocate(d_restart2(3,0:(nres2+1)*nodes))
202 if(.not.allocated(p_c)) allocate(p_c(3,(nres2+2)*nodes))
205 write (*,*) "Processor",me," rest",rest,&
206 "restart1fie",restart1file
207 if(rest.and.restart1file) then
209 inquire(file=mremd_rst_name,exist=file_exist)
210 !d write (*,*) me," Before broadcast: file_exist",file_exist
211 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,&
213 !d write (*,*) me," After broadcast: file_exist",file_exist
215 if(me.eq.king.or..not.out1file) &
216 write (iout,*) 'Master is reading restart1file'
217 call read1restart(i_index)
219 if(me.eq.king.or..not.out1file) &
220 write (iout,*) 'WARNING : no restart1file'
223 if(me.eq.king.or..not.out1file) then
224 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
225 write(iout,*) "i_index"
230 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
239 if (rest.and..not.restart1file) &
240 inquire(file=mremd_rst_name,exist=file_exist)
241 if(.not.file_exist.and.rest.and..not.restart1file) &
242 write(iout,*) 'WARNING : no restart file',mremd_rst_name
243 IF (rest.and.file_exist.and..not.restart1file) THEN
244 write (iout,*) 'Master is reading restart file',&
246 open(irest2,file=mremd_rst_name,status='unknown')
248 read (irest2,*) (i2rep(i),i=0,nodes-1)
250 read (irest2,*) (ifirst(i),i=1,remd_m(1))
253 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
255 read (irest2,*) ndowna(0,il),&
256 (ndowna(i,il),i=1,ndowna(0,il))
262 read (irest2,*) (mset(i),i=1,nset)
264 read (irest2,*) (i2set(i),i=0,nodes-1)
269 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
274 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
275 write(iout,*) "i_index"
280 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
289 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
290 write (iout,'(a6,1000i5)') "ifirst",&
291 (ifirst(i),i=1,remd_m(1))
293 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
294 (nupa(i,il),i=1,nupa(0,il))
295 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
296 (ndowna(i,il),i=1,ndowna(0,il))
298 ELSE IF (.not.(rest.and.file_exist)) THEN
304 if (i2rep(il-1).eq.nrep) then
307 nupa(0,il)=remd_m(i2rep(il-1)+1)
308 k=rep2i(int(i2rep(il-1)))+1
314 if (i2rep(il-1).eq.1) then
317 ndowna(0,il)=remd_m(i2rep(il-1)-1)
318 k=rep2i(i2rep(il-1)-2)+1
326 write (iout,'(a6,100i4)') "ifirst",&
327 (ifirst(i),i=1,remd_m(1))
329 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",&
330 (nupa(i,il),i=1,nupa(0,il))
331 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",&
332 (ndowna(i,il),i=1,ndowna(0,il))
338 ! t_bath=retmin+(retmax-retmin)*me/(nodes-1)
339 if(.not.(rest.and.file_exist.and.restart1file)) then
340 if (me .eq. king) then
343 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
345 !d print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
346 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
351 if(me.eq.king.or..not.out1file) &
352 write(iout,*) me,"iset=",iset,"t_bath=",t_bath
356 stdfp(i)=dsqrt(2*Rb*t_bath/d_time)
360 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
364 ! print *,'irep',me,t_bath
366 if (me.eq.king .or. .not. out1file) &
367 write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
368 call rescale_weights(t_bath)
372 !------copy MD--------------
373 ! The driver for molecular dynamics subroutines
374 !------------------------------------------------
380 if(me.eq.king.or..not.out1file) &
381 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
387 ! Determine the inverse of the inertia matrix.
388 call setup_MD_matrices
392 if (me.eq.king .or. .not. out1file) &
393 write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
395 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
397 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
400 call rescale_weights(t_bath)
404 t_MDsetup = MPI_Wtime()-tt0
406 t_MDsetup = tcpu()-tt0
409 ! Entering the MD loop
415 if (lang.eq.2 .or. lang.eq.3) then
419 call sd_verlet_p_setup
421 call sd_verlet_ciccotti_setup
425 pfric0_mat(i,j,0)=pfric_mat(i,j)
426 afric0_mat(i,j,0)=afric_mat(i,j)
427 vfric0_mat(i,j,0)=vfric_mat(i,j)
428 prand0_mat(i,j,0)=prand_mat(i,j)
429 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
430 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
435 flag_stoch(i)=.false.
439 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
441 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
445 else if (lang.eq.1 .or. lang.eq.4) then
449 if (me.eq.king .or. .not. out1file) &
450 write(iout,*) 'Setup time',time00-walltime
453 t_langsetup=MPI_Wtime()-tt0
456 t_langsetup=tcpu()-tt0
462 do while(.not.end_of_run)
464 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
465 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
467 if (lang.gt.0 .and. surfarea .and. &
468 mod(itime,reset_fricmat).eq.0) then
469 if (lang.eq.2 .or. lang.eq.3) then
473 call sd_verlet_p_setup
475 call sd_verlet_ciccotti_setup
479 pfric0_mat(i,j,0)=pfric_mat(i,j)
480 afric0_mat(i,j,0)=afric_mat(i,j)
481 vfric0_mat(i,j,0)=vfric_mat(i,j)
482 prand0_mat(i,j,0)=prand_mat(i,j)
483 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
484 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
489 flag_stoch(i)=.false.
492 else if (lang.eq.1 .or. lang.eq.4) then
495 write (iout,'(a,i10)') &
496 "Friction matrix reset based on surface area, itime",itime
498 if (reset_vel .and. tbf .and. lang.eq.0 &
499 .and. mod(itime,count_reset_vel).eq.0) then
501 if (me.eq.king .or. .not. out1file) &
502 write(iout,'(a,f20.2)') &
503 "Velocities reset to random values, time",totT
506 d_t_old(j,i)=d_t(j,i)
510 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
514 d_t(j,0)=d_t(j,0)-vcm(j)
517 kinetic_T=2.0d0/(dimen3*Rb)*EK
518 scalfac=dsqrt(T_bath/kinetic_T)
519 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
522 d_t_old(j,i)=scalfac*d_t(j,i)
528 ! Time-reversible RESPA algorithm
529 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
530 call RESPA_step(itime)
532 ! Variable time step algorithm.
533 call velverlet_step(itime)
537 call brown_step(itime)
539 print *,"Brown dynamics not here!"
541 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
547 if (mod(itime,ntwe).eq.0) then
549 call enerprint(potEcomp)
552 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
553 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
555 call hairpin(.true.,nharp,iharp)
556 call secondary2(.true.)
557 call pdbout(potE,tytul,ipdb)
562 if (mod(itime,ntwx).eq.0.and.traj1file) then
563 if(ntwx_cache.lt.max_cache_traj_use) then
564 ntwx_cache=ntwx_cache+1
566 if (max_cache_traj_use.ne.1) &
567 print *,itime,"processor ",me," over cache ",ntwx_cache
570 totT_cache(i)=totT_cache(i+1)
571 EK_cache(i)=EK_cache(i+1)
572 potE_cache(i)=potE_cache(i+1)
573 t_bath_cache(i)=t_bath_cache(i+1)
574 Uconst_cache(i)=Uconst_cache(i+1)
575 iset_cache(i)=iset_cache(i+1)
578 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
581 qpair_cache(ii,i)=qpair_cache(ii,i+1)
584 utheta_cache(ii,i)=utheta_cache(ii,i+1)
585 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
586 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
592 c_cache(j,ii,i)=c_cache(j,ii,i+1)
598 totT_cache(ntwx_cache)=totT
599 EK_cache(ntwx_cache)=EK
600 potE_cache(ntwx_cache)=potE
601 t_bath_cache(ntwx_cache)=t_bath
602 Uconst_cache(ntwx_cache)=Uconst
603 iset_cache(ntwx_cache)=iset
606 qfrag_cache(i,ntwx_cache)=qfrag(i)
609 qpair_cache(i,ntwx_cache)=qpair(i)
612 utheta_cache(i,ntwx_cache)=utheta(i)
613 ugamma_cache(i,ntwx_cache)=ugamma(i)
614 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
619 c_cache(j,i,ntwx_cache)=c(j,i)
624 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
625 .and..not.restart1file) then
628 open(irest1,file=mremd_rst_name,status='unknown')
629 write (irest1,*) "i2rep"
630 write (irest1,*) (i2rep(i),i=0,nodes-1)
631 write (irest1,*) "ifirst"
632 write (irest1,*) (ifirst(i),i=1,remd_m(1))
634 write (irest1,*) "nupa",il
635 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
636 write (irest1,*) "ndowna",il
637 write (irest1,*) ndowna(0,il),&
638 (ndowna(i,il),i=1,ndowna(0,il))
641 write (irest1,*) "nset"
642 write (irest1,*) nset
643 write (irest1,*) "mset"
644 write (irest1,*) (mset(i),i=1,nset)
645 write (irest1,*) "i2set"
646 write (irest1,*) (i2set(i),i=0,nodes-1)
647 write (irest1,*) "i_index"
651 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
659 open(irest2,file=rest2name,status='unknown')
660 write(irest2,*) totT,EK,potE,totE,t_bath
662 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
665 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
668 write (irest2,*) iset
675 ! forced synchronization
676 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
677 .and. .not. mremdsync) then
679 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
681 call mpi_recv(itime_master, 1, MPI_INTEGER,&
682 0,101,CG_COMM, status, ierr)
683 call mpi_barrier(CG_COMM, ierr)
684 !deb if (out1file.or.traj1file) then
685 !deb call mpi_gather(itime,1,mpi_integer,
686 !deb & icache_all,1,mpi_integer,king,
689 call mpi_gather(ntwx_cache,1,mpi_integer,&
690 icache_all,1,mpi_integer,king,&
693 write(iout,*) 'REMD synchro at3',itime_master,itime
694 if (itime_master.ge.n_timestep .or. ovrtim()) &
696 !time call flush(iout)
701 if ((mod(itime,nstex).eq.0.and.me.eq.king &
702 .or.end_of_run.and.me.eq.king ) &
703 .and. .not. mremdsync ) then
706 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
707 CG_COMM, ireqi(i), ierr)
708 !d write(iout,*) 'REMD synchro with',i
711 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
712 call mpi_barrier(CG_COMM, ierr)
714 write(iout,*) 'REMD synchro at2',itime,'time=',time01-time00
715 if (out1file.or.traj1file) then
716 !deb call mpi_gather(itime,1,mpi_integer,
717 !deb & itime_all,1,mpi_integer,king,
719 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
720 !deb & (itime_all(i),i=1,nodes)
722 !deb imin_itime=itime_all(1)
724 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
726 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
727 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
728 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
729 call mpi_gather(ntwx_cache,1,mpi_integer,&
730 icache_all,1,mpi_integer,king,&
732 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
733 ! & (icache_all(i),i=1,nodes)
734 ii_write=icache_all(1)
736 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
738 ! write(iout,*) "MIN ii_write=",ii_write
741 !time call flush(iout)
743 if(mremdsync .and. mod(itime,nstex).eq.0) then
745 if (me.eq.king .or. .not. out1file) &
746 write(iout,*) 'REMD synchro at1',itime,ntwx_cache,Nprocs,nodes
747 !!!!!!! TRIAL OF MINIM SYNC
748 ! if (me.eq.king) then
750 ! call mpi_isend(itime,1,MPI_INTEGER,i,101, &
751 ! CG_COMM, ireqi(i), ierr)
752 !!d write(iout,*) 'REMD synchro with',i
755 ! call mpi_waitall(nodes-1,ireqi,statusi,ierr)
756 ! call mpi_barrier(CG_COMM, ierr)
760 ! if (me.ne.king) then
761 ! call mpi_recv(itime_master, 1, MPI_INTEGER,&
762 ! 0,101,CG_COMM, status, ierr)
763 ! call mpi_barrier(CG_COMM, ierr)
764 !deb if (out1file.or.traj1file) then
767 write(iout,*) icache_all
769 write(iout,*) "before mpi_gather ntwx_cache"
770 call mpi_gather(ntwx_cache,1,mpi_integer,&
771 icache_all(1),1,mpi_integer,king,& ! CONSULT WITH ADAM
773 write(iout,*) "after mpi_gather ntwx_cache"
776 write(iout,'(a19,8000i8)') ' ntwx_cache',&
777 (icache_all(i),i=1,nodes)
778 ii_write=icache_all(1)
780 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
782 write(iout,*) "MIN ii_write=",ii_write
788 ! Update the time safety limiy
789 if (time001-time00.gt.safety) then
790 safety=time001-time00+600
791 write (iout,*) "****** SAFETY increased to",safety," s"
793 if (ovrtim()) end_of_run=.true.
795 if(synflag.and..not.end_of_run) then
799 write(iout,*) 'REMD before',me,t_bath
801 ! call mpi_gather(t_bath,1,mpi_double_precision,
802 ! & remd_t_bath,1,mpi_double_precision,king,
804 potEcomp(n_ene+1)=t_bath
806 potEcomp(n_ene+2)=iset
807 if (iset.lt.nset) then
811 potEcomp(n_ene+3)=Uconst
818 potEcomp(n_ene+4)=Uconst
822 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
823 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
826 call mpi_gather(elow,1,mpi_double_precision,&
827 elowi,1,mpi_double_precision,king,&
829 call mpi_gather(ehigh,1,mpi_double_precision,&
830 ehighi,1,mpi_double_precision,king,&
835 if (me.eq.king .or. .not. out1file) then
836 write(iout,*) 'REMD gather times=',time03-time01 &
840 if (restart1file) call write1rst(i_index)
843 if (me.eq.king .or. .not. out1file) then
844 write(iout,*) 'REMD writing rst time=',time04-time03
847 if (traj1file) call write1traj
849 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
850 !deb & icache_all,1,mpi_integer,king,
852 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
853 !deb & (icache_all(i),i=1,nodes)
858 if (me.eq.king .or. .not. out1file) then
859 write(iout,*) 'REMD writing traj time=',time05-time04
866 remd_t_bath(i)=remd_ene(n_ene+1,i)
867 iremd_iset(i)=remd_ene(n_ene+2,i)
871 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
873 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
877 write(iout,*) 'REMD exchange temp,ene'
879 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
880 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
884 !-------------------------------------
886 write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
889 write (iout,*) "remd_m(1)",remd_m(1)
891 i=ifirst(iran_num(1,remd_m(1)))
898 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
900 if(i.gt.0.and.nupa(0,i).gt.0) then
902 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
904 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
906 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
908 ! do while (iex.eq.i)
909 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
910 iex=nupa(iran_num(1,int(nupa(0,i))),i)
912 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
914 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
916 ! Swap temperatures between conformations i and iex with recalculating the free energies
917 ! following temperature changes.
918 ene_iex_iex=remd_ene(0,iex)
919 ene_i_i=remd_ene(0,i)
920 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
921 ! & " iex",iex," ene_iex_iex",ene_iex_iex
922 ! write (iout,*) "rescaling weights with temperature",
925 call rescale_weights(remd_t_bath(i))
927 ! write (iout,*) "0,iex",remd_t_bath(i)
928 ! call enerprint(remd_ene(0,iex))
930 call sum_energy(remd_ene(0,iex),.false.)
931 ene_iex_i=remd_ene(0,iex)
932 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
934 ! write (iout,*) "0,i",remd_t_bath(i)
935 ! call enerprint(remd_ene(0,i))
937 call sum_energy(remd_ene(0,i),.false.)
938 ! write (iout,*) "ene_i_i",remd_ene(0,i)
940 ! write (iout,*) "rescaling weights with temperature",
942 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
943 write (iout,*) "ERROR: inconsistent energies:",i,&
944 ene_i_i,remd_ene(0,i)
946 call rescale_weights(remd_t_bath(iex))
948 ! write (iout,*) "0,i",remd_t_bath(iex)
949 ! call enerprint(remd_ene(0,i))
951 call sum_energy(remd_ene(0,i),.false.)
952 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
954 ene_i_iex=remd_ene(0,i)
956 ! write (iout,*) "0,iex",remd_t_bath(iex)
957 call enerprint(remd_ene(0,iex))
959 call sum_energy(remd_ene(0,iex),.false.)
960 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
961 write (iout,*) "ERROR: inconsistent energies:",iex,&
962 ene_iex_iex,remd_ene(0,iex)
964 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
965 ! write (iout,*) "i",i," iex",iex
966 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
967 ! & " ene_i_iex",ene_i_iex,
968 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
970 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
971 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
973 ! write(iout,*) 'delta',delta
974 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
975 ! & (remd_ene(i)-remd_ene(iex))/Rb/
976 ! & (remd_t_bath(i)*remd_t_bath(iex))
978 if (delta .gt. 50.0d0) then
984 else if (delta.lt.-50.0d0) then
993 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
994 xxx=ran_number(0.0d0,1.0d0)
995 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
997 if (delta .gt. xxx) then
999 remd_t_bath(i)=remd_t_bath(iex)
1000 remd_t_bath(iex)=tmp
1001 remd_ene(0,i)=ene_i_iex
1002 remd_ene(0,iex)=ene_iex_i
1008 ehighi(i)=ehighi(iex)
1015 nupa(k,i)=nupa(k,iex)
1018 ndowna(k,i)=ndowna(k,iex)
1022 if (ifirst(il).eq.i) ifirst(il)=iex
1024 if (nupa(k,il).eq.i) then
1026 elseif (nupa(k,il).eq.iex) then
1031 if (ndowna(k,il).eq.i) then
1033 elseif (ndowna(k,il).eq.iex) then
1039 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1041 i2rep(i-1)=i2rep(iex-1)
1044 ! write(iout,*) 'exchange',i,iex
1045 ! write (iout,'(a8,100i4)') "@ ifirst",
1046 ! & (ifirst(k),k=1,remd_m(1))
1048 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1049 ! & (nupa(k,il),k=1,nupa(0,il))
1050 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1051 ! & (ndowna(k,il),k=1,ndowna(0,il))
1056 remd_ene(0,iex)=ene_iex_iex
1057 remd_ene(0,i)=ene_i_i
1063 !d write (iout,*) "exchange completed"
1067 !d write(iout,*) "########",ii
1069 i_temp=iran_num(1,nrep)
1070 i_mult=iran_num(1,remd_m(i_temp))
1071 i_iset=iran_num(1,nset)
1072 i_mset=iran_num(1,mset(i_iset))
1073 i=i_index(i_temp,i_mult,i_iset,i_mset)
1075 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1078 !d write(iout,*) "i_dir=",i_dir
1080 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1083 i_mult1=iran_num(1,remd_m(i_temp1))
1085 i_mset1=iran_num(1,mset(i_iset1))
1086 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1088 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1091 i_mult1=iran_num(1,remd_m(i_temp1))
1093 i_mset1=iran_num(1,mset(i_iset1))
1094 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1095 econstr_temp_i=remd_ene(20,i)
1096 econstr_temp_iex=remd_ene(20,iex)
1097 remd_ene(20,i)=remd_ene(n_ene+3,i)
1098 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1100 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1103 i_mult1=iran_num(1,remd_m(i_temp1))
1105 i_mset1=iran_num(1,mset(i_iset1))
1106 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1107 econstr_temp_i=remd_ene(20,i)
1108 econstr_temp_iex=remd_ene(20,iex)
1109 remd_ene(20,i)=remd_ene(n_ene+3,i)
1110 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1116 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1119 ! Swap temperatures between conformations i and iex with recalculating the free energies
1120 ! following temperature changes.
1121 ene_iex_iex=remd_ene(0,iex)
1122 ene_i_i=remd_ene(0,i)
1123 !o write (iout,*) "rescaling weights with temperature",
1125 call rescale_weights(remd_t_bath(i))
1127 call sum_energy(remd_ene(0,iex),.false.)
1128 ene_iex_i=remd_ene(0,iex)
1129 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1130 ! call sum_energy(remd_ene(0,i),.false.)
1131 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1132 ! write (iout,*) "rescaling weights with temperature",
1133 ! & remd_t_bath(iex)
1134 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1135 ! write (iout,*) "ERROR: inconsistent energies:",i,
1136 ! & ene_i_i,remd_ene(0,i)
1138 call rescale_weights(remd_t_bath(iex))
1139 call sum_energy(remd_ene(0,i),.false.)
1140 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1141 ene_i_iex=remd_ene(0,i)
1142 ! call sum_energy(remd_ene(0,iex),.false.)
1143 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1144 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1145 ! & ene_iex_iex,remd_ene(0,iex)
1147 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1148 ! write (iout,*) "i",i," iex",iex
1149 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1150 !d & " ene_i_iex",ene_i_iex,
1151 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1152 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1153 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1155 !d write(iout,*) 'delta',delta
1156 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1157 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1158 ! & (remd_t_bath(i)*remd_t_bath(iex))
1159 if (delta .gt. 50.0d0) then
1164 if (i_dir.eq.1.or.i_dir.eq.3) &
1165 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1166 if (i_dir.eq.2.or.i_dir.eq.3) &
1167 iremd_tot_usa(int(i2set(i-1)))= &
1168 iremd_tot_usa(int(i2set(i-1)))+1
1169 xxx=ran_number(0.0d0,1.0d0)
1170 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1171 if (delta .gt. xxx) then
1173 remd_t_bath(i)=remd_t_bath(iex)
1174 remd_t_bath(iex)=tmp
1177 iremd_iset(i)=iremd_iset(iex)
1178 iremd_iset(iex)=itmp
1180 remd_ene(0,i)=ene_i_iex
1181 remd_ene(0,iex)=ene_iex_i
1183 if (i_dir.eq.1.or.i_dir.eq.3) &
1184 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1187 i2rep(i-1)=i2rep(iex-1)
1190 if (i_dir.eq.2.or.i_dir.eq.3) &
1191 iremd_acc_usa(int(i2set(i-1)))= &
1192 iremd_acc_usa(int(i2set(i-1)))+1
1195 i2set(i-1)=i2set(iex-1)
1198 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1199 i_index(i_temp,i_mult,i_iset,i_mset)= &
1200 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1201 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1204 remd_ene(0,iex)=ene_iex_iex
1205 remd_ene(0,i)=ene_i_i
1206 remd_ene(20,iex)=econstr_temp_iex
1207 remd_ene(20,i)=econstr_temp_i
1211 !d do il1=1,mset(il)
1214 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1227 !-------------------------------------
1228 write (iout,*) "NREP",nrep
1230 if(iremd_tot(i).ne.0) &
1231 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1232 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1237 if(iremd_tot_usa(i).ne.0) &
1238 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1239 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1245 !d write (iout,'(a6,100i4)') "ifirst",
1246 !d & (ifirst(i),i=1,remd_m(1))
1248 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1249 !d & (nupa(i,il),i=1,nupa(0,il))
1250 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1251 !d & (ndowna(i,il),i=1,ndowna(0,il))
1256 !d write (iout,*) "Before scatter"
1259 if (me.eq.king) then
1260 write (iout,*) "t_bath before scatter",remd_t_bath
1264 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1265 t_bath,1,mpi_double_precision,king,&
1267 !d write (iout,*) "After scatter"
1270 call mpi_scatter(iremd_iset,1,mpi_integer,&
1271 iset,1,mpi_integer,king,&
1275 if (me.eq.king .or. .not. out1file) then
1276 write(iout,*) 'REMD scatter time=',time07-time06
1280 call mpi_scatter(elowi,1,mpi_double_precision,&
1281 elow,1,mpi_double_precision,king,&
1283 call mpi_scatter(ehighi,1,mpi_double_precision,&
1284 ehigh,1,mpi_double_precision,king,&
1287 call rescale_weights(t_bath)
1288 !o write (iout,*) "Processor",me,
1289 !o & " rescaling weights with temperature",t_bath
1291 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
1293 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
1297 !de write(iout,*) 'REMD after',me,t_bath
1299 if (me.eq.king .or. .not. out1file) then
1300 write(iout,*) 'REMD exchange time 8-0=',time08-time00
1301 write(iout,*) 'REMD exchange time 8-7=',time08-time07
1302 write(iout,*) 'REMD exchange time 7-6=',time07-time06
1303 write(iout,*) 'REMD exchange time 6-5=',time06-time05
1304 write(iout,*) 'REMD exchange time 5-4=',time05-time04
1305 write(iout,*) 'REMD exchange time 4-3=',time04-time03
1306 write(iout,*) 'REMD exchange time 3-2=',time03-time02
1307 write(iout,*) 'REMD exchange time 2-1=',time02-time01
1308 write(iout,*) 'REMD exchange time 1-0=',time01-time00
1314 if (restart1file) then
1315 if (me.eq.king .or. .not. out1file) &
1316 write(iout,*) 'writing restart at the end of run'
1317 call write1rst(i_index)
1320 if (traj1file) call write1traj
1322 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1323 !deb & icache_all,1,mpi_integer,king,
1324 !deb & CG_COMM,ierr)
1325 !deb write(iout,'(a40,8000i8)')
1326 !deb & ' ntwx_cache after traj1file at the end',
1327 !deb & (icache_all(i),i=1,nodes)
1332 t_MD=MPI_Wtime()-tt0
1336 if (me.eq.king .or. .not. out1file) then
1337 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1339 'MD calculations setup:',t_MDsetup,&
1340 'Energy & gradient evaluation:',t_enegrad,&
1341 'Stochastic MD setup:',t_langsetup,&
1342 'Stochastic MD step setup:',t_sdsetup,&
1344 write (iout,'(/28(1h=),a25,27(1h=))') &
1345 ' End of MD calculation '
1347 !el common /przechowalnia/
1348 ! deallocate(d_restart1)
1349 ! deallocate(d_restart2)
1353 end subroutine MREMD
1354 !-----------------------------------------------------------------------------
1355 subroutine write1rst(i_index)
1358 ! implicit real*8 (a-h,o-z)
1359 ! include 'DIMENSIONS'
1361 ! include 'COMMON.MD'
1362 ! include 'COMMON.IOUNITS'
1363 ! include 'COMMON.REMD'
1364 ! include 'COMMON.SETUP'
1365 ! include 'COMMON.CHAIN'
1366 ! include 'COMMON.SBRIDGE'
1367 ! include 'COMMON.INTERACT'
1369 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1370 !el d_restart2(3,2*nres*maxprocs)
1371 real(kind=4) :: r_d(3,0:2*nres)
1372 real(kind=4) :: t5_restart1(5)
1373 integer :: iret,itmp
1374 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1375 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1377 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1378 !el common /przechowalnia/ d_restart1,d_restart2
1379 integer :: i,j,il,il1,ierr,ixdrf
1384 t5_restart1(4)=t_bath
1385 t5_restart1(5)=Uconst
1387 call mpi_gather(t5_restart1,5,mpi_real,&
1388 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1396 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1397 d_restart1,3*2*nres+3,mpi_real,king,&
1409 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1410 d_restart2,3*2*nres+3,mpi_real,king,&
1415 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1417 call xdrfint_(ixdrf, i2rep(i), iret)
1420 call xdrfint_(ixdrf, ifirst(i), iret)
1424 call xdrfint_(ixdrf, nupa(i,il), iret)
1428 call xdrfint_(ixdrf, ndowna(i,il), iret)
1434 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1441 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1448 call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1454 call xdrfint_(ixdrf, nset, iret)
1456 call xdrfint_(ixdrf,mset(i), iret)
1459 call xdrfint_(ixdrf,i2set(i), iret)
1465 itmp=i_index(i,j,il,il1)
1466 call xdrfint_(ixdrf,itmp, iret)
1473 call xdrfclose_(ixdrf, iret)
1475 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1477 call xdrfint(ixdrf, i2rep(i), iret)
1480 call xdrfint(ixdrf, ifirst(i), iret)
1484 call xdrfint(ixdrf, nupa(i,il), iret)
1488 call xdrfint(ixdrf, ndowna(i,il), iret)
1494 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1501 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1508 call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1515 call xdrfint(ixdrf, nset, iret)
1517 call xdrfint(ixdrf,mset(i), iret)
1520 call xdrfint(ixdrf,i2set(i), iret)
1526 itmp=i_index(i,j,il,il1)
1527 call xdrfint(ixdrf,itmp, iret)
1534 call xdrfclose(ixdrf, iret)
1538 end subroutine write1rst
1539 !-----------------------------------------------------------------------------
1540 subroutine write1traj
1542 ! implicit real*8 (a-h,o-z)
1543 ! include 'DIMENSIONS'
1545 ! include 'COMMON.MD'
1546 ! include 'COMMON.IOUNITS'
1547 ! include 'COMMON.REMD'
1548 ! include 'COMMON.SETUP'
1549 ! include 'COMMON.CHAIN'
1550 ! include 'COMMON.SBRIDGE'
1551 ! include 'COMMON.INTERACT'
1553 real(kind=4) :: t5_restart1(5)
1554 integer :: iret,itmp
1555 real(kind=4) :: xcoord(3,2*nres+2),prec
1556 real(kind=4) :: r_qfrag(50),r_qpair(100)
1557 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1558 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1559 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1560 p_uscdiff(100*Nprocs) !(100*maxprocs)
1561 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1562 real(kind=4) :: r_c(3,2*nres+2)
1563 !el common /przechowalnia/ p_c
1565 integer :: i,j,il,ierr,ii,ixdrf
1567 call mpi_bcast(ii_write,1,mpi_integer,&
1571 ! print *,'traj1file',me,ii_write,ntwx_cache
1575 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1577 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1580 ! write (iout,*) "before gather write1traj: from node",ii
1582 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1584 t5_restart1(1)=totT_cache(ii)
1585 t5_restart1(2)=EK_cache(ii)
1586 t5_restart1(3)=potE_cache(ii)
1587 t5_restart1(4)=t_bath_cache(ii)
1588 t5_restart1(5)=Uconst_cache(ii)
1589 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1591 call mpi_gather(t5_restart1,5,mpi_real,&
1592 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1594 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1597 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1598 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1601 r_qfrag(i)=qfrag_cache(i,ii)
1604 r_qpair(i)=qpair_cache(i,ii)
1607 r_utheta(i)=utheta_cache(i,ii)
1608 r_ugamma(i)=ugamma_cache(i,ii)
1609 r_uscdiff(i)=uscdiff_cache(i,ii)
1612 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1613 p_qfrag,nfrag,mpi_real,king,&
1615 call mpi_gather(r_qpair,npair,mpi_real,&
1616 p_qpair,npair,mpi_real,king,&
1618 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1619 p_utheta,nfrag_back,mpi_real,king,&
1621 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1622 p_ugamma,nfrag_back,mpi_real,king,&
1624 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1625 p_uscdiff,nfrag_back,mpi_real,king,&
1629 write (iout,*) "p_qfrag"
1631 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1633 write (iout,*) "p_qpair"
1635 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1641 r_c(j,i)=c_cache(j,i,ii)
1645 call mpi_gather(r_c,3*2*nres,mpi_real,&
1646 p_c,3*2*nres,mpi_real,king,&
1652 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1653 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1654 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1655 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1656 call xdrfint_(ixdrf, nss, iret)
1659 call xdrfint(ixdrf, idssb(j)+nres, iret)
1660 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1662 call xdrfint_(ixdrf, ihpb(j), iret)
1663 call xdrfint_(ixdrf, jhpb(j), iret)
1666 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1667 call xdrfint_(ixdrf, iset_restart1(il), iret)
1669 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1672 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1675 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1676 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1677 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1682 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1687 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1691 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1695 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1696 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1697 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1698 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1699 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1700 call xdrfint(ixdrf, nss, iret)
1703 call xdrfint(ixdrf, idssb(j)+nres, iret)
1704 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1706 call xdrfint(ixdrf, ihpb(j), iret)
1707 call xdrfint(ixdrf, jhpb(j), iret)
1710 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1711 call xdrfint(ixdrf, iset_restart1(il), iret)
1713 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1716 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1719 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1720 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1721 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1726 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1731 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1735 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1741 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1743 if(me.eq.king) call xdrfclose(ixdrf, iret)
1745 do i=1,ntwx_cache-ii_write
1747 totT_cache(i)=totT_cache(ii_write+i)
1748 EK_cache(i)=EK_cache(ii_write+i)
1749 potE_cache(i)=potE_cache(ii_write+i)
1750 t_bath_cache(i)=t_bath_cache(ii_write+i)
1751 Uconst_cache(i)=Uconst_cache(ii_write+i)
1752 iset_cache(i)=iset_cache(ii_write+i)
1755 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1758 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1761 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1762 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1763 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1768 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1772 ntwx_cache=ntwx_cache-ii_write
1774 end subroutine write1traj
1775 !-----------------------------------------------------------------------------
1776 subroutine read1restart(i_index)
1778 ! implicit real*8 (a-h,o-z)
1779 ! include 'DIMENSIONS'
1781 ! include 'COMMON.MD'
1782 ! include 'COMMON.IOUNITS'
1783 ! include 'COMMON.REMD'
1784 ! include 'COMMON.SETUP'
1785 ! include 'COMMON.CHAIN'
1786 ! include 'COMMON.SBRIDGE'
1787 ! include 'COMMON.INTERACT'
1788 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1789 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1790 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1791 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1793 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1794 !el common /przechowalnia/ d_restart1
1795 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1797 write (*,*) "Processor",me," called read1restart"
1800 open(irest2,file=mremd_rst_name,status='unknown')
1801 read(irest2,*,err=334) i
1802 write(iout,*) "Reading old rst in ASCI format"
1804 call read1restart_old
1808 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1811 call xdrfint_(ixdrf, i2rep(i), iret)
1814 call xdrfint_(ixdrf, ifirst(i), iret)
1817 call xdrfint_(ixdrf, nupa(0,il), iret)
1819 call xdrfint_(ixdrf, nupa(i,il), iret)
1822 call xdrfint_(ixdrf, ndowna(0,il), iret)
1824 call xdrfint_(ixdrf, ndowna(i,il), iret)
1829 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1833 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1836 call xdrfint(ixdrf, i2rep(i), iret)
1839 call xdrfint(ixdrf, ifirst(i), iret)
1842 call xdrfint(ixdrf, nupa(0,il), iret)
1844 call xdrfint(ixdrf, nupa(i,il), iret)
1847 call xdrfint(ixdrf, ndowna(0,il), iret)
1849 call xdrfint(ixdrf, ndowna(i,il), iret)
1854 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1859 call mpi_scatter(t_restart1,5,mpi_real,&
1860 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1864 t_bath=t5_restart1(4)
1869 ! read(irest2,'(3e15.5)')
1870 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1873 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1875 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1881 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1882 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1892 ! read(irest2,'(3e15.5)')
1893 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1896 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1898 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1904 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1905 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1916 call xdrfint_(ixdrf, nset, iret)
1918 call xdrfint_(ixdrf,mset(i), iret)
1921 call xdrfint_(ixdrf,i2set(i), iret)
1927 call xdrfint_(ixdrf,itmp, iret)
1928 i_index(i,j,il,il1)=itmp
1936 call xdrfint(ixdrf, nset, iret)
1938 call xdrfint(ixdrf,mset(i), iret)
1941 call xdrfint(ixdrf,i2set(i), iret)
1947 call xdrfint(ixdrf,itmp, iret)
1948 i_index(i,j,il,il1)=itmp
1955 call mpi_scatter(i2set,1,mpi_integer,&
1956 iset,1,mpi_integer,king,&
1961 if(me.eq.king) close(irest2)
1963 end subroutine read1restart
1964 !-----------------------------------------------------------------------------
1965 subroutine read1restart_old
1967 ! implicit real*8 (a-h,o-z)
1968 ! include 'DIMENSIONS'
1970 ! include 'COMMON.MD'
1971 ! include 'COMMON.IOUNITS'
1972 ! include 'COMMON.REMD'
1973 ! include 'COMMON.SETUP'
1974 ! include 'COMMON.CHAIN'
1975 ! include 'COMMON.SBRIDGE'
1976 ! include 'COMMON.INTERACT'
1977 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1978 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1979 !el common /przechowalnia/ d_restart1
1981 integer :: i,j,il,ierr
1984 open(irest2,file=mremd_rst_name,status='unknown')
1985 read (irest2,*) (i2rep(i),i=0,nodes-1)
1986 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1988 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1989 read (irest2,*) ndowna(0,il),&
1990 (ndowna(i,il),i=1,ndowna(0,il))
1993 read(irest2,*) (t_restart1(j,il),j=1,4)
1996 call mpi_scatter(t_restart1,5,mpi_real,&
1997 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2001 t_bath=t5_restart1(4)
2006 read(irest2,'(3e15.5)') &
2007 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2011 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2012 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2022 read(irest2,'(3e15.5)') &
2023 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2027 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2028 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2034 if(me.eq.king) close(irest2)
2036 end subroutine read1restart_old
2037 !----------------------------------------------------------------
2038 subroutine alloc_MREMD_arrays
2040 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
2041 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1,5)) !(ntyp1))
2043 ! common /remdcommon/ in io: read_REMDpar
2044 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2045 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2046 ! common /remdrestart/
2047 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2049 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2050 allocate(ifirst(0:nodes)) !(maxprocs)
2051 allocate(nupa(0:nodes,0:2*nodes))
2052 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2053 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2054 allocate(iset_restart1(nodes)) !(maxprocs)
2055 ! common /traj1cache/
2056 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2057 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2058 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2059 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2060 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2061 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2062 allocate(utheta_cache(nfrag_back,max_cache_traj))
2063 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2064 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2065 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2068 end subroutine alloc_MREMD_arrays
2069 !-----------------------------------------------------------------------------
2070 !-----------------------------------------------------------------------------