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/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))
414 ! write(iout,*) "stdforcp=",stdforcp(i),itype(i,mnum),i
418 if (itype(i,mnum).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum)&
419 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
420 ! write(iout,*) "stdforcsc=",stdforcsc(i),itype(i,mnum),i
423 call rescale_weights(t_bath)
427 t_MDsetup = MPI_Wtime()-tt0
429 t_MDsetup = tcpu()-tt0
432 ! Entering the MD loop
438 if (lang.eq.2 .or. lang.eq.3) then
442 call sd_verlet_p_setup
444 call sd_verlet_ciccotti_setup
448 pfric0_mat(i,j,0)=pfric_mat(i,j)
449 afric0_mat(i,j,0)=afric_mat(i,j)
450 vfric0_mat(i,j,0)=vfric_mat(i,j)
451 prand0_mat(i,j,0)=prand_mat(i,j)
452 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
453 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
458 flag_stoch(i)=.false.
462 "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
464 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
468 else if (lang.eq.1 .or. lang.eq.4) then
472 if (me.eq.king .or. .not. out1file) &
473 write(iout,*) 'Setup time',time00-walltime
476 t_langsetup=MPI_Wtime()-tt0
479 t_langsetup=tcpu()-tt0
486 do while(.not.end_of_run)
489 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
490 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
492 if (lang.gt.0 .and. surfarea .and. &
493 mod(itime,reset_fricmat).eq.0) then
494 if (lang.eq.2 .or. lang.eq.3) then
498 call sd_verlet_p_setup
500 call sd_verlet_ciccotti_setup
504 pfric0_mat(i,j,0)=pfric_mat(i,j)
505 afric0_mat(i,j,0)=afric_mat(i,j)
506 vfric0_mat(i,j,0)=vfric_mat(i,j)
507 prand0_mat(i,j,0)=prand_mat(i,j)
508 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
509 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
514 flag_stoch(i)=.false.
517 else if (lang.eq.1 .or. lang.eq.4) then
520 write (iout,'(a,i10)') &
521 "Friction matrix reset based on surface area, itime",itime
523 if (reset_vel .and. tbf .and. lang.eq.0 &
524 .and. mod(itime,count_reset_vel).eq.0) then
526 if (me.eq.king .or. .not. out1file) &
527 write(iout,'(a,f20.2)') &
528 "Velocities reset to random values, time",totT
531 d_t_old(j,i)=d_t(j,i)
535 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
539 d_t(j,0)=d_t(j,0)-vcm(j)
542 kinetic_T=2.0d0/(dimen3*Rb)*EK
543 scalfac=dsqrt(T_bath/kinetic_T)
544 !d write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
547 d_t_old(j,i)=scalfac*d_t(j,i)
553 ! Time-reversible RESPA algorithm
554 ! (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
555 call RESPA_step(itime)
557 ! Variable time step algorithm.
558 call velverlet_step(itime)
562 call brown_step(itime)
564 print *,"Brown dynamics not here!"
566 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
573 if (mod(itime,ntwe).eq.0) then
575 ! write(iout,*) "before returnbox"
576 ! call etotal(potEcomp)
577 ! call enerprint(potEcomp)
579 ! write(iout,*) "after returnbox"
581 ! call etotal(potEcomp)
583 call enerprint(potEcomp)
586 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
587 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
589 call hairpin(.true.,nharp,iharp)
590 call secondary2(.true.)
591 call pdbout(potE,tytul,ipdb)
597 if (mod(itime,ntwx).eq.0.and.traj1file) then
598 if(ntwx_cache.lt.max_cache_traj_use) then
599 ntwx_cache=ntwx_cache+1
601 if (max_cache_traj_use.ne.1) &
602 print *,itime,"processor ",me," over cache ",ntwx_cache
605 totT_cache(i)=totT_cache(i+1)
606 EK_cache(i)=EK_cache(i+1)
607 potE_cache(i)=potE_cache(i+1)
608 t_bath_cache(i)=t_bath_cache(i+1)
609 Uconst_cache(i)=Uconst_cache(i+1)
610 iset_cache(i)=iset_cache(i+1)
613 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
616 qpair_cache(ii,i)=qpair_cache(ii,i+1)
619 utheta_cache(ii,i)=utheta_cache(ii,i+1)
620 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
621 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
627 c_cache(j,ii,i)=c_cache(j,ii,i+1)
633 totT_cache(ntwx_cache)=totT
634 EK_cache(ntwx_cache)=EK
635 potE_cache(ntwx_cache)=potE
636 t_bath_cache(ntwx_cache)=t_bath
637 Uconst_cache(ntwx_cache)=Uconst
638 iset_cache(ntwx_cache)=iset
641 qfrag_cache(i,ntwx_cache)=qfrag(i)
644 qpair_cache(i,ntwx_cache)=qpair(i)
647 utheta_cache(i,ntwx_cache)=utheta(i)
648 ugamma_cache(i,ntwx_cache)=ugamma(i)
649 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
654 c_cache(j,i,ntwx_cache)=c(j,i)
659 if ((rstcount.eq.1000.or.itime.eq.n_timestep) &
660 .and..not.restart1file) then
663 open(irest1,file=mremd_rst_name,status='unknown')
664 write (irest1,*) "i2rep"
665 write (irest1,*) (i2rep(i),i=0,nodes-1)
666 write (irest1,*) "ifirst"
667 write (irest1,*) (ifirst(i),i=1,remd_m(1))
669 write (irest1,*) "nupa",il
670 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
671 write (irest1,*) "ndowna",il
672 write (irest1,*) ndowna(0,il),&
673 (ndowna(i,il),i=1,ndowna(0,il))
676 write (irest1,*) "nset"
677 write (irest1,*) nset
678 write (irest1,*) "mset"
679 write (irest1,*) (mset(i),i=1,nset)
680 write (irest1,*) "i2set"
681 write (irest1,*) (i2set(i),i=0,nodes-1)
682 write (irest1,*) "i_index"
686 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
694 open(irest2,file=rest2name,status='unknown')
695 write(irest2,*) totT,EK,potE,totE,t_bath
697 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
700 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
703 write (irest2,*) iset
710 ! forced synchronization
711 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king &
712 .and. .not. mremdsync) then
714 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
716 call mpi_recv(itime_master, 1, MPI_INTEGER,&
717 0,101,CG_COMM, status, ierr)
718 call mpi_barrier(CG_COMM, ierr)
719 !deb if (out1file.or.traj1file) then
720 !deb call mpi_gather(itime,1,mpi_integer,
721 !deb & icache_all,1,mpi_integer,king,
724 call mpi_gather(ntwx_cache,1,mpi_integer,&
725 icache_all,1,mpi_integer,king,&
728 write(iout,*) 'REMD synchro at3',itime_master,itime
729 if (itime_master.ge.n_timestep .or. ovrtim()) &
731 !time call flush(iout)
736 if ((mod(itime,nstex).eq.0.and.me.eq.king &
737 .or.end_of_run.and.me.eq.king ) &
738 .and. .not. mremdsync ) then
741 call mpi_isend(itime,1,MPI_INTEGER,i,101, &
742 CG_COMM, ireqi(i), ierr)
743 !d write(iout,*) 'REMD synchro with',i
746 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
747 call mpi_barrier(CG_COMM, ierr)
749 write(iout,*) 'REMD synchro at2',itime,'time=',time01-time00
750 if (out1file.or.traj1file) then
751 !deb call mpi_gather(itime,1,mpi_integer,
752 !deb & itime_all,1,mpi_integer,king,
754 !deb write(iout,'(a19,8000i8)') ' REMD synchro itime',
755 !deb & (itime_all(i),i=1,nodes)
757 !deb imin_itime=itime_all(1)
759 !deb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
761 !deb ii_write=(imin_itime-imin_itime_old)/ntwx
762 !deb imin_itime_old=int(imin_itime/ntwx)*ntwx
763 !deb write(iout,*) imin_itime,imin_itime_old,ii_write
764 call mpi_gather(ntwx_cache,1,mpi_integer,&
765 icache_all,1,mpi_integer,king,&
767 ! write(iout,'(a19,8000i8)') ' ntwx_cache',
768 ! & (icache_all(i),i=1,nodes)
769 ii_write=icache_all(1)
771 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
773 ! write(iout,*) "MIN ii_write=",ii_write
776 !time call flush(iout)
778 if(mremdsync .and. mod(itime,nstex).eq.0) then
780 if (me.eq.king .or. .not. out1file) &
781 write(iout,*) 'REMD synchro at1',itime,ntwx_cache,Nprocs,nodes
782 !!!!!!! TRIAL OF MINIM SYNC
783 ! if (me.eq.king) then
785 ! call mpi_isend(itime,1,MPI_INTEGER,i,101, &
786 ! CG_COMM, ireqi(i), ierr)
787 !!d write(iout,*) 'REMD synchro with',i
790 ! call mpi_waitall(nodes-1,ireqi,statusi,ierr)
791 ! call mpi_barrier(CG_COMM, ierr)
795 ! if (me.ne.king) then
796 ! call mpi_recv(itime_master, 1, MPI_INTEGER,&
797 ! 0,101,CG_COMM, status, ierr)
798 ! call mpi_barrier(CG_COMM, ierr)
799 !deb if (out1file.or.traj1file) then
802 write(iout,*) icache_all
804 write(iout,*) "before mpi_gather ntwx_cache"
805 call mpi_gather(ntwx_cache,1,mpi_integer,&
806 icache_all(1),1,mpi_integer,king,& ! CONSULT WITH ADAM
808 write(iout,*) "after mpi_gather ntwx_cache"
811 write(iout,'(a19,8000i8)') ' ntwx_cache',&
812 (icache_all(i),i=1,nodes)
813 ii_write=icache_all(1)
815 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
817 write(iout,*) "MIN ii_write=",ii_write
823 ! Update the time safety limiy
824 if (time001-time00.gt.safety) then
825 safety=time001-time00+600
826 write (iout,*) "****** SAFETY increased to",safety," s"
828 if (ovrtim()) end_of_run=.true.
830 if(synflag.and..not.end_of_run) then
834 write(iout,*) 'REMD before',me,t_bath
836 ! call mpi_gather(t_bath,1,mpi_double_precision,
837 ! & remd_t_bath,1,mpi_double_precision,king,
839 potEcomp(n_ene+1)=t_bath
841 potEcomp(n_ene+2)=iset
842 if (iset.lt.nset) then
846 potEcomp(n_ene+3)=Uconst
853 potEcomp(n_ene+4)=Uconst
857 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,&
858 remd_ene(0,1),n_ene+5,mpi_double_precision,king,&
861 call mpi_gather(elow,1,mpi_double_precision,&
862 elowi,1,mpi_double_precision,king,&
864 call mpi_gather(ehigh,1,mpi_double_precision,&
865 ehighi,1,mpi_double_precision,king,&
870 if (me.eq.king .or. .not. out1file) then
871 write(iout,*) 'REMD gather times=',time03-time01 &
875 if (restart1file) call write1rst(i_index)
878 if (me.eq.king .or. .not. out1file) then
879 write(iout,*) 'REMD writing rst time=',time04-time03
882 if (traj1file) call write1traj
884 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
885 !deb & icache_all,1,mpi_integer,king,
887 !deb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
888 !deb & (icache_all(i),i=1,nodes)
893 if (me.eq.king .or. .not. out1file) then
894 write(iout,*) 'REMD writing traj time=',time05-time04
901 remd_t_bath(i)=remd_ene(n_ene+1,i)
902 iremd_iset(i)=remd_ene(n_ene+2,i)
906 !o write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
908 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),&
912 write(iout,*) 'REMD exchange temp,ene'
914 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
915 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
919 !-------------------------------------
921 write (iout,*) "Enter exchnge, remd_m",remd_m(1),&
924 write (iout,*) "remd_m(1)",remd_m(1)
926 i=ifirst(iran_num(1,remd_m(1)))
933 write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
935 if(i.gt.0.and.nupa(0,i).gt.0) then
937 ! if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
939 ! & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
941 ! call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
943 ! do while (iex.eq.i)
944 ! write (iout,*) "upper",nupa(int(nupa(0,i)),i)
945 iex=nupa(iran_num(1,int(nupa(0,i))),i)
947 ! write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
949 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
951 ! Swap temperatures between conformations i and iex with recalculating the free energies
952 ! following temperature changes.
953 ene_iex_iex=remd_ene(0,iex)
954 ene_i_i=remd_ene(0,i)
955 ! write (iout,*) "i",i," ene_i_i",ene_i_i,
956 ! & " iex",iex," ene_iex_iex",ene_iex_iex
957 ! write (iout,*) "rescaling weights with temperature",
960 call rescale_weights(remd_t_bath(i))
962 ! write (iout,*) "0,iex",remd_t_bath(i)
963 ! call enerprint(remd_ene(0,iex))
965 call sum_energy(remd_ene(0,iex),.false.)
966 ene_iex_i=remd_ene(0,iex)
967 ! write (iout,*) "ene_iex_i",remd_ene(0,iex)
969 ! write (iout,*) "0,i",remd_t_bath(i)
970 ! call enerprint(remd_ene(0,i))
972 call sum_energy(remd_ene(0,i),.false.)
973 ! write (iout,*) "ene_i_i",remd_ene(0,i)
975 ! write (iout,*) "rescaling weights with temperature",
977 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
978 write (iout,*) "ERROR: inconsistent energies:",i,&
979 ene_i_i,remd_ene(0,i)
981 call rescale_weights(remd_t_bath(iex))
983 ! write (iout,*) "0,i",remd_t_bath(iex)
984 ! call enerprint(remd_ene(0,i))
986 call sum_energy(remd_ene(0,i),.false.)
987 ! write (iout,*) "ene_i_iex",remd_ene(0,i)
989 ene_i_iex=remd_ene(0,i)
991 ! write (iout,*) "0,iex",remd_t_bath(iex)
992 call enerprint(remd_ene(0,iex))
994 call sum_energy(remd_ene(0,iex),.false.)
995 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
996 write (iout,*) "ERROR: inconsistent energies:",iex,&
997 ene_iex_iex,remd_ene(0,iex)
999 ! write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1000 ! write (iout,*) "i",i," iex",iex
1001 ! write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1002 ! & " ene_i_iex",ene_i_iex,
1003 ! & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1005 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1006 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1008 ! write(iout,*) 'delta',delta
1009 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1010 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1011 ! & (remd_t_bath(i)*remd_t_bath(iex))
1013 if (delta .gt. 50.0d0) then
1017 if(isnan(delta))then
1019 else if (delta.lt.-50.0d0) then
1028 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1029 xxx=ran_number(0.0d0,1.0d0)
1030 ! write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1032 if (delta .gt. xxx) then
1034 remd_t_bath(i)=remd_t_bath(iex)
1035 remd_t_bath(iex)=tmp
1036 remd_ene(0,i)=ene_i_iex
1037 remd_ene(0,iex)=ene_iex_i
1043 ehighi(i)=ehighi(iex)
1050 nupa(k,i)=nupa(k,iex)
1053 ndowna(k,i)=ndowna(k,iex)
1057 if (ifirst(il).eq.i) ifirst(il)=iex
1059 if (nupa(k,il).eq.i) then
1061 elseif (nupa(k,il).eq.iex) then
1066 if (ndowna(k,il).eq.i) then
1068 elseif (ndowna(k,il).eq.iex) then
1074 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1076 i2rep(i-1)=i2rep(iex-1)
1079 ! write(iout,*) 'exchange',i,iex
1080 ! write (iout,'(a8,100i4)') "@ ifirst",
1081 ! & (ifirst(k),k=1,remd_m(1))
1083 ! write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1084 ! & (nupa(k,il),k=1,nupa(0,il))
1085 ! write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1086 ! & (ndowna(k,il),k=1,ndowna(0,il))
1091 remd_ene(0,iex)=ene_iex_iex
1092 remd_ene(0,i)=ene_i_i
1098 !d write (iout,*) "exchange completed"
1102 !d write(iout,*) "########",ii
1104 i_temp=iran_num(1,nrep)
1105 i_mult=iran_num(1,remd_m(i_temp))
1106 i_iset=iran_num(1,nset)
1107 i_mset=iran_num(1,mset(i_iset))
1108 i=i_index(i_temp,i_mult,i_iset,i_mset)
1110 !d write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1113 !d write(iout,*) "i_dir=",i_dir
1115 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1118 i_mult1=iran_num(1,remd_m(i_temp1))
1120 i_mset1=iran_num(1,mset(i_iset1))
1121 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1123 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1126 i_mult1=iran_num(1,remd_m(i_temp1))
1128 i_mset1=iran_num(1,mset(i_iset1))
1129 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1130 econstr_temp_i=remd_ene(20,i)
1131 econstr_temp_iex=remd_ene(20,iex)
1132 remd_ene(20,i)=remd_ene(n_ene+3,i)
1133 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1135 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1138 i_mult1=iran_num(1,remd_m(i_temp1))
1140 i_mset1=iran_num(1,mset(i_iset1))
1141 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1142 econstr_temp_i=remd_ene(20,i)
1143 econstr_temp_iex=remd_ene(20,iex)
1144 remd_ene(20,i)=remd_ene(n_ene+3,i)
1145 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1151 !d write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1154 ! Swap temperatures between conformations i and iex with recalculating the free energies
1155 ! following temperature changes.
1156 ene_iex_iex=remd_ene(0,iex)
1157 ene_i_i=remd_ene(0,i)
1158 !o write (iout,*) "rescaling weights with temperature",
1160 call rescale_weights(remd_t_bath(i))
1162 call sum_energy(remd_ene(0,iex),.false.)
1163 ene_iex_i=remd_ene(0,iex)
1164 !d write (iout,*) "ene_iex_i",remd_ene(0,iex)
1165 ! call sum_energy(remd_ene(0,i),.false.)
1166 !d write (iout,*) "ene_i_i",remd_ene(0,i)
1167 ! write (iout,*) "rescaling weights with temperature",
1168 ! & remd_t_bath(iex)
1169 ! if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1170 ! write (iout,*) "ERROR: inconsistent energies:",i,
1171 ! & ene_i_i,remd_ene(0,i)
1173 call rescale_weights(remd_t_bath(iex))
1174 call sum_energy(remd_ene(0,i),.false.)
1175 !d write (iout,*) "ene_i_iex",remd_ene(0,i)
1176 ene_i_iex=remd_ene(0,i)
1177 ! call sum_energy(remd_ene(0,iex),.false.)
1178 ! if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1179 ! write (iout,*) "ERROR: inconsistent energies:",iex,
1180 ! & ene_iex_iex,remd_ene(0,iex)
1182 !d write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1183 ! write (iout,*) "i",i," iex",iex
1184 !d write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1185 !d & " ene_i_iex",ene_i_iex,
1186 !d & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1187 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))- &
1188 (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1190 !d write(iout,*) 'delta',delta
1191 ! delta=(remd_t_bath(i)-remd_t_bath(iex))*
1192 ! & (remd_ene(i)-remd_ene(iex))/Rb/
1193 ! & (remd_t_bath(i)*remd_t_bath(iex))
1194 if (delta .gt. 50.0d0) then
1199 if (i_dir.eq.1.or.i_dir.eq.3) &
1200 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1201 if (i_dir.eq.2.or.i_dir.eq.3) &
1202 iremd_tot_usa(int(i2set(i-1)))= &
1203 iremd_tot_usa(int(i2set(i-1)))+1
1204 xxx=ran_number(0.0d0,1.0d0)
1205 !d write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1206 if (delta .gt. xxx) then
1208 remd_t_bath(i)=remd_t_bath(iex)
1209 remd_t_bath(iex)=tmp
1212 iremd_iset(i)=iremd_iset(iex)
1213 iremd_iset(iex)=itmp
1215 remd_ene(0,i)=ene_i_iex
1216 remd_ene(0,iex)=ene_iex_i
1218 if (i_dir.eq.1.or.i_dir.eq.3) &
1219 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1222 i2rep(i-1)=i2rep(iex-1)
1225 if (i_dir.eq.2.or.i_dir.eq.3) &
1226 iremd_acc_usa(int(i2set(i-1)))= &
1227 iremd_acc_usa(int(i2set(i-1)))+1
1230 i2set(i-1)=i2set(iex-1)
1233 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1234 i_index(i_temp,i_mult,i_iset,i_mset)= &
1235 i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1236 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1239 remd_ene(0,iex)=ene_iex_iex
1240 remd_ene(0,i)=ene_i_i
1241 remd_ene(20,iex)=econstr_temp_iex
1242 remd_ene(20,i)=econstr_temp_i
1246 !d do il1=1,mset(il)
1249 !d write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1262 !-------------------------------------
1263 write (iout,*) "NREP",nrep
1265 if(iremd_tot(i).ne.0) &
1266 write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i) &
1267 ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1272 if(iremd_tot_usa(i).ne.0) &
1273 write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,&
1274 iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1280 !d write (iout,'(a6,100i4)') "ifirst",
1281 !d & (ifirst(i),i=1,remd_m(1))
1283 !d write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1284 !d & (nupa(i,il),i=1,nupa(0,il))
1285 !d write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1286 !d & (ndowna(i,il),i=1,ndowna(0,il))
1291 !d write (iout,*) "Before scatter"
1294 if (me.eq.king) then
1295 write (iout,*) "t_bath before scatter",remd_t_bath
1299 call mpi_scatter(remd_t_bath,1,mpi_double_precision,&
1300 t_bath,1,mpi_double_precision,king,&
1302 !d write (iout,*) "After scatter"
1305 call mpi_scatter(iremd_iset,1,mpi_integer,&
1306 iset,1,mpi_integer,king,&
1310 if (me.eq.king .or. .not. out1file) then
1311 write(iout,*) 'REMD scatter time=',time07-time06
1315 call mpi_scatter(elowi,1,mpi_double_precision,&
1316 elow,1,mpi_double_precision,king,&
1318 call mpi_scatter(ehighi,1,mpi_double_precision,&
1319 ehigh,1,mpi_double_precision,king,&
1322 call rescale_weights(t_bath)
1323 !o write (iout,*) "Processor",me,
1324 !o & " rescaling weights with temperature",t_bath
1326 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
1328 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
1331 !c Compute the standard deviations of stochastic forces for Langevin dynamics
1332 !c if the friction coefficients do not depend on surface area
1333 if (lang.gt.0 .and. .not.surfarea) then
1336 stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum))
1340 if (itype(i,mnum).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum)&
1341 *dsqrt(gamsc(iabs(itype(i,mnum)),mnum))
1345 !de write(iout,*) 'REMD after',me,t_bath
1347 if (me.eq.king .or. .not. out1file) then
1348 write(iout,*) 'REMD exchange time 8-0=',time08-time00
1349 write(iout,*) 'REMD exchange time 8-7=',time08-time07
1350 write(iout,*) 'REMD exchange time 7-6=',time07-time06
1351 write(iout,*) 'REMD exchange time 6-5=',time06-time05
1352 write(iout,*) 'REMD exchange time 5-4=',time05-time04
1353 write(iout,*) 'REMD exchange time 4-3=',time04-time03
1354 write(iout,*) 'REMD exchange time 3-2=',time03-time02
1355 write(iout,*) 'REMD exchange time 2-1=',time02-time01
1356 write(iout,*) 'REMD exchange time 1-0=',time01-time00
1362 if (restart1file) then
1363 if (me.eq.king .or. .not. out1file) &
1364 write(iout,*) 'writing restart at the end of run'
1365 call write1rst(i_index)
1368 if (traj1file) call write1traj
1370 !deb call mpi_gather(ntwx_cache,1,mpi_integer,
1371 !deb & icache_all,1,mpi_integer,king,
1372 !deb & CG_COMM,ierr)
1373 !deb write(iout,'(a40,8000i8)')
1374 !deb & ' ntwx_cache after traj1file at the end',
1375 !deb & (icache_all(i),i=1,nodes)
1380 t_MD=MPI_Wtime()-tt0
1384 if (me.eq.king .or. .not. out1file) then
1385 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') &
1387 'MD calculations setup:',t_MDsetup,&
1388 'Energy & gradient evaluation:',t_enegrad,&
1389 'Stochastic MD setup:',t_langsetup,&
1390 'Stochastic MD step setup:',t_sdsetup,&
1392 write (iout,'(/28(1h=),a25,27(1h=))') &
1393 ' End of MD calculation '
1395 !el common /przechowalnia/
1396 ! deallocate(d_restart1)
1397 ! deallocate(d_restart2)
1401 end subroutine MREMD
1402 !-----------------------------------------------------------------------------
1403 subroutine write1rst(i_index)
1406 ! implicit real*8 (a-h,o-z)
1407 ! include 'DIMENSIONS'
1409 ! include 'COMMON.MD'
1410 ! include 'COMMON.IOUNITS'
1411 ! include 'COMMON.REMD'
1412 ! include 'COMMON.SETUP'
1413 ! include 'COMMON.CHAIN'
1414 ! include 'COMMON.SBRIDGE'
1415 ! include 'COMMON.INTERACT'
1417 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs),&
1418 !el d_restart2(3,2*nres*maxprocs)
1419 real(kind=4) :: r_d(3,0:2*nres)
1420 real(kind=4) :: t5_restart1(5)
1421 integer :: iret,itmp
1422 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1423 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1425 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1426 !el common /przechowalnia/ d_restart1,d_restart2
1427 integer :: i,j,il,il1,ierr,ixdrf
1432 t5_restart1(4)=t_bath
1433 t5_restart1(5)=Uconst
1435 call mpi_gather(t5_restart1,5,mpi_real,&
1436 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1444 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1445 d_restart1,3*2*nres+3,mpi_real,king,&
1457 call mpi_gather(r_d,3*2*nres+3,mpi_real,&
1458 d_restart2,3*2*nres+3,mpi_real,king,&
1463 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1465 call xdrfint_(ixdrf, i2rep(i), iret)
1468 call xdrfint_(ixdrf, ifirst(i), iret)
1472 call xdrfint_(ixdrf, nupa(i,il), iret)
1476 call xdrfint_(ixdrf, ndowna(i,il), iret)
1482 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1489 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1496 call xdrffloat_(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1502 call xdrfint_(ixdrf, nset, iret)
1504 call xdrfint_(ixdrf,mset(i), iret)
1507 call xdrfint_(ixdrf,i2set(i), iret)
1513 itmp=i_index(i,j,il,il1)
1514 call xdrfint_(ixdrf,itmp, iret)
1521 call xdrfclose_(ixdrf, iret)
1523 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1525 call xdrfint(ixdrf, i2rep(i), iret)
1528 call xdrfint(ixdrf, ifirst(i), iret)
1532 call xdrfint(ixdrf, nupa(i,il), iret)
1536 call xdrfint(ixdrf, ndowna(i,il), iret)
1542 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1549 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1556 call xdrffloat(ixdrf, d_restart2(j,i+(2*nres+1)*il), iret)
1563 call xdrfint(ixdrf, nset, iret)
1565 call xdrfint(ixdrf,mset(i), iret)
1568 call xdrfint(ixdrf,i2set(i), iret)
1574 itmp=i_index(i,j,il,il1)
1575 call xdrfint(ixdrf,itmp, iret)
1582 call xdrfclose(ixdrf, iret)
1586 end subroutine write1rst
1587 !-----------------------------------------------------------------------------
1588 subroutine write1traj
1590 ! implicit real*8 (a-h,o-z)
1591 ! include 'DIMENSIONS'
1593 ! include 'COMMON.MD'
1594 ! include 'COMMON.IOUNITS'
1595 ! include 'COMMON.REMD'
1596 ! include 'COMMON.SETUP'
1597 ! include 'COMMON.CHAIN'
1598 ! include 'COMMON.SBRIDGE'
1599 ! include 'COMMON.INTERACT'
1601 real(kind=4) :: t5_restart1(5)
1602 integer :: iret,itmp
1603 real(kind=4) :: xcoord(3,2*nres+2),prec
1604 real(kind=4) :: r_qfrag(50),r_qpair(100)
1605 real(kind=4) :: r_utheta(50),r_ugamma(100),r_uscdiff(100)
1606 real(kind=4) :: p_qfrag(50*Nprocs),p_qpair(100*Nprocs) !(100*maxprocs)
1607 real(kind=4) :: p_utheta(50*Nprocs),p_ugamma(100*Nprocs),&
1608 p_uscdiff(100*Nprocs) !(100*maxprocs)
1609 !el real(kind=4) :: p_c(3,(nres2+2)*maxprocs)
1610 real(kind=4) :: r_c(3,2*nres+2)
1611 !el common /przechowalnia/ p_c
1613 integer :: i,j,il,ierr,ii,ixdrf
1615 call mpi_bcast(ii_write,1,mpi_integer,&
1619 ! print *,'traj1file',me,ii_write,ntwx_cache
1623 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1625 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1628 ! write (iout,*) "before gather write1traj: from node",ii
1630 ! write (iout,*) totT_cache(ii),EK_cache(ii),potE_cache(ii),t_bath_cache(ii),Uconst_cache(ii)
1632 t5_restart1(1)=totT_cache(ii)
1633 t5_restart1(2)=EK_cache(ii)
1634 t5_restart1(3)=potE_cache(ii)
1635 t5_restart1(4)=t_bath_cache(ii)
1636 t5_restart1(5)=Uconst_cache(ii)
1637 ! write (iout,*) "before gather write1traj: from node",ii,t5_restart1(1),t5_restart1(3),t5_restart1(5),t5_restart1(4)
1639 call mpi_gather(t5_restart1,5,mpi_real,&
1640 t_restart1,5,mpi_real,king,CG_COMM,ierr)
1642 ! write (iout,*) "after gather write1traj: from node",il,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1645 call mpi_gather(iset_cache(ii),1,mpi_integer,&
1646 iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1649 r_qfrag(i)=qfrag_cache(i,ii)
1652 r_qpair(i)=qpair_cache(i,ii)
1655 r_utheta(i)=utheta_cache(i,ii)
1656 r_ugamma(i)=ugamma_cache(i,ii)
1657 r_uscdiff(i)=uscdiff_cache(i,ii)
1660 call mpi_gather(r_qfrag,nfrag,mpi_real,&
1661 p_qfrag,nfrag,mpi_real,king,&
1663 call mpi_gather(r_qpair,npair,mpi_real,&
1664 p_qpair,npair,mpi_real,king,&
1666 call mpi_gather(r_utheta,nfrag_back,mpi_real,&
1667 p_utheta,nfrag_back,mpi_real,king,&
1669 call mpi_gather(r_ugamma,nfrag_back,mpi_real,&
1670 p_ugamma,nfrag_back,mpi_real,king,&
1672 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,&
1673 p_uscdiff,nfrag_back,mpi_real,king,&
1677 write (iout,*) "p_qfrag"
1679 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1681 write (iout,*) "p_qpair"
1683 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1689 r_c(j,i)=c_cache(j,i,ii)
1693 call mpi_gather(r_c,3*2*nres,mpi_real,&
1694 p_c,3*2*nres,mpi_real,king,&
1700 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1701 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1702 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1703 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1704 call xdrfint_(ixdrf, nss, iret)
1707 call xdrfint(ixdrf, idssb(j)+nres, iret)
1708 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1710 call xdrfint_(ixdrf, ihpb(j), iret)
1711 call xdrfint_(ixdrf, jhpb(j), iret)
1714 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1715 call xdrfint_(ixdrf, iset_restart1(il), iret)
1717 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1720 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1723 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1724 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1725 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1730 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1735 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1739 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1743 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1744 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1745 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1746 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1747 ! write (iout,*) "write1traj: from node",ii,t_restart1(1,il),t_restart1(3,il),t_restart1(5,il),t_restart1(4,il)
1748 call xdrfint(ixdrf, nss, iret)
1751 call xdrfint(ixdrf, idssb(j)+nres, iret)
1752 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1754 call xdrfint(ixdrf, ihpb(j), iret)
1755 call xdrfint(ixdrf, jhpb(j), iret)
1758 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1759 call xdrfint(ixdrf, iset_restart1(il), iret)
1761 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1764 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1767 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1768 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1769 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1774 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1779 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1783 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1789 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1791 if(me.eq.king) call xdrfclose(ixdrf, iret)
1793 do i=1,ntwx_cache-ii_write
1795 totT_cache(i)=totT_cache(ii_write+i)
1796 EK_cache(i)=EK_cache(ii_write+i)
1797 potE_cache(i)=potE_cache(ii_write+i)
1798 t_bath_cache(i)=t_bath_cache(ii_write+i)
1799 Uconst_cache(i)=Uconst_cache(ii_write+i)
1800 iset_cache(i)=iset_cache(ii_write+i)
1803 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1806 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1809 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1810 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1811 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1816 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1820 ntwx_cache=ntwx_cache-ii_write
1822 end subroutine write1traj
1823 !-----------------------------------------------------------------------------
1824 subroutine read1restart(i_index)
1826 ! implicit real*8 (a-h,o-z)
1827 ! include 'DIMENSIONS'
1829 ! include 'COMMON.MD'
1830 ! include 'COMMON.IOUNITS'
1831 ! include 'COMMON.REMD'
1832 ! include 'COMMON.SETUP'
1833 ! include 'COMMON.CHAIN'
1834 ! include 'COMMON.SBRIDGE'
1835 ! include 'COMMON.INTERACT'
1836 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
1837 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
1838 ! integer(kind=2) :: i_index(Nprocs/4,Nprocs/20,Nprocs/200,Nprocs/200)
1839 integer(kind=2) :: i_index(Nprocs,Nprocs,Nprocs,Nprocs)
1841 !(maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1842 !el common /przechowalnia/ d_restart1
1843 integer :: i,j,il,il1,ierr,itmp,iret,ixdrf
1845 write (*,*) "Processor",me," called read1restart"
1848 open(irest2,file=mremd_rst_name,status='unknown')
1849 read(irest2,*,err=334) i
1850 write(iout,*) "Reading old rst in ASCI format"
1852 call read1restart_old
1856 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1859 call xdrfint_(ixdrf, i2rep(i), iret)
1862 call xdrfint_(ixdrf, ifirst(i), iret)
1865 call xdrfint_(ixdrf, nupa(0,il), iret)
1867 call xdrfint_(ixdrf, nupa(i,il), iret)
1870 call xdrfint_(ixdrf, ndowna(0,il), iret)
1872 call xdrfint_(ixdrf, ndowna(i,il), iret)
1877 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1881 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1884 call xdrfint(ixdrf, i2rep(i), iret)
1887 call xdrfint(ixdrf, ifirst(i), iret)
1890 call xdrfint(ixdrf, nupa(0,il), iret)
1892 call xdrfint(ixdrf, nupa(i,il), iret)
1895 call xdrfint(ixdrf, ndowna(0,il), iret)
1897 call xdrfint(ixdrf, ndowna(i,il), iret)
1902 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1907 call mpi_scatter(t_restart1,5,mpi_real,&
1908 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1912 t_bath=t5_restart1(4)
1917 ! read(irest2,'(3e15.5)')
1918 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1921 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1923 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1929 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1930 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1940 ! read(irest2,'(3e15.5)')
1941 ! & (d_restart1(j,i+2*nres*il),j=1,3)
1944 call xdrffloat_(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1946 call xdrffloat(ixdrf, d_restart1(j,i+(2*nres+1)*il), iret)
1952 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
1953 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
1964 call xdrfint_(ixdrf, nset, iret)
1966 call xdrfint_(ixdrf,mset(i), iret)
1969 call xdrfint_(ixdrf,i2set(i), iret)
1975 call xdrfint_(ixdrf,itmp, iret)
1976 i_index(i,j,il,il1)=itmp
1984 call xdrfint(ixdrf, nset, iret)
1986 call xdrfint(ixdrf,mset(i), iret)
1989 call xdrfint(ixdrf,i2set(i), iret)
1995 call xdrfint(ixdrf,itmp, iret)
1996 i_index(i,j,il,il1)=itmp
2003 call mpi_scatter(i2set,1,mpi_integer,&
2004 iset,1,mpi_integer,king,&
2009 if(me.eq.king) close(irest2)
2011 end subroutine read1restart
2012 !-----------------------------------------------------------------------------
2013 subroutine read1restart_old
2015 ! implicit real*8 (a-h,o-z)
2016 ! include 'DIMENSIONS'
2018 ! include 'COMMON.MD'
2019 ! include 'COMMON.IOUNITS'
2020 ! include 'COMMON.REMD'
2021 ! include 'COMMON.SETUP'
2022 ! include 'COMMON.CHAIN'
2023 ! include 'COMMON.SBRIDGE'
2024 ! include 'COMMON.INTERACT'
2025 !el real(kind=4) :: d_restart1(3,2*nres*maxprocs)
2026 real(kind=4) :: r_d(3,0:2*nres),t5_restart1(5)
2027 !el common /przechowalnia/ d_restart1
2029 integer :: i,j,il,ierr
2032 open(irest2,file=mremd_rst_name,status='unknown')
2033 read (irest2,*) (i2rep(i),i=0,nodes-1)
2034 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2036 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2037 read (irest2,*) ndowna(0,il),&
2038 (ndowna(i,il),i=1,ndowna(0,il))
2041 read(irest2,*) (t_restart1(j,il),j=1,4)
2044 call mpi_scatter(t_restart1,5,mpi_real,&
2045 t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2049 t_bath=t5_restart1(4)
2054 read(irest2,'(3e15.5)') &
2055 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2059 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2060 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2070 read(irest2,'(3e15.5)') &
2071 (d_restart1(j,i+(2*nres+1)*il),j=1,3)
2075 call mpi_scatter(d_restart1,3*2*nres+3,mpi_real,&
2076 r_d,3*2*nres+3,mpi_real,king,CG_COMM,ierr)
2082 if(me.eq.king) close(irest2)
2084 end subroutine read1restart_old
2085 !----------------------------------------------------------------
2086 subroutine alloc_MREMD_arrays
2088 ! if(.not.allocated(mset)) allocate(mset(max0(nset,1)))
2089 if(.not.allocated(stdfsc)) allocate(stdfsc(ntyp1,5)) !(ntyp1))
2091 ! common /remdcommon/ in io: read_REMDpar
2092 ! real(kind=8),dimension(:),allocatable :: remd_t !(maxprocs)
2093 ! integer,dimension(:),allocatable :: remd_m !(maxprocs)
2094 ! common /remdrestart/
2095 if(.not.allocated(i2rep)) allocate(i2rep(0:2*nodes))
2097 allocate(i2set(0:2*nodes)) !(0:maxprocs)
2098 allocate(ifirst(0:nodes)) !(maxprocs)
2099 allocate(nupa(0:nodes,0:2*nodes))
2100 allocate(ndowna(0:nodes,0:2*nodes)) !(0:maxprocs/4,0:maxprocs)
2101 allocate(t_restart1(5,nodes)) !(5,maxprocs)
2102 allocate(iset_restart1(nodes)) !(maxprocs)
2103 ! common /traj1cache/
2104 allocate(totT_cache(max_cache_traj),EK_cache(max_cache_traj))
2105 allocate(potE_cache(max_cache_traj),t_bath_cache(max_cache_traj))
2106 allocate(Uconst_cache(max_cache_traj)) !(max_cache_traj)
2107 allocate(qfrag_cache(nfrag,max_cache_traj)) !(50,max_cache_traj)
2108 allocate(qpair_cache(npair,max_cache_traj)) !(100,max_cache_traj)
2109 allocate(ugamma_cache(nfrag_back,max_cache_traj))
2110 allocate(utheta_cache(nfrag_back,max_cache_traj))
2111 allocate(uscdiff_cache(nfrag_back,max_cache_traj)) !(maxfrag_back,max_cache_traj)
2112 allocate(c_cache(3,2*nres+2,max_cache_traj)) !(3,maxres2+2,max_cache_traj)
2113 allocate(iset_cache(max_cache_traj)) !(max_cache_traj)
2116 end subroutine alloc_MREMD_arrays
2117 !-----------------------------------------------------------------------------
2118 !-----------------------------------------------------------------------------