2 implicit real*8 (a-h,o-z)
5 include 'COMMON.CONTROL'
9 include 'COMMON.LANGEVIN'
11 include 'COMMON.LANGEVIN.lang0'
13 include 'COMMON.CHAIN'
14 include 'COMMON.DERIV'
16 include 'COMMON.LOCAL'
17 include 'COMMON.INTERACT'
18 include 'COMMON.IOUNITS'
19 include 'COMMON.NAMES'
20 include 'COMMON.TIME1'
22 include 'COMMON.SETUP'
24 double precision cm(3),L(3),vcm(3)
25 double precision energia(0:n_ene)
26 double precision remd_t_bath(maxprocs)
27 integer iremd_iset(maxprocs)
29 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
30 double precision remd_ene(0:n_ene+4,maxprocs)
31 integer iremd_acc(maxprocs),iremd_tot(maxprocs)
32 integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
38 cold integer nup(0:maxprocs),ndown(0:maxprocs)
39 integer rep2i(0:maxprocs),ireqi(maxprocs)
40 integer icache_all(maxprocs)
41 integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
42 logical synflag,end_of_run,file_exist /.false./
48 if(me.eq.king.or..not.out1file) then
49 write (iout,*) 'MREMD',nodes,'time before',time00-walltime
50 write (iout,*) "NREP=",nrep
54 if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
55 call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
57 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
59 cd print *,'MREMD',nodes
60 cd print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
61 cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath
65 do il1=1,max0(mset(il),1)
81 if(me.eq.king.or..not.out1file) then
82 write(iout,*) (i2rep(i),i=0,nodes-1)
83 write(iout,*) (i2set(i),i=0,nodes-1)
88 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
95 c print *,'i2rep',me,i2rep(me)
96 c print *,'rep2i',(rep2i(i),i=0,nrep)
98 cold if (i2rep(me).eq.nrep) then
101 cold nup(0)=remd_m(i2rep(me)+1)
102 cold k=rep2i(int(i2rep(me)))+1
109 cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
111 cold if (i2rep(me).eq.1) then
114 cold ndown(0)=remd_m(i2rep(me)-1)
115 cold k=rep2i(i2rep(me)-2)+1
122 cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
125 write (*,*) "Processor",me," rest",rest,"
126 & restart1fie",restart1file
127 if(rest.and.restart1file) then
129 & inquire(file=mremd_rst_name,exist=file_exist)
130 cd write (*,*) me," Before broadcast: file_exist",file_exist
131 call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
133 cd write (*,*) me," After broadcast: file_exist",file_exist
135 if(me.eq.king.or..not.out1file)
136 & write (iout,*) 'Master is reading restart1file'
137 call read1restart(i_index)
139 if(me.eq.king.or..not.out1file)
140 & write (iout,*) 'WARNING : no restart1file'
143 if(me.eq.king.or..not.out1file) then
144 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
145 write(iout,*) "i_index"
150 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
159 if (rest.and..not.restart1file)
160 & inquire(file=mremd_rst_name,exist=file_exist)
161 if(.not.file_exist.and.rest.and..not.restart1file)
162 & write(iout,*) 'WARNING : no restart file',mremd_rst_name
163 IF (rest.and.file_exist.and..not.restart1file) THEN
164 write (iout,*) 'Master is reading restart file',
166 open(irest2,file=mremd_rst_name,status='unknown')
168 read (irest2,*) (i2rep(i),i=0,nodes-1)
170 read (irest2,*) (ifirst(i),i=1,remd_m(1))
173 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
175 read (irest2,*) ndowna(0,il),
176 & (ndowna(i,il),i=1,ndowna(0,il))
182 read (irest2,*) (mset(i),i=1,nset)
184 read (irest2,*) (i2set(i),i=0,nodes-1)
189 read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
194 write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
195 write(iout,*) "i_index"
200 write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
209 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
210 write (iout,'(a6,1000i5)') "ifirst",
211 & (ifirst(i),i=1,remd_m(1))
213 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
214 & (nupa(i,il),i=1,nupa(0,il))
215 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
216 & (ndowna(i,il),i=1,ndowna(0,il))
218 ELSE IF (.not.(rest.and.file_exist)) THEN
224 if (i2rep(il-1).eq.nrep) then
227 nupa(0,il)=remd_m(i2rep(il-1)+1)
228 k=rep2i(int(i2rep(il-1)))+1
234 if (i2rep(il-1).eq.1) then
237 ndowna(0,il)=remd_m(i2rep(il-1)-1)
238 k=rep2i(i2rep(il-1)-2)+1
246 c write (iout,'(a6,100i4)') "ifirst",
247 c & (ifirst(i),i=1,remd_m(1))
249 c write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
250 c & (nupa(i,il),i=1,nupa(0,il))
251 c write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
252 c & (ndowna(i,il),i=1,ndowna(0,il))
258 c t_bath=retmin+(retmax-retmin)*me/(nodes-1)
259 if(.not.(rest.and.file_exist.and.restart1file)) then
260 if (me .eq. king) then
263 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
265 cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
266 if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
271 if(me.eq.king.or..not.out1file)
272 & write(iout,*) me,"iset=",iset,"t_bath=",t_bath
275 stdfp=dsqrt(2*Rb*t_bath/d_time)
277 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
280 c print *,'irep',me,t_bath
282 if (me.eq.king .or. .not. out1file)
283 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
284 call rescale_weights(t_bath)
288 c------copy MD--------------
289 c The driver for molecular dynamics subroutines
290 c------------------------------------------------
296 if(me.eq.king.or..not.out1file)
297 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
303 c Determine the inverse of the inertia matrix.
304 call setup_MD_matrices
308 if (me.eq.king .or. .not. out1file)
309 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
310 stdfp=dsqrt(2*Rb*t_bath/d_time)
312 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
314 call rescale_weights(t_bath)
318 t_MDsetup = MPI_Wtime()
320 t_MDsetup = tcpu()-tt0
323 c Entering the MD loop
329 if (lang.eq.2 .or. lang.eq.3) then
332 call sd_verlet_p_setup
334 call sd_verlet_ciccotti_setup
339 pfric0_mat(i,j,0)=pfric_mat(i,j)
340 afric0_mat(i,j,0)=afric_mat(i,j)
341 vfric0_mat(i,j,0)=vfric_mat(i,j)
342 prand0_mat(i,j,0)=prand_mat(i,j)
343 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
344 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
350 flag_stoch(i)=.false.
352 else if (lang.eq.1 .or. lang.eq.4) then
356 if (me.eq.king .or. .not. out1file)
357 & write(iout,*) 'Setup time',time00-walltime
360 t_langsetup=MPI_Wtime()
363 t_langsetup=tcpu()-tt0
368 do while(.not.end_of_run)
370 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
371 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
373 if (lang.gt.0 .and. surfarea .and.
374 & mod(itime,reset_fricmat).eq.0) then
375 if (lang.eq.2 .or. lang.eq.3) then
378 call sd_verlet_p_setup
380 call sd_verlet_ciccotti_setup
385 pfric0_mat(i,j,0)=pfric_mat(i,j)
386 afric0_mat(i,j,0)=afric_mat(i,j)
387 vfric0_mat(i,j,0)=vfric_mat(i,j)
388 prand0_mat(i,j,0)=prand_mat(i,j)
389 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
390 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
396 flag_stoch(i)=.false.
398 else if (lang.eq.1 .or. lang.eq.4) then
401 write (iout,'(a,i10)')
402 & "Friction matrix reset based on surface area, itime",itime
404 if (reset_vel .and. tbf .and. lang.eq.0
405 & .and. mod(itime,count_reset_vel).eq.0) then
407 if (me.eq.king .or. .not. out1file)
408 & write(iout,'(a,f20.2)')
409 & "Velocities reset to random values, time",totT
412 d_t_old(j,i)=d_t(j,i)
416 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
420 d_t(j,0)=d_t(j,0)-vcm(j)
423 kinetic_T=2.0d0/(dimen*Rb)*EK
424 scalfac=dsqrt(T_bath/kinetic_T)
425 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
428 d_t_old(j,i)=scalfac*d_t(j,i)
434 c Time-reversible RESPA algorithm
435 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
436 call RESPA_step(itime)
438 c Variable time step algorithm.
439 call velverlet_step(itime)
442 call brown_step(itime)
445 if (mod(itime,ntwe).eq.0) call statout(itime)
447 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
448 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
450 call pdbout(potE,tytul,ipdb)
455 if (mod(itime,ntwx).eq.0.and.traj1file) then
456 if(ntwx_cache.lt.max_cache_traj_use) then
457 ntwx_cache=ntwx_cache+1
459 if (max_cache_traj_use.ne.1)
460 & print *,itime,"processor ",me," over cache ",ntwx_cache
463 totT_cache(i)=totT_cache(i+1)
464 EK_cache(i)=EK_cache(i+1)
465 potE_cache(i)=potE_cache(i+1)
466 t_bath_cache(i)=t_bath_cache(i+1)
467 Uconst_cache(i)=Uconst_cache(i+1)
468 iset_cache(i)=iset_cache(i+1)
471 qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
474 qpair_cache(ii,i)=qpair_cache(ii,i+1)
477 utheta_cache(ii,i)=utheta_cache(ii,i+1)
478 ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
479 uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
485 c_cache(j,ii,i)=c_cache(j,ii,i+1)
491 totT_cache(ntwx_cache)=totT
492 EK_cache(ntwx_cache)=EK
493 potE_cache(ntwx_cache)=potE
494 t_bath_cache(ntwx_cache)=t_bath
495 Uconst_cache(ntwx_cache)=Uconst
496 iset_cache(ntwx_cache)=iset
499 qfrag_cache(i,ntwx_cache)=qfrag(i)
502 qpair_cache(i,ntwx_cache)=qpair(i)
505 utheta_cache(i,ntwx_cache)=utheta(i)
506 ugamma_cache(i,ntwx_cache)=ugamma(i)
507 uscdiff_cache(i,ntwx_cache)=uscdiff(i)
512 c_cache(j,i,ntwx_cache)=c(j,i)
517 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
518 & .and..not.restart1file) then
521 open(irest1,file=mremd_rst_name,status='unknown')
522 write (irest1,*) "i2rep"
523 write (irest1,*) (i2rep(i),i=0,nodes-1)
524 write (irest1,*) "ifirst"
525 write (irest1,*) (ifirst(i),i=1,remd_m(1))
527 write (irest1,*) "nupa",il
528 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
529 write (irest1,*) "ndowna",il
530 write (irest1,*) ndowna(0,il),
531 & (ndowna(i,il),i=1,ndowna(0,il))
534 write (irest1,*) "nset"
535 write (irest1,*) nset
536 write (irest1,*) "mset"
537 write (irest1,*) (mset(i),i=1,nset)
538 write (irest1,*) "i2set"
539 write (irest1,*) (i2set(i),i=0,nodes-1)
540 write (irest1,*) "i_index"
544 write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
552 open(irest2,file=rest2name,status='unknown')
553 write(irest2,*) totT,EK,potE,totE,t_bath
555 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
558 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
561 write (irest2,*) iset
568 c forced synchronization
569 if (mod(itime,i_sync_step).eq.0 .and. me.ne.king
570 & .and. .not. mremdsync) then
572 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
574 call mpi_recv(itime_master, 1, MPI_INTEGER,
575 & 0,101,CG_COMM, status, ierr)
576 call mpi_barrier(CG_COMM, ierr)
577 cdeb if (out1file.or.traj1file) then
578 cdeb call mpi_gather(itime,1,mpi_integer,
579 cdeb & icache_all,1,mpi_integer,king,
582 & call mpi_gather(ntwx_cache,1,mpi_integer,
583 & icache_all,1,mpi_integer,king,
586 & write(iout,*) 'REMD synchro at',itime_master,itime
587 if(itime_master.ge.n_timestep) end_of_run=.true.
588 ctime call flush(iout)
593 if ((mod(itime,nstex).eq.0.and.me.eq.king
594 & .or.end_of_run.and.me.eq.king )
595 & .and. .not. mremdsync ) then
599 call mpi_isend(itime,1,MPI_INTEGER,i,101,
600 & CG_COMM, ireqi(i), ierr)
601 cd write(iout,*) 'REMD synchro with',i
604 call mpi_waitall(nodes-1,ireqi,statusi,ierr)
605 call mpi_barrier(CG_COMM, ierr)
607 write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
608 if (out1file.or.traj1file) then
609 cdeb call mpi_gather(itime,1,mpi_integer,
610 cdeb & itime_all,1,mpi_integer,king,
612 cdeb write(iout,'(a19,8000i8)') ' REMD synchro itime',
613 cdeb & (itime_all(i),i=1,nodes)
615 cdeb imin_itime=itime_all(1)
617 cdeb if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
619 cdeb ii_write=(imin_itime-imin_itime_old)/ntwx
620 cdeb imin_itime_old=int(imin_itime/ntwx)*ntwx
621 cdeb write(iout,*) imin_itime,imin_itime_old,ii_write
622 call mpi_gather(ntwx_cache,1,mpi_integer,
623 & icache_all,1,mpi_integer,king,
625 c write(iout,'(a19,8000i8)') ' ntwx_cache',
626 c & (icache_all(i),i=1,nodes)
627 ii_write=icache_all(1)
629 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
631 c write(iout,*) "MIN ii_write=",ii_write
634 ctime call flush(iout)
636 if(mremdsync .and. mod(itime,nstex).eq.0) then
638 if (me.eq.king .or. .not. out1file)
639 & write(iout,*) 'REMD synchro at',itime
642 call mpi_gather(ntwx_cache,1,mpi_integer,
643 & icache_all,1,mpi_integer,king,
646 write(iout,'(a19,8000i8)') ' ntwx_cache',
647 & (icache_all(i),i=1,nodes)
648 ii_write=icache_all(1)
650 if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
652 write(iout,*) "MIN ii_write=",ii_write
658 if(synflag.and..not.end_of_run) then
662 cd write(iout,*) 'REMD before',me,t_bath
664 c call mpi_gather(t_bath,1,mpi_double_precision,
665 c & remd_t_bath,1,mpi_double_precision,king,
667 potEcomp(n_ene+1)=t_bath
669 potEcomp(n_ene+2)=iset
670 if (iset.lt.nset) then
674 potEcomp(n_ene+3)=Uconst
681 potEcomp(n_ene+4)=Uconst
685 call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
686 & remd_ene(0,1),n_ene+5,mpi_double_precision,king,
689 call mpi_gather(elow,1,mpi_double_precision,
690 & elowi,1,mpi_double_precision,king,
692 call mpi_gather(ehigh,1,mpi_double_precision,
693 & ehighi,1,mpi_double_precision,king,
698 if (me.eq.king .or. .not. out1file) then
699 write(iout,*) 'REMD gather times=',time03-time01
703 if (restart1file) call write1rst(i_index)
706 if (me.eq.king .or. .not. out1file) then
707 write(iout,*) 'REMD writing rst time=',time04-time03
710 if (traj1file) call write1traj
712 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
713 cdeb & icache_all,1,mpi_integer,king,
715 cdeb write(iout,'(a19,8000i8)') ' ntwx_cache after traj1file',
716 cdeb & (icache_all(i),i=1,nodes)
721 if (me.eq.king .or. .not. out1file) then
722 write(iout,*) 'REMD writing traj time=',time05-time04
729 remd_t_bath(i)=remd_ene(n_ene+1,i)
730 iremd_iset(i)=remd_ene(n_ene+2,i)
733 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
735 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
739 cd write(iout,*) 'REMD exchange temp,ene'
741 co write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
742 c write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
745 c-------------------------------------
747 write (iout,*) "Enter exchnge, remd_m",remd_m(1),
751 i=ifirst(iran_num(1,remd_m(1)))
756 if(i.gt.0.and.nupa(0,i).gt.0) then
757 iex=nupa(iran_num(1,int(nupa(0,i))),i)
759 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
761 c Swap temperatures between conformations i and iex with recalculating the free energies
762 c following temperature changes.
763 ene_iex_iex=remd_ene(0,iex)
764 ene_i_i=remd_ene(0,i)
765 cd write (iout,*) "rescaling weights with temperature",
768 call rescale_weights(remd_t_bath(i))
770 cd write (iout,*) "0,iex",remd_t_bath(i)
771 cd call enerprint(remd_ene(0,iex))
773 call sum_energy(remd_ene(0,iex),.false.)
774 ene_iex_i=remd_ene(0,iex)
775 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
777 cd write (iout,*) "0,i",remd_t_bath(i)
778 cd call enerprint(remd_ene(0,i))
780 call sum_energy(remd_ene(0,i),.false.)
781 cd write (iout,*) "ene_i_i",remd_ene(0,i)
783 cd write (iout,*) "rescaling weights with temperature",
784 cd & remd_t_bath(iex)
785 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
786 write (iout,*) "ERROR: inconsistent energies:",i,
787 & ene_i_i,remd_ene(0,i)
789 call rescale_weights(remd_t_bath(iex))
791 cd write (iout,*) "0,i",remd_t_bath(iex)
792 cd call enerprint(remd_ene(0,i))
794 call sum_energy(remd_ene(0,i),.false.)
795 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
797 ene_i_iex=remd_ene(0,i)
799 cd write (iout,*) "0,iex",remd_t_bath(iex)
800 cd call enerprint(remd_ene(0,iex))
802 call sum_energy(remd_ene(0,iex),.false.)
803 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
804 write (iout,*) "ERROR: inconsistent energies:",iex,
805 & ene_iex_iex,remd_ene(0,iex)
807 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
808 cd write (iout,*) "i",i," iex",iex
809 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
810 cd & " ene_i_iex",ene_i_iex,
811 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
813 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
814 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
816 cd write(iout,*) 'delta',delta
817 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
818 c & (remd_ene(i)-remd_ene(iex))/Rb/
819 c & (remd_t_bath(i)*remd_t_bath(iex))
821 if (delta .gt. 50.0d0) then
827 else if (delta.lt.-50.0d0) then
836 iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
837 xxx=ran_number(0.0d0,1.0d0)
838 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
840 if (delta .gt. xxx) then
842 remd_t_bath(i)=remd_t_bath(iex)
844 remd_ene(0,i)=ene_i_iex
845 remd_ene(0,iex)=ene_iex_i
851 ehighi(i)=ehighi(iex)
858 nupa(k,i)=nupa(k,iex)
861 ndowna(k,i)=ndowna(k,iex)
865 if (ifirst(il).eq.i) ifirst(il)=iex
867 if (nupa(k,il).eq.i) then
869 elseif (nupa(k,il).eq.iex) then
874 if (ndowna(k,il).eq.i) then
876 elseif (ndowna(k,il).eq.iex) then
882 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
884 i2rep(i-1)=i2rep(iex-1)
887 cd write(iout,*) 'exchange',i,iex
888 cd write (iout,'(a8,100i4)') "@ ifirst",
889 cd & (ifirst(k),k=1,remd_m(1))
891 cd write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
892 cd & (nupa(k,il),k=1,nupa(0,il))
893 cd write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
894 cd & (ndowna(k,il),k=1,ndowna(0,il))
899 remd_ene(0,iex)=ene_iex_iex
900 remd_ene(0,i)=ene_i_i
906 cd write (iout,*) "exchange completed"
910 cd write(iout,*) "########",ii
912 i_temp=iran_num(1,nrep)
913 i_mult=iran_num(1,remd_m(i_temp))
914 i_iset=iran_num(1,nset)
915 i_mset=iran_num(1,mset(i_iset))
916 i=i_index(i_temp,i_mult,i_iset,i_mset)
918 cd write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
921 cd write(iout,*) "i_dir=",i_dir
923 if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then
926 i_mult1=iran_num(1,remd_m(i_temp1))
928 i_mset1=iran_num(1,mset(i_iset1))
929 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
931 elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
934 i_mult1=iran_num(1,remd_m(i_temp1))
936 i_mset1=iran_num(1,mset(i_iset1))
937 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
938 econstr_temp_i=remd_ene(20,i)
939 econstr_temp_iex=remd_ene(20,iex)
940 remd_ene(20,i)=remd_ene(n_ene+3,i)
941 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
943 elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
946 i_mult1=iran_num(1,remd_m(i_temp1))
948 i_mset1=iran_num(1,mset(i_iset1))
949 iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
950 econstr_temp_i=remd_ene(20,i)
951 econstr_temp_iex=remd_ene(20,iex)
952 remd_ene(20,i)=remd_ene(n_ene+3,i)
953 remd_ene(20,iex)=remd_ene(n_ene+4,iex)
959 cd write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
962 c Swap temperatures between conformations i and iex with recalculating the free energies
963 c following temperature changes.
964 ene_iex_iex=remd_ene(0,iex)
965 ene_i_i=remd_ene(0,i)
966 co write (iout,*) "rescaling weights with temperature",
968 call rescale_weights(remd_t_bath(i))
970 call sum_energy(remd_ene(0,iex),.false.)
971 ene_iex_i=remd_ene(0,iex)
972 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
973 c call sum_energy(remd_ene(0,i),.false.)
974 cd write (iout,*) "ene_i_i",remd_ene(0,i)
975 c write (iout,*) "rescaling weights with temperature",
977 c if (real(ene_i_i).ne.real(remd_ene(0,i))) then
978 c write (iout,*) "ERROR: inconsistent energies:",i,
979 c & ene_i_i,remd_ene(0,i)
981 call rescale_weights(remd_t_bath(iex))
982 call sum_energy(remd_ene(0,i),.false.)
983 cd write (iout,*) "ene_i_iex",remd_ene(0,i)
984 ene_i_iex=remd_ene(0,i)
985 c call sum_energy(remd_ene(0,iex),.false.)
986 c if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
987 c write (iout,*) "ERROR: inconsistent energies:",iex,
988 c & ene_iex_iex,remd_ene(0,iex)
990 cd write (iout,*) "ene_iex_iex",remd_ene(0,iex)
991 c write (iout,*) "i",i," iex",iex
992 cd write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
993 cd & " ene_i_iex",ene_i_iex,
994 cd & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
995 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
996 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
998 cd write(iout,*) 'delta',delta
999 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
1000 c & (remd_ene(i)-remd_ene(iex))/Rb/
1001 c & (remd_t_bath(i)*remd_t_bath(iex))
1002 if (delta .gt. 50.0d0) then
1007 if (i_dir.eq.1.or.i_dir.eq.3)
1008 & iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1009 if (i_dir.eq.2.or.i_dir.eq.3)
1010 & iremd_tot_usa(int(i2set(i-1)))=
1011 & iremd_tot_usa(int(i2set(i-1)))+1
1012 xxx=ran_number(0.0d0,1.0d0)
1013 cd write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1014 if (delta .gt. xxx) then
1016 remd_t_bath(i)=remd_t_bath(iex)
1017 remd_t_bath(iex)=tmp
1020 iremd_iset(i)=iremd_iset(iex)
1021 iremd_iset(iex)=itmp
1023 remd_ene(0,i)=ene_i_iex
1024 remd_ene(0,iex)=ene_iex_i
1026 if (i_dir.eq.1.or.i_dir.eq.3)
1027 & iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1030 i2rep(i-1)=i2rep(iex-1)
1033 if (i_dir.eq.2.or.i_dir.eq.3)
1034 & iremd_acc_usa(int(i2set(i-1)))=
1035 & iremd_acc_usa(int(i2set(i-1)))+1
1038 i2set(i-1)=i2set(iex-1)
1041 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1042 i_index(i_temp,i_mult,i_iset,i_mset)=
1043 & i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1044 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1047 remd_ene(0,iex)=ene_iex_iex
1048 remd_ene(0,i)=ene_i_i
1049 remd_ene(20,iex)=econstr_temp_iex
1050 remd_ene(20,i)=econstr_temp_i
1054 cd do il1=1,mset(il)
1057 cd write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1070 c-------------------------------------
1071 write (iout,*) "NREP",nrep
1073 if(iremd_tot(i).ne.0)
1074 & write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1075 & ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1080 if(iremd_tot_usa(i).ne.0)
1081 & write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1082 & iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1088 cd write (iout,'(a6,100i4)') "ifirst",
1089 cd & (ifirst(i),i=1,remd_m(1))
1091 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1092 cd & (nupa(i,il),i=1,nupa(0,il))
1093 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1094 cd & (ndowna(i,il),i=1,ndowna(0,il))
1099 cd write (iout,*) "Before scatter"
1101 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1102 & t_bath,1,mpi_double_precision,king,
1104 cd write (iout,*) "After scatter"
1107 & call mpi_scatter(iremd_iset,1,mpi_integer,
1108 & iset,1,mpi_integer,king,
1112 if (me.eq.king .or. .not. out1file) then
1113 write(iout,*) 'REMD scatter time=',time07-time06
1117 call mpi_scatter(elowi,1,mpi_double_precision,
1118 & elow,1,mpi_double_precision,king,
1120 call mpi_scatter(ehighi,1,mpi_double_precision,
1121 & ehigh,1,mpi_double_precision,king,
1124 call rescale_weights(t_bath)
1125 co write (iout,*) "Processor",me,
1126 co & " rescaling weights with temperature",t_bath
1128 stdfp=dsqrt(2*Rb*t_bath/d_time)
1130 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1133 cde write(iout,*) 'REMD after',me,t_bath
1135 if (me.eq.king .or. .not. out1file) then
1136 write(iout,*) 'REMD exchange time=',time08-time00
1142 if (restart1file) then
1143 if (me.eq.king .or. .not. out1file)
1144 & write(iout,*) 'writing restart at the end of run'
1145 call write1rst(i_index)
1148 if (traj1file) call write1traj
1150 cdeb call mpi_gather(ntwx_cache,1,mpi_integer,
1151 cdeb & icache_all,1,mpi_integer,king,
1152 cdeb & CG_COMM,ierr)
1153 cdeb write(iout,'(a40,8000i8)')
1154 cdeb & ' ntwx_cache after traj1file at the end',
1155 cdeb & (icache_all(i),i=1,nodes)
1164 if (me.eq.king .or. .not. out1file) then
1165 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
1167 & 'MD calculations setup:',t_MDsetup,
1168 & 'Energy & gradient evaluation:',t_enegrad,
1169 & 'Stochastic MD setup:',t_langsetup,
1170 & 'Stochastic MD step setup:',t_sdsetup,
1172 write (iout,'(/28(1h=),a25,27(1h=))')
1173 & ' End of MD calculation '
1178 c-----------------------------------------------------------------------
1180 subroutine write1rst_oldtxt
1181 implicit real*8 (a-h,o-z)
1182 include 'DIMENSIONS'
1185 include 'COMMON.IOUNITS'
1186 include 'COMMON.REMD'
1187 include 'COMMON.SETUP'
1188 include 'COMMON.CHAIN'
1189 include 'COMMON.SBRIDGE'
1190 include 'COMMON.INTERACT'
1192 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1193 & d_restart2(3,2*maxres*maxprocs)
1201 t5_restart1(4)=t_bath
1202 t5_restart1(5)=Uconst
1204 call mpi_gather(t5_restart1,5,mpi_real,
1205 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1213 call mpi_gather(r_d,3*2*nres,mpi_real,
1214 & d_restart1,3*2*nres,mpi_real,king,
1223 call mpi_gather(r_d,3*2*nres,mpi_real,
1224 & d_restart2,3*2*nres,mpi_real,king,
1228 open(irest1,file=mremd_rst_name,status='unknown')
1229 write (irest1,*) (i2rep(i),i=0,nodes-1)
1230 write (irest1,*) (ifirst(i),i=1,remd_m(1))
1232 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1233 write (irest1,*) ndowna(0,il),
1234 & (ndowna(i,il),i=1,ndowna(0,il))
1238 write(irest1,*) (t_restart1(j,il),j=1,4)
1243 write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
1248 write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
1257 subroutine write1rst(i_index)
1258 implicit real*8 (a-h,o-z)
1259 include 'DIMENSIONS'
1262 include 'COMMON.IOUNITS'
1263 include 'COMMON.REMD'
1264 include 'COMMON.SETUP'
1265 include 'COMMON.CHAIN'
1266 include 'COMMON.SBRIDGE'
1267 include 'COMMON.INTERACT'
1269 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1270 & d_restart2(3,2*maxres*maxprocs)
1274 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1280 t5_restart1(4)=t_bath
1281 t5_restart1(5)=Uconst
1283 call mpi_gather(t5_restart1,5,mpi_real,
1284 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1292 call mpi_gather(r_d,3*2*nres,mpi_real,
1293 & d_restart1,3*2*nres,mpi_real,king,
1302 call mpi_gather(r_d,3*2*nres,mpi_real,
1303 & d_restart2,3*2*nres,mpi_real,king,
1307 c open(irest1,file=mremd_rst_name,status='unknown')
1308 call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1309 c write (irest1,*) (i2rep(i),i=0,nodes-1)
1311 call xdrfint(ixdrf, i2rep(i), iret)
1313 c write (irest1,*) (ifirst(i),i=1,remd_m(1))
1315 call xdrfint(ixdrf, ifirst(i), iret)
1318 c write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1320 call xdrfint(ixdrf, nupa(i,il), iret)
1323 c write (irest1,*) ndowna(0,il),
1324 c & (ndowna(i,il),i=1,ndowna(0,il))
1326 call xdrfint(ixdrf, ndowna(i,il), iret)
1331 c write(irest1,*) (t_restart1(j,il),j=1,4)
1333 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1339 c write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
1341 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1347 c write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
1349 call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1356 c write (irest1,*) nset
1357 call xdrfint(ixdrf, nset, iret)
1358 c write (irest1,*) (mset(i),i=1,nset)
1360 call xdrfint(ixdrf,mset(i), iret)
1362 c write (irest1,*) (i2set(i),i=0,nodes-1)
1364 call xdrfint(ixdrf,i2set(i), iret)
1366 c write (irest1,*) "i_index"
1370 c write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
1372 itmp=i_index(i,j,il,il1)
1373 call xdrfint(ixdrf,itmp, iret)
1382 call xdrfclose(ixdrf, iret)
1389 subroutine write1traj
1390 implicit real*8 (a-h,o-z)
1391 include 'DIMENSIONS'
1394 include 'COMMON.IOUNITS'
1395 include 'COMMON.REMD'
1396 include 'COMMON.SETUP'
1397 include 'COMMON.CHAIN'
1398 include 'COMMON.SBRIDGE'
1399 include 'COMMON.INTERACT'
1403 real xcoord(3,maxres2+2),prec
1404 real r_qfrag(50),r_qpair(100)
1405 real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1406 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1407 real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1408 & p_uscdiff(100*maxprocs)
1409 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1411 call mpi_bcast(ii_write,1,mpi_integer,
1412 & king,CG_COMM,ierr)
1415 print *,'traj1file',me,ii_write,ntwx_cache
1418 if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1420 t5_restart1(1)=totT_cache(ii)
1421 t5_restart1(2)=EK_cache(ii)
1422 t5_restart1(3)=potE_cache(ii)
1423 t5_restart1(4)=t_bath_cache(ii)
1424 t5_restart1(5)=Uconst_cache(ii)
1425 call mpi_gather(t5_restart1,5,mpi_real,
1426 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
1428 call mpi_gather(iset_cache(ii),1,mpi_integer,
1429 & iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1432 r_qfrag(i)=qfrag_cache(i,ii)
1435 r_qpair(i)=qpair_cache(i,ii)
1438 r_utheta(i)=utheta_cache(i,ii)
1439 r_ugamma(i)=ugamma_cache(i,ii)
1440 r_uscdiff(i)=uscdiff_cache(i,ii)
1443 call mpi_gather(r_qfrag,nfrag,mpi_real,
1444 & p_qfrag,nfrag,mpi_real,king,
1446 call mpi_gather(r_qpair,npair,mpi_real,
1447 & p_qpair,npair,mpi_real,king,
1449 call mpi_gather(r_utheta,nfrag_back,mpi_real,
1450 & p_utheta,nfrag_back,mpi_real,king,
1452 call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1453 & p_ugamma,nfrag_back,mpi_real,king,
1455 call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1456 & p_uscdiff,nfrag_back,mpi_real,king,
1460 write (iout,*) "p_qfrag"
1462 write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1464 write (iout,*) "p_qpair"
1466 write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1472 r_c(j,i)=c_cache(j,i,ii)
1476 call mpi_gather(r_c,3*2*nres,mpi_real,
1477 & p_c,3*2*nres,mpi_real,king,
1483 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1484 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1485 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1486 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1487 call xdrfint(ixdrf, nss, iret)
1489 call xdrfint(ixdrf, ihpb(j), iret)
1490 call xdrfint(ixdrf, jhpb(j), iret)
1492 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1493 call xdrfint(ixdrf, iset_restart1(il), iret)
1495 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1498 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1501 call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1502 call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1503 call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1508 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1513 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1517 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1521 if(me.eq.king) call xdrfclose(ixdrf, iret)
1522 do i=1,ntwx_cache-ii_write
1524 totT_cache(i)=totT_cache(ii_write+i)
1525 EK_cache(i)=EK_cache(ii_write+i)
1526 potE_cache(i)=potE_cache(ii_write+i)
1527 t_bath_cache(i)=t_bath_cache(ii_write+i)
1528 Uconst_cache(i)=Uconst_cache(ii_write+i)
1529 iset_cache(i)=iset_cache(ii_write+i)
1532 qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1535 qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1538 utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1539 ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1540 uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1545 c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1549 ntwx_cache=ntwx_cache-ii_write
1554 subroutine read1restart(i_index)
1555 implicit real*8 (a-h,o-z)
1556 include 'DIMENSIONS'
1559 include 'COMMON.IOUNITS'
1560 include 'COMMON.REMD'
1561 include 'COMMON.SETUP'
1562 include 'COMMON.CHAIN'
1563 include 'COMMON.SBRIDGE'
1564 include 'COMMON.INTERACT'
1565 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1568 & (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1569 write (*,*) "Processor",me," called read1restart"
1572 open(irest2,file=mremd_rst_name,status='unknown')
1573 read(irest2,*,err=334) i
1574 write(iout,*) "Reading old rst in ASCI format"
1576 call read1restart_old
1579 call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1581 c read (irest2,*) (i2rep(i),i=0,nodes-1)
1583 call xdrfint(ixdrf, i2rep(i), iret)
1585 c read (irest2,*) (ifirst(i),i=1,remd_m(1))
1587 call xdrfint(ixdrf, ifirst(i), iret)
1590 c read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1591 call xdrfint(ixdrf, nupa(0,il), iret)
1593 call xdrfint(ixdrf, nupa(i,il), iret)
1596 c read (irest2,*) ndowna(0,il),
1597 c & (ndowna(i,il),i=1,ndowna(0,il))
1598 call xdrfint(ixdrf, ndowna(0,il), iret)
1600 call xdrfint(ixdrf, ndowna(i,il), iret)
1604 c read(irest2,*) (t_restart1(j,il),j=1,4)
1606 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1610 call mpi_scatter(t_restart1,5,mpi_real,
1611 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1615 t_bath=t5_restart1(4)
1620 c read(irest2,'(3e15.5)')
1621 c & (d_restart1(j,i+2*nres*il),j=1,3)
1623 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1628 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1629 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1639 c read(irest2,'(3e15.5)')
1640 c & (d_restart1(j,i+2*nres*il),j=1,3)
1642 call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1647 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1648 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1658 call xdrfint(ixdrf, nset, iret)
1660 call xdrfint(ixdrf,mset(i), iret)
1663 call xdrfint(ixdrf,i2set(i), iret)
1669 call xdrfint(ixdrf,itmp, iret)
1670 i_index(i,j,il,il1)=itmp
1677 call mpi_scatter(i2set,1,mpi_integer,
1678 & iset,1,mpi_integer,king,
1684 if(me.eq.king) close(irest2)
1688 subroutine read1restart_old
1689 implicit real*8 (a-h,o-z)
1690 include 'DIMENSIONS'
1693 include 'COMMON.IOUNITS'
1694 include 'COMMON.REMD'
1695 include 'COMMON.SETUP'
1696 include 'COMMON.CHAIN'
1697 include 'COMMON.SBRIDGE'
1698 include 'COMMON.INTERACT'
1699 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1703 open(irest2,file=mremd_rst_name,status='unknown')
1704 read (irest2,*) (i2rep(i),i=0,nodes-1)
1705 read (irest2,*) (ifirst(i),i=1,remd_m(1))
1707 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1708 read (irest2,*) ndowna(0,il),
1709 & (ndowna(i,il),i=1,ndowna(0,il))
1712 read(irest2,*) (t_restart1(j,il),j=1,4)
1715 call mpi_scatter(t_restart1,5,mpi_real,
1716 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1720 t_bath=t5_restart1(4)
1725 read(irest2,'(3e15.5)')
1726 & (d_restart1(j,i+2*nres*il),j=1,3)
1730 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1731 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1741 read(irest2,'(3e15.5)')
1742 & (d_restart1(j,i+2*nres*il),j=1,3)
1746 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1747 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1753 if(me.eq.king) close(irest2)