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/2,Nprocs/2,Nprocs/10,Nprocs/10)
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,itime
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
463 do while(.not.end_of_run)
466 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
467 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
469 if (lang.gt.0 .and. surfarea .and. &
470 mod(itime,reset_fricmat).eq.0) then
471 if (lang.eq.2 .or. lang.eq.3) then
475 call sd_verlet_p_setup
477 call sd_verlet_ciccotti_setup
481 pfric0_mat(i,j,0)=pfric_mat(i,j)
482 afric0_mat(i,j,0)=afric_mat(i,j)
483 vfric0_mat(i,j,0)=vfric_mat(i,j)
484 prand0_mat(i,j,0)=prand_mat(i,j)
485 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
486 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
491 flag_stoch(i)=.false.
494 else if (lang.eq.1 .or. lang.eq.4) then
497 write (iout,'(a,i10)') &
498 "Friction matrix reset based on surface area, itime",itime
500 if (reset_vel .and. tbf .and. lang.eq.0 &
501 .and. mod(itime,count_reset_vel).eq.0) then
503 if (me.eq.king .or. .not. out1file) &
504 write(iout,'(a,f20.2)') &
505 "Velocities reset to random values, time",totT
508 d_t_old(j,i)=d_t(j,i)
512 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
516 d_t(j,0)=d_t(j,0)-vcm(j)
519 kinetic_T=2.0d0/(dimen3*Rb)*EK
520 scalfac=dsqrt(T_bath/kinetic_T)
521 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
524 d_t_old(j,i)=scalfac*d_t(j,i)
530 ! Time-reversible RESPA algorithm
531 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
532 call RESPA_step(itime)
534 ! Variable time step algorithm.
535 call velverlet_step(itime)
539 call brown_step(itime)
541 print *,"Brown dynamics not here!"
543 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
550 if (mod(itime,ntwe).eq.0) then
552 call enerprint(potEcomp)
555 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
556 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
558 call hairpin(.true.,nharp,iharp)
559 call secondary2(.true.)
560 call pdbout(potE,tytul,ipdb)
565 if (mod(itime,ntwx).eq.0.and.traj1file) then
566 if(ntwx_cache.lt.max_cache_traj_use) then
567 ntwx_cache=ntwx_cache+1
569 if (max_cache_traj_use.ne.1) &
570 print *,itime,"processor ",me," over cache ",ntwx_cache
573 totT_cache(i)=totT_cache(i+1)
574 EK_cache(i)=EK_cache(i+1)
575 potE_cache(i)=potE_cache(i+1)
576 t_bath_cache(i)=t_bath_cache(i+1)
577 Uconst_cache(i)=Uconst_cache(i+1)
578 iset_cache(i)=iset_cache(i+1)
581 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
584 qpair_cache(ii,i)=qpair_cache(ii,i+1)
587 utheta_cache(ii,i)=utheta_cache(ii,i+1)
588 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
589 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
595 c_cache(j,ii,i)=c_cache(j,ii,i+1)
601 totT_cache(ntwx_cache)=totT
602 EK_cache(ntwx_cache)=EK
603 potE_cache(ntwx_cache)=potE
604 t_bath_cache(ntwx_cache)=t_bath
605 Uconst_cache(ntwx_cache)=Uconst
606 iset_cache(ntwx_cache)=iset
609 qfrag_cache(i,ntwx_cache)=qfrag(i)
612 qpair_cache(i,ntwx_cache)=qpair(i)
615 utheta_cache(i,ntwx_cache)=utheta(i)
616 ugamma_cache(i,ntwx_cache)=ugamma(i)
617 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
622 c_cache(j,i,ntwx_cache)=c(j,i)
627 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
628 .and..not.restart1file) then
631 open(irest1,file=mremd_rst_name,status='unknown')
632 write (irest1,*) "i2rep"
633 write (irest1,*) (i2rep(i),i=0,nodes-1)
634 write (irest1,*) "ifirst"
635 write (irest1,*) (ifirst(i),i=1,remd_m(1))
637 write (irest1,*) "nupa",il
638 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
639 write (irest1,*) "ndowna",il
640 write (irest1,*) ndowna(0,il),&
641 (ndowna(i,il),i=1,ndowna(0,il))
644 write (irest1,*) "nset"
645 write (irest1,*) nset
646 write (irest1,*) "mset"
647 write (irest1,*) (mset(i),i=1,nset)
648 write (irest1,*) "i2set"
649 write (irest1,*) (i2set(i),i=0,nodes-1)
650 write (irest1,*) "i_index"
654 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
662 open(irest2,file=rest2name,status='unknown')
663 write(irest2,*) totT,EK,potE,totE,t_bath
665 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
668 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
671 write (irest2,*) iset
678 ! forced synchronization
679 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
680 .and. .not. mremdsync) then
682 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
684 call mpi_recv(itime_master, 1, MPI_INTEGER,&
685 0,101,CG_COMM, status, ierr)
686 call mpi_barrier(CG_COMM, ierr)
687 !deb if (out1file.or.traj1file) then
688 !deb call mpi_gather(itime,1,mpi_integer,
689 !deb & icache_all,1,mpi_integer,king,
692 call mpi_gather(ntwx_cache,1,mpi_integer,&
693 icache_all,1,mpi_integer,king,&
696 write(iout,*) 'REMD synchro at3',itime_master,itime
697 if (itime_master.ge.n_timestep .or. ovrtim()) &
699 !time call flush(iout)
704 if ((mod(itime,nstex).eq.0.and.me.eq.king &
705 .or.end_of_run.and.me.eq.king ) &
706 .and. .not. mremdsync ) then
709 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
710 CG_COMM, ireqi(i), ierr)
711 !d write(iout,*) 'REMD synchro with',i
714 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
715 call mpi_barrier(CG_COMM, ierr)
717 write(iout,*) 'REMD synchro at2',itime,'time=',time01-time00
718 if (out1file.or.traj1file) then
719 !deb call mpi_gather(itime,1,mpi_integer,
720 !deb & itime_all,1,mpi_integer,king,
722 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
723 !deb & (itime_all(i),i=1,nodes)
725 !deb imin_itime=itime_all(1)
727 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
729 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
730 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
731 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
732 call mpi_gather(ntwx_cache,1,mpi_integer,&
733 icache_all,1,mpi_integer,king,&
735 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
736 ! & (icache_all(i),i=1,nodes)
737 ii_write=icache_all(1)
739 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
741 ! write(iout,*) "MIN ii_write=",ii_write
744 !time call flush(iout)
746 if(mremdsync .and. mod(itime,nstex).eq.0) then
748 if (me.eq.king .or. .not. out1file) &
749 write(iout,*) 'REMD synchro at1',itime,ntwx_cache,Nprocs,nodes
750 !!!!!!! TRIAL OF MINIM SYNC
751 ! if (me.eq.king) then
753 ! call mpi_isend(itime,1,MPI_INTEGER,i,101, &
754 ! CG_COMM, ireqi(i), ierr)
755 !!d write(iout,*) 'REMD synchro with',i
758 ! call mpi_waitall(nodes-1,ireqi,statusi,ierr)
759 ! call mpi_barrier(CG_COMM, ierr)
763 ! if (me.ne.king) then
764 ! call mpi_recv(itime_master, 1, MPI_INTEGER,&
765 ! 0,101,CG_COMM, status, ierr)
766 ! call mpi_barrier(CG_COMM, ierr)
767 !deb if (out1file.or.traj1file) then
770 write(iout,*) icache_all
772 write(iout,*) "before mpi_gather ntwx_cache"
773 call mpi_gather(ntwx_cache,1,mpi_integer,&
774 icache_all(1),1,mpi_integer,king,& ! CONSULT WITH ADAM
776 write(iout,*) "after mpi_gather ntwx_cache"
779 write(iout,'(a19,8000i8)') ' ntwx_cache',&
780 (icache_all(i),i=1,nodes)
781 ii_write=icache_all(1)
783 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
785 write(iout,*) "MIN ii_write=",ii_write
791 ! Update the time safety limiy
792 if (time001-time00.gt.safety) then
793 safety=time001-time00+600
794 write (iout,*) "****** SAFETY increased to",safety," s"
796 if (ovrtim()) end_of_run=.true.
798 if(synflag.and..not.end_of_run) then
802 write(iout,*) 'REMD before',me,t_bath
804 ! call mpi_gather(t_bath,1,mpi_double_precision,
805 ! & remd_t_bath,1,mpi_double_precision,king,
807 potEcomp(n_ene+1)=t_bath
809 potEcomp(n_ene+2)=iset
810 if (iset.lt.nset) then
814 potEcomp(n_ene+3)=Uconst
821 potEcomp(n_ene+4)=Uconst
825 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
826 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
829 call mpi_gather(elow,1,mpi_double_precision,&
830 elowi,1,mpi_double_precision,king,&
832 call mpi_gather(ehigh,1,mpi_double_precision,&
833 ehighi,1,mpi_double_precision,king,&
838 if (me.eq.king .or. .not. out1file) then
839 write(iout,*) 'REMD gather times=',time03-time01 &
843 if (restart1file) call write1rst(i_index)
846 if (me.eq.king .or. .not. out1file) then
847 write(iout,*) 'REMD writing rst time=',time04-time03
850 if (traj1file) call write1traj
852 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
853 !deb & icache_all,1,mpi_integer,king,
855 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
856 !deb & (icache_all(i),i=1,nodes)
861 if (me.eq.king .or. .not. out1file) then
862 write(iout,*) 'REMD writing traj time=',time05-time04
869 remd_t_bath(i)=remd_ene(n_ene+1,i)
870 iremd_iset(i)=remd_ene(n_ene+2,i)
874 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
876 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
880 write(iout,*) 'REMD exchange temp,ene'
882 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
883 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
887 !-------------------------------------
889 write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
892 write (iout,*) "remd_m(1)",remd_m(1)
894 i=ifirst(iran_num(1,remd_m(1)))
901 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
903 if(i.gt.0.and.nupa(0,i).gt.0) then
905 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
907 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
909 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
911 ! do while (iex.eq.i)
912 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
913 iex=nupa(iran_num(1,int(nupa(0,i))),i)
915 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
917 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
919 ! Swap temperatures between conformations i and iex with recalculating the free energies
920 ! following temperature changes.
921 ene_iex_iex=remd_ene(0,iex)
922 ene_i_i=remd_ene(0,i)
923 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
924 ! & " iex",iex," ene_iex_iex",ene_iex_iex
925 ! write (iout,*) "rescaling weights with temperature",
928 call rescale_weights(remd_t_bath(i))
930 ! write (iout,*) "0,iex",remd_t_bath(i)
931 ! call enerprint(remd_ene(0,iex))
933 call sum_energy(remd_ene(0,iex),.false.)
934 ene_iex_i=remd_ene(0,iex)
935 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
937 ! write (iout,*) "0,i",remd_t_bath(i)
938 ! call enerprint(remd_ene(0,i))
940 call sum_energy(remd_ene(0,i),.false.)
941 ! write (iout,*) "ene_i_i",remd_ene(0,i)
943 ! write (iout,*) "rescaling weights with temperature",
945 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
946 write (iout,*) "ERROR: inconsistent energies:",i,&
947 ene_i_i,remd_ene(0,i)
949 call rescale_weights(remd_t_bath(iex))
951 ! write (iout,*) "0,i",remd_t_bath(iex)
952 ! call enerprint(remd_ene(0,i))
954 call sum_energy(remd_ene(0,i),.false.)
955 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
957 ene_i_iex=remd_ene(0,i)
959 ! write (iout,*) "0,iex",remd_t_bath(iex)
960 call enerprint(remd_ene(0,iex))
962 call sum_energy(remd_ene(0,iex),.false.)
963 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
964 write (iout,*) "ERROR: inconsistent energies:",iex,&
965 ene_iex_iex,remd_ene(0,iex)
967 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
968 ! write (iout,*) "i",i," iex",iex
969 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
970 ! & " ene_i_iex",ene_i_iex,
971 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
973 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
974 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
976 ! write(iout,*) 'delta',delta
977 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
978 ! & (remd_ene(i)-remd_ene(iex))/Rb/
979 ! & (remd_t_bath(i)*remd_t_bath(iex))
981 if (delta .gt. 50.0d0) then
987 else if (delta.lt.-50.0d0) then
996 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
997 xxx=ran_number(0.0d0,1.0d0)
998 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1000 if (delta .gt. xxx) then
1002 remd_t_bath(i)=remd_t_bath(iex)
1003 remd_t_bath(iex)=tmp
1004 remd_ene(0,i)=ene_i_iex
1005 remd_ene(0,iex)=ene_iex_i
1011 ehighi(i)=ehighi(iex)
1018 nupa(k,i)=nupa(k,iex)
1021 ndowna(k,i)=ndowna(k,iex)
1025 if (ifirst(il).eq.i) ifirst(il)=iex
1027 if (nupa(k,il).eq.i) then
1029 elseif (nupa(k,il).eq.iex) then
1034 if (ndowna(k,il).eq.i) then
1036 elseif (ndowna(k,il).eq.iex) then
1042 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1044 i2rep(i-1)=i2rep(iex-1)
1047 ! write(iout,*) 'exchange',i,iex
1048 ! write (iout,'(a8,100i4)') "@ ifirst",
1049 ! & (ifirst(k),k=1,remd_m(1))
1051 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1052 ! & (nupa(k,il),k=1,nupa(0,il))
1053 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1054 ! & (ndowna(k,il),k=1,ndowna(0,il))
1059 remd_ene(0,iex)=ene_iex_iex
1060 remd_ene(0,i)=ene_i_i
1066 !d write (iout,*) "exchange completed"
1070 !d write(iout,*) "########",ii
1072 i_temp=iran_num(1,nrep)
1073 i_mult=iran_num(1,remd_m(i_temp))
1074 i_iset=iran_num(1,nset)
1075 i_mset=iran_num(1,mset(i_iset))
1076 i=i_index(i_temp,i_mult,i_iset,i_mset)
1078 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1081 !d write(iout,*) "i_dir=",i_dir
1083 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1086 i_mult1=iran_num(1,remd_m(i_temp1))
1088 i_mset1=iran_num(1,mset(i_iset1))
1089 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1091 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1094 i_mult1=iran_num(1,remd_m(i_temp1))
1096 i_mset1=iran_num(1,mset(i_iset1))
1097 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1098 econstr_temp_i=remd_ene(20,i)
1099 econstr_temp_iex=remd_ene(20,iex)
1100 remd_ene(20,i)=remd_ene(n_ene+3,i)
1101 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1103 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1106 i_mult1=iran_num(1,remd_m(i_temp1))
1108 i_mset1=iran_num(1,mset(i_iset1))
1109 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1110 econstr_temp_i=remd_ene(20,i)
1111 econstr_temp_iex=remd_ene(20,iex)
1112 remd_ene(20,i)=remd_ene(n_ene+3,i)
1113 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1119 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1122 ! Swap temperatures between conformations i and iex with recalculating the free energies
1123 ! following temperature changes.
1124 ene_iex_iex=remd_ene(0,iex)
1125 ene_i_i=remd_ene(0,i)
1126 !o write (iout,*) "rescaling weights with temperature",
1128 call rescale_weights(remd_t_bath(i))
1130 call sum_energy(remd_ene(0,iex),.false.)
1131 ene_iex_i=remd_ene(0,iex)
1132 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1133 ! call sum_energy(remd_ene(0,i),.false.)
1134 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1135 ! write (iout,*) "rescaling weights with temperature",
1136 ! & remd_t_bath(iex)
1137 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1138 ! write (iout,*) "ERROR: inconsistent energies:",i,
1139 ! & ene_i_i,remd_ene(0,i)
1141 call rescale_weights(remd_t_bath(iex))
1142 call sum_energy(remd_ene(0,i),.false.)
1143 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1144 ene_i_iex=remd_ene(0,i)
1145 ! call sum_energy(remd_ene(0,iex),.false.)
1146 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1147 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1148 ! & ene_iex_iex,remd_ene(0,iex)
1150 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1151 ! write (iout,*) "i",i," iex",iex
1152 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1153 !d & " ene_i_iex",ene_i_iex,
1154 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1155 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1156 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1158 !d write(iout,*) 'delta',delta
1159 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1160 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1161 ! & (remd_t_bath(i)*remd_t_bath(iex))
1162 if (delta .gt. 50.0d0) then
1167 if (i_dir.eq.1.or.i_dir.eq.3) &
1168 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1169 if (i_dir.eq.2.or.i_dir.eq.3) &
1170 iremd_tot_usa(int(i2set(i-1)))= &
1171 iremd_tot_usa(int(i2set(i-1)))+1
1172 xxx=ran_number(0.0d0,1.0d0)
1173 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1174 if (delta .gt. xxx) then
1176 remd_t_bath(i)=remd_t_bath(iex)
1177 remd_t_bath(iex)=tmp
1180 iremd_iset(i)=iremd_iset(iex)
1181 iremd_iset(iex)=itmp
1183 remd_ene(0,i)=ene_i_iex
1184 remd_ene(0,iex)=ene_iex_i
1186 if (i_dir.eq.1.or.i_dir.eq.3) &
1187 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1190 i2rep(i-1)=i2rep(iex-1)
1193 if (i_dir.eq.2.or.i_dir.eq.3) &
1194 iremd_acc_usa(int(i2set(i-1)))= &
1195 iremd_acc_usa(int(i2set(i-1)))+1
1198 i2set(i-1)=i2set(iex-1)
1201 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1202 i_index(i_temp,i_mult,i_iset,i_mset)= &
1203 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1204 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1207 remd_ene(0,iex)=ene_iex_iex
1208 remd_ene(0,i)=ene_i_i
1209 remd_ene(20,iex)=econstr_temp_iex
1210 remd_ene(20,i)=econstr_temp_i
1214 !d do il1=1,mset(il)
1217 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1230 !-------------------------------------
1231 write (iout,*) "NREP",nrep
1233 if(iremd_tot(i).ne.0) &
1234 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1235 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1240 if(iremd_tot_usa(i).ne.0) &
1241 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1242 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1248 !d write (iout,'(a6,100i4)') "ifirst",
1249 !d & (ifirst(i),i=1,remd_m(1))
1251 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1252 !d & (nupa(i,il),i=1,nupa(0,il))
1253 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1254 !d & (ndowna(i,il),i=1,ndowna(0,il))
1259 !d write (iout,*) "Before scatter"
1262 if (me.eq.king) then
1263 write (iout,*) "t_bath before scatter",remd_t_bath
1267 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1268 t_bath,1,mpi_double_precision,king,&
1270 !d write (iout,*) "After scatter"
1273 call mpi_scatter(iremd_iset,1,mpi_integer,&
1274 iset,1,mpi_integer,king,&
1278 if (me.eq.king .or. .not. out1file) then
1279 write(iout,*) 'REMD scatter time=',time07-time06
1283 call mpi_scatter(elowi,1,mpi_double_precision,&
1284 elow,1,mpi_double_precision,king,&
1286 call mpi_scatter(ehighi,1,mpi_double_precision,&
1287 ehigh,1,mpi_double_precision,king,&
1290 call rescale_weights(t_bath)
1291 !o write (iout,*) "Processor",me,
1292 !o & " rescaling weights with temperature",t_bath
1294 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
1296 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
1300 !de write(iout,*) 'REMD after',me,t_bath
1302 if (me.eq.king .or. .not. out1file) then
1303 write(iout,*) 'REMD exchange time 8-0=',time08-time00
1304 write(iout,*) 'REMD exchange time 8-7=',time08-time07
1305 write(iout,*) 'REMD exchange time 7-6=',time07-time06
1306 write(iout,*) 'REMD exchange time 6-5=',time06-time05
1307 write(iout,*) 'REMD exchange time 5-4=',time05-time04
1308 write(iout,*) 'REMD exchange time 4-3=',time04-time03
1309 write(iout,*) 'REMD exchange time 3-2=',time03-time02
1310 write(iout,*) 'REMD exchange time 2-1=',time02-time01
1311 write(iout,*) 'REMD exchange time 1-0=',time01-time00
1317 if (restart1file) then
1318 if (me.eq.king .or. .not. out1file) &
1319 write(iout,*) 'writing restart at the end of run'
1320 call write1rst(i_index)
1323 if (traj1file) call write1traj
1325 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1326 !deb & icache_all,1,mpi_integer,king,
1327 !deb & CG_COMM,ierr)
1328 !deb write(iout,'(a40,8000i8)')
1329 !deb & ' ntwx_cache after traj1file at the end',
1330 !deb & (icache_all(i),i=1,nodes)
1335 t_MD=MPI_Wtime()-tt0
1339 if (me.eq.king .or. .not. out1file) then
1340 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1342 'MD calculations setup:',t_MDsetup,&
1343 'Energy & gradient evaluation:',t_enegrad,&
1344 'Stochastic MD setup:',t_langsetup,&
1345 'Stochastic MD step setup:',t_sdsetup,&
1347 write (iout,'(/28(1h=),a25,27(1h=))') &
1348 ' End of MD calculation '
1350 !el common /przechowalnia/
1351 ! deallocate(d_restart1)
1352 ! deallocate(d_restart2)
1356 end subroutine MREMD
1357 !-----------------------------------------------------------------------------
1358 subroutine write1rst(i_index)
1361 ! implicit real*8 (a-h,o-z)
1362 ! include 'DIMENSIONS'
1364 ! include 'COMMON.MD'
1365 ! include 'COMMON.IOUNITS'
1366 ! include 'COMMON.REMD'
1367 ! include 'COMMON.SETUP'
1368 ! include 'COMMON.CHAIN'
1369 ! include 'COMMON.SBRIDGE'
1370 ! include 'COMMON.INTERACT'
1372 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1373 !el d_restart2(3,2*nres*maxprocs)
1374 real(kind=4) :: r_d(3,0:2*nres)
1375 real(kind=4) :: t5_restart1(5)
1376 integer :: iret,itmp
1377 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1378 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1380 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1381 !el common /przechowalnia/ d_restart1,d_restart2
1382 integer :: i,j,il,il1,ierr,ixdrf
1387 t5_restart1(4)=t_bath
1388 t5_restart1(5)=Uconst
1390 call mpi_gather(t5_restart1,5,mpi_real,&
1391 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1399 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1400 d_restart1,3*2*nres+3,mpi_real,king,&
1412 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1413 d_restart2,3*2*nres+3,mpi_real,king,&
1418 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1420 call xdrfint_(ixdrf, i2rep(i), iret)
1423 call xdrfint_(ixdrf, ifirst(i), iret)
1427 call xdrfint_(ixdrf, nupa(i,il), iret)
1431 call xdrfint_(ixdrf, ndowna(i,il), iret)
1437 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1444 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1451 call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1457 call xdrfint_(ixdrf, nset, iret)
1459 call xdrfint_(ixdrf,mset(i), iret)
1462 call xdrfint_(ixdrf,i2set(i), iret)
1468 itmp=i_index(i,j,il,il1)
1469 call xdrfint_(ixdrf,itmp, iret)
1476 call xdrfclose_(ixdrf, iret)
1478 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1480 call xdrfint(ixdrf, i2rep(i), iret)
1483 call xdrfint(ixdrf, ifirst(i), iret)
1487 call xdrfint(ixdrf, nupa(i,il), iret)
1491 call xdrfint(ixdrf, ndowna(i,il), iret)
1497 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1504 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1511 call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1518 call xdrfint(ixdrf, nset, iret)
1520 call xdrfint(ixdrf,mset(i), iret)
1523 call xdrfint(ixdrf,i2set(i), iret)
1529 itmp=i_index(i,j,il,il1)
1530 call xdrfint(ixdrf,itmp, iret)
1537 call xdrfclose(ixdrf, iret)
1541 end subroutine write1rst
1542 !-----------------------------------------------------------------------------
1543 subroutine write1traj
1545 ! implicit real*8 (a-h,o-z)
1546 ! include 'DIMENSIONS'
1548 ! include 'COMMON.MD'
1549 ! include 'COMMON.IOUNITS'
1550 ! include 'COMMON.REMD'
1551 ! include 'COMMON.SETUP'
1552 ! include 'COMMON.CHAIN'
1553 ! include 'COMMON.SBRIDGE'
1554 ! include 'COMMON.INTERACT'
1556 real(kind=4) :: t5_restart1(5)
1557 integer :: iret,itmp
1558 real(kind=4) :: xcoord(3,2*nres+2),prec
1559 real(kind=4) :: r_qfrag(50),r_qpair(100)
1560 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1561 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1562 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1563 p_uscdiff(100*Nprocs) !(100*maxprocs)
1564 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1565 real(kind=4) :: r_c(3,2*nres+2)
1566 !el common /przechowalnia/ p_c
1568 integer :: i,j,il,ierr,ii,ixdrf
1570 call mpi_bcast(ii_write,1,mpi_integer,&
1574 ! print *,'traj1file',me,ii_write,ntwx_cache
1578 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1580 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1583 ! write (iout,*) "before gather write1traj: from node",ii
1585 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1587 t5_restart1(1)=totT_cache(ii)
1588 t5_restart1(2)=EK_cache(ii)
1589 t5_restart1(3)=potE_cache(ii)
1590 t5_restart1(4)=t_bath_cache(ii)
1591 t5_restart1(5)=Uconst_cache(ii)
1592 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1594 call mpi_gather(t5_restart1,5,mpi_real,&
1595 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1597 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1600 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1601 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1604 r_qfrag(i)=qfrag_cache(i,ii)
1607 r_qpair(i)=qpair_cache(i,ii)
1610 r_utheta(i)=utheta_cache(i,ii)
1611 r_ugamma(i)=ugamma_cache(i,ii)
1612 r_uscdiff(i)=uscdiff_cache(i,ii)
1615 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1616 p_qfrag,nfrag,mpi_real,king,&
1618 call mpi_gather(r_qpair,npair,mpi_real,&
1619 p_qpair,npair,mpi_real,king,&
1621 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1622 p_utheta,nfrag_back,mpi_real,king,&
1624 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1625 p_ugamma,nfrag_back,mpi_real,king,&
1627 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1628 p_uscdiff,nfrag_back,mpi_real,king,&
1632 write (iout,*) "p_qfrag"
1634 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1636 write (iout,*) "p_qpair"
1638 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1644 r_c(j,i)=c_cache(j,i,ii)
1648 call mpi_gather(r_c,3*2*nres,mpi_real,&
1649 p_c,3*2*nres,mpi_real,king,&
1655 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1656 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1657 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1658 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1659 call xdrfint_(ixdrf, nss, iret)
1662 call xdrfint(ixdrf, idssb(j)+nres, iret)
1663 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1665 call xdrfint_(ixdrf, ihpb(j), iret)
1666 call xdrfint_(ixdrf, jhpb(j), iret)
1669 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1670 call xdrfint_(ixdrf, iset_restart1(il), iret)
1672 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1675 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1678 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1679 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1680 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1685 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1690 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1694 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1698 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1699 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1700 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1701 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1702 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1703 call xdrfint(ixdrf, nss, iret)
1706 call xdrfint(ixdrf, idssb(j)+nres, iret)
1707 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1709 call xdrfint(ixdrf, ihpb(j), iret)
1710 call xdrfint(ixdrf, jhpb(j), iret)
1713 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1714 call xdrfint(ixdrf, iset_restart1(il), iret)
1716 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1719 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1722 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1723 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1724 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1729 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1734 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1738 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1744 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1746 if(me.eq.king) call xdrfclose(ixdrf, iret)
1748 do i=1,ntwx_cache-ii_write
1750 totT_cache(i)=totT_cache(ii_write+i)
1751 EK_cache(i)=EK_cache(ii_write+i)
1752 potE_cache(i)=potE_cache(ii_write+i)
1753 t_bath_cache(i)=t_bath_cache(ii_write+i)
1754 Uconst_cache(i)=Uconst_cache(ii_write+i)
1755 iset_cache(i)=iset_cache(ii_write+i)
1758 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1761 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1764 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1765 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1766 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1771 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1775 ntwx_cache=ntwx_cache-ii_write
1777 end subroutine write1traj
1778 !-----------------------------------------------------------------------------
1779 subroutine read1restart(i_index)
1781 ! implicit real*8 (a-h,o-z)
1782 ! include 'DIMENSIONS'
1784 ! include 'COMMON.MD'
1785 ! include 'COMMON.IOUNITS'
1786 ! include 'COMMON.REMD'
1787 ! include 'COMMON.SETUP'
1788 ! include 'COMMON.CHAIN'
1789 ! include 'COMMON.SBRIDGE'
1790 ! include 'COMMON.INTERACT'
1791 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1792 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1793 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1794 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1796 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1797 !el common /przechowalnia/ d_restart1
1798 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1800 write (*,*) "Processor",me," called read1restart"
1803 open(irest2,file=mremd_rst_name,status='unknown')
1804 read(irest2,*,err=334) i
1805 write(iout,*) "Reading old rst in ASCI format"
1807 call read1restart_old
1811 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1814 call xdrfint_(ixdrf, i2rep(i), iret)
1817 call xdrfint_(ixdrf, ifirst(i), iret)
1820 call xdrfint_(ixdrf, nupa(0,il), iret)
1822 call xdrfint_(ixdrf, nupa(i,il), iret)
1825 call xdrfint_(ixdrf, ndowna(0,il), iret)
1827 call xdrfint_(ixdrf, ndowna(i,il), iret)
1832 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1836 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1839 call xdrfint(ixdrf, i2rep(i), iret)
1842 call xdrfint(ixdrf, ifirst(i), iret)
1845 call xdrfint(ixdrf, nupa(0,il), iret)
1847 call xdrfint(ixdrf, nupa(i,il), iret)
1850 call xdrfint(ixdrf, ndowna(0,il), iret)
1852 call xdrfint(ixdrf, ndowna(i,il), iret)
1857 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1862 call mpi_scatter(t_restart1,5,mpi_real,&
1863 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1867 t_bath=t5_restart1(4)
1872 ! read(irest2,'(3e15.5)')
1873 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1876 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1878 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1884 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1885 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1895 ! read(irest2,'(3e15.5)')
1896 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1899 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1901 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1907 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1908 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1919 call xdrfint_(ixdrf, nset, iret)
1921 call xdrfint_(ixdrf,mset(i), iret)
1924 call xdrfint_(ixdrf,i2set(i), iret)
1930 call xdrfint_(ixdrf,itmp, iret)
1931 i_index(i,j,il,il1)=itmp
1939 call xdrfint(ixdrf, nset, iret)
1941 call xdrfint(ixdrf,mset(i), iret)
1944 call xdrfint(ixdrf,i2set(i), iret)
1950 call xdrfint(ixdrf,itmp, iret)
1951 i_index(i,j,il,il1)=itmp
1958 call mpi_scatter(i2set,1,mpi_integer,&
1959 iset,1,mpi_integer,king,&
1964 if(me.eq.king) close(irest2)
1966 end subroutine read1restart
1967 !-----------------------------------------------------------------------------
1968 subroutine read1restart_old
1970 ! implicit real*8 (a-h,o-z)
1971 ! include 'DIMENSIONS'
1973 ! include 'COMMON.MD'
1974 ! include 'COMMON.IOUNITS'
1975 ! include 'COMMON.REMD'
1976 ! include 'COMMON.SETUP'
1977 ! include 'COMMON.CHAIN'
1978 ! include 'COMMON.SBRIDGE'
1979 ! include 'COMMON.INTERACT'
1980 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1981 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1982 !el common /przechowalnia/ d_restart1
1984 integer :: i,j,il,ierr
1987 open(irest2,file=mremd_rst_name,status='unknown')
1988 read (irest2,*) (i2rep(i),i=0,nodes-1)
1989 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1991 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1992 read (irest2,*) ndowna(0,il),&
1993 (ndowna(i,il),i=1,ndowna(0,il))
1996 read(irest2,*) (t_restart1(j,il),j=1,4)
1999 call mpi_scatter(t_restart1,5,mpi_real,&
2000 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2004 t_bath=t5_restart1(4)
2009 read(irest2,'(3e15.5)') &
2010 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2014 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2015 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2025 read(irest2,'(3e15.5)') &
2026 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2030 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2031 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2037 if(me.eq.king) close(irest2)
2039 end subroutine read1restart_old
2040 !----------------------------------------------------------------
2041 subroutine alloc_MREMD_arrays
2043 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
2044 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1,5)) !(ntyp1))
2046 ! common /remdcommon/ in io: read_REMDpar
2047 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2048 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2049 ! common /remdrestart/
2050 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2052 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2053 allocate(ifirst(0:nodes)) !(maxprocs)
2054 allocate(nupa(0:nodes,0:2*nodes))
2055 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2056 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2057 allocate(iset_restart1(nodes)) !(maxprocs)
2058 ! common /traj1cache/
2059 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2060 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2061 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2062 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2063 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2064 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2065 allocate(utheta_cache(nfrag_back,max_cache_traj))
2066 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2067 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2068 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2071 end subroutine alloc_MREMD_arrays
2072 !-----------------------------------------------------------------------------
2073 !-----------------------------------------------------------------------------