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
70 if(me.eq.king.or..not.out1file) then
71 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
72 write (iout,*) "NREP=",nrep
76 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
77 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
79 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
81 cd print *,'MREMD',nodes
82 cd print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
83 cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath
87 do il1=1,max0(mset(il),1)
103 if(me.eq.king.or..not.out1file) then
104 write(iout,*) (i2rep(i),i=0,nodes-1)
105 write(iout,*) (i2set(i),i=0,nodes-1)
110 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
117 c print *,'i2rep',me,i2rep(me)
118 c print *,'rep2i',(rep2i(i),i=0,nrep)
120 cold if (i2rep(me).eq.nrep) then
123 cold nup(0)=remd_m(i2rep(me)+1)
124 cold k=rep2i(int(i2rep(me)))+1
131 cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
133 cold if (i2rep(me).eq.1) then
136 cold ndown(0)=remd_m(i2rep(me)-1)
137 cold k=rep2i(i2rep(me)-2)+1
144 cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
147 write (*,*) "Processor",me," rest",rest,"
148 & restart1fie",restart1file
149 if(rest.and.restart1file) then
151 & inquire(file=mremd_rst_name,exist=file_exist)
152 cd write (*,*) me," Before broadcast: file_exist",file_exist
153 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
155 cd write (*,*) me," After broadcast: file_exist",file_exist
157 if(me.eq.king.or..not.out1file)
158 & write (iout,*) 'Master is reading restart1file'
159 call read1restart(i_index)
161 if(me.eq.king.or..not.out1file)
162 & write (iout,*) 'WARNING : no restart1file'
165 if(me.eq.king.or..not.out1file) then
166 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
167 write(iout,*) "i_index"
172 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
178 stdfp=dsqrt(2*Rb*t_bath/d_time)
180 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
182 if (lang.gt.0 .and. .not.surfarea) then
184 stdforcp(i)=stdfp*dsqrt(gamp)
187 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
188 & *dsqrt(gamsc(iabs(itype(i))))
194 if (rest.and..not.restart1file)
195 & inquire(file=mremd_rst_name,exist=file_exist)
196 if(.not.file_exist.and.rest.and..not.restart1file)
197 & write(iout,*) 'WARNING : no restart file',mremd_rst_name
198 IF (rest.and.file_exist.and..not.restart1file) THEN
199 write (iout,*) 'Master is reading restart file',
201 open(irest2,file=mremd_rst_name,status='unknown')
203 read (irest2,*) (i2rep(i),i=0,nodes-1)
205 read (irest2,*) (ifirst(i),i=1,remd_m(1))
208 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
210 read (irest2,*) ndowna(0,il),
211 & (ndowna(i,il),i=1,ndowna(0,il))
217 read (irest2,*) (mset(i),i=1,nset)
219 read (irest2,*) (i2set(i),i=0,nodes-1)
224 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
229 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
230 write(iout,*) "i_index"
235 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
244 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
245 write (iout,'(a6,1000i5)') "ifirst",
246 & (ifirst(i),i=1,remd_m(1))
248 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
249 & (nupa(i,il),i=1,nupa(0,il))
250 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
251 & (ndowna(i,il),i=1,ndowna(0,il))
253 stdfp=dsqrt(2*Rb*t_bath/d_time)
255 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
257 if (lang.gt.0 .and. .not.surfarea) then
259 stdforcp(i)=stdfp*dsqrt(gamp)
262 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
263 & *dsqrt(gamsc(iabs(itype(i))))
267 ELSE IF (.not.(rest.and.file_exist)) THEN
273 if (i2rep(il-1).eq.nrep) then
276 nupa(0,il)=remd_m(i2rep(il-1)+1)
277 k=rep2i(int(i2rep(il-1)))+1
283 if (i2rep(il-1).eq.1) then
286 ndowna(0,il)=remd_m(i2rep(il-1)-1)
287 k=rep2i(i2rep(il-1)-2)+1
295 write (iout,'(a6,100i4)') "ifirst",
296 & (ifirst(i),i=1,remd_m(1))
298 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
299 & (nupa(i,il),i=1,nupa(0,il))
300 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
301 & (ndowna(i,il),i=1,ndowna(0,il))
307 c t_bath=retmin+(retmax-retmin)*me/(nodes-1)
308 if(.not.(rest.and.file_exist.and.restart1file)) then
309 if (me .eq. king) then
312 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
314 cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
315 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
320 if(me.eq.king.or..not.out1file)
321 & write(iout,*) me,"iset=",iset,"t_bath=",t_bath
325 stdfp=dsqrt(2*Rb*t_bath/d_time)
327 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
330 c print *,'irep',me,t_bath
332 if (me.eq.king .or. .not. out1file)
333 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
334 call rescale_weights(t_bath)
338 c------copy MD--------------
339 c The driver for molecular dynamics subroutines
340 c------------------------------------------------
346 if(me.eq.king.or..not.out1file)
347 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
353 c Determine the inverse of the inertia matrix.
354 call setup_MD_matrices
358 if (me.eq.king .or. .not. out1file)
359 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
360 stdfp=dsqrt(2*Rb*t_bath/d_time)
362 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
364 call rescale_weights(t_bath)
368 t_MDsetup = MPI_Wtime()-tt0
370 t_MDsetup = tcpu()-tt0
373 c Entering the MD loop
379 if (lang.eq.2 .or. lang.eq.3) then
383 call sd_verlet_p_setup
385 call sd_verlet_ciccotti_setup
389 pfric0_mat(i,j,0)=pfric_mat(i,j)
390 afric0_mat(i,j,0)=afric_mat(i,j)
391 vfric0_mat(i,j,0)=vfric_mat(i,j)
392 prand0_mat(i,j,0)=prand_mat(i,j)
393 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
394 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
399 flag_stoch(i)=.false.
403 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
405 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
409 else if (lang.eq.1 .or. lang.eq.4) then
413 if (me.eq.king .or. .not. out1file)
414 & write(iout,*) 'Setup time',time00-walltime
417 t_langsetup=MPI_Wtime()-tt0
420 t_langsetup=tcpu()-tt0
425 do while(.not.end_of_run)
427 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
428 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
430 if (lang.gt.0 .and. surfarea .and.
431 & mod(itime,reset_fricmat).eq.0) then
432 if (lang.eq.2 .or. lang.eq.3) then
436 call sd_verlet_p_setup
438 call sd_verlet_ciccotti_setup
442 pfric0_mat(i,j,0)=pfric_mat(i,j)
443 afric0_mat(i,j,0)=afric_mat(i,j)
444 vfric0_mat(i,j,0)=vfric_mat(i,j)
445 prand0_mat(i,j,0)=prand_mat(i,j)
446 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
447 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
452 flag_stoch(i)=.false.
455 else if (lang.eq.1 .or. lang.eq.4) then
458 write (iout,'(a,i10)')
459 & "Friction matrix reset based on surface area, itime",itime
461 if (reset_vel .and. tbf .and. lang.eq.0
462 & .and. mod(itime,count_reset_vel).eq.0) then
464 if (me.eq.king .or. .not. out1file)
465 & write(iout,'(a,f20.2)')
466 & "Velocities reset to random values, time",totT
469 d_t_old(j,i)=d_t(j,i)
473 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
477 d_t(j,0)=d_t(j,0)-vcm(j)
480 kinetic_T=2.0d0/(dimen3*Rb)*EK
481 scalfac=dsqrt(T_bath/kinetic_T)
482 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
485 d_t_old(j,i)=scalfac*d_t(j,i)
491 c Time-reversible RESPA algorithm
492 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
493 call RESPA_step(itime)
495 c Variable time step algorithm.
496 call velverlet_step(itime)
500 call brown_step(itime)
502 print *,"Brown dynamics not here!"
504 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
510 if (mod(itime,ntwe).eq.0) call statout(itime)
512 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
513 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
515 call hairpin(.true.,nharp,iharp)
516 call secondary2(.true.)
517 call pdbout(potE,tytul,ipdb)
522 if (mod(itime,ntwx).eq.0.and.traj1file) then
523 if(ntwx_cache.lt.max_cache_traj_use) then
524 ntwx_cache=ntwx_cache+1
526 if (max_cache_traj_use.ne.1)
527 & print *,itime,"processor ",me," over cache ",ntwx_cache
530 totT_cache(i)=totT_cache(i+1)
531 EK_cache(i)=EK_cache(i+1)
532 potE_cache(i)=potE_cache(i+1)
533 t_bath_cache(i)=t_bath_cache(i+1)
534 Uconst_cache(i)=Uconst_cache(i+1)
535 iset_cache(i)=iset_cache(i+1)
538 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
541 qpair_cache(ii,i)=qpair_cache(ii,i+1)
544 utheta_cache(ii,i)=utheta_cache(ii,i+1)
545 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
546 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
552 c_cache(j,ii,i)=c_cache(j,ii,i+1)
558 totT_cache(ntwx_cache)=totT
559 EK_cache(ntwx_cache)=EK
560 potE_cache(ntwx_cache)=potE
561 t_bath_cache(ntwx_cache)=t_bath
562 Uconst_cache(ntwx_cache)=Uconst
563 iset_cache(ntwx_cache)=iset
566 qfrag_cache(i,ntwx_cache)=qfrag(i)
569 qpair_cache(i,ntwx_cache)=qpair(i)
572 utheta_cache(i,ntwx_cache)=utheta(i)
573 ugamma_cache(i,ntwx_cache)=ugamma(i)
574 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
576 C print *,'przed returnbox'
578 C call enerprint(remd_ene(0,i))
581 c_cache(j,i,ntwx_cache)=c(j,i)
586 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
587 & .and..not.restart1file) then
590 open(irest1,file=mremd_rst_name,status='unknown')
591 write (irest1,*) "i2rep"
592 write (irest1,*) (i2rep(i),i=0,nodes-1)
593 write (irest1,*) "ifirst"
594 write (irest1,*) (ifirst(i),i=1,remd_m(1))
596 write (irest1,*) "nupa",il
597 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
598 write (irest1,*) "ndowna",il
599 write (irest1,*) ndowna(0,il),
600 & (ndowna(i,il),i=1,ndowna(0,il))
603 write (irest1,*) "nset"
604 write (irest1,*) nset
605 write (irest1,*) "mset"
606 write (irest1,*) (mset(i),i=1,nset)
607 write (irest1,*) "i2set"
608 write (irest1,*) (i2set(i),i=0,nodes-1)
609 write (irest1,*) "i_index"
613 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
621 open(irest2,file=rest2name,status='unknown')
622 write(irest2,*) totT,EK,potE,totE,t_bath
624 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
627 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
630 write (irest2,*) iset
637 c forced synchronization
638 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
639 & .and. .not. mremdsync) then
641 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
643 call mpi_recv(itime_master, 1, MPI_INTEGER,
644 & 0,101,CG_COMM, status, ierr)
645 call mpi_barrier(CG_COMM, ierr)
646 cdeb if (out1file.or.traj1file) then
647 cdeb call mpi_gather(itime,1,mpi_integer,
648 cdeb & icache_all,1,mpi_integer,king,
651 & call mpi_gather(ntwx_cache,1,mpi_integer,
652 & icache_all,1,mpi_integer,king,
655 & write(iout,*) 'REMD synchro at',itime_master,itime
656 if (itime_master.ge.n_timestep .or. ovrtim())
658 ctime call flush(iout)
663 if ((mod(itime,nstex).eq.0.and.me.eq.king
664 & .or.end_of_run.and.me.eq.king )
665 & .and. .not. mremdsync ) then
668 call mpi_isend(itime,1,MPI_INTEGER,i,101,
669 & CG_COMM, ireqi(i), ierr)
670 cd write(iout,*) 'REMD synchro with',i
673 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
674 call mpi_barrier(CG_COMM, ierr)
676 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
677 if (out1file.or.traj1file) then
678 cdeb call mpi_gather(itime,1,mpi_integer,
679 cdeb & itime_all,1,mpi_integer,king,
681 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
682 cdeb & (itime_all(i),i=1,nodes)
684 cdeb imin_itime=itime_all(1)
686 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
688 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
689 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
690 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
691 call mpi_gather(ntwx_cache,1,mpi_integer,
692 & icache_all,1,mpi_integer,king,
694 c write(iout,'(a19,8000i8)') ' ntwx_cache',
695 c & (icache_all(i),i=1,nodes)
696 ii_write=icache_all(1)
698 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
700 c write(iout,*) "MIN ii_write=",ii_write
703 ctime call flush(iout)
705 if(mremdsync .and. mod(itime,nstex).eq.0) then
707 if (me.eq.king .or. .not. out1file)
708 & write(iout,*) 'REMD synchro at',itime
711 call mpi_gather(ntwx_cache,1,mpi_integer,
712 & icache_all,1,mpi_integer,king,
715 write(iout,'(a19,8000i8)') ' ntwx_cache',
716 & (icache_all(i),i=1,nodes)
717 ii_write=icache_all(1)
719 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
721 write(iout,*) "MIN ii_write=",ii_write
727 c Update the time safety limiy
728 if (time001-time00.gt.safety) then
729 safety=time001-time00+600
730 write (iout,*) "****** SAFETY increased to",safety," s"
732 if (ovrtim()) end_of_run=.true.
734 if(synflag.and..not.end_of_run) then
738 write(iout,*) 'REMD before',me,t_bath
740 c call mpi_gather(t_bath,1,mpi_double_precision,
741 c & remd_t_bath,1,mpi_double_precision,king,
743 potEcomp(n_ene+1)=t_bath
745 potEcomp(n_ene+2)=iset
746 if (iset.lt.nset) then
751 call Econstr_back_qlike
755 potEcomp(n_ene+3)=Uconst+Uconst_back
763 call Econstr_back_qlike
767 potEcomp(n_ene+4)=Uconst+Uconst_back
771 c 9/11/17 Adaptive US
774 write (iout,*) "me ",me," itt",itt
776 if (itt.lt.nrep) then
780 potEcomp(n_ene+5)=Uconst
782 write (iout,*) "t_bath",t_bath_temp,t_bath,
785 if (iset.lt.nset) then
789 potEcomp(n_ene+7)=Uconst
791 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
801 potEcomp(n_ene+6)=Uconst
803 write (iout,*) "t_bath",t_bath_temp,t_bath,
810 potEcomp(n_ene+8)=Uconst
812 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
820 call mpi_gather(potEcomp(0),n_ene+9,mpi_double_precision,
821 & remd_ene(0,1),n_ene+9,mpi_double_precision,king,
824 call mpi_gather(elow,1,mpi_double_precision,
825 & elowi,1,mpi_double_precision,king,
827 call mpi_gather(ehigh,1,mpi_double_precision,
828 & ehighi,1,mpi_double_precision,king,
833 if (me.eq.king .or. .not. out1file) then
834 write(iout,*) 'REMD gather times=',time03-time01
838 if (restart1file) call write1rst(i_index)
841 if (me.eq.king .or. .not. out1file) then
842 write(iout,*) 'REMD writing rst time=',time04-time03
845 if (traj1file) call write1traj
847 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
848 cdeb & icache_all,1,mpi_integer,king,
850 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
851 cdeb & (icache_all(i),i=1,nodes)
856 if (me.eq.king .or. .not. out1file) then
857 write(iout,*) 'REMD writing traj time=',time05-time04
864 remd_t_bath(i)=remd_ene(n_ene+1,i)
865 iremd_iset(i)=remd_ene(n_ene+2,i)
869 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
871 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
875 write(iout,*) 'REMD exchange temp,ene'
877 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
878 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
882 c-------------------------------------
885 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
887 write (iout,*) "remd_m(1)",remd_m(1)
890 i=ifirst(iran_num(1,remd_m(1)))
895 & write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
896 if(i.gt.0.and.nupa(0,i).gt.0) then
898 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
900 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
902 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
904 c do while (iex.eq.i)
905 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
906 iex=nupa(iran_num(1,int(nupa(0,i))),i)
908 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
911 write (iout,*) "i",i,"iex",iex," temperatures",
912 & remd_t_bath(i),remd_t_bath(iex)
915 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
917 c Swap temperatures between conformations i and iex with recalculating the free energies
918 c following temperature changes.
919 ene_iex_iex=remd_ene(0,iex)
920 ene_i_i=remd_ene(0,i)
921 c write (iout,*) "i",i," ene_i_i",ene_i_i,
922 c & " iex",iex," ene_iex_iex",ene_iex_iex
923 c write (iout,*) "rescaling weights with temperature",
926 call rescale_weights(remd_t_bath(i))
928 c write (iout,*) "0,iex",remd_t_bath(i)
929 c call enerprint(remd_ene(0,iex))
931 call sum_energy(remd_ene(0,iex),.false.)
932 ene_iex_i=remd_ene(0,iex)
933 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
935 c write (iout,*) "0,i",remd_t_bath(i)
936 c call enerprint(remd_ene(0,i))
938 call sum_energy(remd_ene(0,i),.false.)
939 c write (iout,*) "ene_i_i",remd_ene(0,i)
941 c write (iout,*) "rescaling weights with temperature",
943 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
944 write (iout,*) "ERROR: inconsistent energies:",i,
945 & ene_i_i,remd_ene(0,i)
947 call rescale_weights(remd_t_bath(iex))
949 c write (iout,*) "0,i",remd_t_bath(iex)
950 c call enerprint(remd_ene(0,i))
952 call sum_energy(remd_ene(0,i),.false.)
953 c write (iout,*) "ene_i_iex",remd_ene(0,i)
955 ene_i_iex=remd_ene(0,i)
957 c write (iout,*) "0,iex",remd_t_bath(iex)
958 c call enerprint(remd_ene(0,iex))
960 call sum_energy(remd_ene(0,iex),.false.)
961 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
962 write (iout,*) "ERROR: inconsistent energies:",iex,
963 & ene_iex_iex,remd_ene(0,iex)
966 write (iout,*) "ene_iex_iex",remd_ene(0,iex)
967 write (iout,*) "i",i," iex",iex
968 write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
969 & " ene_i_iex",ene_i_iex,
970 & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
973 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
974 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
976 c write(iout,*) 'delta',delta
977 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
978 c & (remd_ene(i)-remd_ene(iex))/Rb/
979 c & (remd_t_bath(i)*remd_t_bath(iex))
981 if (delta .gt. 50.0d0) then
987 else if (delta.lt.-50.0d0) then
996 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
997 xxx=ran_number(0.0d0,1.0d0)
998 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1000 if (delta .gt. xxx) then
1002 remd_t_bath(i)=remd_t_bath(iex)
1003 remd_t_bath(iex)=tmp
1004 remd_ene(0,i)=ene_i_iex
1005 remd_ene(0,iex)=ene_iex_i
1011 ehighi(i)=ehighi(iex)
1018 nupa(k,i)=nupa(k,iex)
1021 ndowna(k,i)=ndowna(k,iex)
1025 if (ifirst(il).eq.i) ifirst(il)=iex
1027 if (nupa(k,il).eq.i) then
1029 elseif (nupa(k,il).eq.iex) then
1034 if (ndowna(k,il).eq.i) then
1036 elseif (ndowna(k,il).eq.iex) then
1042 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1044 i2rep(i-1)=i2rep(iex-1)
1047 write(iout,*) 'exchange',i,iex
1048 write (iout,'(a8,100i4)') "@ ifirst",
1049 & (ifirst(k),k=1,remd_m(1))
1051 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1052 & (nupa(k,il),k=1,nupa(0,il))
1053 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1054 & (ndowna(k,il),k=1,ndowna(0,il))
1059 remd_ene(0,iex)=ene_iex_iex
1060 remd_ene(0,i)=ene_i_i
1066 cd write (iout,*) "exchange completed"
1070 cd write(iout,*) "########",ii
1072 i_temp=iran_num(1,nrep)
1073 i_mult=iran_num(1,remd_m(i_temp))
1074 i_iset=iran_num(1,nset)
1075 i_mset=iran_num(1,mset(i_iset))
1076 i=i_index(i_temp,i_mult,i_iset,i_mset)
1078 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1081 cd write(iout,*) "i_dir=",i_dir
1083 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1086 i_mult1=iran_num(1,remd_m(i_temp1))
1088 i_mset1=iran_num(1,mset(i_iset1))
1089 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1090 c 9/1/17 AL: Correction; otherwise the restraint energies are mis-assigned
1091 c on failed replica exchange attempt
1092 econstr_temp_i=remd_ene(20,i)
1093 econstr_temp_iex=remd_ene(20,iex)
1094 c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
1096 remd_ene(20,i)=remd_ene(n_ene+5,i)
1097 remd_ene(20,iex)=remd_ene(n_ene+6,iex)
1099 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1102 i_mult1=iran_num(1,remd_m(i_temp1))
1104 i_mset1=iran_num(1,mset(i_iset1))
1105 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1106 econstr_temp_i=remd_ene(20,i)
1107 econstr_temp_iex=remd_ene(20,iex)
1108 remd_ene(20,i)=remd_ene(n_ene+3,i)
1109 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1111 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1114 i_mult1=iran_num(1,remd_m(i_temp1))
1116 i_mset1=iran_num(1,mset(i_iset1))
1117 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1118 econstr_temp_i=remd_ene(20,i)
1119 econstr_temp_iex=remd_ene(20,iex)
1121 remd_ene(20,i)=remd_ene(n_ene+7,i)
1122 remd_ene(20,iex)=remd_ene(n_ene+8,iex)
1124 remd_ene(20,i)=remd_ene(n_ene+3,i)
1125 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1131 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1134 c Swap temperatures between conformations i and iex with recalculating the free energies
1135 c following temperature changes.
1137 write (iout,*) "i_dir",i_dir," i",i," iex",iex,
1138 & " t_bath",remd_t_bath(i),remd_t_bath(iex)
1140 ene_iex_iex=remd_ene(0,iex)
1141 ene_i_i=remd_ene(0,i)
1142 co write (iout,*) "rescaling weights with temperature",
1144 call rescale_weights(remd_t_bath(i))
1146 call sum_energy(remd_ene(0,iex),.false.)
1147 ene_iex_i=remd_ene(0,iex)
1148 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
1149 c call sum_energy(remd_ene(0,i),.false.)
1150 cd write (iout,*) "ene_i_i",remd_ene(0,i)
1151 c write (iout,*) "rescaling weights with temperature",
1152 c & remd_t_bath(iex)
1153 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1154 c write (iout,*) "ERROR: inconsistent energies:",i,
1155 c & ene_i_i,remd_ene(0,i)
1157 call rescale_weights(remd_t_bath(iex))
1158 call sum_energy(remd_ene(0,i),.false.)
1159 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1160 ene_i_iex=remd_ene(0,i)
1161 c call sum_energy(remd_ene(0,iex),.false.)
1162 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1163 c write (iout,*) "ERROR: inconsistent energies:",iex,
1164 c & ene_iex_iex,remd_ene(0,iex)
1166 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1167 c write (iout,*) "i",i," iex",iex
1168 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1169 cd & " ene_i_iex",ene_i_iex,
1170 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1171 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1172 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1174 cd write(iout,*) 'delta',delta
1175 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1176 c & (remd_ene(i)-remd_ene(iex))/Rb/
1177 c & (remd_t_bath(i)*remd_t_bath(iex))
1178 if (delta .gt. 50.0d0) then
1183 if (i_dir.eq.1.or.i_dir.eq.3)
1184 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1185 if (i_dir.eq.2.or.i_dir.eq.3)
1186 & iremd_tot_usa(int(i2set(i-1)))=
1187 & iremd_tot_usa(int(i2set(i-1)))+1
1188 xxx=ran_number(0.0d0,1.0d0)
1189 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1190 if (delta .gt. xxx) then
1192 remd_t_bath(i)=remd_t_bath(iex)
1193 remd_t_bath(iex)=tmp
1196 iremd_iset(i)=iremd_iset(iex)
1197 iremd_iset(iex)=itmp
1199 remd_ene(0,i)=ene_i_iex
1200 remd_ene(0,iex)=ene_iex_i
1202 if (i_dir.eq.1.or.i_dir.eq.3)
1203 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1206 i2rep(i-1)=i2rep(iex-1)
1209 if (i_dir.eq.2.or.i_dir.eq.3)
1210 & iremd_acc_usa(int(i2set(i-1)))=
1211 & iremd_acc_usa(int(i2set(i-1)))+1
1214 i2set(i-1)=i2set(iex-1)
1217 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1218 i_index(i_temp,i_mult,i_iset,i_mset)=
1219 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1220 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1223 remd_ene(0,iex)=ene_iex_iex
1224 remd_ene(0,i)=ene_i_i
1225 remd_ene(20,iex)=econstr_temp_iex
1226 remd_ene(20,i)=econstr_temp_i
1230 cd do il1=1,mset(il)
1233 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1246 c-------------------------------------
1247 write (iout,*) "NREP",nrep
1249 if(iremd_tot(i).ne.0)
1250 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1251 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1256 if(iremd_tot_usa(i).ne.0)
1257 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1258 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1264 cd write (iout,'(a6,100i4)') "ifirst",
1265 cd & (ifirst(i),i=1,remd_m(1))
1267 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1268 cd & (nupa(i,il),i=1,nupa(0,il))
1269 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1270 cd & (ndowna(i,il),i=1,ndowna(0,il))
1275 cd write (iout,*) "Before scatter"
1277 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1278 & t_bath,1,mpi_double_precision,king,
1280 cd write (iout,*) "After scatter"
1283 & call mpi_scatter(iremd_iset,1,mpi_integer,
1284 & iset,1,mpi_integer,king,
1288 if (me.eq.king .or. .not. out1file) then
1289 write(iout,*) 'REMD scatter time=',time07-time06
1293 call mpi_scatter(elowi,1,mpi_double_precision,
1294 & elow,1,mpi_double_precision,king,
1296 call mpi_scatter(ehighi,1,mpi_double_precision,
1297 & ehigh,1,mpi_double_precision,king,
1300 call rescale_weights(t_bath)
1301 co write (iout,*) "Processor",me,
1302 co & " rescaling weights with temperature",t_bath
1304 stdfp=dsqrt(2*Rb*t_bath/d_time)
1306 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1309 c Compute the standard deviations of stochastic forces for Langevin dynamics
1310 c if the friction coefficients do not depend on surface area
1311 if (lang.gt.0 .and. .not.surfarea) then
1313 stdforcp(i)=stdfp*dsqrt(gamp)
1316 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1317 & *dsqrt(gamsc(iabs(itype(i))))
1320 cde write(iout,*) 'REMD after',me,t_bath
1322 if (me.eq.king .or. .not. out1file) then
1323 write(iout,*) 'REMD exchange time=',time08-time00
1329 if (restart1file) then
1330 if (me.eq.king .or. .not. out1file)
1331 & write(iout,*) 'writing restart at the end of run'
1332 call write1rst(i_index)
1335 if (traj1file) call write1traj
1337 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1338 cdeb & icache_all,1,mpi_integer,king,
1339 cdeb & CG_COMM,ierr)
1340 cdeb write(iout,'(a40,8000i8)')
1341 cdeb & ' ntwx_cache after traj1file at the end',
1342 cdeb & (icache_all(i),i=1,nodes)
1347 t_MD=MPI_Wtime()-tt0
1351 if (me.eq.king .or. .not. out1file) then
1352 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1354 & 'MD calculations setup:',t_MDsetup,
1355 & 'Energy & gradient evaluation:',t_enegrad,
1356 & 'Stochastic MD setup:',t_langsetup,
1357 & 'Stochastic MD step setup:',t_sdsetup,
1359 write (iout,'(/28(1h=),a25,27(1h=))')
1360 & ' End of MD calculation '
1365 c-----------------------------------------------------------------------
1366 subroutine write1rst(i_index)
1368 include 'DIMENSIONS'
1370 include 'COMMON.CONTROL'
1373 include 'COMMON.LAGRANGE.5diag'
1375 include 'COMMON.LAGRANGE'
1377 include 'COMMON.QRESTR'
1378 include 'COMMON.IOUNITS'
1379 include 'COMMON.REMD'
1380 include 'COMMON.SETUP'
1381 include 'COMMON.CHAIN'
1382 include 'COMMON.SBRIDGE'
1383 include 'COMMON.INTERACT'
1385 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1386 & d_restart2(3,2*maxres*maxprocs)
1390 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1391 common /przechowalnia/ d_restart1,d_restart2
1392 integer i,j,il1,il,ixdrf
1398 t5_restart1(4)=t_bath
1399 t5_restart1(5)=Uconst
1401 call mpi_gather(t5_restart1,5,mpi_real,
1402 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1410 call mpi_gather(r_d,3*2*nres,mpi_real,
1411 & d_restart1,3*2*nres,mpi_real,king,
1420 call mpi_gather(r_d,3*2*nres,mpi_real,
1421 & d_restart2,3*2*nres,mpi_real,king,
1426 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1428 call xdrfint_(ixdrf, i2rep(i), iret)
1431 call xdrfint_(ixdrf, ifirst(i), iret)
1435 call xdrfint_(ixdrf, nupa(i,il), iret)
1439 call xdrfint_(ixdrf, ndowna(i,il), iret)
1445 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1452 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1459 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1465 call xdrfint_(ixdrf, nset, iret)
1467 call xdrfint_(ixdrf,mset(i), iret)
1470 call xdrfint_(ixdrf,i2set(i), iret)
1476 itmp=i_index(i,j,il,il1)
1477 call xdrfint_(ixdrf,itmp, iret)
1484 call xdrfclose_(ixdrf, iret)
1486 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1488 call xdrfint(ixdrf, i2rep(i), iret)
1491 call xdrfint(ixdrf, ifirst(i), iret)
1495 call xdrfint(ixdrf, nupa(i,il), iret)
1499 call xdrfint(ixdrf, ndowna(i,il), iret)
1505 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1512 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1519 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1526 call xdrfint(ixdrf, nset, iret)
1528 call xdrfint(ixdrf,mset(i), iret)
1531 call xdrfint(ixdrf,i2set(i), iret)
1537 itmp=i_index(i,j,il,il1)
1538 call xdrfint(ixdrf,itmp, iret)
1545 call xdrfclose(ixdrf, iret)
1552 subroutine write1traj
1554 include 'DIMENSIONS'
1557 include 'COMMON.QRESTR'
1558 include 'COMMON.IOUNITS'
1559 include 'COMMON.REMD'
1560 include 'COMMON.SETUP'
1561 include 'COMMON.CHAIN'
1562 include 'COMMON.SBRIDGE'
1563 include 'COMMON.INTERACT'
1567 real xcoord(3,maxres2+2),prec
1568 real r_qfrag(50),r_qpair(100)
1569 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1570 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1571 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1572 & p_uscdiff(100*maxprocs)
1573 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1574 common /przechowalnia/ p_c
1575 integer ii,i,il,j,ixdrf
1578 call mpi_bcast(ii_write,1,mpi_integer,
1579 & king,CG_COMM,ierr)
1582 c print *,'traj1file',me,ii_write,ntwx_cache
1586 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1588 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1591 t5_restart1(1)=totT_cache(ii)
1592 t5_restart1(2)=EK_cache(ii)
1593 t5_restart1(3)=potE_cache(ii)
1594 t5_restart1(4)=t_bath_cache(ii)
1595 t5_restart1(5)=Uconst_cache(ii)
1596 call mpi_gather(t5_restart1,5,mpi_real,
1597 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1599 call mpi_gather(iset_cache(ii),1,mpi_integer,
1600 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1603 r_qfrag(i)=qfrag_cache(i,ii)
1606 r_qpair(i)=qpair_cache(i,ii)
1609 r_utheta(i)=utheta_cache(i,ii)
1610 r_ugamma(i)=ugamma_cache(i,ii)
1611 r_uscdiff(i)=uscdiff_cache(i,ii)
1614 call mpi_gather(r_qfrag,nfrag,mpi_real,
1615 & p_qfrag,nfrag,mpi_real,king,
1617 call mpi_gather(r_qpair,npair,mpi_real,
1618 & p_qpair,npair,mpi_real,king,
1620 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1621 & p_utheta,nfrag_back,mpi_real,king,
1623 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1624 & p_ugamma,nfrag_back,mpi_real,king,
1626 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1627 & p_uscdiff,nfrag_back,mpi_real,king,
1631 write (iout,*) "p_qfrag"
1633 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1635 write (iout,*) "p_qpair"
1637 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1643 r_c(j,i)=c_cache(j,i,ii)
1647 call mpi_gather(r_c,3*2*nres,mpi_real,
1648 & p_c,3*2*nres,mpi_real,king,
1654 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1655 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1656 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1657 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1658 call xdrfint_(ixdrf, nss, iret)
1661 call xdrfint(ixdrf, idssb(j)+nres, iret)
1662 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1664 call xdrfint_(ixdrf, ihpb(j), iret)
1665 call xdrfint_(ixdrf, jhpb(j), iret)
1668 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1669 call xdrfint_(ixdrf, iset_restart1(il), iret)
1671 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1674 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1677 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1678 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1679 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1684 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1689 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1693 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1697 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1698 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1699 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1700 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1701 call xdrfint(ixdrf, nss, iret)
1704 call xdrfint(ixdrf, idssb(j)+nres, iret)
1705 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1707 call xdrfint(ixdrf, ihpb(j), iret)
1708 call xdrfint(ixdrf, jhpb(j), iret)
1711 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1712 call xdrfint(ixdrf, iset_restart1(il), iret)
1714 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1717 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1720 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1721 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1722 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1727 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1732 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1736 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1742 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1744 if(me.eq.king) call xdrfclose(ixdrf, iret)
1746 do i=1,ntwx_cache-ii_write
1748 totT_cache(i)=totT_cache(ii_write+i)
1749 EK_cache(i)=EK_cache(ii_write+i)
1750 potE_cache(i)=potE_cache(ii_write+i)
1751 t_bath_cache(i)=t_bath_cache(ii_write+i)
1752 Uconst_cache(i)=Uconst_cache(ii_write+i)
1753 iset_cache(i)=iset_cache(ii_write+i)
1756 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1759 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1762 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1763 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1764 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1769 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1773 ntwx_cache=ntwx_cache-ii_write
1778 subroutine read1restart(i_index)
1780 include 'DIMENSIONS'
1782 include 'COMMON.CONTROL'
1785 include 'COMMON.LAGRANGE.5diag'
1787 include 'COMMON.LAGRANGE'
1789 include 'COMMON.QRESTR'
1790 include 'COMMON.IOUNITS'
1791 include 'COMMON.REMD'
1792 include 'COMMON.SETUP'
1793 include 'COMMON.CHAIN'
1794 include 'COMMON.SBRIDGE'
1795 include 'COMMON.INTERACT'
1796 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1799 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1800 common /przechowalnia/ d_restart1
1801 integer i,j,il,il1,ixdrf,iret,itmp
1803 write (*,*) "Processor",me," called read1restart"
1806 open(irest2,file=mremd_rst_name,status='unknown')
1807 read(irest2,*,err=334) i
1808 write(iout,*) "Reading old rst in ASCI format"
1810 call read1restart_old
1814 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1817 call xdrfint_(ixdrf, i2rep(i), iret)
1820 call xdrfint_(ixdrf, ifirst(i), iret)
1823 call xdrfint_(ixdrf, nupa(0,il), iret)
1825 call xdrfint_(ixdrf, nupa(i,il), iret)
1828 call xdrfint_(ixdrf, ndowna(0,il), iret)
1830 call xdrfint_(ixdrf, ndowna(i,il), iret)
1835 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1839 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1842 call xdrfint(ixdrf, i2rep(i), iret)
1845 call xdrfint(ixdrf, ifirst(i), iret)
1848 call xdrfint(ixdrf, nupa(0,il), iret)
1850 call xdrfint(ixdrf, nupa(i,il), iret)
1853 call xdrfint(ixdrf, ndowna(0,il), iret)
1855 call xdrfint(ixdrf, ndowna(i,il), iret)
1860 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1865 call mpi_scatter(t_restart1,5,mpi_real,
1866 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1870 t_bath=t5_restart1(4)
1875 c read(irest2,'(3e15.5)')
1876 c & (d_restart1(j,i+2*nres*il),j=1,3)
1879 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1881 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1887 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1888 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1898 c read(irest2,'(3e15.5)')
1899 c & (d_restart1(j,i+2*nres*il),j=1,3)
1902 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1904 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1910 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1911 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1922 call xdrfint_(ixdrf, nset, iret)
1924 call xdrfint_(ixdrf,mset(i), iret)
1927 call xdrfint_(ixdrf,i2set(i), iret)
1933 call xdrfint_(ixdrf,itmp, iret)
1934 i_index(i,j,il,il1)=itmp
1942 call xdrfint(ixdrf, nset, iret)
1944 call xdrfint(ixdrf,mset(i), iret)
1947 call xdrfint(ixdrf,i2set(i), iret)
1953 call xdrfint(ixdrf,itmp, iret)
1954 i_index(i,j,il,il1)=itmp
1961 Corrected AL 8/19/2014: each processor needs whole iset array not only its
1963 c call mpi_scatter(i2set,1,mpi_integer,
1964 c & iset,1,mpi_integer,king,
1966 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1973 if(me.eq.king) close(irest2)
1977 subroutine read1restart_old
1979 include 'DIMENSIONS'
1983 include 'COMMON.LAGRANGE.5diag'
1985 include 'COMMON.LAGRANGE'
1987 include 'COMMON.IOUNITS'
1988 include 'COMMON.REMD'
1989 include 'COMMON.SETUP'
1990 include 'COMMON.CHAIN'
1991 include 'COMMON.SBRIDGE'
1992 include 'COMMON.INTERACT'
1993 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1995 common /przechowalnia/ d_restart1
1999 open(irest2,file=mremd_rst_name,status='unknown')
2000 read (irest2,*) (i2rep(i),i=0,nodes-1)
2001 read (irest2,*) (ifirst(i),i=1,remd_m(1))
2003 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
2004 read (irest2,*) ndowna(0,il),
2005 & (ndowna(i,il),i=1,ndowna(0,il))
2008 read(irest2,*) (t_restart1(j,il),j=1,4)
2011 call mpi_scatter(t_restart1,5,mpi_real,
2012 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
2016 t_bath=t5_restart1(4)
2021 read(irest2,'(3e15.5)')
2022 & (d_restart1(j,i+2*nres*il),j=1,3)
2026 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2027 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2037 read(irest2,'(3e15.5)')
2038 & (d_restart1(j,i+2*nres*il),j=1,3)
2042 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2043 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2049 if(me.eq.king) close(irest2)
2052 c------------------------------------------