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
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 write (*,*) "Processor",me," rest",rest,"
149 & 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
427 do while(.not.end_of_run)
429 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
430 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
432 if (lang.gt.0 .and. surfarea .and.
433 & mod(itime,reset_fricmat).eq.0) then
434 if (lang.eq.2 .or. lang.eq.3) then
438 call sd_verlet_p_setup
440 call sd_verlet_ciccotti_setup
444 pfric0_mat(i,j,0)=pfric_mat(i,j)
445 afric0_mat(i,j,0)=afric_mat(i,j)
446 vfric0_mat(i,j,0)=vfric_mat(i,j)
447 prand0_mat(i,j,0)=prand_mat(i,j)
448 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
449 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
454 flag_stoch(i)=.false.
457 else if (lang.eq.1 .or. lang.eq.4) then
460 write (iout,'(a,i10)')
461 & "Friction matrix reset based on surface area, itime",itime
463 if (reset_vel .and. tbf .and. lang.eq.0
464 & .and. mod(itime,count_reset_vel).eq.0) then
466 if (me.eq.king .or. .not. out1file)
467 & write(iout,'(a,f20.2)')
468 & "Velocities reset to random values, time",totT
471 d_t_old(j,i)=d_t(j,i)
475 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
479 d_t(j,0)=d_t(j,0)-vcm(j)
482 kinetic_T=2.0d0/(dimen3*Rb)*EK
483 scalfac=dsqrt(T_bath/kinetic_T)
484 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
487 d_t_old(j,i)=scalfac*d_t(j,i)
493 c Time-reversible RESPA algorithm
494 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
495 call RESPA_step(itime)
497 c Variable time step algorithm.
498 call velverlet_step(itime)
502 call brown_step(itime)
504 print *,"Brown dynamics not here!"
506 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
512 if (mod(itime,ntwe).eq.0) call statout(itime)
514 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
515 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
517 call hairpin(.true.,nharp,iharp)
518 call secondary2(.true.)
519 call pdbout(potE,tytul,ipdb)
524 if (mod(itime,ntwx).eq.0.and.traj1file) then
525 if(ntwx_cache.lt.max_cache_traj_use) then
526 ntwx_cache=ntwx_cache+1
528 if (max_cache_traj_use.ne.1)
529 & print *,itime,"processor ",me," over cache ",ntwx_cache
532 totT_cache(i)=totT_cache(i+1)
533 EK_cache(i)=EK_cache(i+1)
534 potE_cache(i)=potE_cache(i+1)
535 t_bath_cache(i)=t_bath_cache(i+1)
536 Uconst_cache(i)=Uconst_cache(i+1)
537 iset_cache(i)=iset_cache(i+1)
540 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
543 qpair_cache(ii,i)=qpair_cache(ii,i+1)
546 utheta_cache(ii,i)=utheta_cache(ii,i+1)
547 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
548 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
554 c_cache(j,ii,i)=c_cache(j,ii,i+1)
560 totT_cache(ntwx_cache)=totT
561 EK_cache(ntwx_cache)=EK
562 potE_cache(ntwx_cache)=potE
563 t_bath_cache(ntwx_cache)=t_bath
564 Uconst_cache(ntwx_cache)=Uconst
565 iset_cache(ntwx_cache)=iset
568 qfrag_cache(i,ntwx_cache)=qfrag(i)
571 qpair_cache(i,ntwx_cache)=qpair(i)
574 utheta_cache(i,ntwx_cache)=utheta(i)
575 ugamma_cache(i,ntwx_cache)=ugamma(i)
576 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
578 C print *,'przed returnbox'
580 C call enerprint(remd_ene(0,i))
583 c_cache(j,i,ntwx_cache)=c(j,i)
588 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
589 & .and..not.restart1file) then
592 open(irest1,file=mremd_rst_name,status='unknown')
593 write (irest1,*) "i2rep"
594 write (irest1,*) (i2rep(i),i=0,nodes-1)
595 write (irest1,*) "ifirst"
596 write (irest1,*) (ifirst(i),i=1,remd_m(1))
598 write (irest1,*) "nupa",il
599 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
600 write (irest1,*) "ndowna",il
601 write (irest1,*) ndowna(0,il),
602 & (ndowna(i,il),i=1,ndowna(0,il))
605 write (irest1,*) "nset"
606 write (irest1,*) nset
607 write (irest1,*) "mset"
608 write (irest1,*) (mset(i),i=1,nset)
609 write (irest1,*) "i2set"
610 write (irest1,*) (i2set(i),i=0,nodes-1)
611 write (irest1,*) "i_index"
615 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
623 open(irest2,file=rest2name,status='unknown')
624 write(irest2,*) totT,EK,potE,totE,t_bath
626 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
629 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
632 write (irest2,*) iset
639 c forced synchronization
640 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
641 & .and. .not. mremdsync) then
643 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
645 call mpi_recv(itime_master, 1, MPI_INTEGER,
646 & 0,101,CG_COMM, status, ierr)
647 call mpi_barrier(CG_COMM, ierr)
648 cdeb if (out1file.or.traj1file) then
649 cdeb call mpi_gather(itime,1,mpi_integer,
650 cdeb & icache_all,1,mpi_integer,king,
653 & call mpi_gather(ntwx_cache,1,mpi_integer,
654 & icache_all,1,mpi_integer,king,
657 & write(iout,*) 'REMD synchro at',itime_master,itime
658 if (itime_master.ge.n_timestep .or. ovrtim())
660 ctime call flush(iout)
665 if ((mod(itime,nstex).eq.0.and.me.eq.king
666 & .or.end_of_run.and.me.eq.king )
667 & .and. .not. mremdsync ) then
670 call mpi_isend(itime,1,MPI_INTEGER,i,101,
671 & CG_COMM, ireqi(i), ierr)
672 cd write(iout,*) 'REMD synchro with',i
675 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
676 call mpi_barrier(CG_COMM, ierr)
678 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
679 if (out1file.or.traj1file) then
680 cdeb call mpi_gather(itime,1,mpi_integer,
681 cdeb & itime_all,1,mpi_integer,king,
683 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
684 cdeb & (itime_all(i),i=1,nodes)
686 cdeb imin_itime=itime_all(1)
688 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
690 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
691 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
692 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
693 call mpi_gather(ntwx_cache,1,mpi_integer,
694 & icache_all,1,mpi_integer,king,
696 c write(iout,'(a19,8000i8)') ' ntwx_cache',
697 c & (icache_all(i),i=1,nodes)
698 ii_write=icache_all(1)
700 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
702 c write(iout,*) "MIN ii_write=",ii_write
705 ctime call flush(iout)
707 if(mremdsync .and. mod(itime,nstex).eq.0) then
709 if (me.eq.king .or. .not. out1file)
710 & write(iout,*) 'REMD synchro at',itime
713 call mpi_gather(ntwx_cache,1,mpi_integer,
714 & icache_all,1,mpi_integer,king,
717 write(iout,'(a19,8000i8)') ' ntwx_cache',
718 & (icache_all(i),i=1,nodes)
719 ii_write=icache_all(1)
721 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
723 write(iout,*) "MIN ii_write=",ii_write
729 c Update the time safety limiy
730 if (time001-time00.gt.safety) then
731 safety=time001-time00+600
732 write (iout,*) "****** SAFETY increased to",safety," s"
734 if (ovrtim()) end_of_run=.true.
736 if(synflag.and..not.end_of_run) then
740 write(iout,*) 'REMD before',me,t_bath
742 c call mpi_gather(t_bath,1,mpi_double_precision,
743 c & remd_t_bath,1,mpi_double_precision,king,
745 potEcomp(n_ene+1)=t_bath
747 potEcomp(n_ene+2)=iset
748 if (iset.lt.nset) then
753 call Econstr_back_qlike
757 potEcomp(n_ene+3)=Uconst+Uconst_back
765 call Econstr_back_qlike
769 potEcomp(n_ene+4)=Uconst+Uconst_back
773 c 9/11/17 Adaptive US
776 write (iout,*) "me ",me," itt",itt
778 if (itt.lt.nrep) then
782 potEcomp(n_ene+5)=Uconst
784 write (iout,*) "t_bath",t_bath_temp,t_bath,
787 if (iset.lt.nset) then
791 potEcomp(n_ene+7)=Uconst
793 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
803 potEcomp(n_ene+6)=Uconst
805 write (iout,*) "t_bath",t_bath_temp,t_bath,
812 potEcomp(n_ene+8)=Uconst
814 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
822 call mpi_gather(potEcomp(0),n_ene+9,mpi_double_precision,
823 & remd_ene(0,1),n_ene+9,mpi_double_precision,king,
826 call mpi_gather(elow,1,mpi_double_precision,
827 & elowi,1,mpi_double_precision,king,
829 call mpi_gather(ehigh,1,mpi_double_precision,
830 & ehighi,1,mpi_double_precision,king,
835 if (me.eq.king .or. .not. out1file) then
836 write(iout,*) 'REMD gather times=',time03-time01
840 if (restart1file) call write1rst(i_index)
843 if (me.eq.king .or. .not. out1file) then
844 write(iout,*) 'REMD writing rst time=',time04-time03
847 if (traj1file) call write1traj
849 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
850 cdeb & icache_all,1,mpi_integer,king,
852 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
853 cdeb & (icache_all(i),i=1,nodes)
858 if (me.eq.king .or. .not. out1file) then
859 write(iout,*) 'REMD writing traj time=',time05-time04
866 remd_t_bath(i)=remd_ene(n_ene+1,i)
867 iremd_iset(i)=remd_ene(n_ene+2,i)
871 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
873 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
877 write(iout,*) 'REMD exchange temp,ene'
879 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
880 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
884 c-------------------------------------
887 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
889 write (iout,*) "remd_m(1)",remd_m(1)
892 i=ifirst(iran_num(1,remd_m(1)))
897 & write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
898 if(i.gt.0.and.nupa(0,i).gt.0) then
900 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
902 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
904 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
906 c do while (iex.eq.i)
907 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
908 iex=nupa(iran_num(1,int(nupa(0,i))),i)
910 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
913 write (iout,*) "i",i,"iex",iex," temperatures",
914 & remd_t_bath(i),remd_t_bath(iex)
917 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
919 c Swap temperatures between conformations i and iex with recalculating the free energies
920 c following temperature changes.
921 ene_iex_iex=remd_ene(0,iex)
922 ene_i_i=remd_ene(0,i)
923 c write (iout,*) "i",i," ene_i_i",ene_i_i,
924 c & " iex",iex," ene_iex_iex",ene_iex_iex
925 c write (iout,*) "rescaling weights with temperature",
928 call rescale_weights(remd_t_bath(i))
930 c write (iout,*) "0,iex",remd_t_bath(i)
931 c call enerprint(remd_ene(0,iex))
933 call sum_energy(remd_ene(0,iex),.false.)
934 ene_iex_i=remd_ene(0,iex)
935 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
937 c write (iout,*) "0,i",remd_t_bath(i)
938 c call enerprint(remd_ene(0,i))
940 call sum_energy(remd_ene(0,i),.false.)
941 c write (iout,*) "ene_i_i",remd_ene(0,i)
943 c write (iout,*) "rescaling weights with temperature",
945 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
946 write (iout,*) "ERROR: inconsistent energies:",i,
947 & ene_i_i,remd_ene(0,i)
949 call rescale_weights(remd_t_bath(iex))
951 c write (iout,*) "0,i",remd_t_bath(iex)
952 c call enerprint(remd_ene(0,i))
954 call sum_energy(remd_ene(0,i),.false.)
955 c write (iout,*) "ene_i_iex",remd_ene(0,i)
957 ene_i_iex=remd_ene(0,i)
959 c write (iout,*) "0,iex",remd_t_bath(iex)
960 c call enerprint(remd_ene(0,iex))
962 call sum_energy(remd_ene(0,iex),.false.)
963 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
964 write (iout,*) "ERROR: inconsistent energies:",iex,
965 & ene_iex_iex,remd_ene(0,iex)
968 write (iout,*) "ene_iex_iex",remd_ene(0,iex)
969 write (iout,*) "i",i," iex",iex
970 write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
971 & " ene_i_iex",ene_i_iex,
972 & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
975 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
976 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
978 c write(iout,*) 'delta',delta
979 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
980 c & (remd_ene(i)-remd_ene(iex))/Rb/
981 c & (remd_t_bath(i)*remd_t_bath(iex))
983 if (delta .gt. 50.0d0) then
989 else if (delta.lt.-50.0d0) then
998 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
999 xxx=ran_number(0.0d0,1.0d0)
1000 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1002 if (delta .gt. xxx) then
1004 remd_t_bath(i)=remd_t_bath(iex)
1005 remd_t_bath(iex)=tmp
1006 remd_ene(0,i)=ene_i_iex
1007 remd_ene(0,iex)=ene_iex_i
1013 ehighi(i)=ehighi(iex)
1020 nupa(k,i)=nupa(k,iex)
1023 ndowna(k,i)=ndowna(k,iex)
1027 if (ifirst(il).eq.i) ifirst(il)=iex
1029 if (nupa(k,il).eq.i) then
1031 elseif (nupa(k,il).eq.iex) then
1036 if (ndowna(k,il).eq.i) then
1038 elseif (ndowna(k,il).eq.iex) then
1044 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1046 i2rep(i-1)=i2rep(iex-1)
1049 write(iout,*) 'exchange',i,iex
1050 write (iout,'(a8,100i4)') "@ ifirst",
1051 & (ifirst(k),k=1,remd_m(1))
1053 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1054 & (nupa(k,il),k=1,nupa(0,il))
1055 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1056 & (ndowna(k,il),k=1,ndowna(0,il))
1061 remd_ene(0,iex)=ene_iex_iex
1062 remd_ene(0,i)=ene_i_i
1068 cd write (iout,*) "exchange completed"
1072 cd write(iout,*) "########",ii
1074 i_temp=iran_num(1,nrep)
1075 i_mult=iran_num(1,remd_m(i_temp))
1076 i_iset=iran_num(1,nset)
1077 i_mset=iran_num(1,mset(i_iset))
1078 i=i_index(i_temp,i_mult,i_iset,i_mset)
1080 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1083 cd write(iout,*) "i_dir=",i_dir
1085 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1088 i_mult1=iran_num(1,remd_m(i_temp1))
1090 i_mset1=iran_num(1,mset(i_iset1))
1091 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1092 c 9/1/17 AL: Correction; otherwise the restraint energies are mis-assigned
1093 c on failed replica exchange attempt
1094 econstr_temp_i=remd_ene(i_econstr,i)
1095 econstr_temp_iex=remd_ene(i_econstr,iex)
1096 c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
1098 remd_ene(i_econstr,i)=remd_ene(n_ene+5,i)
1099 remd_ene(i_econstr,iex)=remd_ene(n_ene+6,iex)
1101 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1104 i_mult1=iran_num(1,remd_m(i_temp1))
1106 i_mset1=iran_num(1,mset(i_iset1))
1107 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1108 econstr_temp_i=remd_ene(i_econstr,i)
1109 econstr_temp_iex=remd_ene(i_econstr,iex)
1110 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1111 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1113 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1116 i_mult1=iran_num(1,remd_m(i_temp1))
1118 i_mset1=iran_num(1,mset(i_iset1))
1119 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1120 econstr_temp_i=remd_ene(i_econstr,i)
1121 econstr_temp_iex=remd_ene(i_econstr,iex)
1123 remd_ene(i_econstr,i)=remd_ene(n_ene+7,i)
1124 remd_ene(i_econstr,iex)=remd_ene(n_ene+8,iex)
1126 remd_ene(i_econstr,i)=remd_ene(n_ene+3,i)
1127 remd_ene(i_econstr,iex)=remd_ene(n_ene+4,iex)
1133 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1136 c Swap temperatures between conformations i and iex with recalculating the free energies
1137 c following temperature changes.
1139 write (iout,*) "i_dir",i_dir," i",i," iex",iex,
1140 & " t_bath",remd_t_bath(i),remd_t_bath(iex)
1142 ene_iex_iex=remd_ene(0,iex)
1143 ene_i_i=remd_ene(0,i)
1144 co write (iout,*) "rescaling weights with temperature",
1146 call rescale_weights(remd_t_bath(i))
1148 call sum_energy(remd_ene(0,iex),.false.)
1149 ene_iex_i=remd_ene(0,iex)
1150 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
1151 c call sum_energy(remd_ene(0,i),.false.)
1152 cd write (iout,*) "ene_i_i",remd_ene(0,i)
1153 c write (iout,*) "rescaling weights with temperature",
1154 c & remd_t_bath(iex)
1155 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1156 c write (iout,*) "ERROR: inconsistent energies:",i,
1157 c & ene_i_i,remd_ene(0,i)
1159 call rescale_weights(remd_t_bath(iex))
1160 call sum_energy(remd_ene(0,i),.false.)
1161 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1162 ene_i_iex=remd_ene(0,i)
1163 c call sum_energy(remd_ene(0,iex),.false.)
1164 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1165 c write (iout,*) "ERROR: inconsistent energies:",iex,
1166 c & ene_iex_iex,remd_ene(0,iex)
1168 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1169 c write (iout,*) "i",i," iex",iex
1170 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1171 cd & " ene_i_iex",ene_i_iex,
1172 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1173 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1174 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1176 cd write(iout,*) 'delta',delta
1177 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1178 c & (remd_ene(i)-remd_ene(iex))/Rb/
1179 c & (remd_t_bath(i)*remd_t_bath(iex))
1180 if (delta .gt. 50.0d0) then
1185 if (i_dir.eq.1.or.i_dir.eq.3)
1186 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1187 if (i_dir.eq.2.or.i_dir.eq.3)
1188 & iremd_tot_usa(int(i2set(i-1)))=
1189 & iremd_tot_usa(int(i2set(i-1)))+1
1190 xxx=ran_number(0.0d0,1.0d0)
1191 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1192 if (delta .gt. xxx) then
1194 remd_t_bath(i)=remd_t_bath(iex)
1195 remd_t_bath(iex)=tmp
1198 iremd_iset(i)=iremd_iset(iex)
1199 iremd_iset(iex)=itmp
1201 remd_ene(0,i)=ene_i_iex
1202 remd_ene(0,iex)=ene_iex_i
1204 if (i_dir.eq.1.or.i_dir.eq.3)
1205 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1208 i2rep(i-1)=i2rep(iex-1)
1211 if (i_dir.eq.2.or.i_dir.eq.3)
1212 & iremd_acc_usa(int(i2set(i-1)))=
1213 & iremd_acc_usa(int(i2set(i-1)))+1
1216 i2set(i-1)=i2set(iex-1)
1219 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1220 i_index(i_temp,i_mult,i_iset,i_mset)=
1221 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1222 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1225 remd_ene(0,iex)=ene_iex_iex
1226 remd_ene(0,i)=ene_i_i
1227 remd_ene(i_econstr,iex)=econstr_temp_iex
1228 remd_ene(i_econstr,i)=econstr_temp_i
1232 cd do il1=1,mset(il)
1235 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1248 c-------------------------------------
1249 write (iout,*) "NREP",nrep
1251 if(iremd_tot(i).ne.0)
1252 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1253 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1258 if(iremd_tot_usa(i).ne.0)
1259 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1260 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1266 cd write (iout,'(a6,100i4)') "ifirst",
1267 cd & (ifirst(i),i=1,remd_m(1))
1269 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1270 cd & (nupa(i,il),i=1,nupa(0,il))
1271 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1272 cd & (ndowna(i,il),i=1,ndowna(0,il))
1277 cd write (iout,*) "Before scatter"
1279 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1280 & t_bath,1,mpi_double_precision,king,
1282 cd write (iout,*) "After scatter"
1285 & call mpi_scatter(iremd_iset,1,mpi_integer,
1286 & iset,1,mpi_integer,king,
1290 if (me.eq.king .or. .not. out1file) then
1291 write(iout,*) 'REMD scatter time=',time07-time06
1295 call mpi_scatter(elowi,1,mpi_double_precision,
1296 & elow,1,mpi_double_precision,king,
1298 call mpi_scatter(ehighi,1,mpi_double_precision,
1299 & ehigh,1,mpi_double_precision,king,
1302 call rescale_weights(t_bath)
1303 co write (iout,*) "Processor",me,
1304 co & " rescaling weights with temperature",t_bath
1306 stdfp=dsqrt(2*Rb*t_bath/d_time)
1308 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1310 c Compute the standard deviations of stochastic forces for Langevin dynamics
1311 c if the friction coefficients do not depend on surface area
1312 if (lang.gt.0 .and. .not.surfarea) then
1314 stdforcp(i)=stdfp*dsqrt(gamp)
1317 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1318 & *dsqrt(gamsc(iabs(itype(i))))
1321 cde write(iout,*) 'REMD after',me,t_bath
1323 if (me.eq.king .or. .not. out1file) then
1324 write(iout,*) 'REMD exchange time=',time08-time00
1330 if (restart1file) then
1331 if (me.eq.king .or. .not. out1file)
1332 & write(iout,*) 'writing restart at the end of run'
1333 call write1rst(i_index)
1336 if (traj1file) call write1traj
1338 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1339 cdeb & icache_all,1,mpi_integer,king,
1340 cdeb & CG_COMM,ierr)
1341 cdeb write(iout,'(a40,8000i8)')
1342 cdeb & ' ntwx_cache after traj1file at the end',
1343 cdeb & (icache_all(i),i=1,nodes)
1348 t_MD=MPI_Wtime()-tt0
1352 if (me.eq.king .or. .not. out1file) then
1353 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1355 & 'MD calculations setup:',t_MDsetup,
1356 & 'Energy & gradient evaluation:',t_enegrad,
1357 & 'Stochastic MD setup:',t_langsetup,
1358 & 'Stochastic MD step setup:',t_sdsetup,
1360 write (iout,'(/28(1h=),a25,27(1h=))')
1361 & ' End of MD calculation '
1366 c-----------------------------------------------------------------------
1367 subroutine write1rst(i_index)
1369 include 'DIMENSIONS'
1371 include 'COMMON.CONTROL'
1374 include 'COMMON.LAGRANGE.5diag'
1376 include 'COMMON.LAGRANGE'
1378 include 'COMMON.QRESTR'
1379 include 'COMMON.IOUNITS'
1380 include 'COMMON.REMD'
1381 include 'COMMON.SETUP'
1382 include 'COMMON.CHAIN'
1383 include 'COMMON.SBRIDGE'
1384 include 'COMMON.INTERACT'
1386 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1387 & d_restart2(3,2*maxres*maxprocs)
1391 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1392 common /przechowalnia/ d_restart1,d_restart2
1393 integer i,j,il1,il,ixdrf
1399 t5_restart1(4)=t_bath
1400 t5_restart1(5)=Uconst
1402 call mpi_gather(t5_restart1,5,mpi_real,
1403 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1411 call mpi_gather(r_d,3*2*nres,mpi_real,
1412 & d_restart1,3*2*nres,mpi_real,king,
1421 call mpi_gather(r_d,3*2*nres,mpi_real,
1422 & d_restart2,3*2*nres,mpi_real,king,
1427 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1429 call xdrfint_(ixdrf, i2rep(i), iret)
1432 call xdrfint_(ixdrf, ifirst(i), iret)
1436 call xdrfint_(ixdrf, nupa(i,il), iret)
1440 call xdrfint_(ixdrf, ndowna(i,il), iret)
1446 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1453 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1460 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1466 call xdrfint_(ixdrf, nset, iret)
1468 call xdrfint_(ixdrf,mset(i), iret)
1471 call xdrfint_(ixdrf,i2set(i), iret)
1477 itmp=i_index(i,j,il,il1)
1478 call xdrfint_(ixdrf,itmp, iret)
1485 call xdrfclose_(ixdrf, iret)
1487 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1489 call xdrfint(ixdrf, i2rep(i), iret)
1492 call xdrfint(ixdrf, ifirst(i), iret)
1496 call xdrfint(ixdrf, nupa(i,il), iret)
1500 call xdrfint(ixdrf, ndowna(i,il), iret)
1506 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1513 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1520 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1527 call xdrfint(ixdrf, nset, iret)
1529 call xdrfint(ixdrf,mset(i), iret)
1532 call xdrfint(ixdrf,i2set(i), iret)
1538 itmp=i_index(i,j,il,il1)
1539 call xdrfint(ixdrf,itmp, iret)
1546 call xdrfclose(ixdrf, iret)
1553 subroutine write1traj
1555 include 'DIMENSIONS'
1558 include 'COMMON.QRESTR'
1559 include 'COMMON.IOUNITS'
1560 include 'COMMON.REMD'
1561 include 'COMMON.SETUP'
1562 include 'COMMON.CHAIN'
1563 include 'COMMON.SBRIDGE'
1564 include 'COMMON.INTERACT'
1568 real xcoord(3,maxres2+2),prec
1569 real r_qfrag(50),r_qpair(100)
1570 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1571 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1572 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1573 & p_uscdiff(100*maxprocs)
1574 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1575 common /przechowalnia/ p_c
1576 integer ii,i,il,j,ixdrf
1579 call mpi_bcast(ii_write,1,mpi_integer,
1580 & king,CG_COMM,ierr)
1583 c print *,'traj1file',me,ii_write,ntwx_cache
1587 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1589 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1592 t5_restart1(1)=totT_cache(ii)
1593 t5_restart1(2)=EK_cache(ii)
1594 t5_restart1(3)=potE_cache(ii)
1595 t5_restart1(4)=t_bath_cache(ii)
1596 t5_restart1(5)=Uconst_cache(ii)
1597 call mpi_gather(t5_restart1,5,mpi_real,
1598 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1600 call mpi_gather(iset_cache(ii),1,mpi_integer,
1601 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1604 r_qfrag(i)=qfrag_cache(i,ii)
1607 r_qpair(i)=qpair_cache(i,ii)
1610 r_utheta(i)=utheta_cache(i,ii)
1611 r_ugamma(i)=ugamma_cache(i,ii)
1612 r_uscdiff(i)=uscdiff_cache(i,ii)
1615 call mpi_gather(r_qfrag,nfrag,mpi_real,
1616 & p_qfrag,nfrag,mpi_real,king,
1618 call mpi_gather(r_qpair,npair,mpi_real,
1619 & p_qpair,npair,mpi_real,king,
1621 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1622 & p_utheta,nfrag_back,mpi_real,king,
1624 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1625 & p_ugamma,nfrag_back,mpi_real,king,
1627 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1628 & p_uscdiff,nfrag_back,mpi_real,king,
1632 write (iout,*) "p_qfrag"
1634 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1636 write (iout,*) "p_qpair"
1638 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1644 r_c(j,i)=c_cache(j,i,ii)
1648 call mpi_gather(r_c,3*2*nres,mpi_real,
1649 & p_c,3*2*nres,mpi_real,king,
1655 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1656 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1657 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1658 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1659 call xdrfint_(ixdrf, nss, iret)
1662 call xdrfint(ixdrf, idssb(j)+nres, iret)
1663 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1665 call xdrfint_(ixdrf, ihpb(j), iret)
1666 call xdrfint_(ixdrf, jhpb(j), iret)
1669 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1670 call xdrfint_(ixdrf, iset_restart1(il), iret)
1672 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1675 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1678 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1679 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1680 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1685 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1690 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1694 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1698 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1699 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1700 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1701 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1702 call xdrfint(ixdrf, nss, iret)
1705 call xdrfint(ixdrf, idssb(j)+nres, iret)
1706 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1708 call xdrfint(ixdrf, ihpb(j), iret)
1709 call xdrfint(ixdrf, jhpb(j), iret)
1712 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1713 call xdrfint(ixdrf, iset_restart1(il), iret)
1715 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1718 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1721 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1722 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1723 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1728 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1733 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1737 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1743 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1745 if(me.eq.king) call xdrfclose(ixdrf, iret)
1747 do i=1,ntwx_cache-ii_write
1749 totT_cache(i)=totT_cache(ii_write+i)
1750 EK_cache(i)=EK_cache(ii_write+i)
1751 potE_cache(i)=potE_cache(ii_write+i)
1752 t_bath_cache(i)=t_bath_cache(ii_write+i)
1753 Uconst_cache(i)=Uconst_cache(ii_write+i)
1754 iset_cache(i)=iset_cache(ii_write+i)
1757 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1760 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1763 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1764 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1765 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1770 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1774 ntwx_cache=ntwx_cache-ii_write
1779 subroutine read1restart(i_index)
1781 include 'DIMENSIONS'
1783 include 'COMMON.CONTROL'
1786 include 'COMMON.LAGRANGE.5diag'
1788 include 'COMMON.LAGRANGE'
1790 include 'COMMON.QRESTR'
1791 include 'COMMON.IOUNITS'
1792 include 'COMMON.REMD'
1793 include 'COMMON.SETUP'
1794 include 'COMMON.CHAIN'
1795 include 'COMMON.SBRIDGE'
1796 include 'COMMON.INTERACT'
1797 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1800 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1801 common /przechowalnia/ d_restart1
1802 integer i,j,il,il1,ixdrf,iret,itmp
1804 write (*,*) "Processor",me," called read1restart"
1807 open(irest2,file=mremd_rst_name,status='unknown')
1808 read(irest2,*,err=334) i
1809 write(iout,*) "Reading old rst in ASCI format"
1811 call read1restart_old
1815 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1818 call xdrfint_(ixdrf, i2rep(i), iret)
1821 call xdrfint_(ixdrf, ifirst(i), iret)
1824 call xdrfint_(ixdrf, nupa(0,il), iret)
1826 call xdrfint_(ixdrf, nupa(i,il), iret)
1829 call xdrfint_(ixdrf, ndowna(0,il), iret)
1831 call xdrfint_(ixdrf, ndowna(i,il), iret)
1836 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1840 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1843 call xdrfint(ixdrf, i2rep(i), iret)
1846 call xdrfint(ixdrf, ifirst(i), iret)
1849 call xdrfint(ixdrf, nupa(0,il), iret)
1851 call xdrfint(ixdrf, nupa(i,il), iret)
1854 call xdrfint(ixdrf, ndowna(0,il), iret)
1856 call xdrfint(ixdrf, ndowna(i,il), iret)
1861 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1866 call mpi_scatter(t_restart1,5,mpi_real,
1867 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1871 t_bath=t5_restart1(4)
1876 c read(irest2,'(3e15.5)')
1877 c & (d_restart1(j,i+2*nres*il),j=1,3)
1880 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1882 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1888 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1889 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1899 c read(irest2,'(3e15.5)')
1900 c & (d_restart1(j,i+2*nres*il),j=1,3)
1903 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1905 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1911 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1912 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1923 call xdrfint_(ixdrf, nset, iret)
1925 call xdrfint_(ixdrf,mset(i), iret)
1928 call xdrfint_(ixdrf,i2set(i), iret)
1934 call xdrfint_(ixdrf,itmp, iret)
1935 i_index(i,j,il,il1)=itmp
1943 call xdrfint(ixdrf, nset, iret)
1945 call xdrfint(ixdrf,mset(i), iret)
1948 call xdrfint(ixdrf,i2set(i), iret)
1954 call xdrfint(ixdrf,itmp, iret)
1955 i_index(i,j,il,il1)=itmp
1962 Corrected AL 8/19/2014: each processor needs whole iset array not only its
1964 c call mpi_scatter(i2set,1,mpi_integer,
1965 c & iset,1,mpi_integer,king,
1967 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1974 if(me.eq.king) close(irest2)
1978 subroutine read1restart_old
1980 include 'DIMENSIONS'
1984 include 'COMMON.LAGRANGE.5diag'
1986 include 'COMMON.LAGRANGE'
1988 include 'COMMON.IOUNITS'
1989 include 'COMMON.REMD'
1990 include 'COMMON.SETUP'
1991 include 'COMMON.CHAIN'
1992 include 'COMMON.SBRIDGE'
1993 include 'COMMON.INTERACT'
1994 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1996 common /przechowalnia/ d_restart1
2000 open(irest2,file=mremd_rst_name,status='unknown')
2001 read (irest2,*) (i2rep(i),i=0,nodes-1)
2002 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2004 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2005 read (irest2,*) ndowna(0,il),
2006 & (ndowna(i,il),i=1,ndowna(0,il))
2009 read(irest2,*) (t_restart1(j,il),j=1,4)
2012 call mpi_scatter(t_restart1,5,mpi_real,
2013 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2017 t_bath=t5_restart1(4)
2022 read(irest2,'(3e15.5)')
2023 & (d_restart1(j,i+2*nres*il),j=1,3)
2027 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2028 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2038 read(irest2,'(3e15.5)')
2039 & (d_restart1(j,i+2*nres*il),j=1,3)
2043 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2044 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2050 if(me.eq.king) close(irest2)
2053 c------------------------------------------