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
416 if (me.eq.king .or. .not. out1file)
417 & write(iout,*) 'Setup time',time00-walltime
420 t_langsetup=MPI_Wtime()-tt0
423 t_langsetup=tcpu()-tt0
429 c write (iout,*) "first_pass",first_pass
430 do while(.not.end_of_run)
432 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
433 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
435 if (lang.gt.0 .and. surfarea .and.
436 & mod(itime,reset_fricmat).eq.0) then
437 if (lang.eq.2 .or. lang.eq.3) then
441 call sd_verlet_p_setup
443 call sd_verlet_ciccotti_setup
447 pfric0_mat(i,j,0)=pfric_mat(i,j)
448 afric0_mat(i,j,0)=afric_mat(i,j)
449 vfric0_mat(i,j,0)=vfric_mat(i,j)
450 prand0_mat(i,j,0)=prand_mat(i,j)
451 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
452 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
457 flag_stoch(i)=.false.
460 else if (lang.eq.1 .or. lang.eq.4) then
463 write (iout,'(a,i10)')
464 & "Friction matrix reset based on surface area, itime",itime
466 if (reset_vel .and. tbf .and. lang.eq.0
467 & .and. mod(itime,count_reset_vel).eq.0) then
469 if (me.eq.king .or. .not. out1file)
470 & write(iout,'(a,f20.2)')
471 & "Velocities reset to random values, time",totT
474 d_t_old(j,i)=d_t(j,i)
478 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
482 d_t(j,0)=d_t(j,0)-vcm(j)
485 kinetic_T=2.0d0/(dimen3*Rb)*EK
486 scalfac=dsqrt(T_bath/kinetic_T)
487 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
490 d_t_old(j,i)=scalfac*d_t(j,i)
496 c Time-reversible RESPA algorithm
497 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
498 call RESPA_step(itime)
500 c Variable time step algorithm.
501 call velverlet_step(itime)
505 call brown_step(itime)
507 print *,"Brown dynamics not here!"
509 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
515 if (mod(itime,ntwe).eq.0) call statout(itime)
517 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
518 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
520 call hairpin(.true.,nharp,iharp)
521 call secondary2(.true.)
522 call pdbout(potE,tytul,ipdb)
527 if (mod(itime,ntwx).eq.0.and.traj1file) then
528 if(ntwx_cache.lt.max_cache_traj_use) then
529 ntwx_cache=ntwx_cache+1
531 if (max_cache_traj_use.ne.1)
532 & print *,itime,"processor ",me," over cache ",ntwx_cache
535 totT_cache(i)=totT_cache(i+1)
536 EK_cache(i)=EK_cache(i+1)
537 potE_cache(i)=potE_cache(i+1)
538 t_bath_cache(i)=t_bath_cache(i+1)
539 Uconst_cache(i)=Uconst_cache(i+1)
540 iset_cache(i)=iset_cache(i+1)
543 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
546 qpair_cache(ii,i)=qpair_cache(ii,i+1)
549 utheta_cache(ii,i)=utheta_cache(ii,i+1)
550 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
551 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
557 c_cache(j,ii,i)=c_cache(j,ii,i+1)
563 totT_cache(ntwx_cache)=totT
564 EK_cache(ntwx_cache)=EK
565 potE_cache(ntwx_cache)=potE
566 t_bath_cache(ntwx_cache)=t_bath
567 Uconst_cache(ntwx_cache)=Uconst
568 iset_cache(ntwx_cache)=iset
571 qfrag_cache(i,ntwx_cache)=qfrag(i)
574 qpair_cache(i,ntwx_cache)=qpair(i)
577 utheta_cache(i,ntwx_cache)=utheta(i)
578 ugamma_cache(i,ntwx_cache)=ugamma(i)
579 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
581 C print *,'przed returnbox'
583 C call enerprint(remd_ene(0,i))
586 c_cache(j,i,ntwx_cache)=c(j,i)
591 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
592 & .and..not.restart1file) then
595 open(irest1,file=mremd_rst_name,status='unknown')
596 write (irest1,*) "i2rep"
597 write (irest1,*) (i2rep(i),i=0,nodes-1)
598 write (irest1,*) "ifirst"
599 write (irest1,*) (ifirst(i),i=1,remd_m(1))
601 write (irest1,*) "nupa",il
602 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
603 write (irest1,*) "ndowna",il
604 write (irest1,*) ndowna(0,il),
605 & (ndowna(i,il),i=1,ndowna(0,il))
608 write (irest1,*) "nset"
609 write (irest1,*) nset
610 write (irest1,*) "mset"
611 write (irest1,*) (mset(i),i=1,nset)
612 write (irest1,*) "i2set"
613 write (irest1,*) (i2set(i),i=0,nodes-1)
614 write (irest1,*) "i_index"
618 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
626 open(irest2,file=rest2name,status='unknown')
627 write(irest2,*) totT,EK,potE,totE,t_bath
629 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
632 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
635 write (irest2,*) iset
642 c forced synchronization
643 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
644 & .and. .not. mremdsync) then
646 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
648 call mpi_recv(itime_master, 1, MPI_INTEGER,
649 & 0,101,CG_COMM, status, ierr)
650 call mpi_barrier(CG_COMM, ierr)
651 cdeb if (out1file.or.traj1file) then
652 cdeb call mpi_gather(itime,1,mpi_integer,
653 cdeb & icache_all,1,mpi_integer,king,
656 & call mpi_gather(ntwx_cache,1,mpi_integer,
657 & icache_all,1,mpi_integer,king,
660 & write(iout,*) 'REMD synchro at',itime_master,itime
661 if (itime_master.ge.n_timestep .or. ovrtim())
663 ctime call flush(iout)
668 if ((mod(itime,nstex).eq.0.and.me.eq.king
669 & .or.end_of_run.and.me.eq.king )
670 & .and. .not. mremdsync ) then
673 call mpi_isend(itime,1,MPI_INTEGER,i,101,
674 & CG_COMM, ireqi(i), ierr)
675 cd write(iout,*) 'REMD synchro with',i
678 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
679 call mpi_barrier(CG_COMM, ierr)
681 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
682 if (out1file.or.traj1file) then
683 cdeb call mpi_gather(itime,1,mpi_integer,
684 cdeb & itime_all,1,mpi_integer,king,
686 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
687 cdeb & (itime_all(i),i=1,nodes)
689 cdeb imin_itime=itime_all(1)
691 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
693 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
694 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
695 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
696 call mpi_gather(ntwx_cache,1,mpi_integer,
697 & icache_all,1,mpi_integer,king,
699 c write(iout,'(a19,8000i8)') ' ntwx_cache',
700 c & (icache_all(i),i=1,nodes)
701 ii_write=icache_all(1)
703 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
705 c write(iout,*) "MIN ii_write=",ii_write
708 ctime call flush(iout)
710 if(mremdsync .and. mod(itime,nstex).eq.0) then
712 if (me.eq.king .or. .not. out1file)
713 & write(iout,*) 'REMD synchro at',itime
716 call mpi_gather(ntwx_cache,1,mpi_integer,
717 & icache_all,1,mpi_integer,king,
720 write(iout,'(a19,8000i8)') ' ntwx_cache',
721 & (icache_all(i),i=1,nodes)
722 ii_write=icache_all(1)
724 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
726 write(iout,*) "MIN ii_write=",ii_write
732 c Update the time safety limiy
733 if (time001-time00.gt.safety) then
734 safety=time001-time00+600
735 write (iout,*) "****** SAFETY increased to",safety," s"
737 if (ovrtim()) end_of_run=.true.
739 if(synflag.and..not.end_of_run) then
743 write(iout,*) 'REMD before',me,t_bath
745 c call mpi_gather(t_bath,1,mpi_double_precision,
746 c & remd_t_bath,1,mpi_double_precision,king,
748 potEcomp(n_ene+1)=t_bath
750 potEcomp(n_ene+2)=iset
751 if (iset.lt.nset) then
756 call Econstr_back_qlike
760 potEcomp(n_ene+3)=Uconst+Uconst_back
768 call Econstr_back_qlike
772 potEcomp(n_ene+4)=Uconst+Uconst_back
776 c 9/11/17 Adaptive US
779 write (iout,*) "me ",me," itt",itt
781 if (itt.lt.nrep) then
785 potEcomp(n_ene+5)=Uconst
787 write (iout,*) "t_bath",t_bath_temp,t_bath,
790 if (iset.lt.nset) then
794 potEcomp(n_ene+7)=Uconst
796 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
806 potEcomp(n_ene+6)=Uconst
808 write (iout,*) "t_bath",t_bath_temp,t_bath,
815 potEcomp(n_ene+8)=Uconst
817 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
825 call mpi_gather(potEcomp(0),n_ene+9,mpi_double_precision,
826 & remd_ene(0,1),n_ene+9,mpi_double_precision,king,
829 call mpi_gather(elow,1,mpi_double_precision,
830 & elowi,1,mpi_double_precision,king,
832 call mpi_gather(ehigh,1,mpi_double_precision,
833 & ehighi,1,mpi_double_precision,king,
838 if (me.eq.king .or. .not. out1file) then
839 write(iout,*) 'REMD gather times=',time03-time01
843 if (restart1file) call write1rst(i_index)
846 if (me.eq.king .or. .not. out1file) then
847 write(iout,*) 'REMD writing rst time=',time04-time03
850 if (traj1file) call write1traj
852 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
853 cdeb & icache_all,1,mpi_integer,king,
855 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
856 cdeb & (icache_all(i),i=1,nodes)
861 if (me.eq.king .or. .not. out1file) then
862 write(iout,*) 'REMD writing traj time=',time05-time04
869 remd_t_bath(i)=remd_ene(n_ene+1,i)
870 iremd_iset(i)=remd_ene(n_ene+2,i)
874 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
876 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
880 write(iout,*) 'REMD exchange temp,ene'
882 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
883 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
887 c-------------------------------------
890 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
892 write (iout,*) "remd_m(1)",remd_m(1)
895 i=ifirst(iran_num(1,remd_m(1)))
900 & write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
901 if(i.gt.0.and.nupa(0,i).gt.0) then
903 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
905 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
907 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
909 c do while (iex.eq.i)
910 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
911 iex=nupa(iran_num(1,int(nupa(0,i))),i)
913 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
916 write (iout,*) "i",i,"iex",iex," temperatures",
917 & remd_t_bath(i),remd_t_bath(iex)
920 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
922 c Swap temperatures between conformations i and iex with recalculating the free energies
923 c following temperature changes.
924 ene_iex_iex=remd_ene(0,iex)
925 ene_i_i=remd_ene(0,i)
926 c write (iout,*) "i",i," ene_i_i",ene_i_i,
927 c & " iex",iex," ene_iex_iex",ene_iex_iex
928 c write (iout,*) "rescaling weights with temperature",
931 call rescale_weights(remd_t_bath(i))
933 c write (iout,*) "0,iex",remd_t_bath(i)
934 c call enerprint(remd_ene(0,iex))
936 call sum_energy(remd_ene(0,iex),.false.)
937 ene_iex_i=remd_ene(0,iex)
938 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
940 c write (iout,*) "0,i",remd_t_bath(i)
941 c call enerprint(remd_ene(0,i))
943 call sum_energy(remd_ene(0,i),.false.)
944 c write (iout,*) "ene_i_i",remd_ene(0,i)
946 c write (iout,*) "rescaling weights with temperature",
948 c write (iout,*) "first_pass",first_pass
949 if (.not.first_pass.and.
950 & real(ene_i_i).ne.real(remd_ene(0,i)))
952 write (iout,*) "ERROR: inconsistent energies:",i,
953 & ene_i_i,remd_ene(0,i)
955 call rescale_weights(remd_t_bath(iex))
957 c write (iout,*) "0,i",remd_t_bath(iex)
958 c call enerprint(remd_ene(0,i))
960 call sum_energy(remd_ene(0,i),.false.)
961 c write (iout,*) "ene_i_iex",remd_ene(0,i)
963 ene_i_iex=remd_ene(0,i)
965 c write (iout,*) "0,iex",remd_t_bath(iex)
966 c call enerprint(remd_ene(0,iex))
968 call sum_energy(remd_ene(0,iex),.false.)
969 if (.not.first_pass.and.
970 & real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
971 write (iout,*) "ERROR: inconsistent energies:",iex,
972 & ene_iex_iex,remd_ene(0,iex)
975 write (iout,*) "ene_iex_iex",remd_ene(0,iex)
976 write (iout,*) "i",i," iex",iex
977 write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
978 & " ene_i_iex",ene_i_iex,
979 & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
982 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
983 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
985 c write(iout,*) 'delta',delta
986 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
987 c & (remd_ene(i)-remd_ene(iex))/Rb/
988 c & (remd_t_bath(i)*remd_t_bath(iex))
990 if (delta .gt. 50.0d0) then
996 else if (delta.lt.-50.0d0) then
1005 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1006 xxx=ran_number(0.0d0,1.0d0)
1007 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1009 if (delta .gt. xxx) then
1011 remd_t_bath(i)=remd_t_bath(iex)
1012 remd_t_bath(iex)=tmp
1013 remd_ene(0,i)=ene_i_iex
1014 remd_ene(0,iex)=ene_iex_i
1020 ehighi(i)=ehighi(iex)
1027 nupa(k,i)=nupa(k,iex)
1030 ndowna(k,i)=ndowna(k,iex)
1034 if (ifirst(il).eq.i) ifirst(il)=iex
1036 if (nupa(k,il).eq.i) then
1038 elseif (nupa(k,il).eq.iex) then
1043 if (ndowna(k,il).eq.i) then
1045 elseif (ndowna(k,il).eq.iex) then
1051 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1053 i2rep(i-1)=i2rep(iex-1)
1056 write(iout,*) 'exchange',i,iex
1057 write (iout,'(a8,100i4)') "@ ifirst",
1058 & (ifirst(k),k=1,remd_m(1))
1060 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1061 & (nupa(k,il),k=1,nupa(0,il))
1062 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1063 & (ndowna(k,il),k=1,ndowna(0,il))
1068 remd_ene(0,iex)=ene_iex_iex
1069 remd_ene(0,i)=ene_i_i
1076 cd write (iout,*) "exchange completed"
1080 cd write(iout,*) "########",ii
1082 i_temp=iran_num(1,nrep)
1083 i_mult=iran_num(1,remd_m(i_temp))
1084 i_iset=iran_num(1,nset)
1085 i_mset=iran_num(1,mset(i_iset))
1086 i=i_index(i_temp,i_mult,i_iset,i_mset)
1088 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1091 cd write(iout,*) "i_dir=",i_dir
1093 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1096 i_mult1=iran_num(1,remd_m(i_temp1))
1098 i_mset1=iran_num(1,mset(i_iset1))
1099 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1100 c 9/1/17 AL: Correction; otherwise the restraint energies are mis-assigned
1101 c on failed replica exchange attempt
1102 econstr_temp_i=remd_ene(i_econstr,i)
1103 econstr_temp_iex=remd_ene(i_econstr,iex)
1104 c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
1106 remd_ene(i_econstr,i)=remd_ene(n_ene+5,i)
1107 remd_ene(i_econstr,iex)=remd_ene(n_ene+6,iex)
1109 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1112 i_mult1=iran_num(1,remd_m(i_temp1))
1114 i_mset1=iran_num(1,mset(i_iset1))
1115 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1116 econstr_temp_i=remd_ene(i_econstr,i)
1117 econstr_temp_iex=remd_ene(i_econstr,iex)
1118 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1119 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1121 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1124 i_mult1=iran_num(1,remd_m(i_temp1))
1126 i_mset1=iran_num(1,mset(i_iset1))
1127 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1128 econstr_temp_i=remd_ene(i_econstr,i)
1129 econstr_temp_iex=remd_ene(i_econstr,iex)
1131 remd_ene(i_econstr,i)=remd_ene(n_ene+7,i)
1132 remd_ene(i_econstr,iex)=remd_ene(n_ene+8,iex)
1134 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1135 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1141 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1144 c Swap temperatures between conformations i and iex with recalculating the free energies
1145 c following temperature changes.
1147 write (iout,*) "i_dir",i_dir," i",i," iex",iex,
1148 & " t_bath",remd_t_bath(i),remd_t_bath(iex)
1150 ene_iex_iex=remd_ene(0,iex)
1151 ene_i_i=remd_ene(0,i)
1152 co write (iout,*) "rescaling weights with temperature",
1154 call rescale_weights(remd_t_bath(i))
1156 call sum_energy(remd_ene(0,iex),.false.)
1157 ene_iex_i=remd_ene(0,iex)
1158 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
1159 c call sum_energy(remd_ene(0,i),.false.)
1160 cd write (iout,*) "ene_i_i",remd_ene(0,i)
1161 c write (iout,*) "rescaling weights with temperature",
1162 c & remd_t_bath(iex)
1163 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1164 c write (iout,*) "ERROR: inconsistent energies:",i,
1165 c & ene_i_i,remd_ene(0,i)
1167 call rescale_weights(remd_t_bath(iex))
1168 call sum_energy(remd_ene(0,i),.false.)
1169 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1170 ene_i_iex=remd_ene(0,i)
1171 c call sum_energy(remd_ene(0,iex),.false.)
1172 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1173 c write (iout,*) "ERROR: inconsistent energies:",iex,
1174 c & ene_iex_iex,remd_ene(0,iex)
1176 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1177 c write (iout,*) "i",i," iex",iex
1178 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1179 cd & " ene_i_iex",ene_i_iex,
1180 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1181 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1182 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1184 cd write(iout,*) 'delta',delta
1185 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1186 c & (remd_ene(i)-remd_ene(iex))/Rb/
1187 c & (remd_t_bath(i)*remd_t_bath(iex))
1188 if (delta .gt. 50.0d0) then
1193 if (i_dir.eq.1.or.i_dir.eq.3)
1194 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1195 if (i_dir.eq.2.or.i_dir.eq.3)
1196 & iremd_tot_usa(int(i2set(i-1)))=
1197 & iremd_tot_usa(int(i2set(i-1)))+1
1198 xxx=ran_number(0.0d0,1.0d0)
1199 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1200 if (delta .gt. xxx) then
1202 remd_t_bath(i)=remd_t_bath(iex)
1203 remd_t_bath(iex)=tmp
1206 iremd_iset(i)=iremd_iset(iex)
1207 iremd_iset(iex)=itmp
1209 remd_ene(0,i)=ene_i_iex
1210 remd_ene(0,iex)=ene_iex_i
1212 if (i_dir.eq.1.or.i_dir.eq.3)
1213 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1216 i2rep(i-1)=i2rep(iex-1)
1219 if (i_dir.eq.2.or.i_dir.eq.3)
1220 & iremd_acc_usa(int(i2set(i-1)))=
1221 & iremd_acc_usa(int(i2set(i-1)))+1
1224 i2set(i-1)=i2set(iex-1)
1227 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1228 i_index(i_temp,i_mult,i_iset,i_mset)=
1229 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1230 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1233 remd_ene(0,iex)=ene_iex_iex
1234 remd_ene(0,i)=ene_i_i
1235 remd_ene(i_econstr,iex)=econstr_temp_iex
1236 remd_ene(i_econstr,i)=econstr_temp_i
1240 cd do il1=1,mset(il)
1243 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1256 c-------------------------------------
1257 write (iout,*) "NREP",nrep
1259 if(iremd_tot(i).ne.0)
1260 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1261 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1266 if(iremd_tot_usa(i).ne.0)
1267 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1268 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1274 cd write (iout,'(a6,100i4)') "ifirst",
1275 cd & (ifirst(i),i=1,remd_m(1))
1277 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1278 cd & (nupa(i,il),i=1,nupa(0,il))
1279 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1280 cd & (ndowna(i,il),i=1,ndowna(0,il))
1285 cd write (iout,*) "Before scatter"
1287 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1288 & t_bath,1,mpi_double_precision,king,
1290 cd write (iout,*) "After scatter"
1293 & call mpi_scatter(iremd_iset,1,mpi_integer,
1294 & iset,1,mpi_integer,king,
1298 if (me.eq.king .or. .not. out1file) then
1299 write(iout,*) 'REMD scatter time=',time07-time06
1303 call mpi_scatter(elowi,1,mpi_double_precision,
1304 & elow,1,mpi_double_precision,king,
1306 call mpi_scatter(ehighi,1,mpi_double_precision,
1307 & ehigh,1,mpi_double_precision,king,
1310 call rescale_weights(t_bath)
1311 co write (iout,*) "Processor",me,
1312 co & " rescaling weights with temperature",t_bath
1314 stdfp=dsqrt(2*Rb*t_bath/d_time)
1316 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1318 c Compute the standard deviations of stochastic forces for Langevin dynamics
1319 c if the friction coefficients do not depend on surface area
1320 if (lang.gt.0 .and. .not.surfarea) then
1322 stdforcp(i)=stdfp*dsqrt(gamp)
1325 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1326 & *dsqrt(gamsc(iabs(itype(i))))
1329 cde write(iout,*) 'REMD after',me,t_bath
1331 if (me.eq.king .or. .not. out1file) then
1332 write(iout,*) 'REMD exchange time=',time08-time00
1338 if (restart1file) then
1339 if (me.eq.king .or. .not. out1file)
1340 & write(iout,*) 'writing restart at the end of run'
1341 call write1rst(i_index)
1344 if (traj1file) call write1traj
1346 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1347 cdeb & icache_all,1,mpi_integer,king,
1348 cdeb & CG_COMM,ierr)
1349 cdeb write(iout,'(a40,8000i8)')
1350 cdeb & ' ntwx_cache after traj1file at the end',
1351 cdeb & (icache_all(i),i=1,nodes)
1356 t_MD=MPI_Wtime()-tt0
1360 if (me.eq.king .or. .not. out1file) then
1361 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1363 & 'MD calculations setup:',t_MDsetup,
1364 & 'Energy & gradient evaluation:',t_enegrad,
1365 & 'Stochastic MD setup:',t_langsetup,
1366 & 'Stochastic MD step setup:',t_sdsetup,
1368 write (iout,'(/28(1h=),a25,27(1h=))')
1369 & ' End of MD calculation '
1374 c-----------------------------------------------------------------------
1375 subroutine write1rst(i_index)
1377 include 'DIMENSIONS'
1379 include 'COMMON.CONTROL'
1382 include 'COMMON.LAGRANGE.5diag'
1384 include 'COMMON.LAGRANGE'
1386 include 'COMMON.QRESTR'
1387 include 'COMMON.IOUNITS'
1388 include 'COMMON.REMD'
1389 include 'COMMON.SETUP'
1390 include 'COMMON.CHAIN'
1391 include 'COMMON.SBRIDGE'
1392 include 'COMMON.INTERACT'
1394 real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1),
1395 & d_restart2(3,2*maxres*maxprocs)
1399 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1400 common /przechowalnia/ d_restart1,d_restart2
1401 integer i,j,il1,il,ixdrf
1407 t5_restart1(4)=t_bath
1408 t5_restart1(5)=Uconst
1410 call mpi_gather(t5_restart1,5,mpi_real,
1411 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1419 call mpi_gather(r_d,3*2*nres,mpi_real,
1420 & d_restart1,3*2*nres,mpi_real,king,
1429 call mpi_gather(r_d,3*2*nres,mpi_real,
1430 & d_restart2,3*2*nres,mpi_real,king,
1435 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1437 call xdrfint_(ixdrf, i2rep(i), iret)
1440 call xdrfint_(ixdrf, ifirst(i), iret)
1444 call xdrfint_(ixdrf, nupa(i,il), iret)
1448 call xdrfint_(ixdrf, ndowna(i,il), iret)
1454 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1461 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1468 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1474 call xdrfint_(ixdrf, nset, iret)
1476 call xdrfint_(ixdrf,mset(i), iret)
1479 call xdrfint_(ixdrf,i2set(i), iret)
1485 itmp=i_index(i,j,il,il1)
1486 call xdrfint_(ixdrf,itmp, iret)
1493 call xdrfclose_(ixdrf, iret)
1495 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1497 call xdrfint(ixdrf, i2rep(i), iret)
1500 call xdrfint(ixdrf, ifirst(i), iret)
1504 call xdrfint(ixdrf, nupa(i,il), iret)
1508 call xdrfint(ixdrf, ndowna(i,il), iret)
1514 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1521 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1528 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1535 call xdrfint(ixdrf, nset, iret)
1537 call xdrfint(ixdrf,mset(i), iret)
1540 call xdrfint(ixdrf,i2set(i), iret)
1546 itmp=i_index(i,j,il,il1)
1547 call xdrfint(ixdrf,itmp, iret)
1554 call xdrfclose(ixdrf, iret)
1561 subroutine write1traj
1563 include 'DIMENSIONS'
1566 include 'COMMON.QRESTR'
1567 include 'COMMON.IOUNITS'
1568 include 'COMMON.REMD'
1569 include 'COMMON.SETUP'
1570 include 'COMMON.CHAIN'
1571 include 'COMMON.SBRIDGE'
1572 include 'COMMON.INTERACT'
1576 real xcoord(3,maxres2+2),prec
1577 real r_qfrag(50),r_qpair(100)
1578 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1579 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1580 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1581 & p_uscdiff(100*maxprocs)
1582 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1583 common /przechowalnia/ p_c
1584 integer ii,i,il,j,ixdrf
1587 call mpi_bcast(ii_write,1,mpi_integer,
1588 & king,CG_COMM,ierr)
1591 c print *,'traj1file',me,ii_write,ntwx_cache
1595 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1597 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1600 t5_restart1(1)=totT_cache(ii)
1601 t5_restart1(2)=EK_cache(ii)
1602 t5_restart1(3)=potE_cache(ii)
1603 t5_restart1(4)=t_bath_cache(ii)
1604 t5_restart1(5)=Uconst_cache(ii)
1605 call mpi_gather(t5_restart1,5,mpi_real,
1606 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1608 call mpi_gather(iset_cache(ii),1,mpi_integer,
1609 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1612 r_qfrag(i)=qfrag_cache(i,ii)
1615 r_qpair(i)=qpair_cache(i,ii)
1618 r_utheta(i)=utheta_cache(i,ii)
1619 r_ugamma(i)=ugamma_cache(i,ii)
1620 r_uscdiff(i)=uscdiff_cache(i,ii)
1623 call mpi_gather(r_qfrag,nfrag,mpi_real,
1624 & p_qfrag,nfrag,mpi_real,king,
1626 call mpi_gather(r_qpair,npair,mpi_real,
1627 & p_qpair,npair,mpi_real,king,
1629 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1630 & p_utheta,nfrag_back,mpi_real,king,
1632 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1633 & p_ugamma,nfrag_back,mpi_real,king,
1635 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1636 & p_uscdiff,nfrag_back,mpi_real,king,
1640 write (iout,*) "p_qfrag"
1642 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1644 write (iout,*) "p_qpair"
1646 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1652 r_c(j,i)=c_cache(j,i,ii)
1656 call mpi_gather(r_c,3*2*nres,mpi_real,
1657 & p_c,3*2*nres,mpi_real,king,
1663 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1664 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1665 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1666 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1667 call xdrfint_(ixdrf, nss, iret)
1670 call xdrfint(ixdrf, idssb(j)+nres, iret)
1671 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1673 call xdrfint_(ixdrf, ihpb(j), iret)
1674 call xdrfint_(ixdrf, jhpb(j), iret)
1677 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1678 call xdrfint_(ixdrf, iset_restart1(il), iret)
1680 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1683 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1686 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1687 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1688 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1693 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1698 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1702 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1706 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1707 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1708 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1709 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1710 call xdrfint(ixdrf, nss, iret)
1713 call xdrfint(ixdrf, idssb(j)+nres, iret)
1714 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1716 call xdrfint(ixdrf, ihpb(j), iret)
1717 call xdrfint(ixdrf, jhpb(j), iret)
1720 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1721 call xdrfint(ixdrf, iset_restart1(il), iret)
1723 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1726 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1729 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1730 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1731 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1736 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1741 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1745 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1751 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1753 if(me.eq.king) call xdrfclose(ixdrf, iret)
1755 do i=1,ntwx_cache-ii_write
1757 totT_cache(i)=totT_cache(ii_write+i)
1758 EK_cache(i)=EK_cache(ii_write+i)
1759 potE_cache(i)=potE_cache(ii_write+i)
1760 t_bath_cache(i)=t_bath_cache(ii_write+i)
1761 Uconst_cache(i)=Uconst_cache(ii_write+i)
1762 iset_cache(i)=iset_cache(ii_write+i)
1765 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1768 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1771 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1772 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1773 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1778 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1782 ntwx_cache=ntwx_cache-ii_write
1787 subroutine read1restart(i_index)
1789 include 'DIMENSIONS'
1791 include 'COMMON.CONTROL'
1794 include 'COMMON.LAGRANGE.5diag'
1796 include 'COMMON.LAGRANGE'
1798 include 'COMMON.QRESTR'
1799 include 'COMMON.IOUNITS'
1800 include 'COMMON.REMD'
1801 include 'COMMON.SETUP'
1802 include 'COMMON.CHAIN'
1803 include 'COMMON.SBRIDGE'
1804 include 'COMMON.INTERACT'
1805 real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1),
1808 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1809 common /przechowalnia/ d_restart1
1810 integer i,j,il,il1,ixdrf,iret,itmp
1812 c write (*,*) "Processor",me," called read1restart"
1815 open(irest2,file=mremd_rst_name,status='unknown')
1816 read(irest2,*,err=334) i
1817 write(iout,*) "Reading old rst in ASCI format"
1819 call read1restart_old
1823 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1826 call xdrfint_(ixdrf, i2rep(i), iret)
1829 call xdrfint_(ixdrf, ifirst(i), iret)
1832 call xdrfint_(ixdrf, nupa(0,il), iret)
1834 call xdrfint_(ixdrf, nupa(i,il), iret)
1837 call xdrfint_(ixdrf, ndowna(0,il), iret)
1839 call xdrfint_(ixdrf, ndowna(i,il), iret)
1844 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1848 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1851 call xdrfint(ixdrf, i2rep(i), iret)
1854 call xdrfint(ixdrf, ifirst(i), iret)
1857 call xdrfint(ixdrf, nupa(0,il), iret)
1859 call xdrfint(ixdrf, nupa(i,il), iret)
1862 call xdrfint(ixdrf, ndowna(0,il), iret)
1864 call xdrfint(ixdrf, ndowna(i,il), iret)
1869 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1874 call mpi_scatter(t_restart1,5,mpi_real,
1875 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1879 t_bath=t5_restart1(4)
1884 c read(irest2,'(3e15.5)')
1885 c & (d_restart1(j,i+2*nres*il),j=1,3)
1888 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1890 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1896 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1897 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1907 c read(irest2,'(3e15.5)')
1908 c & (d_restart1(j,i+2*nres*il),j=1,3)
1911 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1913 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1919 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1920 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1931 call xdrfint_(ixdrf, nset, iret)
1933 call xdrfint_(ixdrf,mset(i), iret)
1936 call xdrfint_(ixdrf,i2set(i), iret)
1942 call xdrfint_(ixdrf,itmp, iret)
1943 i_index(i,j,il,il1)=itmp
1951 call xdrfint(ixdrf, nset, iret)
1953 call xdrfint(ixdrf,mset(i), iret)
1956 call xdrfint(ixdrf,i2set(i), iret)
1962 call xdrfint(ixdrf,itmp, iret)
1963 i_index(i,j,il,il1)=itmp
1970 Corrected AL 8/19/2014: each processor needs whole iset array not only its
1972 c call mpi_scatter(i2set,1,mpi_integer,
1973 c & iset,1,mpi_integer,king,
1975 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1982 if(me.eq.king) close(irest2)
1986 subroutine read1restart_old
1988 include 'DIMENSIONS'
1992 include 'COMMON.LAGRANGE.5diag'
1994 include 'COMMON.LAGRANGE'
1996 include 'COMMON.IOUNITS'
1997 include 'COMMON.REMD'
1998 include 'COMMON.SETUP'
1999 include 'COMMON.CHAIN'
2000 include 'COMMON.SBRIDGE'
2001 include 'COMMON.INTERACT'
2002 real d_restart1(3,2*maxres*maxprocs),r_d(3,0:2*maxres-1),
2004 common /przechowalnia/ d_restart1
2008 open(irest2,file=mremd_rst_name,status='unknown')
2009 read (irest2,*) (i2rep(i),i=0,nodes-1)
2010 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2012 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2013 read (irest2,*) ndowna(0,il),
2014 & (ndowna(i,il),i=1,ndowna(0,il))
2017 read(irest2,*) (t_restart1(j,il),j=1,4)
2020 call mpi_scatter(t_restart1,5,mpi_real,
2021 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2025 t_bath=t5_restart1(4)
2030 read(irest2,'(3e15.5)')
2031 & (d_restart1(j,i+2*nres*il),j=1,3)
2035 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2036 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2046 read(irest2,'(3e15.5)')
2047 & (d_restart1(j,i+2*nres*il),j=1,3)
2051 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2052 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2058 if(me.eq.king) close(irest2)
2061 c------------------------------------------