5 include 'COMMON.CONTROL'
9 include 'COMMON.LAGRANGE.5diag'
11 include 'COMMON.LAGRANGE'
13 include 'COMMON.QRESTR'
15 include 'COMMON.LANGEVIN'
18 include 'COMMON.LANGEVIN.lang0.5diag'
20 include 'COMMON.LANGEVIN.lang0'
23 include 'COMMON.CHAIN'
24 include 'COMMON.DERIV'
26 include 'COMMON.LOCAL'
27 include 'COMMON.INTERACT'
28 include 'COMMON.IOUNITS'
29 include 'COMMON.NAMES'
30 include 'COMMON.TIME1'
32 include 'COMMON.SETUP'
34 include 'COMMON.HAIRPIN'
35 double precision time00,time01,time02,time03,time04,time05,
36 & time06,time07,time08,time001,tt0
37 double precision scalfac
38 integer i,j,k,il,il1,ii,iex,itmp,i_temp,i_mult,i_iset,i_mset,
39 & i_dir,i_temp1,i_mult1,i_mset1
40 integer ERRCODE,ierr,ierror
41 double precision cm(3),L(3),vcm(3)
42 double precision energia(0:n_ene)
43 double precision remd_t_bath(maxprocs)
44 integer iremd_iset(maxprocs)
46 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
47 double precision remd_ene(0:n_ene+8,maxprocs)
48 integer iremd_acc(maxprocs),iremd_tot(maxprocs)
49 integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
54 integer itime,i_set_temp,itt,itime_master,irr,i_iset1
55 integer nharp,iharp(4,maxres/3)
56 cold integer nup(0:maxprocs),ndown(0:maxprocs)
57 integer rep2i(0:maxprocs),ireqi(maxprocs)
58 integer icache_all(maxprocs)
59 integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
60 logical synflag,end_of_run,file_exist /.false./,ovrtim,first_pass
61 double precision t_bath_temp,delta,ene_iex_iex,ene_i_i,ene_iex_i,
62 & ene_i_iex,xxx,tmp,econstr_temp_iex,econstr_temp_i
64 double precision ran_number
71 if(me.eq.king.or..not.out1file) then
72 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
73 write (iout,*) "NREP=",nrep
77 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
78 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
80 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
82 cd print *,'MREMD',nodes
83 cd print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
84 cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath
88 do il1=1,max0(mset(il),1)
104 if(me.eq.king.or..not.out1file) then
105 write(iout,*) (i2rep(i),i=0,nodes-1)
106 write(iout,*) (i2set(i),i=0,nodes-1)
111 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
118 c print *,'i2rep',me,i2rep(me)
119 c print *,'rep2i',(rep2i(i),i=0,nrep)
121 cold if (i2rep(me).eq.nrep) then
124 cold nup(0)=remd_m(i2rep(me)+1)
125 cold k=rep2i(int(i2rep(me)))+1
132 cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
134 cold if (i2rep(me).eq.1) then
137 cold ndown(0)=remd_m(i2rep(me)-1)
138 cold k=rep2i(i2rep(me)-2)+1
145 cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
148 c write (*,*) "Processor",me," rest",rest,"
149 c & restart1fie",restart1file
150 if(rest.and.restart1file) then
152 & inquire(file=mremd_rst_name,exist=file_exist)
153 cd write (*,*) me," Before broadcast: file_exist",file_exist
154 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
156 cd write (*,*) me," After broadcast: file_exist",file_exist
158 if(me.eq.king.or..not.out1file)
159 & write (iout,*) 'Master is reading restart1file'
160 call read1restart(i_index)
162 if(me.eq.king.or..not.out1file)
163 & write (iout,*) 'WARNING : no restart1file'
166 if(me.eq.king.or..not.out1file) then
167 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
168 write(iout,*) "i_index"
173 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
179 stdfp=dsqrt(2*Rb*t_bath/d_time)
181 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
183 if (lang.gt.0 .and. .not.surfarea) then
185 stdforcp(i)=stdfp*dsqrt(gamp)
188 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
189 & *dsqrt(gamsc(iabs(itype(i))))
196 if (rest.and..not.restart1file)
197 & inquire(file=mremd_rst_name,exist=file_exist)
198 if(.not.file_exist.and.rest.and..not.restart1file)
199 & write(iout,*) 'WARNING : no restart file',mremd_rst_name
200 IF (rest.and.file_exist.and..not.restart1file) THEN
201 write (iout,*) 'Master is reading restart file',
203 open(irest2,file=mremd_rst_name,status='unknown')
205 read (irest2,*) (i2rep(i),i=0,nodes-1)
207 read (irest2,*) (ifirst(i),i=1,remd_m(1))
210 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
212 read (irest2,*) ndowna(0,il),
213 & (ndowna(i,il),i=1,ndowna(0,il))
219 read (irest2,*) (mset(i),i=1,nset)
221 read (irest2,*) (i2set(i),i=0,nodes-1)
226 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
231 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
232 write(iout,*) "i_index"
237 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
246 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
247 write (iout,'(a6,1000i5)') "ifirst",
248 & (ifirst(i),i=1,remd_m(1))
250 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
251 & (nupa(i,il),i=1,nupa(0,il))
252 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
253 & (ndowna(i,il),i=1,ndowna(0,il))
256 stdfp=dsqrt(2*Rb*t_bath/d_time)
258 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
260 if (lang.gt.0 .and. .not.surfarea) then
262 stdforcp(i)=stdfp*dsqrt(gamp)
265 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
266 & *dsqrt(gamsc(iabs(itype(i))))
269 ELSE IF (.not.(rest.and.file_exist)) THEN
275 if (i2rep(il-1).eq.nrep) then
278 nupa(0,il)=remd_m(i2rep(il-1)+1)
279 k=rep2i(int(i2rep(il-1)))+1
285 if (i2rep(il-1).eq.1) then
288 ndowna(0,il)=remd_m(i2rep(il-1)-1)
289 k=rep2i(i2rep(il-1)-2)+1
297 write (iout,'(a6,100i4)') "ifirst",
298 & (ifirst(i),i=1,remd_m(1))
300 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
301 & (nupa(i,il),i=1,nupa(0,il))
302 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
303 & (ndowna(i,il),i=1,ndowna(0,il))
309 c t_bath=retmin+(retmax-retmin)*me/(nodes-1)
310 if(.not.(rest.and.file_exist.and.restart1file)) then
311 if (me .eq. king) then
314 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
316 cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
317 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
322 if(me.eq.king.or..not.out1file)
323 & write(iout,*) me,"iset=",iset,"t_bath=",t_bath
327 stdfp=dsqrt(2*Rb*t_bath/d_time)
329 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
332 c print *,'irep',me,t_bath
334 if (me.eq.king .or. .not. out1file)
335 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
336 call rescale_weights(t_bath)
340 c------copy MD--------------
341 c The driver for molecular dynamics subroutines
342 c------------------------------------------------
348 if(me.eq.king.or..not.out1file)
349 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
355 c Determine the inverse of the inertia matrix.
356 call setup_MD_matrices
360 if (me.eq.king .or. .not. out1file)
361 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
362 stdfp=dsqrt(2*Rb*t_bath/d_time)
364 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
366 call rescale_weights(t_bath)
370 t_MDsetup = MPI_Wtime()-tt0
372 t_MDsetup = tcpu()-tt0
375 c Entering the MD loop
381 if (lang.eq.2 .or. lang.eq.3) then
385 call sd_verlet_p_setup
387 call sd_verlet_ciccotti_setup
391 pfric0_mat(i,j,0)=pfric_mat(i,j)
392 afric0_mat(i,j,0)=afric_mat(i,j)
393 vfric0_mat(i,j,0)=vfric_mat(i,j)
394 prand0_mat(i,j,0)=prand_mat(i,j)
395 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
396 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
401 flag_stoch(i)=.false.
405 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
407 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
411 else if (lang.eq.1 .or. lang.eq.4) then
415 if (me.eq.king .or. .not. out1file)
416 & write(iout,*) 'Setup time',time00-walltime
419 t_langsetup=MPI_Wtime()-tt0
422 t_langsetup=tcpu()-tt0
428 c write (iout,*) "first_pass",first_pass
429 do while(.not.end_of_run)
431 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
432 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
434 if (lang.gt.0 .and. surfarea .and.
435 & mod(itime,reset_fricmat).eq.0) then
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.
459 else if (lang.eq.1 .or. lang.eq.4) then
462 write (iout,'(a,i10)')
463 & "Friction matrix reset based on surface area, itime",itime
465 if (reset_vel .and. tbf .and. lang.eq.0
466 & .and. mod(itime,count_reset_vel).eq.0) then
468 if (me.eq.king .or. .not. out1file)
469 & write(iout,'(a,f20.2)')
470 & "Velocities reset to random values, time",totT
473 d_t_old(j,i)=d_t(j,i)
477 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
481 d_t(j,0)=d_t(j,0)-vcm(j)
484 kinetic_T=2.0d0/(dimen3*Rb)*EK
485 scalfac=dsqrt(T_bath/kinetic_T)
486 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
489 d_t_old(j,i)=scalfac*d_t(j,i)
495 c Time-reversible RESPA algorithm
496 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
497 call RESPA_step(itime)
499 c Variable time step algorithm.
500 call velverlet_step(itime)
504 call brown_step(itime)
506 print *,"Brown dynamics not here!"
508 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
514 if (mod(itime,ntwe).eq.0) call statout(itime)
516 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
517 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
519 call hairpin(.true.,nharp,iharp)
520 call secondary2(.true.)
521 call pdbout(potE,tytul,ipdb)
526 if (mod(itime,ntwx).eq.0.and.traj1file) then
527 if(ntwx_cache.lt.max_cache_traj_use) then
528 ntwx_cache=ntwx_cache+1
530 if (max_cache_traj_use.ne.1)
531 & print *,itime,"processor ",me," over cache ",ntwx_cache
534 totT_cache(i)=totT_cache(i+1)
535 EK_cache(i)=EK_cache(i+1)
536 potE_cache(i)=potE_cache(i+1)
537 t_bath_cache(i)=t_bath_cache(i+1)
538 Uconst_cache(i)=Uconst_cache(i+1)
539 iset_cache(i)=iset_cache(i+1)
542 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
545 qpair_cache(ii,i)=qpair_cache(ii,i+1)
548 utheta_cache(ii,i)=utheta_cache(ii,i+1)
549 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
550 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
556 c_cache(j,ii,i)=c_cache(j,ii,i+1)
562 totT_cache(ntwx_cache)=totT
563 EK_cache(ntwx_cache)=EK
564 potE_cache(ntwx_cache)=potE
565 t_bath_cache(ntwx_cache)=t_bath
566 Uconst_cache(ntwx_cache)=Uconst
567 iset_cache(ntwx_cache)=iset
570 qfrag_cache(i,ntwx_cache)=qfrag(i)
573 qpair_cache(i,ntwx_cache)=qpair(i)
576 utheta_cache(i,ntwx_cache)=utheta(i)
577 ugamma_cache(i,ntwx_cache)=ugamma(i)
578 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
580 C print *,'przed returnbox'
582 C call enerprint(remd_ene(0,i))
585 c_cache(j,i,ntwx_cache)=c(j,i)
590 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
591 & .and..not.restart1file) then
594 open(irest1,file=mremd_rst_name,status='unknown')
595 write (irest1,*) "i2rep"
596 write (irest1,*) (i2rep(i),i=0,nodes-1)
597 write (irest1,*) "ifirst"
598 write (irest1,*) (ifirst(i),i=1,remd_m(1))
600 write (irest1,*) "nupa",il
601 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
602 write (irest1,*) "ndowna",il
603 write (irest1,*) ndowna(0,il),
604 & (ndowna(i,il),i=1,ndowna(0,il))
607 write (irest1,*) "nset"
608 write (irest1,*) nset
609 write (irest1,*) "mset"
610 write (irest1,*) (mset(i),i=1,nset)
611 write (irest1,*) "i2set"
612 write (irest1,*) (i2set(i),i=0,nodes-1)
613 write (irest1,*) "i_index"
617 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
625 open(irest2,file=rest2name,status='unknown')
626 write(irest2,*) totT,EK,potE,totE,t_bath
628 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
631 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
634 write (irest2,*) iset
641 c forced synchronization
642 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
643 & .and. .not. mremdsync) then
645 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
647 call mpi_recv(itime_master, 1, MPI_INTEGER,
648 & 0,101,CG_COMM, status, ierr)
649 call mpi_barrier(CG_COMM, ierr)
650 cdeb if (out1file.or.traj1file) then
651 cdeb call mpi_gather(itime,1,mpi_integer,
652 cdeb & icache_all,1,mpi_integer,king,
655 & call mpi_gather(ntwx_cache,1,mpi_integer,
656 & icache_all,1,mpi_integer,king,
659 & write(iout,*) 'REMD synchro at',itime_master,itime
660 if (itime_master.ge.n_timestep .or. ovrtim())
662 ctime call flush(iout)
667 if ((mod(itime,nstex).eq.0.and.me.eq.king
668 & .or.end_of_run.and.me.eq.king )
669 & .and. .not. mremdsync ) then
672 call mpi_isend(itime,1,MPI_INTEGER,i,101,
673 & CG_COMM, ireqi(i), ierr)
674 cd write(iout,*) 'REMD synchro with',i
677 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
678 call mpi_barrier(CG_COMM, ierr)
680 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
681 if (out1file.or.traj1file) then
682 cdeb call mpi_gather(itime,1,mpi_integer,
683 cdeb & itime_all,1,mpi_integer,king,
685 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
686 cdeb & (itime_all(i),i=1,nodes)
688 cdeb imin_itime=itime_all(1)
690 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
692 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
693 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
694 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
695 call mpi_gather(ntwx_cache,1,mpi_integer,
696 & icache_all,1,mpi_integer,king,
698 c write(iout,'(a19,8000i8)') ' ntwx_cache',
699 c & (icache_all(i),i=1,nodes)
700 ii_write=icache_all(1)
702 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
704 c write(iout,*) "MIN ii_write=",ii_write
707 ctime call flush(iout)
709 if(mremdsync .and. mod(itime,nstex).eq.0) then
711 if (me.eq.king .or. .not. out1file)
712 & write(iout,*) 'REMD synchro at',itime
715 call mpi_gather(ntwx_cache,1,mpi_integer,
716 & icache_all,1,mpi_integer,king,
719 write(iout,'(a19,8000i8)') ' ntwx_cache',
720 & (icache_all(i),i=1,nodes)
721 ii_write=icache_all(1)
723 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
725 write(iout,*) "MIN ii_write=",ii_write
731 c Update the time safety limiy
732 if (time001-time00.gt.safety) then
733 safety=time001-time00+600
734 write (iout,*) "****** SAFETY increased to",safety," s"
736 if (ovrtim()) end_of_run=.true.
738 if(synflag.and..not.end_of_run) then
742 write(iout,*) 'REMD before',me,t_bath
744 c call mpi_gather(t_bath,1,mpi_double_precision,
745 c & remd_t_bath,1,mpi_double_precision,king,
747 potEcomp(n_ene+1)=t_bath
749 potEcomp(n_ene+2)=iset
750 if (iset.lt.nset) then
755 call Econstr_back_qlike
759 potEcomp(n_ene+3)=Uconst+Uconst_back
767 call Econstr_back_qlike
771 potEcomp(n_ene+4)=Uconst+Uconst_back
775 c 9/11/17 Adaptive US
778 write (iout,*) "me ",me," itt",itt
780 if (itt.lt.nrep) then
784 potEcomp(n_ene+5)=Uconst
786 write (iout,*) "t_bath",t_bath_temp,t_bath,
789 if (iset.lt.nset) then
793 potEcomp(n_ene+7)=Uconst
795 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
805 potEcomp(n_ene+6)=Uconst
807 write (iout,*) "t_bath",t_bath_temp,t_bath,
814 potEcomp(n_ene+8)=Uconst
816 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
824 call mpi_gather(potEcomp(0),n_ene+9,mpi_double_precision,
825 & remd_ene(0,1),n_ene+9,mpi_double_precision,king,
828 call mpi_gather(elow,1,mpi_double_precision,
829 & elowi,1,mpi_double_precision,king,
831 call mpi_gather(ehigh,1,mpi_double_precision,
832 & ehighi,1,mpi_double_precision,king,
837 if (me.eq.king .or. .not. out1file) then
838 write(iout,*) 'REMD gather times=',time03-time01
842 if (restart1file) call write1rst(i_index)
845 if (me.eq.king .or. .not. out1file) then
846 write(iout,*) 'REMD writing rst time=',time04-time03
849 if (traj1file) call write1traj
851 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
852 cdeb & icache_all,1,mpi_integer,king,
854 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
855 cdeb & (icache_all(i),i=1,nodes)
860 if (me.eq.king .or. .not. out1file) then
861 write(iout,*) 'REMD writing traj time=',time05-time04
868 remd_t_bath(i)=remd_ene(n_ene+1,i)
869 iremd_iset(i)=remd_ene(n_ene+2,i)
873 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
875 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
879 write(iout,*) 'REMD exchange temp,ene'
881 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
882 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
886 c-------------------------------------
889 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
891 write (iout,*) "remd_m(1)",remd_m(1)
894 i=ifirst(iran_num(1,remd_m(1)))
899 & write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
900 if(i.gt.0.and.nupa(0,i).gt.0) then
902 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
904 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
906 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
908 c do while (iex.eq.i)
909 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
910 iex=nupa(iran_num(1,int(nupa(0,i))),i)
912 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
915 write (iout,*) "i",i,"iex",iex," temperatures",
916 & remd_t_bath(i),remd_t_bath(iex)
919 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
921 c Swap temperatures between conformations i and iex with recalculating the free energies
922 c following temperature changes.
923 ene_iex_iex=remd_ene(0,iex)
924 ene_i_i=remd_ene(0,i)
925 c write (iout,*) "i",i," ene_i_i",ene_i_i,
926 c & " iex",iex," ene_iex_iex",ene_iex_iex
927 c write (iout,*) "rescaling weights with temperature",
930 call rescale_weights(remd_t_bath(i))
932 c write (iout,*) "0,iex",remd_t_bath(i)
933 c call enerprint(remd_ene(0,iex))
935 call sum_energy(remd_ene(0,iex),.false.)
936 ene_iex_i=remd_ene(0,iex)
937 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
939 c write (iout,*) "0,i",remd_t_bath(i)
940 c call enerprint(remd_ene(0,i))
942 call sum_energy(remd_ene(0,i),.false.)
943 c write (iout,*) "ene_i_i",remd_ene(0,i)
945 c write (iout,*) "rescaling weights with temperature",
947 c write (iout,*) "first_pass",first_pass
948 if (.not.first_pass.and.
949 & real(ene_i_i).ne.real(remd_ene(0,i)))
951 write (iout,*) "ERROR: inconsistent energies:",i,
952 & ene_i_i,remd_ene(0,i)
954 call rescale_weights(remd_t_bath(iex))
956 c write (iout,*) "0,i",remd_t_bath(iex)
957 c call enerprint(remd_ene(0,i))
959 call sum_energy(remd_ene(0,i),.false.)
960 c write (iout,*) "ene_i_iex",remd_ene(0,i)
962 ene_i_iex=remd_ene(0,i)
964 c write (iout,*) "0,iex",remd_t_bath(iex)
965 c call enerprint(remd_ene(0,iex))
967 call sum_energy(remd_ene(0,iex),.false.)
968 if (.not.first_pass.and.
969 & real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
970 write (iout,*) "ERROR: inconsistent energies:",iex,
971 & ene_iex_iex,remd_ene(0,iex)
974 write (iout,*) "ene_iex_iex",remd_ene(0,iex)
975 write (iout,*) "i",i," iex",iex
976 write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
977 & " ene_i_iex",ene_i_iex,
978 & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
981 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
982 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
984 c write(iout,*) 'delta',delta
985 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
986 c & (remd_ene(i)-remd_ene(iex))/Rb/
987 c & (remd_t_bath(i)*remd_t_bath(iex))
989 if (delta .gt. 50.0d0) then
995 else if (delta.lt.-50.0d0) then
1004 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1005 xxx=ran_number(0.0d0,1.0d0)
1006 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1008 if (delta .gt. xxx) then
1010 remd_t_bath(i)=remd_t_bath(iex)
1011 remd_t_bath(iex)=tmp
1012 remd_ene(0,i)=ene_i_iex
1013 remd_ene(0,iex)=ene_iex_i
1019 ehighi(i)=ehighi(iex)
1026 nupa(k,i)=nupa(k,iex)
1029 ndowna(k,i)=ndowna(k,iex)
1033 if (ifirst(il).eq.i) ifirst(il)=iex
1035 if (nupa(k,il).eq.i) then
1037 elseif (nupa(k,il).eq.iex) then
1042 if (ndowna(k,il).eq.i) then
1044 elseif (ndowna(k,il).eq.iex) then
1050 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1052 i2rep(i-1)=i2rep(iex-1)
1055 write(iout,*) 'exchange',i,iex
1056 write (iout,'(a8,100i4)') "@ ifirst",
1057 & (ifirst(k),k=1,remd_m(1))
1059 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1060 & (nupa(k,il),k=1,nupa(0,il))
1061 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1062 & (ndowna(k,il),k=1,ndowna(0,il))
1067 remd_ene(0,iex)=ene_iex_iex
1068 remd_ene(0,i)=ene_i_i
1075 cd write (iout,*) "exchange completed"
1079 cd write(iout,*) "########",ii
1081 i_temp=iran_num(1,nrep)
1082 i_mult=iran_num(1,remd_m(i_temp))
1083 i_iset=iran_num(1,nset)
1084 i_mset=iran_num(1,mset(i_iset))
1085 i=i_index(i_temp,i_mult,i_iset,i_mset)
1087 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1090 cd write(iout,*) "i_dir=",i_dir
1092 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1095 i_mult1=iran_num(1,remd_m(i_temp1))
1097 i_mset1=iran_num(1,mset(i_iset1))
1098 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1099 c 9/1/17 AL: Correction; otherwise the restraint energies are mis-assigned
1100 c on failed replica exchange attempt
1101 econstr_temp_i=remd_ene(i_econstr,i)
1102 econstr_temp_iex=remd_ene(i_econstr,iex)
1103 c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
1105 remd_ene(i_econstr,i)=remd_ene(n_ene+5,i)
1106 remd_ene(i_econstr,iex)=remd_ene(n_ene+6,iex)
1108 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1111 i_mult1=iran_num(1,remd_m(i_temp1))
1113 i_mset1=iran_num(1,mset(i_iset1))
1114 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1115 econstr_temp_i=remd_ene(i_econstr,i)
1116 econstr_temp_iex=remd_ene(i_econstr,iex)
1117 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1118 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1120 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1123 i_mult1=iran_num(1,remd_m(i_temp1))
1125 i_mset1=iran_num(1,mset(i_iset1))
1126 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1127 econstr_temp_i=remd_ene(i_econstr,i)
1128 econstr_temp_iex=remd_ene(i_econstr,iex)
1130 remd_ene(i_econstr,i)=remd_ene(n_ene+7,i)
1131 remd_ene(i_econstr,iex)=remd_ene(n_ene+8,iex)
1133 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1134 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1140 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1143 c Swap temperatures between conformations i and iex with recalculating the free energies
1144 c following temperature changes.
1146 write (iout,*) "i_dir",i_dir," i",i," iex",iex,
1147 & " t_bath",remd_t_bath(i),remd_t_bath(iex)
1149 ene_iex_iex=remd_ene(0,iex)
1150 ene_i_i=remd_ene(0,i)
1151 co write (iout,*) "rescaling weights with temperature",
1153 call rescale_weights(remd_t_bath(i))
1155 call sum_energy(remd_ene(0,iex),.false.)
1156 ene_iex_i=remd_ene(0,iex)
1157 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
1158 c call sum_energy(remd_ene(0,i),.false.)
1159 cd write (iout,*) "ene_i_i",remd_ene(0,i)
1160 c write (iout,*) "rescaling weights with temperature",
1161 c & remd_t_bath(iex)
1162 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1163 c write (iout,*) "ERROR: inconsistent energies:",i,
1164 c & ene_i_i,remd_ene(0,i)
1166 call rescale_weights(remd_t_bath(iex))
1167 call sum_energy(remd_ene(0,i),.false.)
1168 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1169 ene_i_iex=remd_ene(0,i)
1170 c call sum_energy(remd_ene(0,iex),.false.)
1171 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1172 c write (iout,*) "ERROR: inconsistent energies:",iex,
1173 c & ene_iex_iex,remd_ene(0,iex)
1175 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1176 c write (iout,*) "i",i," iex",iex
1177 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1178 cd & " ene_i_iex",ene_i_iex,
1179 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1180 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1181 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1183 cd write(iout,*) 'delta',delta
1184 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1185 c & (remd_ene(i)-remd_ene(iex))/Rb/
1186 c & (remd_t_bath(i)*remd_t_bath(iex))
1187 if (delta .gt. 50.0d0) then
1192 if (i_dir.eq.1.or.i_dir.eq.3)
1193 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1194 if (i_dir.eq.2.or.i_dir.eq.3)
1195 & iremd_tot_usa(int(i2set(i-1)))=
1196 & iremd_tot_usa(int(i2set(i-1)))+1
1197 xxx=ran_number(0.0d0,1.0d0)
1198 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1199 if (delta .gt. xxx) then
1201 remd_t_bath(i)=remd_t_bath(iex)
1202 remd_t_bath(iex)=tmp
1205 iremd_iset(i)=iremd_iset(iex)
1206 iremd_iset(iex)=itmp
1208 remd_ene(0,i)=ene_i_iex
1209 remd_ene(0,iex)=ene_iex_i
1211 if (i_dir.eq.1.or.i_dir.eq.3)
1212 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1215 i2rep(i-1)=i2rep(iex-1)
1218 if (i_dir.eq.2.or.i_dir.eq.3)
1219 & iremd_acc_usa(int(i2set(i-1)))=
1220 & iremd_acc_usa(int(i2set(i-1)))+1
1223 i2set(i-1)=i2set(iex-1)
1226 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1227 i_index(i_temp,i_mult,i_iset,i_mset)=
1228 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1229 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1232 remd_ene(0,iex)=ene_iex_iex
1233 remd_ene(0,i)=ene_i_i
1234 remd_ene(i_econstr,iex)=econstr_temp_iex
1235 remd_ene(i_econstr,i)=econstr_temp_i
1239 cd do il1=1,mset(il)
1242 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1255 c-------------------------------------
1256 write (iout,*) "NREP",nrep
1258 if(iremd_tot(i).ne.0)
1259 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1260 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1265 if(iremd_tot_usa(i).ne.0)
1266 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1267 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1273 cd write (iout,'(a6,100i4)') "ifirst",
1274 cd & (ifirst(i),i=1,remd_m(1))
1276 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1277 cd & (nupa(i,il),i=1,nupa(0,il))
1278 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1279 cd & (ndowna(i,il),i=1,ndowna(0,il))
1284 cd write (iout,*) "Before scatter"
1286 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1287 & t_bath,1,mpi_double_precision,king,
1289 cd write (iout,*) "After scatter"
1292 & call mpi_scatter(iremd_iset,1,mpi_integer,
1293 & iset,1,mpi_integer,king,
1297 if (me.eq.king .or. .not. out1file) then
1298 write(iout,*) 'REMD scatter time=',time07-time06
1302 call mpi_scatter(elowi,1,mpi_double_precision,
1303 & elow,1,mpi_double_precision,king,
1305 call mpi_scatter(ehighi,1,mpi_double_precision,
1306 & ehigh,1,mpi_double_precision,king,
1309 call rescale_weights(t_bath)
1310 co write (iout,*) "Processor",me,
1311 co & " rescaling weights with temperature",t_bath
1313 stdfp=dsqrt(2*Rb*t_bath/d_time)
1315 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1317 c Compute the standard deviations of stochastic forces for Langevin dynamics
1318 c if the friction coefficients do not depend on surface area
1319 if (lang.gt.0 .and. .not.surfarea) then
1321 stdforcp(i)=stdfp*dsqrt(gamp)
1324 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1325 & *dsqrt(gamsc(iabs(itype(i))))
1328 cde write(iout,*) 'REMD after',me,t_bath
1330 if (me.eq.king .or. .not. out1file) then
1331 write(iout,*) 'REMD exchange time=',time08-time00
1337 if (restart1file) then
1338 if (me.eq.king .or. .not. out1file)
1339 & write(iout,*) 'writing restart at the end of run'
1340 call write1rst(i_index)
1343 if (traj1file) call write1traj
1345 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1346 cdeb & icache_all,1,mpi_integer,king,
1347 cdeb & CG_COMM,ierr)
1348 cdeb write(iout,'(a40,8000i8)')
1349 cdeb & ' ntwx_cache after traj1file at the end',
1350 cdeb & (icache_all(i),i=1,nodes)
1355 t_MD=MPI_Wtime()-tt0
1359 if (me.eq.king .or. .not. out1file) then
1360 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1362 & 'MD calculations setup:',t_MDsetup,
1363 & 'Energy & gradient evaluation:',t_enegrad,
1364 & 'Stochastic MD setup:',t_langsetup,
1365 & 'Stochastic MD step setup:',t_sdsetup,
1367 write (iout,'(/28(1h=),a25,27(1h=))')
1368 & ' End of MD calculation '
1373 c-----------------------------------------------------------------------
1374 subroutine write1rst(i_index)
1376 include 'DIMENSIONS'
1378 include 'COMMON.CONTROL'
1381 include 'COMMON.LAGRANGE.5diag'
1383 include 'COMMON.LAGRANGE'
1385 include 'COMMON.QRESTR'
1386 include 'COMMON.IOUNITS'
1387 include 'COMMON.REMD'
1388 include 'COMMON.SETUP'
1389 include 'COMMON.CHAIN'
1390 include 'COMMON.SBRIDGE'
1391 include 'COMMON.INTERACT'
1393 real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1),
1394 & d_restart2(3,2*maxres*maxprocs)
1398 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1399 common /przechowalnia/ d_restart1,d_restart2
1400 integer i,j,il1,il,ixdrf
1406 t5_restart1(4)=t_bath
1407 t5_restart1(5)=Uconst
1409 call mpi_gather(t5_restart1,5,mpi_real,
1410 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1418 call mpi_gather(r_d,3*2*nres,mpi_real,
1419 & d_restart1,3*2*nres,mpi_real,king,
1428 call mpi_gather(r_d,3*2*nres,mpi_real,
1429 & d_restart2,3*2*nres,mpi_real,king,
1434 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1436 call xdrfint_(ixdrf, i2rep(i), iret)
1439 call xdrfint_(ixdrf, ifirst(i), iret)
1443 call xdrfint_(ixdrf, nupa(i,il), iret)
1447 call xdrfint_(ixdrf, ndowna(i,il), iret)
1453 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1460 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1467 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1473 call xdrfint_(ixdrf, nset, iret)
1475 call xdrfint_(ixdrf,mset(i), iret)
1478 call xdrfint_(ixdrf,i2set(i), iret)
1484 itmp=i_index(i,j,il,il1)
1485 call xdrfint_(ixdrf,itmp, iret)
1492 call xdrfclose_(ixdrf, iret)
1494 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1496 call xdrfint(ixdrf, i2rep(i), iret)
1499 call xdrfint(ixdrf, ifirst(i), iret)
1503 call xdrfint(ixdrf, nupa(i,il), iret)
1507 call xdrfint(ixdrf, ndowna(i,il), iret)
1513 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1520 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1527 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1534 call xdrfint(ixdrf, nset, iret)
1536 call xdrfint(ixdrf,mset(i), iret)
1539 call xdrfint(ixdrf,i2set(i), iret)
1545 itmp=i_index(i,j,il,il1)
1546 call xdrfint(ixdrf,itmp, iret)
1553 call xdrfclose(ixdrf, iret)
1560 subroutine write1traj
1562 include 'DIMENSIONS'
1565 include 'COMMON.QRESTR'
1566 include 'COMMON.IOUNITS'
1567 include 'COMMON.REMD'
1568 include 'COMMON.SETUP'
1569 include 'COMMON.CHAIN'
1570 include 'COMMON.SBRIDGE'
1571 include 'COMMON.INTERACT'
1575 real xcoord(3,maxres2+2),prec
1576 real r_qfrag(50),r_qpair(100)
1577 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1578 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1579 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1580 & p_uscdiff(100*maxprocs)
1581 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1582 common /przechowalnia/ p_c
1583 integer ii,i,il,j,ixdrf
1586 call mpi_bcast(ii_write,1,mpi_integer,
1587 & king,CG_COMM,ierr)
1590 c print *,'traj1file',me,ii_write,ntwx_cache
1594 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1596 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1599 t5_restart1(1)=totT_cache(ii)
1600 t5_restart1(2)=EK_cache(ii)
1601 t5_restart1(3)=potE_cache(ii)
1602 t5_restart1(4)=t_bath_cache(ii)
1603 t5_restart1(5)=Uconst_cache(ii)
1604 call mpi_gather(t5_restart1,5,mpi_real,
1605 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1607 call mpi_gather(iset_cache(ii),1,mpi_integer,
1608 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1611 r_qfrag(i)=qfrag_cache(i,ii)
1614 r_qpair(i)=qpair_cache(i,ii)
1617 r_utheta(i)=utheta_cache(i,ii)
1618 r_ugamma(i)=ugamma_cache(i,ii)
1619 r_uscdiff(i)=uscdiff_cache(i,ii)
1622 call mpi_gather(r_qfrag,nfrag,mpi_real,
1623 & p_qfrag,nfrag,mpi_real,king,
1625 call mpi_gather(r_qpair,npair,mpi_real,
1626 & p_qpair,npair,mpi_real,king,
1628 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1629 & p_utheta,nfrag_back,mpi_real,king,
1631 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1632 & p_ugamma,nfrag_back,mpi_real,king,
1634 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1635 & p_uscdiff,nfrag_back,mpi_real,king,
1639 write (iout,*) "p_qfrag"
1641 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1643 write (iout,*) "p_qpair"
1645 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1651 r_c(j,i)=c_cache(j,i,ii)
1655 call mpi_gather(r_c,3*2*nres,mpi_real,
1656 & p_c,3*2*nres,mpi_real,king,
1662 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1663 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1664 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1665 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1666 call xdrfint_(ixdrf, nss, iret)
1669 call xdrfint(ixdrf, idssb(j)+nres, iret)
1670 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1672 call xdrfint_(ixdrf, ihpb(j), iret)
1673 call xdrfint_(ixdrf, jhpb(j), iret)
1676 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1677 call xdrfint_(ixdrf, iset_restart1(il), iret)
1679 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1682 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1685 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1686 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1687 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1692 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1697 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1701 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1705 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1706 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1707 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1708 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1709 call xdrfint(ixdrf, nss, iret)
1712 call xdrfint(ixdrf, idssb(j)+nres, iret)
1713 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1715 call xdrfint(ixdrf, ihpb(j), iret)
1716 call xdrfint(ixdrf, jhpb(j), iret)
1719 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1720 call xdrfint(ixdrf, iset_restart1(il), iret)
1722 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1725 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1728 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1729 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1730 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1735 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1740 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1744 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1750 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1752 if(me.eq.king) call xdrfclose(ixdrf, iret)
1754 do i=1,ntwx_cache-ii_write
1756 totT_cache(i)=totT_cache(ii_write+i)
1757 EK_cache(i)=EK_cache(ii_write+i)
1758 potE_cache(i)=potE_cache(ii_write+i)
1759 t_bath_cache(i)=t_bath_cache(ii_write+i)
1760 Uconst_cache(i)=Uconst_cache(ii_write+i)
1761 iset_cache(i)=iset_cache(ii_write+i)
1764 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1767 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1770 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1771 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1772 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1777 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1781 ntwx_cache=ntwx_cache-ii_write
1786 subroutine read1restart(i_index)
1788 include 'DIMENSIONS'
1790 include 'COMMON.CONTROL'
1793 include 'COMMON.LAGRANGE.5diag'
1795 include 'COMMON.LAGRANGE'
1797 include 'COMMON.QRESTR'
1798 include 'COMMON.IOUNITS'
1799 include 'COMMON.REMD'
1800 include 'COMMON.SETUP'
1801 include 'COMMON.CHAIN'
1802 include 'COMMON.SBRIDGE'
1803 include 'COMMON.INTERACT'
1804 real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1),
1807 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1808 common /przechowalnia/ d_restart1
1809 integer i,j,il,il1,ixdrf,iret,itmp
1811 c write (*,*) "Processor",me," called read1restart"
1814 open(irest2,file=mremd_rst_name,status='unknown')
1815 read(irest2,*,err=334) i
1816 write(iout,*) "Reading old rst in ASCI format"
1818 call read1restart_old
1822 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1825 call xdrfint_(ixdrf, i2rep(i), iret)
1828 call xdrfint_(ixdrf, ifirst(i), iret)
1831 call xdrfint_(ixdrf, nupa(0,il), iret)
1833 call xdrfint_(ixdrf, nupa(i,il), iret)
1836 call xdrfint_(ixdrf, ndowna(0,il), iret)
1838 call xdrfint_(ixdrf, ndowna(i,il), iret)
1843 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1847 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1850 call xdrfint(ixdrf, i2rep(i), iret)
1853 call xdrfint(ixdrf, ifirst(i), iret)
1856 call xdrfint(ixdrf, nupa(0,il), iret)
1858 call xdrfint(ixdrf, nupa(i,il), iret)
1861 call xdrfint(ixdrf, ndowna(0,il), iret)
1863 call xdrfint(ixdrf, ndowna(i,il), iret)
1868 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1873 call mpi_scatter(t_restart1,5,mpi_real,
1874 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1878 t_bath=t5_restart1(4)
1883 c read(irest2,'(3e15.5)')
1884 c & (d_restart1(j,i+2*nres*il),j=1,3)
1887 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1889 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1895 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1896 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1906 c read(irest2,'(3e15.5)')
1907 c & (d_restart1(j,i+2*nres*il),j=1,3)
1910 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1912 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1918 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1919 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1930 call xdrfint_(ixdrf, nset, iret)
1932 call xdrfint_(ixdrf,mset(i), iret)
1935 call xdrfint_(ixdrf,i2set(i), iret)
1941 call xdrfint_(ixdrf,itmp, iret)
1942 i_index(i,j,il,il1)=itmp
1950 call xdrfint(ixdrf, nset, iret)
1952 call xdrfint(ixdrf,mset(i), iret)
1955 call xdrfint(ixdrf,i2set(i), iret)
1961 call xdrfint(ixdrf,itmp, iret)
1962 i_index(i,j,il,il1)=itmp
1969 Corrected AL 8/19/2014: each processor needs whole iset array not only its
1971 c call mpi_scatter(i2set,1,mpi_integer,
1972 c & iset,1,mpi_integer,king,
1974 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1981 if(me.eq.king) close(irest2)
1985 subroutine read1restart_old
1987 include 'DIMENSIONS'
1991 include 'COMMON.LAGRANGE.5diag'
1993 include 'COMMON.LAGRANGE'
1995 include 'COMMON.IOUNITS'
1996 include 'COMMON.REMD'
1997 include 'COMMON.SETUP'
1998 include 'COMMON.CHAIN'
1999 include 'COMMON.SBRIDGE'
2000 include 'COMMON.INTERACT'
2001 real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1),
2003 common /przechowalnia/ d_restart1
2007 open(irest2,file=mremd_rst_name,status='unknown')
2008 read (irest2,*) (i2rep(i),i=0,nodes-1)
2009 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2011 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2012 read (irest2,*) ndowna(0,il),
2013 & (ndowna(i,il),i=1,ndowna(0,il))
2016 read(irest2,*) (t_restart1(j,il),j=1,4)
2019 call mpi_scatter(t_restart1,5,mpi_real,
2020 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2024 t_bath=t5_restart1(4)
2029 read(irest2,'(3e15.5)')
2030 & (d_restart1(j,i+2*nres*il),j=1,3)
2034 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2035 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2045 read(irest2,'(3e15.5)')
2046 & (d_restart1(j,i+2*nres*il),j=1,3)
2050 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2051 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2057 if(me.eq.king) close(irest2)
2060 c------------------------------------------