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 c Compute the standard deviations of stochastic forces for Langevin dynamics
1283 c if the friction coefficients do not depend on surface area
1284 if (lang.gt.0 .and. .not.surfarea) then
1286 stdforcp(i)=stdfp*dsqrt(gamp)
1289 if (itype(i).ne.ntyp1) stdforcsc(i)=stdfsc(iabs(itype(i)))
1290 & *dsqrt(gamsc(iabs(itype(i))))
1295 cde write(iout,*) 'REMD after',me,t_bath
1297 if (me.eq.king .or. .not. out1file) then
1298 write(iout,*) 'REMD exchange time=',time08-time00
1304 if (restart1file) then
1305 if (me.eq.king .or. .not. out1file)
1306 & write(iout,*) 'writing restart at the end of run'
1307 call write1rst(i_index)
1310 if (traj1file) call write1traj
1312 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1313 cdeb & icache_all,1,mpi_integer,king,
1314 cdeb & CG_COMM,ierr)
1315 cdeb write(iout,'(a40,8000i8)')
1316 cdeb & ' ntwx_cache after traj1file at the end',
1317 cdeb & (icache_all(i),i=1,nodes)
1322 t_MD=MPI_Wtime()-tt0
1326 if (me.eq.king .or. .not. out1file) then
1327 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1329 & 'MD calculations setup:',t_MDsetup,
1330 & 'Energy & gradient evaluation:',t_enegrad,
1331 & 'Stochastic MD setup:',t_langsetup,
1332 & 'Stochastic MD step setup:',t_sdsetup,
1334 write (iout,'(/28(1h=),a25,27(1h=))')
1335 & ' End of MD calculation '
1340 c-----------------------------------------------------------------------
1341 subroutine write1rst(i_index)
1343 include 'DIMENSIONS'
1345 include 'COMMON.CONTROL'
1348 include 'COMMON.LAGRANGE.5diag'
1350 include 'COMMON.LAGRANGE'
1352 include 'COMMON.QRESTR'
1353 include 'COMMON.IOUNITS'
1354 include 'COMMON.REMD'
1355 include 'COMMON.SETUP'
1356 include 'COMMON.CHAIN'
1357 include 'COMMON.SBRIDGE'
1358 include 'COMMON.INTERACT'
1360 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1361 & d_restart2(3,2*maxres*maxprocs)
1365 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1366 common /przechowalnia/ d_restart1,d_restart2
1367 integer i,j,il1,il,ixdrf
1373 t5_restart1(4)=t_bath
1374 t5_restart1(5)=Uconst
1376 call mpi_gather(t5_restart1,5,mpi_real,
1377 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1385 call mpi_gather(r_d,3*2*nres,mpi_real,
1386 & d_restart1,3*2*nres,mpi_real,king,
1395 call mpi_gather(r_d,3*2*nres,mpi_real,
1396 & d_restart2,3*2*nres,mpi_real,king,
1401 call xdrfopen_(ixdrf,mremd_rst_name, "w", iret)
1403 call xdrfint_(ixdrf, i2rep(i), iret)
1406 call xdrfint_(ixdrf, ifirst(i), iret)
1410 call xdrfint_(ixdrf, nupa(i,il), iret)
1414 call xdrfint_(ixdrf, ndowna(i,il), iret)
1420 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
1427 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1434 call xdrffloat_(ixdrf, d_restart2(j,i+2*nres*il), iret)
1440 call xdrfint_(ixdrf, nset, iret)
1442 call xdrfint_(ixdrf,mset(i), iret)
1445 call xdrfint_(ixdrf,i2set(i), iret)
1451 itmp=i_index(i,j,il,il1)
1452 call xdrfint_(ixdrf,itmp, iret)
1459 call xdrfclose_(ixdrf, iret)
1461 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1463 call xdrfint(ixdrf, i2rep(i), iret)
1466 call xdrfint(ixdrf, ifirst(i), iret)
1470 call xdrfint(ixdrf, nupa(i,il), iret)
1474 call xdrfint(ixdrf, ndowna(i,il), iret)
1480 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1487 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1494 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1501 call xdrfint(ixdrf, nset, iret)
1503 call xdrfint(ixdrf,mset(i), iret)
1506 call xdrfint(ixdrf,i2set(i), iret)
1512 itmp=i_index(i,j,il,il1)
1513 call xdrfint(ixdrf,itmp, iret)
1520 call xdrfclose(ixdrf, iret)
1527 subroutine write1traj
1529 include 'DIMENSIONS'
1532 include 'COMMON.QRESTR'
1533 include 'COMMON.IOUNITS'
1534 include 'COMMON.REMD'
1535 include 'COMMON.SETUP'
1536 include 'COMMON.CHAIN'
1537 include 'COMMON.SBRIDGE'
1538 include 'COMMON.INTERACT'
1542 real xcoord(3,maxres2+2),prec
1543 real r_qfrag(50),r_qpair(100)
1544 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1545 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1546 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1547 & p_uscdiff(100*maxprocs)
1548 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1549 common /przechowalnia/ p_c
1550 integer ii,i,il,j,ixdrf
1553 call mpi_bcast(ii_write,1,mpi_integer,
1554 & king,CG_COMM,ierr)
1557 c print *,'traj1file',me,ii_write,ntwx_cache
1561 if(me.eq.king) call xdrfopen_(ixdrf,cartname, "a", iret)
1563 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1566 t5_restart1(1)=totT_cache(ii)
1567 t5_restart1(2)=EK_cache(ii)
1568 t5_restart1(3)=potE_cache(ii)
1569 t5_restart1(4)=t_bath_cache(ii)
1570 t5_restart1(5)=Uconst_cache(ii)
1571 call mpi_gather(t5_restart1,5,mpi_real,
1572 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1574 call mpi_gather(iset_cache(ii),1,mpi_integer,
1575 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1578 r_qfrag(i)=qfrag_cache(i,ii)
1581 r_qpair(i)=qpair_cache(i,ii)
1584 r_utheta(i)=utheta_cache(i,ii)
1585 r_ugamma(i)=ugamma_cache(i,ii)
1586 r_uscdiff(i)=uscdiff_cache(i,ii)
1589 call mpi_gather(r_qfrag,nfrag,mpi_real,
1590 & p_qfrag,nfrag,mpi_real,king,
1592 call mpi_gather(r_qpair,npair,mpi_real,
1593 & p_qpair,npair,mpi_real,king,
1595 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1596 & p_utheta,nfrag_back,mpi_real,king,
1598 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1599 & p_ugamma,nfrag_back,mpi_real,king,
1601 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1602 & p_uscdiff,nfrag_back,mpi_real,king,
1606 write (iout,*) "p_qfrag"
1608 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1610 write (iout,*) "p_qpair"
1612 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1618 r_c(j,i)=c_cache(j,i,ii)
1622 call mpi_gather(r_c,3*2*nres,mpi_real,
1623 & p_c,3*2*nres,mpi_real,king,
1629 call xdrffloat_(ixdrf, real(t_restart1(1,il)), iret)
1630 call xdrffloat_(ixdrf, real(t_restart1(3,il)), iret)
1631 call xdrffloat_(ixdrf, real(t_restart1(5,il)), iret)
1632 call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret)
1633 call xdrfint_(ixdrf, nss, iret)
1636 call xdrfint(ixdrf, idssb(j)+nres, iret)
1637 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1639 call xdrfint_(ixdrf, ihpb(j), iret)
1640 call xdrfint_(ixdrf, jhpb(j), iret)
1643 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
1644 call xdrfint_(ixdrf, iset_restart1(il), iret)
1646 call xdrffloat_(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1649 call xdrffloat_(ixdrf, p_qpair(i+(il-1)*npair), iret)
1652 call xdrffloat_(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1653 call xdrffloat_(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1654 call xdrffloat_(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1659 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1664 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1668 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
1672 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1673 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1674 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1675 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1676 call xdrfint(ixdrf, nss, iret)
1679 call xdrfint(ixdrf, idssb(j)+nres, iret)
1680 call xdrfint(ixdrf, jdssb(j)+nres, iret)
1682 call xdrfint(ixdrf, ihpb(j), iret)
1683 call xdrfint(ixdrf, jhpb(j), iret)
1686 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1687 call xdrfint(ixdrf, iset_restart1(il), iret)
1689 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1692 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1695 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1696 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1697 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1702 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1707 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1711 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1717 if(me.eq.king) call xdrfclose_(ixdrf, iret)
1719 if(me.eq.king) call xdrfclose(ixdrf, iret)
1721 do i=1,ntwx_cache-ii_write
1723 totT_cache(i)=totT_cache(ii_write+i)
1724 EK_cache(i)=EK_cache(ii_write+i)
1725 potE_cache(i)=potE_cache(ii_write+i)
1726 t_bath_cache(i)=t_bath_cache(ii_write+i)
1727 Uconst_cache(i)=Uconst_cache(ii_write+i)
1728 iset_cache(i)=iset_cache(ii_write+i)
1731 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1734 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1737 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1738 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1739 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1744 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1748 ntwx_cache=ntwx_cache-ii_write
1753 subroutine read1restart(i_index)
1755 include 'DIMENSIONS'
1757 include 'COMMON.CONTROL'
1760 include 'COMMON.LAGRANGE.5diag'
1762 include 'COMMON.LAGRANGE'
1764 include 'COMMON.QRESTR'
1765 include 'COMMON.IOUNITS'
1766 include 'COMMON.REMD'
1767 include 'COMMON.SETUP'
1768 include 'COMMON.CHAIN'
1769 include 'COMMON.SBRIDGE'
1770 include 'COMMON.INTERACT'
1771 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1774 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1775 common /przechowalnia/ d_restart1
1776 integer i,j,il,il1,ixdrf,iret,itmp
1778 write (*,*) "Processor",me," called read1restart"
1781 open(irest2,file=mremd_rst_name,status='unknown')
1782 read(irest2,*,err=334) i
1783 write(iout,*) "Reading old rst in ASCI format"
1785 call read1restart_old
1789 call xdrfopen_(ixdrf,mremd_rst_name, "r", iret)
1792 call xdrfint_(ixdrf, i2rep(i), iret)
1795 call xdrfint_(ixdrf, ifirst(i), iret)
1798 call xdrfint_(ixdrf, nupa(0,il), iret)
1800 call xdrfint_(ixdrf, nupa(i,il), iret)
1803 call xdrfint_(ixdrf, ndowna(0,il), iret)
1805 call xdrfint_(ixdrf, ndowna(i,il), iret)
1810 call xdrffloat_(ixdrf, t_restart1(j,il), iret)
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)
1840 call mpi_scatter(t_restart1,5,mpi_real,
1841 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1845 t_bath=t5_restart1(4)
1850 c read(irest2,'(3e15.5)')
1851 c & (d_restart1(j,i+2*nres*il),j=1,3)
1854 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1856 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1862 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1863 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1873 c read(irest2,'(3e15.5)')
1874 c & (d_restart1(j,i+2*nres*il),j=1,3)
1877 call xdrffloat_(ixdrf, d_restart1(j,i+2*nres*il), iret)
1879 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1885 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1886 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1897 call xdrfint_(ixdrf, nset, iret)
1899 call xdrfint_(ixdrf,mset(i), iret)
1902 call xdrfint_(ixdrf,i2set(i), iret)
1908 call xdrfint_(ixdrf,itmp, iret)
1909 i_index(i,j,il,il1)=itmp
1917 call xdrfint(ixdrf, nset, iret)
1919 call xdrfint(ixdrf,mset(i), iret)
1922 call xdrfint(ixdrf,i2set(i), iret)
1928 call xdrfint(ixdrf,itmp, iret)
1929 i_index(i,j,il,il1)=itmp
1936 Corrected AL 8/19/2014: each processor needs whole iset array not only its
1938 c call mpi_scatter(i2set,1,mpi_integer,
1939 c & iset,1,mpi_integer,king,
1941 call mpi_bcast(i2set(0),nodes,mpi_integer,king,
1948 if(me.eq.king) close(irest2)
1952 subroutine read1restart_old
1954 include 'DIMENSIONS'
1958 include 'COMMON.LAGRANGE.5diag'
1960 include 'COMMON.LAGRANGE'
1962 include 'COMMON.IOUNITS'
1963 include 'COMMON.REMD'
1964 include 'COMMON.SETUP'
1965 include 'COMMON.CHAIN'
1966 include 'COMMON.SBRIDGE'
1967 include 'COMMON.INTERACT'
1968 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1970 common /przechowalnia/ d_restart1
1974 open(irest2,file=mremd_rst_name,status='unknown')
1975 read (irest2,*) (i2rep(i),i=0,nodes-1)
1976 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1978 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1979 read (irest2,*) ndowna(0,il),
1980 & (ndowna(i,il),i=1,ndowna(0,il))
1983 read(irest2,*) (t_restart1(j,il),j=1,4)
1986 call mpi_scatter(t_restart1,5,mpi_real,
1987 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1991 t_bath=t5_restart1(4)
1996 read(irest2,'(3e15.5)')
1997 & (d_restart1(j,i+2*nres*il),j=1,3)
2001 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2002 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2012 read(irest2,'(3e15.5)')
2013 & (d_restart1(j,i+2*nres*il),j=1,3)
2017 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
2018 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
2024 if(me.eq.king) close(irest2)
2027 c------------------------------------------