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,mnum
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)
363 if (lang.gt.0 .and. .not.surfarea) then
366 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
370 if (itype(i,mnum).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum)&
371 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
374 ! print *,'irep',me,t_bath
376 if (me.eq.king .or. .not. out1file) &
377 write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
378 call rescale_weights(t_bath)
382 !------copy MD--------------
383 ! The driver for molecular dynamics subroutines
384 !------------------------------------------------
390 if(me.eq.king.or..not.out1file) &
391 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
397 ! Determine the inverse of the inertia matrix.
398 call setup_MD_matrices
402 if (me.eq.king .or. .not. out1file) &
403 write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
405 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
407 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
410 if (lang.gt.0 .and. .not.surfarea) then
413 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
417 if (itype(i,mnum).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum)&
418 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
421 call rescale_weights(t_bath)
425 t_MDsetup = MPI_Wtime()-tt0
427 t_MDsetup = tcpu()-tt0
430 ! Entering the MD loop
436 if (lang.eq.2 .or. lang.eq.3) then
440 call sd_verlet_p_setup
442 call sd_verlet_ciccotti_setup
446 pfric0_mat(i,j,0)=pfric_mat(i,j)
447 afric0_mat(i,j,0)=afric_mat(i,j)
448 vfric0_mat(i,j,0)=vfric_mat(i,j)
449 prand0_mat(i,j,0)=prand_mat(i,j)
450 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
451 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
456 flag_stoch(i)=.false.
460 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
462 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
466 else if (lang.eq.1 .or. lang.eq.4) then
470 if (me.eq.king .or. .not. out1file) &
471 write(iout,*) 'Setup time',time00-walltime
474 t_langsetup=MPI_Wtime()-tt0
477 t_langsetup=tcpu()-tt0
484 do while(.not.end_of_run)
487 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
488 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
490 if (lang.gt.0 .and. surfarea .and. &
491 mod(itime,reset_fricmat).eq.0) then
492 if (lang.eq.2 .or. lang.eq.3) then
496 call sd_verlet_p_setup
498 call sd_verlet_ciccotti_setup
502 pfric0_mat(i,j,0)=pfric_mat(i,j)
503 afric0_mat(i,j,0)=afric_mat(i,j)
504 vfric0_mat(i,j,0)=vfric_mat(i,j)
505 prand0_mat(i,j,0)=prand_mat(i,j)
506 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
507 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
512 flag_stoch(i)=.false.
515 else if (lang.eq.1 .or. lang.eq.4) then
518 write (iout,'(a,i10)') &
519 "Friction matrix reset based on surface area, itime",itime
521 if (reset_vel .and. tbf .and. lang.eq.0 &
522 .and. mod(itime,count_reset_vel).eq.0) then
524 if (me.eq.king .or. .not. out1file) &
525 write(iout,'(a,f20.2)') &
526 "Velocities reset to random values, time",totT
529 d_t_old(j,i)=d_t(j,i)
533 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
537 d_t(j,0)=d_t(j,0)-vcm(j)
540 kinetic_T=2.0d0/(dimen3*Rb)*EK
541 scalfac=dsqrt(T_bath/kinetic_T)
542 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
545 d_t_old(j,i)=scalfac*d_t(j,i)
551 ! Time-reversible RESPA algorithm
552 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
553 call RESPA_step(itime)
555 ! Variable time step algorithm.
556 call velverlet_step(itime)
560 call brown_step(itime)
562 print *,"Brown dynamics not here!"
564 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
571 if (mod(itime,ntwe).eq.0) then
573 call enerprint(potEcomp)
576 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
577 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
579 call hairpin(.true.,nharp,iharp)
580 call secondary2(.true.)
581 call pdbout(potE,tytul,ipdb)
586 if (mod(itime,ntwx).eq.0.and.traj1file) then
587 if(ntwx_cache.lt.max_cache_traj_use) then
588 ntwx_cache=ntwx_cache+1
590 if (max_cache_traj_use.ne.1) &
591 print *,itime,"processor ",me," over cache ",ntwx_cache
594 totT_cache(i)=totT_cache(i+1)
595 EK_cache(i)=EK_cache(i+1)
596 potE_cache(i)=potE_cache(i+1)
597 t_bath_cache(i)=t_bath_cache(i+1)
598 Uconst_cache(i)=Uconst_cache(i+1)
599 iset_cache(i)=iset_cache(i+1)
602 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
605 qpair_cache(ii,i)=qpair_cache(ii,i+1)
608 utheta_cache(ii,i)=utheta_cache(ii,i+1)
609 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
610 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
616 c_cache(j,ii,i)=c_cache(j,ii,i+1)
622 totT_cache(ntwx_cache)=totT
623 EK_cache(ntwx_cache)=EK
624 potE_cache(ntwx_cache)=potE
625 t_bath_cache(ntwx_cache)=t_bath
626 Uconst_cache(ntwx_cache)=Uconst
627 iset_cache(ntwx_cache)=iset
630 qfrag_cache(i,ntwx_cache)=qfrag(i)
633 qpair_cache(i,ntwx_cache)=qpair(i)
636 utheta_cache(i,ntwx_cache)=utheta(i)
637 ugamma_cache(i,ntwx_cache)=ugamma(i)
638 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
643 c_cache(j,i,ntwx_cache)=c(j,i)
648 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
649 .and..not.restart1file) then
652 open(irest1,file=mremd_rst_name,status='unknown')
653 write (irest1,*) "i2rep"
654 write (irest1,*) (i2rep(i),i=0,nodes-1)
655 write (irest1,*) "ifirst"
656 write (irest1,*) (ifirst(i),i=1,remd_m(1))
658 write (irest1,*) "nupa",il
659 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
660 write (irest1,*) "ndowna",il
661 write (irest1,*) ndowna(0,il),&
662 (ndowna(i,il),i=1,ndowna(0,il))
665 write (irest1,*) "nset"
666 write (irest1,*) nset
667 write (irest1,*) "mset"
668 write (irest1,*) (mset(i),i=1,nset)
669 write (irest1,*) "i2set"
670 write (irest1,*) (i2set(i),i=0,nodes-1)
671 write (irest1,*) "i_index"
675 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
683 open(irest2,file=rest2name,status='unknown')
684 write(irest2,*) totT,EK,potE,totE,t_bath
686 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
689 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
692 write (irest2,*) iset
699 ! forced synchronization
700 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
701 .and. .not. mremdsync) then
703 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
705 call mpi_recv(itime_master, 1, MPI_INTEGER,&
706 0,101,CG_COMM, status, ierr)
707 call mpi_barrier(CG_COMM, ierr)
708 !deb if (out1file.or.traj1file) then
709 !deb call mpi_gather(itime,1,mpi_integer,
710 !deb & icache_all,1,mpi_integer,king,
713 call mpi_gather(ntwx_cache,1,mpi_integer,&
714 icache_all,1,mpi_integer,king,&
717 write(iout,*) 'REMD synchro at3',itime_master,itime
718 if (itime_master.ge.n_timestep .or. ovrtim()) &
720 !time call flush(iout)
725 if ((mod(itime,nstex).eq.0.and.me.eq.king &
726 .or.end_of_run.and.me.eq.king ) &
727 .and. .not. mremdsync ) then
730 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
731 CG_COMM, ireqi(i), ierr)
732 !d write(iout,*) 'REMD synchro with',i
735 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
736 call mpi_barrier(CG_COMM, ierr)
738 write(iout,*) 'REMD synchro at2',itime,'time=',time01-time00
739 if (out1file.or.traj1file) then
740 !deb call mpi_gather(itime,1,mpi_integer,
741 !deb & itime_all,1,mpi_integer,king,
743 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
744 !deb & (itime_all(i),i=1,nodes)
746 !deb imin_itime=itime_all(1)
748 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
750 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
751 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
752 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
753 call mpi_gather(ntwx_cache,1,mpi_integer,&
754 icache_all,1,mpi_integer,king,&
756 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
757 ! & (icache_all(i),i=1,nodes)
758 ii_write=icache_all(1)
760 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
762 ! write(iout,*) "MIN ii_write=",ii_write
765 !time call flush(iout)
767 if(mremdsync .and. mod(itime,nstex).eq.0) then
769 if (me.eq.king .or. .not. out1file) &
770 write(iout,*) 'REMD synchro at1',itime,ntwx_cache,Nprocs,nodes
771 !!!!!!! TRIAL OF MINIM SYNC
772 ! if (me.eq.king) then
774 ! call mpi_isend(itime,1,MPI_INTEGER,i,101, &
775 ! CG_COMM, ireqi(i), ierr)
776 !!d write(iout,*) 'REMD synchro with',i
779 ! call mpi_waitall(nodes-1,ireqi,statusi,ierr)
780 ! call mpi_barrier(CG_COMM, ierr)
784 ! if (me.ne.king) then
785 ! call mpi_recv(itime_master, 1, MPI_INTEGER,&
786 ! 0,101,CG_COMM, status, ierr)
787 ! call mpi_barrier(CG_COMM, ierr)
788 !deb if (out1file.or.traj1file) then
791 write(iout,*) icache_all
793 write(iout,*) "before mpi_gather ntwx_cache"
794 call mpi_gather(ntwx_cache,1,mpi_integer,&
795 icache_all(1),1,mpi_integer,king,& ! CONSULT WITH ADAM
797 write(iout,*) "after mpi_gather ntwx_cache"
800 write(iout,'(a19,8000i8)') ' ntwx_cache',&
801 (icache_all(i),i=1,nodes)
802 ii_write=icache_all(1)
804 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
806 write(iout,*) "MIN ii_write=",ii_write
812 ! Update the time safety limiy
813 if (time001-time00.gt.safety) then
814 safety=time001-time00+600
815 write (iout,*) "****** SAFETY increased to",safety," s"
817 if (ovrtim()) end_of_run=.true.
819 if(synflag.and..not.end_of_run) then
823 write(iout,*) 'REMD before',me,t_bath
825 ! call mpi_gather(t_bath,1,mpi_double_precision,
826 ! & remd_t_bath,1,mpi_double_precision,king,
828 potEcomp(n_ene+1)=t_bath
830 potEcomp(n_ene+2)=iset
831 if (iset.lt.nset) then
835 potEcomp(n_ene+3)=Uconst
842 potEcomp(n_ene+4)=Uconst
846 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
847 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
850 call mpi_gather(elow,1,mpi_double_precision,&
851 elowi,1,mpi_double_precision,king,&
853 call mpi_gather(ehigh,1,mpi_double_precision,&
854 ehighi,1,mpi_double_precision,king,&
859 if (me.eq.king .or. .not. out1file) then
860 write(iout,*) 'REMD gather times=',time03-time01 &
864 if (restart1file) call write1rst(i_index)
867 if (me.eq.king .or. .not. out1file) then
868 write(iout,*) 'REMD writing rst time=',time04-time03
871 if (traj1file) call write1traj
873 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
874 !deb & icache_all,1,mpi_integer,king,
876 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
877 !deb & (icache_all(i),i=1,nodes)
882 if (me.eq.king .or. .not. out1file) then
883 write(iout,*) 'REMD writing traj time=',time05-time04
890 remd_t_bath(i)=remd_ene(n_ene+1,i)
891 iremd_iset(i)=remd_ene(n_ene+2,i)
895 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
897 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
901 write(iout,*) 'REMD exchange temp,ene'
903 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
904 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
908 !-------------------------------------
910 write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
913 write (iout,*) "remd_m(1)",remd_m(1)
915 i=ifirst(iran_num(1,remd_m(1)))
922 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
924 if(i.gt.0.and.nupa(0,i).gt.0) then
926 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
928 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
930 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
932 ! do while (iex.eq.i)
933 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
934 iex=nupa(iran_num(1,int(nupa(0,i))),i)
936 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
938 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
940 ! Swap temperatures between conformations i and iex with recalculating the free energies
941 ! following temperature changes.
942 ene_iex_iex=remd_ene(0,iex)
943 ene_i_i=remd_ene(0,i)
944 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
945 ! & " iex",iex," ene_iex_iex",ene_iex_iex
946 ! write (iout,*) "rescaling weights with temperature",
949 call rescale_weights(remd_t_bath(i))
951 ! write (iout,*) "0,iex",remd_t_bath(i)
952 ! call enerprint(remd_ene(0,iex))
954 call sum_energy(remd_ene(0,iex),.false.)
955 ene_iex_i=remd_ene(0,iex)
956 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
958 ! write (iout,*) "0,i",remd_t_bath(i)
959 ! call enerprint(remd_ene(0,i))
961 call sum_energy(remd_ene(0,i),.false.)
962 ! write (iout,*) "ene_i_i",remd_ene(0,i)
964 ! write (iout,*) "rescaling weights with temperature",
966 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
967 write (iout,*) "ERROR: inconsistent energies:",i,&
968 ene_i_i,remd_ene(0,i)
970 call rescale_weights(remd_t_bath(iex))
972 ! write (iout,*) "0,i",remd_t_bath(iex)
973 ! call enerprint(remd_ene(0,i))
975 call sum_energy(remd_ene(0,i),.false.)
976 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
978 ene_i_iex=remd_ene(0,i)
980 ! write (iout,*) "0,iex",remd_t_bath(iex)
981 call enerprint(remd_ene(0,iex))
983 call sum_energy(remd_ene(0,iex),.false.)
984 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
985 write (iout,*) "ERROR: inconsistent energies:",iex,&
986 ene_iex_iex,remd_ene(0,iex)
988 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
989 ! write (iout,*) "i",i," iex",iex
990 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
991 ! & " ene_i_iex",ene_i_iex,
992 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
994 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
995 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
997 ! write(iout,*) 'delta',delta
998 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
999 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1000 ! & (remd_t_bath(i)*remd_t_bath(iex))
1002 if (delta .gt. 50.0d0) then
1006 if(isnan(delta))then
1008 else if (delta.lt.-50.0d0) then
1017 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1018 xxx=ran_number(0.0d0,1.0d0)
1019 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1021 if (delta .gt. xxx) then
1023 remd_t_bath(i)=remd_t_bath(iex)
1024 remd_t_bath(iex)=tmp
1025 remd_ene(0,i)=ene_i_iex
1026 remd_ene(0,iex)=ene_iex_i
1032 ehighi(i)=ehighi(iex)
1039 nupa(k,i)=nupa(k,iex)
1042 ndowna(k,i)=ndowna(k,iex)
1046 if (ifirst(il).eq.i) ifirst(il)=iex
1048 if (nupa(k,il).eq.i) then
1050 elseif (nupa(k,il).eq.iex) then
1055 if (ndowna(k,il).eq.i) then
1057 elseif (ndowna(k,il).eq.iex) then
1063 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1065 i2rep(i-1)=i2rep(iex-1)
1068 ! write(iout,*) 'exchange',i,iex
1069 ! write (iout,'(a8,100i4)') "@ ifirst",
1070 ! & (ifirst(k),k=1,remd_m(1))
1072 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1073 ! & (nupa(k,il),k=1,nupa(0,il))
1074 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1075 ! & (ndowna(k,il),k=1,ndowna(0,il))
1080 remd_ene(0,iex)=ene_iex_iex
1081 remd_ene(0,i)=ene_i_i
1087 !d write (iout,*) "exchange completed"
1091 !d write(iout,*) "########",ii
1093 i_temp=iran_num(1,nrep)
1094 i_mult=iran_num(1,remd_m(i_temp))
1095 i_iset=iran_num(1,nset)
1096 i_mset=iran_num(1,mset(i_iset))
1097 i=i_index(i_temp,i_mult,i_iset,i_mset)
1099 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1102 !d write(iout,*) "i_dir=",i_dir
1104 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1107 i_mult1=iran_num(1,remd_m(i_temp1))
1109 i_mset1=iran_num(1,mset(i_iset1))
1110 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1112 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1115 i_mult1=iran_num(1,remd_m(i_temp1))
1117 i_mset1=iran_num(1,mset(i_iset1))
1118 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1119 econstr_temp_i=remd_ene(20,i)
1120 econstr_temp_iex=remd_ene(20,iex)
1121 remd_ene(20,i)=remd_ene(n_ene+3,i)
1122 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1124 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1127 i_mult1=iran_num(1,remd_m(i_temp1))
1129 i_mset1=iran_num(1,mset(i_iset1))
1130 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1131 econstr_temp_i=remd_ene(20,i)
1132 econstr_temp_iex=remd_ene(20,iex)
1133 remd_ene(20,i)=remd_ene(n_ene+3,i)
1134 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1140 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1143 ! Swap temperatures between conformations i and iex with recalculating the free energies
1144 ! following temperature changes.
1145 ene_iex_iex=remd_ene(0,iex)
1146 ene_i_i=remd_ene(0,i)
1147 !o write (iout,*) "rescaling weights with temperature",
1149 call rescale_weights(remd_t_bath(i))
1151 call sum_energy(remd_ene(0,iex),.false.)
1152 ene_iex_i=remd_ene(0,iex)
1153 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1154 ! call sum_energy(remd_ene(0,i),.false.)
1155 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1156 ! write (iout,*) "rescaling weights with temperature",
1157 ! & remd_t_bath(iex)
1158 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1159 ! write (iout,*) "ERROR: inconsistent energies:",i,
1160 ! & ene_i_i,remd_ene(0,i)
1162 call rescale_weights(remd_t_bath(iex))
1163 call sum_energy(remd_ene(0,i),.false.)
1164 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1165 ene_i_iex=remd_ene(0,i)
1166 ! call sum_energy(remd_ene(0,iex),.false.)
1167 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1168 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1169 ! & ene_iex_iex,remd_ene(0,iex)
1171 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1172 ! write (iout,*) "i",i," iex",iex
1173 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1174 !d & " ene_i_iex",ene_i_iex,
1175 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1176 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1177 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1179 !d write(iout,*) 'delta',delta
1180 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1181 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1182 ! & (remd_t_bath(i)*remd_t_bath(iex))
1183 if (delta .gt. 50.0d0) then
1188 if (i_dir.eq.1.or.i_dir.eq.3) &
1189 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1190 if (i_dir.eq.2.or.i_dir.eq.3) &
1191 iremd_tot_usa(int(i2set(i-1)))= &
1192 iremd_tot_usa(int(i2set(i-1)))+1
1193 xxx=ran_number(0.0d0,1.0d0)
1194 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1195 if (delta .gt. xxx) then
1197 remd_t_bath(i)=remd_t_bath(iex)
1198 remd_t_bath(iex)=tmp
1201 iremd_iset(i)=iremd_iset(iex)
1202 iremd_iset(iex)=itmp
1204 remd_ene(0,i)=ene_i_iex
1205 remd_ene(0,iex)=ene_iex_i
1207 if (i_dir.eq.1.or.i_dir.eq.3) &
1208 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1211 i2rep(i-1)=i2rep(iex-1)
1214 if (i_dir.eq.2.or.i_dir.eq.3) &
1215 iremd_acc_usa(int(i2set(i-1)))= &
1216 iremd_acc_usa(int(i2set(i-1)))+1
1219 i2set(i-1)=i2set(iex-1)
1222 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1223 i_index(i_temp,i_mult,i_iset,i_mset)= &
1224 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1225 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1228 remd_ene(0,iex)=ene_iex_iex
1229 remd_ene(0,i)=ene_i_i
1230 remd_ene(20,iex)=econstr_temp_iex
1231 remd_ene(20,i)=econstr_temp_i
1235 !d do il1=1,mset(il)
1238 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1251 !-------------------------------------
1252 write (iout,*) "NREP",nrep
1254 if(iremd_tot(i).ne.0) &
1255 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1256 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1261 if(iremd_tot_usa(i).ne.0) &
1262 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1263 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1269 !d write (iout,'(a6,100i4)') "ifirst",
1270 !d & (ifirst(i),i=1,remd_m(1))
1272 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1273 !d & (nupa(i,il),i=1,nupa(0,il))
1274 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1275 !d & (ndowna(i,il),i=1,ndowna(0,il))
1280 !d write (iout,*) "Before scatter"
1283 if (me.eq.king) then
1284 write (iout,*) "t_bath before scatter",remd_t_bath
1288 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1289 t_bath,1,mpi_double_precision,king,&
1291 !d write (iout,*) "After scatter"
1294 call mpi_scatter(iremd_iset,1,mpi_integer,&
1295 iset,1,mpi_integer,king,&
1299 if (me.eq.king .or. .not. out1file) then
1300 write(iout,*) 'REMD scatter time=',time07-time06
1304 call mpi_scatter(elowi,1,mpi_double_precision,&
1305 elow,1,mpi_double_precision,king,&
1307 call mpi_scatter(ehighi,1,mpi_double_precision,&
1308 ehigh,1,mpi_double_precision,king,&
1311 call rescale_weights(t_bath)
1312 !o write (iout,*) "Processor",me,
1313 !o & " rescaling weights with temperature",t_bath
1315 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
1317 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
1320 !c Compute the standard deviations of stochastic forces for Langevin dynamics
1321 !c if the friction coefficients do not depend on surface area
1322 if (lang.gt.0 .and. .not.surfarea) then
1325 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
1329 if (itype(i,mnum).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum)&
1330 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
1334 !de write(iout,*) 'REMD after',me,t_bath
1336 if (me.eq.king .or. .not. out1file) then
1337 write(iout,*) 'REMD exchange time 8-0=',time08-time00
1338 write(iout,*) 'REMD exchange time 8-7=',time08-time07
1339 write(iout,*) 'REMD exchange time 7-6=',time07-time06
1340 write(iout,*) 'REMD exchange time 6-5=',time06-time05
1341 write(iout,*) 'REMD exchange time 5-4=',time05-time04
1342 write(iout,*) 'REMD exchange time 4-3=',time04-time03
1343 write(iout,*) 'REMD exchange time 3-2=',time03-time02
1344 write(iout,*) 'REMD exchange time 2-1=',time02-time01
1345 write(iout,*) 'REMD exchange time 1-0=',time01-time00
1351 if (restart1file) then
1352 if (me.eq.king .or. .not. out1file) &
1353 write(iout,*) 'writing restart at the end of run'
1354 call write1rst(i_index)
1357 if (traj1file) call write1traj
1359 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1360 !deb & icache_all,1,mpi_integer,king,
1361 !deb & CG_COMM,ierr)
1362 !deb write(iout,'(a40,8000i8)')
1363 !deb & ' ntwx_cache after traj1file at the end',
1364 !deb & (icache_all(i),i=1,nodes)
1369 t_MD=MPI_Wtime()-tt0
1373 if (me.eq.king .or. .not. out1file) then
1374 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1376 'MD calculations setup:',t_MDsetup,&
1377 'Energy & gradient evaluation:',t_enegrad,&
1378 'Stochastic MD setup:',t_langsetup,&
1379 'Stochastic MD step setup:',t_sdsetup,&
1381 write (iout,'(/28(1h=),a25,27(1h=))') &
1382 ' End of MD calculation '
1384 !el common /przechowalnia/
1385 ! deallocate(d_restart1)
1386 ! deallocate(d_restart2)
1390 end subroutine MREMD
1391 !-----------------------------------------------------------------------------
1392 subroutine write1rst(i_index)
1395 ! implicit real*8 (a-h,o-z)
1396 ! include 'DIMENSIONS'
1398 ! include 'COMMON.MD'
1399 ! include 'COMMON.IOUNITS'
1400 ! include 'COMMON.REMD'
1401 ! include 'COMMON.SETUP'
1402 ! include 'COMMON.CHAIN'
1403 ! include 'COMMON.SBRIDGE'
1404 ! include 'COMMON.INTERACT'
1406 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1407 !el d_restart2(3,2*nres*maxprocs)
1408 real(kind=4) :: r_d(3,0:2*nres)
1409 real(kind=4) :: t5_restart1(5)
1410 integer :: iret,itmp
1411 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1412 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1414 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1415 !el common /przechowalnia/ d_restart1,d_restart2
1416 integer :: i,j,il,il1,ierr,ixdrf
1421 t5_restart1(4)=t_bath
1422 t5_restart1(5)=Uconst
1424 call mpi_gather(t5_restart1,5,mpi_real,&
1425 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1433 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1434 d_restart1,3*2*nres+3,mpi_real,king,&
1446 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1447 d_restart2,3*2*nres+3,mpi_real,king,&
1452 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1454 call xdrfint_(ixdrf, i2rep(i), iret)
1457 call xdrfint_(ixdrf, ifirst(i), iret)
1461 call xdrfint_(ixdrf, nupa(i,il), iret)
1465 call xdrfint_(ixdrf, ndowna(i,il), iret)
1471 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1478 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1485 call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1491 call xdrfint_(ixdrf, nset, iret)
1493 call xdrfint_(ixdrf,mset(i), iret)
1496 call xdrfint_(ixdrf,i2set(i), iret)
1502 itmp=i_index(i,j,il,il1)
1503 call xdrfint_(ixdrf,itmp, iret)
1510 call xdrfclose_(ixdrf, iret)
1512 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1514 call xdrfint(ixdrf, i2rep(i), iret)
1517 call xdrfint(ixdrf, ifirst(i), iret)
1521 call xdrfint(ixdrf, nupa(i,il), iret)
1525 call xdrfint(ixdrf, ndowna(i,il), iret)
1531 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1538 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1545 call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1552 call xdrfint(ixdrf, nset, iret)
1554 call xdrfint(ixdrf,mset(i), iret)
1557 call xdrfint(ixdrf,i2set(i), iret)
1563 itmp=i_index(i,j,il,il1)
1564 call xdrfint(ixdrf,itmp, iret)
1571 call xdrfclose(ixdrf, iret)
1575 end subroutine write1rst
1576 !-----------------------------------------------------------------------------
1577 subroutine write1traj
1579 ! implicit real*8 (a-h,o-z)
1580 ! include 'DIMENSIONS'
1582 ! include 'COMMON.MD'
1583 ! include 'COMMON.IOUNITS'
1584 ! include 'COMMON.REMD'
1585 ! include 'COMMON.SETUP'
1586 ! include 'COMMON.CHAIN'
1587 ! include 'COMMON.SBRIDGE'
1588 ! include 'COMMON.INTERACT'
1590 real(kind=4) :: t5_restart1(5)
1591 integer :: iret,itmp
1592 real(kind=4) :: xcoord(3,2*nres+2),prec
1593 real(kind=4) :: r_qfrag(50),r_qpair(100)
1594 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1595 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1596 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1597 p_uscdiff(100*Nprocs) !(100*maxprocs)
1598 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1599 real(kind=4) :: r_c(3,2*nres+2)
1600 !el common /przechowalnia/ p_c
1602 integer :: i,j,il,ierr,ii,ixdrf
1604 call mpi_bcast(ii_write,1,mpi_integer,&
1608 ! print *,'traj1file',me,ii_write,ntwx_cache
1612 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1614 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1617 ! write (iout,*) "before gather write1traj: from node",ii
1619 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1621 t5_restart1(1)=totT_cache(ii)
1622 t5_restart1(2)=EK_cache(ii)
1623 t5_restart1(3)=potE_cache(ii)
1624 t5_restart1(4)=t_bath_cache(ii)
1625 t5_restart1(5)=Uconst_cache(ii)
1626 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1628 call mpi_gather(t5_restart1,5,mpi_real,&
1629 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1631 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1634 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1635 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1638 r_qfrag(i)=qfrag_cache(i,ii)
1641 r_qpair(i)=qpair_cache(i,ii)
1644 r_utheta(i)=utheta_cache(i,ii)
1645 r_ugamma(i)=ugamma_cache(i,ii)
1646 r_uscdiff(i)=uscdiff_cache(i,ii)
1649 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1650 p_qfrag,nfrag,mpi_real,king,&
1652 call mpi_gather(r_qpair,npair,mpi_real,&
1653 p_qpair,npair,mpi_real,king,&
1655 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1656 p_utheta,nfrag_back,mpi_real,king,&
1658 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1659 p_ugamma,nfrag_back,mpi_real,king,&
1661 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1662 p_uscdiff,nfrag_back,mpi_real,king,&
1666 write (iout,*) "p_qfrag"
1668 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1670 write (iout,*) "p_qpair"
1672 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1678 r_c(j,i)=c_cache(j,i,ii)
1682 call mpi_gather(r_c,3*2*nres,mpi_real,&
1683 p_c,3*2*nres,mpi_real,king,&
1689 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1690 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1691 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1692 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1693 call xdrfint_(ixdrf, nss, iret)
1696 call xdrfint(ixdrf, idssb(j)+nres, iret)
1697 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1699 call xdrfint_(ixdrf, ihpb(j), iret)
1700 call xdrfint_(ixdrf, jhpb(j), iret)
1703 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1704 call xdrfint_(ixdrf, iset_restart1(il), iret)
1706 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1709 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1712 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1713 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1714 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1719 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1724 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1728 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1732 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1733 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1734 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1735 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1736 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1737 call xdrfint(ixdrf, nss, iret)
1740 call xdrfint(ixdrf, idssb(j)+nres, iret)
1741 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1743 call xdrfint(ixdrf, ihpb(j), iret)
1744 call xdrfint(ixdrf, jhpb(j), iret)
1747 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1748 call xdrfint(ixdrf, iset_restart1(il), iret)
1750 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1753 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1756 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1757 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1758 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1763 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1768 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1772 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1778 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1780 if(me.eq.king) call xdrfclose(ixdrf, iret)
1782 do i=1,ntwx_cache-ii_write
1784 totT_cache(i)=totT_cache(ii_write+i)
1785 EK_cache(i)=EK_cache(ii_write+i)
1786 potE_cache(i)=potE_cache(ii_write+i)
1787 t_bath_cache(i)=t_bath_cache(ii_write+i)
1788 Uconst_cache(i)=Uconst_cache(ii_write+i)
1789 iset_cache(i)=iset_cache(ii_write+i)
1792 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1795 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1798 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1799 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1800 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1805 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1809 ntwx_cache=ntwx_cache-ii_write
1811 end subroutine write1traj
1812 !-----------------------------------------------------------------------------
1813 subroutine read1restart(i_index)
1815 ! implicit real*8 (a-h,o-z)
1816 ! include 'DIMENSIONS'
1818 ! include 'COMMON.MD'
1819 ! include 'COMMON.IOUNITS'
1820 ! include 'COMMON.REMD'
1821 ! include 'COMMON.SETUP'
1822 ! include 'COMMON.CHAIN'
1823 ! include 'COMMON.SBRIDGE'
1824 ! include 'COMMON.INTERACT'
1825 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1826 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1827 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1828 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1830 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1831 !el common /przechowalnia/ d_restart1
1832 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1834 write (*,*) "Processor",me," called read1restart"
1837 open(irest2,file=mremd_rst_name,status='unknown')
1838 read(irest2,*,err=334) i
1839 write(iout,*) "Reading old rst in ASCI format"
1841 call read1restart_old
1845 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1848 call xdrfint_(ixdrf, i2rep(i), iret)
1851 call xdrfint_(ixdrf, ifirst(i), iret)
1854 call xdrfint_(ixdrf, nupa(0,il), iret)
1856 call xdrfint_(ixdrf, nupa(i,il), iret)
1859 call xdrfint_(ixdrf, ndowna(0,il), iret)
1861 call xdrfint_(ixdrf, ndowna(i,il), iret)
1866 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1870 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1873 call xdrfint(ixdrf, i2rep(i), iret)
1876 call xdrfint(ixdrf, ifirst(i), iret)
1879 call xdrfint(ixdrf, nupa(0,il), iret)
1881 call xdrfint(ixdrf, nupa(i,il), iret)
1884 call xdrfint(ixdrf, ndowna(0,il), iret)
1886 call xdrfint(ixdrf, ndowna(i,il), iret)
1891 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1896 call mpi_scatter(t_restart1,5,mpi_real,&
1897 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1901 t_bath=t5_restart1(4)
1906 ! read(irest2,'(3e15.5)')
1907 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1910 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1912 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1918 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1919 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1929 ! read(irest2,'(3e15.5)')
1930 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1933 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1935 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1941 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1942 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1953 call xdrfint_(ixdrf, nset, iret)
1955 call xdrfint_(ixdrf,mset(i), iret)
1958 call xdrfint_(ixdrf,i2set(i), iret)
1964 call xdrfint_(ixdrf,itmp, iret)
1965 i_index(i,j,il,il1)=itmp
1973 call xdrfint(ixdrf, nset, iret)
1975 call xdrfint(ixdrf,mset(i), iret)
1978 call xdrfint(ixdrf,i2set(i), iret)
1984 call xdrfint(ixdrf,itmp, iret)
1985 i_index(i,j,il,il1)=itmp
1992 call mpi_scatter(i2set,1,mpi_integer,&
1993 iset,1,mpi_integer,king,&
1998 if(me.eq.king) close(irest2)
2000 end subroutine read1restart
2001 !-----------------------------------------------------------------------------
2002 subroutine read1restart_old
2004 ! implicit real*8 (a-h,o-z)
2005 ! include 'DIMENSIONS'
2007 ! include 'COMMON.MD'
2008 ! include 'COMMON.IOUNITS'
2009 ! include 'COMMON.REMD'
2010 ! include 'COMMON.SETUP'
2011 ! include 'COMMON.CHAIN'
2012 ! include 'COMMON.SBRIDGE'
2013 ! include 'COMMON.INTERACT'
2014 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
2015 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
2016 !el common /przechowalnia/ d_restart1
2018 integer :: i,j,il,ierr
2021 open(irest2,file=mremd_rst_name,status='unknown')
2022 read (irest2,*) (i2rep(i),i=0,nodes-1)
2023 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2025 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2026 read (irest2,*) ndowna(0,il),&
2027 (ndowna(i,il),i=1,ndowna(0,il))
2030 read(irest2,*) (t_restart1(j,il),j=1,4)
2033 call mpi_scatter(t_restart1,5,mpi_real,&
2034 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2038 t_bath=t5_restart1(4)
2043 read(irest2,'(3e15.5)') &
2044 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2048 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2049 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2059 read(irest2,'(3e15.5)') &
2060 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2064 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2065 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2071 if(me.eq.king) close(irest2)
2073 end subroutine read1restart_old
2074 !----------------------------------------------------------------
2075 subroutine alloc_MREMD_arrays
2077 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
2078 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1,5)) !(ntyp1))
2080 ! common /remdcommon/ in io: read_REMDpar
2081 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2082 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2083 ! common /remdrestart/
2084 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2086 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2087 allocate(ifirst(0:nodes)) !(maxprocs)
2088 allocate(nupa(0:nodes,0:2*nodes))
2089 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2090 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2091 allocate(iset_restart1(nodes)) !(maxprocs)
2092 ! common /traj1cache/
2093 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2094 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2095 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2096 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2097 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2098 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2099 allocate(utheta_cache(nfrag_back,max_cache_traj))
2100 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2101 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2102 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2105 end subroutine alloc_MREMD_arrays
2106 !-----------------------------------------------------------------------------
2107 !-----------------------------------------------------------------------------