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)
181 if (rest.and..not.restart1file)
182 & inquire(file=mremd_rst_name,exist=file_exist)
183 if(.not.file_exist.and.rest.and..not.restart1file)
184 & write(iout,*) 'WARNING : no restart file',mremd_rst_name
185 IF (rest.and.file_exist.and..not.restart1file) THEN
186 write (iout,*) 'Master is reading restart file',
188 open(irest2,file=mremd_rst_name,status='unknown')
190 read (irest2,*) (i2rep(i),i=0,nodes-1)
192 read (irest2,*) (ifirst(i),i=1,remd_m(1))
195 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
197 read (irest2,*) ndowna(0,il),
198 & (ndowna(i,il),i=1,ndowna(0,il))
204 read (irest2,*) (mset(i),i=1,nset)
206 read (irest2,*) (i2set(i),i=0,nodes-1)
211 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
216 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
217 write(iout,*) "i_index"
222 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
231 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
232 write (iout,'(a6,1000i5)') "ifirst",
233 & (ifirst(i),i=1,remd_m(1))
235 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
236 & (nupa(i,il),i=1,nupa(0,il))
237 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
238 & (ndowna(i,il),i=1,ndowna(0,il))
240 ELSE IF (.not.(rest.and.file_exist)) THEN
246 if (i2rep(il-1).eq.nrep) then
249 nupa(0,il)=remd_m(i2rep(il-1)+1)
250 k=rep2i(int(i2rep(il-1)))+1
256 if (i2rep(il-1).eq.1) then
259 ndowna(0,il)=remd_m(i2rep(il-1)-1)
260 k=rep2i(i2rep(il-1)-2)+1
268 write (iout,'(a6,100i4)') "ifirst",
269 & (ifirst(i),i=1,remd_m(1))
271 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
272 & (nupa(i,il),i=1,nupa(0,il))
273 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
274 & (ndowna(i,il),i=1,ndowna(0,il))
280 c t_bath=retmin+(retmax-retmin)*me/(nodes-1)
281 if(.not.(rest.and.file_exist.and.restart1file)) then
282 if (me .eq. king) then
285 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
287 cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
288 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
293 if(me.eq.king.or..not.out1file)
294 & write(iout,*) me,"iset=",iset,"t_bath=",t_bath
298 stdfp=dsqrt(2*Rb*t_bath/d_time)
300 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
303 c print *,'irep',me,t_bath
305 if (me.eq.king .or. .not. out1file)
306 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
307 call rescale_weights(t_bath)
311 c------copy MD--------------
312 c The driver for molecular dynamics subroutines
313 c------------------------------------------------
319 if(me.eq.king.or..not.out1file)
320 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
326 c Determine the inverse of the inertia matrix.
327 call setup_MD_matrices
331 if (me.eq.king .or. .not. out1file)
332 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
333 stdfp=dsqrt(2*Rb*t_bath/d_time)
335 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
337 call rescale_weights(t_bath)
341 t_MDsetup = MPI_Wtime()-tt0
343 t_MDsetup = tcpu()-tt0
346 c Entering the MD loop
352 if (lang.eq.2 .or. lang.eq.3) then
356 call sd_verlet_p_setup
358 call sd_verlet_ciccotti_setup
362 pfric0_mat(i,j,0)=pfric_mat(i,j)
363 afric0_mat(i,j,0)=afric_mat(i,j)
364 vfric0_mat(i,j,0)=vfric_mat(i,j)
365 prand0_mat(i,j,0)=prand_mat(i,j)
366 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
367 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
372 flag_stoch(i)=.false.
376 & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
378 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
382 else if (lang.eq.1 .or. lang.eq.4) then
386 if (me.eq.king .or. .not. out1file)
387 & write(iout,*) 'Setup time',time00-walltime
390 t_langsetup=MPI_Wtime()-tt0
393 t_langsetup=tcpu()-tt0
398 do while(.not.end_of_run)
400 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
401 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
403 if (lang.gt.0 .and. surfarea .and.
404 & mod(itime,reset_fricmat).eq.0) then
405 if (lang.eq.2 .or. lang.eq.3) then
409 call sd_verlet_p_setup
411 call sd_verlet_ciccotti_setup
415 pfric0_mat(i,j,0)=pfric_mat(i,j)
416 afric0_mat(i,j,0)=afric_mat(i,j)
417 vfric0_mat(i,j,0)=vfric_mat(i,j)
418 prand0_mat(i,j,0)=prand_mat(i,j)
419 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
420 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
425 flag_stoch(i)=.false.
428 else if (lang.eq.1 .or. lang.eq.4) then
431 write (iout,'(a,i10)')
432 & "Friction matrix reset based on surface area, itime",itime
434 if (reset_vel .and. tbf .and. lang.eq.0
435 & .and. mod(itime,count_reset_vel).eq.0) then
437 if (me.eq.king .or. .not. out1file)
438 & write(iout,'(a,f20.2)')
439 & "Velocities reset to random values, time",totT
442 d_t_old(j,i)=d_t(j,i)
446 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
450 d_t(j,0)=d_t(j,0)-vcm(j)
453 kinetic_T=2.0d0/(dimen3*Rb)*EK
454 scalfac=dsqrt(T_bath/kinetic_T)
455 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
458 d_t_old(j,i)=scalfac*d_t(j,i)
464 c Time-reversible RESPA algorithm
465 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
466 call RESPA_step(itime)
468 c Variable time step algorithm.
469 call velverlet_step(itime)
473 call brown_step(itime)
475 print *,"Brown dynamics not here!"
477 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
483 if (mod(itime,ntwe).eq.0) call statout(itime)
485 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
486 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
488 call hairpin(.true.,nharp,iharp)
489 call secondary2(.true.)
490 call pdbout(potE,tytul,ipdb)
495 if (mod(itime,ntwx).eq.0.and.traj1file) then
496 if(ntwx_cache.lt.max_cache_traj_use) then
497 ntwx_cache=ntwx_cache+1
499 if (max_cache_traj_use.ne.1)
500 & print *,itime,"processor ",me," over cache ",ntwx_cache
503 totT_cache(i)=totT_cache(i+1)
504 EK_cache(i)=EK_cache(i+1)
505 potE_cache(i)=potE_cache(i+1)
506 t_bath_cache(i)=t_bath_cache(i+1)
507 Uconst_cache(i)=Uconst_cache(i+1)
508 iset_cache(i)=iset_cache(i+1)
511 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
514 qpair_cache(ii,i)=qpair_cache(ii,i+1)
517 utheta_cache(ii,i)=utheta_cache(ii,i+1)
518 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
519 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
525 c_cache(j,ii,i)=c_cache(j,ii,i+1)
531 totT_cache(ntwx_cache)=totT
532 EK_cache(ntwx_cache)=EK
533 potE_cache(ntwx_cache)=potE
534 t_bath_cache(ntwx_cache)=t_bath
535 Uconst_cache(ntwx_cache)=Uconst
536 iset_cache(ntwx_cache)=iset
539 qfrag_cache(i,ntwx_cache)=qfrag(i)
542 qpair_cache(i,ntwx_cache)=qpair(i)
545 utheta_cache(i,ntwx_cache)=utheta(i)
546 ugamma_cache(i,ntwx_cache)=ugamma(i)
547 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
549 C print *,'przed returnbox'
551 C call enerprint(remd_ene(0,i))
554 c_cache(j,i,ntwx_cache)=c(j,i)
559 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
560 & .and..not.restart1file) then
563 open(irest1,file=mremd_rst_name,status='unknown')
564 write (irest1,*) "i2rep"
565 write (irest1,*) (i2rep(i),i=0,nodes-1)
566 write (irest1,*) "ifirst"
567 write (irest1,*) (ifirst(i),i=1,remd_m(1))
569 write (irest1,*) "nupa",il
570 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
571 write (irest1,*) "ndowna",il
572 write (irest1,*) ndowna(0,il),
573 & (ndowna(i,il),i=1,ndowna(0,il))
576 write (irest1,*) "nset"
577 write (irest1,*) nset
578 write (irest1,*) "mset"
579 write (irest1,*) (mset(i),i=1,nset)
580 write (irest1,*) "i2set"
581 write (irest1,*) (i2set(i),i=0,nodes-1)
582 write (irest1,*) "i_index"
586 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
594 open(irest2,file=rest2name,status='unknown')
595 write(irest2,*) totT,EK,potE,totE,t_bath
597 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
600 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
603 write (irest2,*) iset
610 c forced synchronization
611 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
612 & .and. .not. mremdsync) then
614 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
616 call mpi_recv(itime_master, 1, MPI_INTEGER,
617 & 0,101,CG_COMM, status, ierr)
618 call mpi_barrier(CG_COMM, ierr)
619 cdeb if (out1file.or.traj1file) then
620 cdeb call mpi_gather(itime,1,mpi_integer,
621 cdeb & icache_all,1,mpi_integer,king,
624 & call mpi_gather(ntwx_cache,1,mpi_integer,
625 & icache_all,1,mpi_integer,king,
628 & write(iout,*) 'REMD synchro at',itime_master,itime
629 if (itime_master.ge.n_timestep .or. ovrtim())
631 ctime call flush(iout)
636 if ((mod(itime,nstex).eq.0.and.me.eq.king
637 & .or.end_of_run.and.me.eq.king )
638 & .and. .not. mremdsync ) then
641 call mpi_isend(itime,1,MPI_INTEGER,i,101,
642 & CG_COMM, ireqi(i), ierr)
643 cd write(iout,*) 'REMD synchro with',i
646 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
647 call mpi_barrier(CG_COMM, ierr)
649 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
650 if (out1file.or.traj1file) then
651 cdeb call mpi_gather(itime,1,mpi_integer,
652 cdeb & itime_all,1,mpi_integer,king,
654 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
655 cdeb & (itime_all(i),i=1,nodes)
657 cdeb imin_itime=itime_all(1)
659 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
661 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
662 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
663 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
664 call mpi_gather(ntwx_cache,1,mpi_integer,
665 & icache_all,1,mpi_integer,king,
667 c write(iout,'(a19,8000i8)') ' ntwx_cache',
668 c & (icache_all(i),i=1,nodes)
669 ii_write=icache_all(1)
671 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
673 c write(iout,*) "MIN ii_write=",ii_write
676 ctime call flush(iout)
678 if(mremdsync .and. mod(itime,nstex).eq.0) then
680 if (me.eq.king .or. .not. out1file)
681 & write(iout,*) 'REMD synchro at',itime
684 call mpi_gather(ntwx_cache,1,mpi_integer,
685 & icache_all,1,mpi_integer,king,
688 write(iout,'(a19,8000i8)') ' ntwx_cache',
689 & (icache_all(i),i=1,nodes)
690 ii_write=icache_all(1)
692 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
694 write(iout,*) "MIN ii_write=",ii_write
700 c Update the time safety limiy
701 if (time001-time00.gt.safety) then
702 safety=time001-time00+600
703 write (iout,*) "****** SAFETY increased to",safety," s"
705 if (ovrtim()) end_of_run=.true.
707 if(synflag.and..not.end_of_run) then
711 write(iout,*) 'REMD before',me,t_bath
713 c call mpi_gather(t_bath,1,mpi_double_precision,
714 c & remd_t_bath,1,mpi_double_precision,king,
716 potEcomp(n_ene+1)=t_bath
718 potEcomp(n_ene+2)=iset
719 if (iset.lt.nset) then
724 call Econstr_back_qlike
728 potEcomp(n_ene+3)=Uconst+Uconst_back
736 call Econstr_back_qlike
740 potEcomp(n_ene+4)=Uconst+Uconst_back
744 c 9/11/17 Adaptive US
747 write (iout,*) "me ",me," itt",itt
749 if (itt.lt.nrep) then
753 potEcomp(n_ene+5)=Uconst
755 write (iout,*) "t_bath",t_bath_temp,t_bath,
758 if (iset.lt.nset) then
762 potEcomp(n_ene+7)=Uconst
764 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
774 potEcomp(n_ene+6)=Uconst
776 write (iout,*) "t_bath",t_bath_temp,t_bath,
783 potEcomp(n_ene+8)=Uconst
785 write (iout,*)"iset",i_set_temp,iset," Uconst",Uconst
793 call mpi_gather(potEcomp(0),n_ene+9,mpi_double_precision,
794 & remd_ene(0,1),n_ene+9,mpi_double_precision,king,
797 call mpi_gather(elow,1,mpi_double_precision,
798 & elowi,1,mpi_double_precision,king,
800 call mpi_gather(ehigh,1,mpi_double_precision,
801 & ehighi,1,mpi_double_precision,king,
806 if (me.eq.king .or. .not. out1file) then
807 write(iout,*) 'REMD gather times=',time03-time01
811 if (restart1file) call write1rst(i_index)
814 if (me.eq.king .or. .not. out1file) then
815 write(iout,*) 'REMD writing rst time=',time04-time03
818 if (traj1file) call write1traj
820 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
821 cdeb & icache_all,1,mpi_integer,king,
823 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
824 cdeb & (icache_all(i),i=1,nodes)
829 if (me.eq.king .or. .not. out1file) then
830 write(iout,*) 'REMD writing traj time=',time05-time04
837 remd_t_bath(i)=remd_ene(n_ene+1,i)
838 iremd_iset(i)=remd_ene(n_ene+2,i)
842 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
844 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
848 write(iout,*) 'REMD exchange temp,ene'
850 write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
851 write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
855 c-------------------------------------
858 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
860 write (iout,*) "remd_m(1)",remd_m(1)
863 i=ifirst(iran_num(1,remd_m(1)))
868 & write (iout,*) "i",i," nupa(0,i)",int(nupa(0,i))
869 if(i.gt.0.and.nupa(0,i).gt.0) then
871 c if (i.eq.1 .and. int(nupa(0,i)).eq.1) then
873 c & "CHUJ ABSOLUTNY!!! No way to sample a distinct replica in MREMD"
875 c call MPI_Abort(MPI_COMM_WORLD,ERRCODE,ierr)
877 c do while (iex.eq.i)
878 c write (iout,*) "upper",nupa(int(nupa(0,i)),i)
879 iex=nupa(iran_num(1,int(nupa(0,i))),i)
881 c write (iout,*) "nupa(0,i)",nupa(0,i)," iex",iex
884 write (iout,*) "i",i,"iex",iex," temperatures",
885 & remd_t_bath(i),remd_t_bath(iex)
888 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
890 c Swap temperatures between conformations i and iex with recalculating the free energies
891 c following temperature changes.
892 ene_iex_iex=remd_ene(0,iex)
893 ene_i_i=remd_ene(0,i)
894 c write (iout,*) "i",i," ene_i_i",ene_i_i,
895 c & " iex",iex," ene_iex_iex",ene_iex_iex
896 c write (iout,*) "rescaling weights with temperature",
899 call rescale_weights(remd_t_bath(i))
901 c write (iout,*) "0,iex",remd_t_bath(i)
902 c call enerprint(remd_ene(0,iex))
904 call sum_energy(remd_ene(0,iex),.false.)
905 ene_iex_i=remd_ene(0,iex)
906 c write (iout,*) "ene_iex_i",remd_ene(0,iex)
908 c write (iout,*) "0,i",remd_t_bath(i)
909 c call enerprint(remd_ene(0,i))
911 call sum_energy(remd_ene(0,i),.false.)
912 c write (iout,*) "ene_i_i",remd_ene(0,i)
914 c write (iout,*) "rescaling weights with temperature",
916 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
917 write (iout,*) "ERROR: inconsistent energies:",i,
918 & ene_i_i,remd_ene(0,i)
920 call rescale_weights(remd_t_bath(iex))
922 c write (iout,*) "0,i",remd_t_bath(iex)
923 c call enerprint(remd_ene(0,i))
925 call sum_energy(remd_ene(0,i),.false.)
926 c write (iout,*) "ene_i_iex",remd_ene(0,i)
928 ene_i_iex=remd_ene(0,i)
930 c write (iout,*) "0,iex",remd_t_bath(iex)
931 c call enerprint(remd_ene(0,iex))
933 call sum_energy(remd_ene(0,iex),.false.)
934 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
935 write (iout,*) "ERROR: inconsistent energies:",iex,
936 & ene_iex_iex,remd_ene(0,iex)
939 write (iout,*) "ene_iex_iex",remd_ene(0,iex)
940 write (iout,*) "i",i," iex",iex
941 write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
942 & " ene_i_iex",ene_i_iex,
943 & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
946 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
947 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
949 c write(iout,*) 'delta',delta
950 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
951 c & (remd_ene(i)-remd_ene(iex))/Rb/
952 c & (remd_t_bath(i)*remd_t_bath(iex))
954 if (delta .gt. 50.0d0) then
960 else if (delta.lt.-50.0d0) then
969 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
970 xxx=ran_number(0.0d0,1.0d0)
971 c write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
973 if (delta .gt. xxx) then
975 remd_t_bath(i)=remd_t_bath(iex)
977 remd_ene(0,i)=ene_i_iex
978 remd_ene(0,iex)=ene_iex_i
984 ehighi(i)=ehighi(iex)
991 nupa(k,i)=nupa(k,iex)
994 ndowna(k,i)=ndowna(k,iex)
998 if (ifirst(il).eq.i) ifirst(il)=iex
1000 if (nupa(k,il).eq.i) then
1002 elseif (nupa(k,il).eq.iex) then
1007 if (ndowna(k,il).eq.i) then
1009 elseif (ndowna(k,il).eq.iex) then
1015 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1017 i2rep(i-1)=i2rep(iex-1)
1020 write(iout,*) 'exchange',i,iex
1021 write (iout,'(a8,100i4)') "@ ifirst",
1022 & (ifirst(k),k=1,remd_m(1))
1024 write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
1025 & (nupa(k,il),k=1,nupa(0,il))
1026 write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
1027 & (ndowna(k,il),k=1,ndowna(0,il))
1032 remd_ene(0,iex)=ene_iex_iex
1033 remd_ene(0,i)=ene_i_i
1039 cd write (iout,*) "exchange completed"
1043 cd write(iout,*) "########",ii
1045 i_temp=iran_num(1,nrep)
1046 i_mult=iran_num(1,remd_m(i_temp))
1047 i_iset=iran_num(1,nset)
1048 i_mset=iran_num(1,mset(i_iset))
1049 i=i_index(i_temp,i_mult,i_iset,i_mset)
1051 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
1054 cd write(iout,*) "i_dir=",i_dir
1056 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
1059 i_mult1=iran_num(1,remd_m(i_temp1))
1061 i_mset1=iran_num(1,mset(i_iset1))
1062 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1063 c 9/1/17 AL: Correction; otherwise the restraint energies are mis-assigned
1064 c on failed replica exchange attempt
1065 econstr_temp_i=remd_ene(20,i)
1066 econstr_temp_iex=remd_ene(20,iex)
1067 c 9/11/17 AL: Adaptive sampling (temperature dependent restraints potentials)
1069 remd_ene(20,i)=remd_ene(n_ene+5,i)
1070 remd_ene(20,iex)=remd_ene(n_ene+6,iex)
1072 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
1075 i_mult1=iran_num(1,remd_m(i_temp1))
1077 i_mset1=iran_num(1,mset(i_iset1))
1078 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1079 econstr_temp_i=remd_ene(20,i)
1080 econstr_temp_iex=remd_ene(20,iex)
1081 remd_ene(20,i)=remd_ene(n_ene+3,i)
1082 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1084 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
1087 i_mult1=iran_num(1,remd_m(i_temp1))
1089 i_mset1=iran_num(1,mset(i_iset1))
1090 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1091 econstr_temp_i=remd_ene(20,i)
1092 econstr_temp_iex=remd_ene(20,iex)
1094 remd_ene(20,i)=remd_ene(n_ene+7,i)
1095 remd_ene(20,iex)=remd_ene(n_ene+8,iex)
1097 remd_ene(20,i)=remd_ene(n_ene+3,i)
1098 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
1104 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
1107 c Swap temperatures between conformations i and iex with recalculating the free energies
1108 c following temperature changes.
1110 write (iout,*) "i_dir",i_dir," i",i," iex",iex,
1111 & " t_bath",remd_t_bath(i),remd_t_bath(iex)
1113 ene_iex_iex=remd_ene(0,iex)
1114 ene_i_i=remd_ene(0,i)
1115 co write (iout,*) "rescaling weights with temperature",
1117 call rescale_weights(remd_t_bath(i))
1119 call sum_energy(remd_ene(0,iex),.false.)
1120 ene_iex_i=remd_ene(0,iex)
1121 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
1122 c call sum_energy(remd_ene(0,i),.false.)
1123 cd write (iout,*) "ene_i_i",remd_ene(0,i)
1124 c write (iout,*) "rescaling weights with temperature",
1125 c & remd_t_bath(iex)
1126 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
1127 c write (iout,*) "ERROR: inconsistent energies:",i,
1128 c & ene_i_i,remd_ene(0,i)
1130 call rescale_weights(remd_t_bath(iex))
1131 call sum_energy(remd_ene(0,i),.false.)
1132 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
1133 ene_i_iex=remd_ene(0,i)
1134 c call sum_energy(remd_ene(0,iex),.false.)
1135 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
1136 c write (iout,*) "ERROR: inconsistent energies:",iex,
1137 c & ene_iex_iex,remd_ene(0,iex)
1139 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
1140 c write (iout,*) "i",i," iex",iex
1141 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
1142 cd & " ene_i_iex",ene_i_iex,
1143 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
1144 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
1145 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
1147 cd write(iout,*) 'delta',delta
1148 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1149 c & (remd_ene(i)-remd_ene(iex))/Rb/
1150 c & (remd_t_bath(i)*remd_t_bath(iex))
1151 if (delta .gt. 50.0d0) then
1156 if (i_dir.eq.1.or.i_dir.eq.3)
1157 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1158 if (i_dir.eq.2.or.i_dir.eq.3)
1159 & iremd_tot_usa(int(i2set(i-1)))=
1160 & iremd_tot_usa(int(i2set(i-1)))+1
1161 xxx=ran_number(0.0d0,1.0d0)
1162 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1163 if (delta .gt. xxx) then
1165 remd_t_bath(i)=remd_t_bath(iex)
1166 remd_t_bath(iex)=tmp
1169 iremd_iset(i)=iremd_iset(iex)
1170 iremd_iset(iex)=itmp
1172 remd_ene(0,i)=ene_i_iex
1173 remd_ene(0,iex)=ene_iex_i
1175 if (i_dir.eq.1.or.i_dir.eq.3)
1176 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1179 i2rep(i-1)=i2rep(iex-1)
1182 if (i_dir.eq.2.or.i_dir.eq.3)
1183 & iremd_acc_usa(int(i2set(i-1)))=
1184 & iremd_acc_usa(int(i2set(i-1)))+1
1187 i2set(i-1)=i2set(iex-1)
1190 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1191 i_index(i_temp,i_mult,i_iset,i_mset)=
1192 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1193 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1196 remd_ene(0,iex)=ene_iex_iex
1197 remd_ene(0,i)=ene_i_i
1198 remd_ene(20,iex)=econstr_temp_iex
1199 remd_ene(20,i)=econstr_temp_i
1203 cd do il1=1,mset(il)
1206 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1219 c-------------------------------------
1220 write (iout,*) "NREP",nrep
1222 if(iremd_tot(i).ne.0)
1223 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1224 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1229 if(iremd_tot_usa(i).ne.0)
1230 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1231 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1237 cd write (iout,'(a6,100i4)') "ifirst",
1238 cd & (ifirst(i),i=1,remd_m(1))
1240 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1241 cd & (nupa(i,il),i=1,nupa(0,il))
1242 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1243 cd & (ndowna(i,il),i=1,ndowna(0,il))
1248 cd write (iout,*) "Before scatter"
1250 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1251 & t_bath,1,mpi_double_precision,king,
1253 cd write (iout,*) "After scatter"
1256 & call mpi_scatter(iremd_iset,1,mpi_integer,
1257 & iset,1,mpi_integer,king,
1261 if (me.eq.king .or. .not. out1file) then
1262 write(iout,*) 'REMD scatter time=',time07-time06
1266 call mpi_scatter(elowi,1,mpi_double_precision,
1267 & elow,1,mpi_double_precision,king,
1269 call mpi_scatter(ehighi,1,mpi_double_precision,
1270 & ehigh,1,mpi_double_precision,king,
1273 call rescale_weights(t_bath)
1274 co write (iout,*) "Processor",me,
1275 co & " rescaling weights with temperature",t_bath
1277 stdfp=dsqrt(2*Rb*t_bath/d_time)
1279 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1282 cde write(iout,*) 'REMD after',me,t_bath
1284 if (me.eq.king .or. .not. out1file) then
1285 write(iout,*) 'REMD exchange time=',time08-time00
1291 if (restart1file) then
1292 if (me.eq.king .or. .not. out1file)
1293 & write(iout,*) 'writing restart at the end of run'
1294 call write1rst(i_index)
1297 if (traj1file) call write1traj
1299 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1300 cdeb & icache_all,1,mpi_integer,king,
1301 cdeb & CG_COMM,ierr)
1302 cdeb write(iout,'(a40,8000i8)')
1303 cdeb & ' ntwx_cache after traj1file at the end',
1304 cdeb & (icache_all(i),i=1,nodes)
1309 t_MD=MPI_Wtime()-tt0
1313 if (me.eq.king .or. .not. out1file) then
1314 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1316 & 'MD calculations setup:',t_MDsetup,
1317 & 'Energy & gradient evaluation:',t_enegrad,
1318 & 'Stochastic MD setup:',t_langsetup,
1319 & 'Stochastic MD step setup:',t_sdsetup,
1321 write (iout,'(/28(1h=),a25,27(1h=))')
1322 & ' End of MD calculation '
1327 c-----------------------------------------------------------------------
1328 subroutine write1rst(i_index)
1330 include 'DIMENSIONS'
1332 include 'COMMON.CONTROL'
1335 include 'COMMON.LAGRANGE.5diag'
1337 include 'COMMON.LAGRANGE'
1339 include 'COMMON.QRESTR'
1340 include 'COMMON.IOUNITS'
1341 include 'COMMON.REMD'
1342 include 'COMMON.SETUP'
1343 include 'COMMON.CHAIN'
1344 include 'COMMON.SBRIDGE'
1345 include 'COMMON.INTERACT'
1347 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1348 & d_restart2(3,2*maxres*maxprocs)
1352 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1353 common /przechowalnia/ d_restart1,d_restart2
1354 integer i,j,il1,il,ixdrf
1360 t5_restart1(4)=t_bath
1361 t5_restart1(5)=Uconst
1363 call mpi_gather(t5_restart1,5,mpi_real,
1364 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1372 call mpi_gather(r_d,3*2*nres,mpi_real,
1373 & d_restart1,3*2*nres,mpi_real,king,
1382 call mpi_gather(r_d,3*2*nres,mpi_real,
1383 & d_restart2,3*2*nres,mpi_real,king,
1388 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1390 call xdrfint_(ixdrf, i2rep(i), iret)
1393 call xdrfint_(ixdrf, ifirst(i), iret)
1397 call xdrfint_(ixdrf, nupa(i,il), iret)
1401 call xdrfint_(ixdrf, ndowna(i,il), iret)
1407 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1414 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1421 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1427 call xdrfint_(ixdrf, nset, iret)
1429 call xdrfint_(ixdrf,mset(i), iret)
1432 call xdrfint_(ixdrf,i2set(i), iret)
1438 itmp=i_index(i,j,il,il1)
1439 call xdrfint_(ixdrf,itmp, iret)
1446 call xdrfclose_(ixdrf, iret)
1448 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1450 call xdrfint(ixdrf, i2rep(i), iret)
1453 call xdrfint(ixdrf, ifirst(i), iret)
1457 call xdrfint(ixdrf, nupa(i,il), iret)
1461 call xdrfint(ixdrf, ndowna(i,il), iret)
1467 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1474 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1481 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1488 call xdrfint(ixdrf, nset, iret)
1490 call xdrfint(ixdrf,mset(i), iret)
1493 call xdrfint(ixdrf,i2set(i), iret)
1499 itmp=i_index(i,j,il,il1)
1500 call xdrfint(ixdrf,itmp, iret)
1507 call xdrfclose(ixdrf, iret)
1514 subroutine write1traj
1516 include 'DIMENSIONS'
1519 include 'COMMON.QRESTR'
1520 include 'COMMON.IOUNITS'
1521 include 'COMMON.REMD'
1522 include 'COMMON.SETUP'
1523 include 'COMMON.CHAIN'
1524 include 'COMMON.SBRIDGE'
1525 include 'COMMON.INTERACT'
1529 real xcoord(3,maxres2+2),prec
1530 real r_qfrag(50),r_qpair(100)
1531 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1532 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1533 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1534 & p_uscdiff(100*maxprocs)
1535 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1536 common /przechowalnia/ p_c
1537 integer ii,i,il,j,ixdrf
1540 call mpi_bcast(ii_write,1,mpi_integer,
1541 & king,CG_COMM,ierr)
1544 c print *,'traj1file',me,ii_write,ntwx_cache
1548 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1550 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1553 t5_restart1(1)=totT_cache(ii)
1554 t5_restart1(2)=EK_cache(ii)
1555 t5_restart1(3)=potE_cache(ii)
1556 t5_restart1(4)=t_bath_cache(ii)
1557 t5_restart1(5)=Uconst_cache(ii)
1558 call mpi_gather(t5_restart1,5,mpi_real,
1559 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1561 call mpi_gather(iset_cache(ii),1,mpi_integer,
1562 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1565 r_qfrag(i)=qfrag_cache(i,ii)
1568 r_qpair(i)=qpair_cache(i,ii)
1571 r_utheta(i)=utheta_cache(i,ii)
1572 r_ugamma(i)=ugamma_cache(i,ii)
1573 r_uscdiff(i)=uscdiff_cache(i,ii)
1576 call mpi_gather(r_qfrag,nfrag,mpi_real,
1577 & p_qfrag,nfrag,mpi_real,king,
1579 call mpi_gather(r_qpair,npair,mpi_real,
1580 & p_qpair,npair,mpi_real,king,
1582 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1583 & p_utheta,nfrag_back,mpi_real,king,
1585 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1586 & p_ugamma,nfrag_back,mpi_real,king,
1588 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1589 & p_uscdiff,nfrag_back,mpi_real,king,
1593 write (iout,*) "p_qfrag"
1595 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1597 write (iout,*) "p_qpair"
1599 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1605 r_c(j,i)=c_cache(j,i,ii)
1609 call mpi_gather(r_c,3*2*nres,mpi_real,
1610 & p_c,3*2*nres,mpi_real,king,
1616 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1617 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1618 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1619 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1620 call xdrfint_(ixdrf, nss, iret)
1623 call xdrfint(ixdrf, idssb(j)+nres, iret)
1624 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1626 call xdrfint_(ixdrf, ihpb(j), iret)
1627 call xdrfint_(ixdrf, jhpb(j), iret)
1630 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1631 call xdrfint_(ixdrf, iset_restart1(il), iret)
1633 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1636 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1639 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1640 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1641 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1646 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1651 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1655 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1659 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1660 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1661 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1662 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1663 call xdrfint(ixdrf, nss, iret)
1666 call xdrfint(ixdrf, idssb(j)+nres, iret)
1667 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1669 call xdrfint(ixdrf, ihpb(j), iret)
1670 call xdrfint(ixdrf, jhpb(j), iret)
1673 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1674 call xdrfint(ixdrf, iset_restart1(il), iret)
1676 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1679 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1682 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1683 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1684 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1689 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1694 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1698 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1704 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1706 if(me.eq.king) call xdrfclose(ixdrf, iret)
1708 do i=1,ntwx_cache-ii_write
1710 totT_cache(i)=totT_cache(ii_write+i)
1711 EK_cache(i)=EK_cache(ii_write+i)
1712 potE_cache(i)=potE_cache(ii_write+i)
1713 t_bath_cache(i)=t_bath_cache(ii_write+i)
1714 Uconst_cache(i)=Uconst_cache(ii_write+i)
1715 iset_cache(i)=iset_cache(ii_write+i)
1718 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1721 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1724 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1725 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1726 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1731 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1735 ntwx_cache=ntwx_cache-ii_write
1740 subroutine read1restart(i_index)
1742 include 'DIMENSIONS'
1744 include 'COMMON.CONTROL'
1747 include 'COMMON.LAGRANGE.5diag'
1749 include 'COMMON.LAGRANGE'
1751 include 'COMMON.QRESTR'
1752 include 'COMMON.IOUNITS'
1753 include 'COMMON.REMD'
1754 include 'COMMON.SETUP'
1755 include 'COMMON.CHAIN'
1756 include 'COMMON.SBRIDGE'
1757 include 'COMMON.INTERACT'
1758 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1761 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1762 common /przechowalnia/ d_restart1
1763 integer i,j,il,il1,ixdrf,iret,itmp
1765 write (*,*) "Processor",me," called read1restart"
1768 open(irest2,file=mremd_rst_name,status='unknown')
1769 read(irest2,*,err=334) i
1770 write(iout,*) "Reading old rst in ASCI format"
1772 call read1restart_old
1776 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1779 call xdrfint_(ixdrf, i2rep(i), iret)
1782 call xdrfint_(ixdrf, ifirst(i), iret)
1785 call xdrfint_(ixdrf, nupa(0,il), iret)
1787 call xdrfint_(ixdrf, nupa(i,il), iret)
1790 call xdrfint_(ixdrf, ndowna(0,il), iret)
1792 call xdrfint_(ixdrf, ndowna(i,il), iret)
1797 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1801 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1804 call xdrfint(ixdrf, i2rep(i), iret)
1807 call xdrfint(ixdrf, ifirst(i), iret)
1810 call xdrfint(ixdrf, nupa(0,il), iret)
1812 call xdrfint(ixdrf, nupa(i,il), iret)
1815 call xdrfint(ixdrf, ndowna(0,il), iret)
1817 call xdrfint(ixdrf, ndowna(i,il), iret)
1822 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1827 call mpi_scatter(t_restart1,5,mpi_real,
1828 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1832 t_bath=t5_restart1(4)
1837 c read(irest2,'(3e15.5)')
1838 c & (d_restart1(j,i+2*nres*il),j=1,3)
1841 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1843 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1849 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1850 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1860 c read(irest2,'(3e15.5)')
1861 c & (d_restart1(j,i+2*nres*il),j=1,3)
1864 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1866 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1872 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1873 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1884 call xdrfint_(ixdrf, nset, iret)
1886 call xdrfint_(ixdrf,mset(i), iret)
1889 call xdrfint_(ixdrf,i2set(i), iret)
1895 call xdrfint_(ixdrf,itmp, iret)
1896 i_index(i,j,il,il1)=itmp
1904 call xdrfint(ixdrf, nset, iret)
1906 call xdrfint(ixdrf,mset(i), iret)
1909 call xdrfint(ixdrf,i2set(i), iret)
1915 call xdrfint(ixdrf,itmp, iret)
1916 i_index(i,j,il,il1)=itmp
1923 Corrected AL 8/19/2014: each processor needs whole iset array not only its
1925 c call mpi_scatter(i2set,1,mpi_integer,
1926 c & iset,1,mpi_integer,king,
1928 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1935 if(me.eq.king) close(irest2)
1939 subroutine read1restart_old
1941 include 'DIMENSIONS'
1945 include 'COMMON.LAGRANGE.5diag'
1947 include 'COMMON.LAGRANGE'
1949 include 'COMMON.IOUNITS'
1950 include 'COMMON.REMD'
1951 include 'COMMON.SETUP'
1952 include 'COMMON.CHAIN'
1953 include 'COMMON.SBRIDGE'
1954 include 'COMMON.INTERACT'
1955 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1957 common /przechowalnia/ d_restart1
1961 open(irest2,file=mremd_rst_name,status='unknown')
1962 read (irest2,*) (i2rep(i),i=0,nodes-1)
1963 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1965 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1966 read (irest2,*) ndowna(0,il),
1967 & (ndowna(i,il),i=1,ndowna(0,il))
1970 read(irest2,*) (t_restart1(j,il),j=1,4)
1973 call mpi_scatter(t_restart1,5,mpi_real,
1974 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1978 t_bath=t5_restart1(4)
1983 read(irest2,'(3e15.5)')
1984 & (d_restart1(j,i+2*nres*il),j=1,3)
1988 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1989 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1999 read(irest2,'(3e15.5)')
2000 & (d_restart1(j,i+2*nres*il),j=1,3)
2004 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2005 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2011 if(me.eq.king) close(irest2)
2014 c------------------------------------------