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 double precision remd_ene(0:n_ene+1,maxprocs)
28 integer iremd_acc(maxprocs),iremd_tot(maxprocs)
34 integer nup(0:maxprocs),ndown(0:maxprocs)
35 integer rep2i(0:maxprocs)
36 integer itime_all(maxprocs)
37 integer status(MPI_STATUS_SIZE)
38 logical synflag,end_of_run,file_exist
41 if(me.eq.king.or..not.out1file)
42 & write (iout,*) 'MREMD',nodes,'time before',time00-walltime
45 mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
47 cd print *,'MREMD',nodes
49 cd print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
51 cde write (iout,*) "Start MREMD: me",me," t_bath",t_bath
64 c print *,'i2rep',me,i2rep(me)
65 c print *,'rep2i',(rep2i(i),i=0,nrep)
67 if (i2rep(me).eq.nrep) then
70 nup(0)=remd_m(i2rep(me)+1)
78 cd print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
80 if (i2rep(me).eq.1) then
83 ndown(0)=remd_m(i2rep(me)-1)
84 k=rep2i(i2rep(me)-2)+1
91 cd print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
94 if(rest.and.restart1file) then
95 inquire(file=mremd_rst_name,exist=file_exist)
96 if(file_exist) call read1restart
100 if (rest.and..not.restart1file)
101 & inquire(file=mremd_rst_name,exist=file_exist)
102 IF (rest.and.file_exist.and..not.restart1file) THEN
103 open(irest2,file=mremd_rst_name,status='unknown')
105 read (irest2,*) (i2rep(i),i=0,nodes-1)
107 read (irest2,*) (ifirst(i),i=1,remd_m(1))
110 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
112 read (irest2,*) ndowna(0,il),
113 & (ndowna(i,il),i=1,ndowna(0,il))
117 write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
118 write (iout,'(a6,1000i5)') "ifirst",
119 & (ifirst(i),i=1,remd_m(1))
121 write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
122 & (nupa(i,il),i=1,nupa(0,il))
123 write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
124 & (ndowna(i,il),i=1,ndowna(0,il))
133 if (i2rep(il-1).eq.nrep) then
136 nupa(0,il)=remd_m(i2rep(il-1)+1)
137 k=rep2i(i2rep(il-1))+1
144 if (i2rep(il-1).eq.1) then
147 ndowna(0,il)=remd_m(i2rep(il-1)-1)
148 k=rep2i(i2rep(il-1)-2)+1
157 c write (iout,'(a6,100i4)') "ifirst",
158 c & (ifirst(i),i=1,remd_m(1))
160 c write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
161 c & (nupa(i,il),i=1,nupa(0,il))
162 c write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
163 c & (ndowna(i,il),i=1,ndowna(0,il))
169 c t_bath=retmin+(retmax-retmin)*me/(nodes-1)
170 if(.not.(rest.and.file_exist.and.restart1file)) then
171 if (me .eq. king) then
174 t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
176 cd print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
177 if (remd_tlist) t_bath=remd_t(i2rep(me))
180 stdfp=dsqrt(2*Rb*t_bath/d_time)
182 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
185 c print *,'irep',me,t_bath
187 if (me.eq.king .or. .not. out1file)
188 & write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
189 call rescale_weights(t_bath)
193 c------copy MD--------------
194 c The driver for molecular dynamics subroutines
195 c------------------------------------------------
201 if(me.eq.king.or..not.out1file)
202 & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
204 c Determine the inverse of the inertia matrix.
205 call setup_MD_matrices
209 if (me.eq.king .or. .not. out1file)
210 & write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
211 stdfp=dsqrt(2*Rb*t_bath/d_time)
213 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
215 call rescale_weights(t_bath)
218 t_MDsetup = tcpu()-tt0
220 c Entering the MD loop
222 if (lang.eq.2 .or. lang.eq.3) then
225 call sd_verlet_p_setup
227 call sd_verlet_ciccotti_setup
232 pfric0_mat(i,j,0)=pfric_mat(i,j)
233 afric0_mat(i,j,0)=afric_mat(i,j)
234 vfric0_mat(i,j,0)=vfric_mat(i,j)
235 prand0_mat(i,j,0)=prand_mat(i,j)
236 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
237 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
243 flag_stoch(i)=.false.
245 else if (lang.eq.1 .or. lang.eq.4) then
249 if (me.eq.king .or. .not. out1file)
250 & write(iout,*) 'Setup time',time00-walltime
252 t_langsetup=tcpu()-tt0
256 do while(.not.end_of_run)
258 if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
259 if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
261 if (lang.gt.0 .and. surfarea .and.
262 & mod(itime,reset_fricmat).eq.0) then
263 if (lang.eq.2 .or. lang.eq.3) then
266 call sd_verlet_p_setup
268 call sd_verlet_ciccotti_setup
273 pfric0_mat(i,j,0)=pfric_mat(i,j)
274 afric0_mat(i,j,0)=afric_mat(i,j)
275 vfric0_mat(i,j,0)=vfric_mat(i,j)
276 prand0_mat(i,j,0)=prand_mat(i,j)
277 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
278 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
284 flag_stoch(i)=.false.
286 else if (lang.eq.1 .or. lang.eq.4) then
289 write (iout,'(a,i10)')
290 & "Friction matrix reset based on surface area, itime",itime
292 if (reset_vel .and. tbf .and. lang.eq.0
293 & .and. mod(itime,count_reset_vel).eq.0) then
295 if (me.eq.king .or. .not. out1file)
296 & write(iout,'(a,f20.2)')
297 & "Velocities reset to random values, time",totT
300 d_t_old(j,i)=d_t(j,i)
304 if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
308 d_t(j,0)=d_t(j,0)-vcm(j)
311 kinetic_T=2.0d0/(dimen*Rb)*EK
312 scalfac=dsqrt(T_bath/kinetic_T)
313 cd write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT
316 d_t_old(j,i)=scalfac*d_t(j,i)
322 c Time-reversible RESPA algorithm
323 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
324 call RESPA_step(itime)
326 c Variable time step algorithm.
327 call velverlet_step(itime)
330 call brown_step(itime)
333 if (mod(itime,ntwe).eq.0) call statout(itime)
335 if (mod(itime,ntwx).eq.0.and..not.traj1file) then
336 write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
338 call pdbout(potE,tytul,ipdb)
343 if ((rstcount.eq.1000.or.itime.eq.n_timestep)
344 & .and..not.restart1file) then
347 open(irest1,file=mremd_rst_name,status='unknown')
348 write (irest1,*) "i2rep"
349 write (irest1,*) (i2rep(i),i=0,nodes-1)
350 write (irest1,*) "ifirst"
351 write (irest1,*) (ifirst(i),i=1,remd_m(1))
353 write (irest1,*) "nupa",il
354 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
355 write (irest1,*) "ndowna",il
356 write (irest1,*) ndowna(0,il),
357 & (ndowna(i,il),i=1,ndowna(0,il))
361 open(irest2,file=rest2name,status='unknown')
362 write(irest2,*) totT,EK,potE,totE,t_bath
364 write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
367 write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
374 c forced synchronization
375 if (mod(itime,100).eq.0 .and. me.ne.king
376 & .and. .not. mremdsync) then
378 call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
380 call mpi_recv(itime_master, 1, MPI_INTEGER,
381 & 0,101,CG_COMM, status, ierr)
382 if (.not. out1file) then
383 write(iout,*) 'REMD synchro at',itime_master,itime
385 call mpi_gather(itime,1,mpi_integer,
386 & itime_all,1,mpi_integer,king,
389 if(itime_master.ge.n_timestep) end_of_run=.true.
395 if ((mod(itime,nstex).eq.0.and.me.eq.king
396 & .or.end_of_run.and.me.eq.king )
397 & .and. .not. mremdsync ) then
401 call mpi_isend(itime,1,MPI_INTEGER,i,101,
402 & CG_COMM, ireq, ierr)
403 cd write(iout,*) 'REMD synchro with',i
407 write(iout,*) 'REMD synchro at',itime,'time=',time02-time00
409 call mpi_gather(itime,1,mpi_integer,
410 & itime_all,1,mpi_integer,king,
412 write(iout,'(a18,8000i8)') 'REMD synchro itime',
413 & (itime_all(i),i=1,nodes)
417 if(mremdsync .and. mod(itime,nstex).eq.0) then
419 if (me.eq.king .or. .not. out1file)
420 & write(iout,*) 'REMD synchro at',itime
423 if(synflag.and..not.end_of_run) then
427 cd write(iout,*) 'REMD before',me,t_bath
429 c call mpi_gather(t_bath,1,mpi_double_precision,
430 c & remd_t_bath,1,mpi_double_precision,king,
432 potEcomp(n_ene+1)=t_bath
433 call mpi_gather(potEcomp(0),n_ene+2,mpi_double_precision,
434 & remd_ene(0,1),n_ene+2,mpi_double_precision,king,
437 call mpi_gather(elow,1,mpi_double_precision,
438 & elowi,1,mpi_double_precision,king,
440 call mpi_gather(ehigh,1,mpi_double_precision,
441 & ehighi,1,mpi_double_precision,king,
445 if (restart1file) call write1rst
446 if (traj1file) call write1traj
450 remd_t_bath(i)=remd_ene(n_ene+1,i)
453 co write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
455 write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
459 cd write(iout,*) 'REMD exchange temp,ene'
461 co write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
462 c write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
465 c-------------------------------------
467 i=ifirst(iran_num(1,remd_m(1)))
470 if(i.gt.0.and.nupa(0,i).gt.0) then
471 iex=nupa(iran_num(1,int(nupa(0,i))),i)
473 call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
475 c Swap temperatures between conformations i and iex with recalculating the free energies
476 c following temperature changes.
477 ene_iex_iex=remd_ene(0,iex)
478 ene_i_i=remd_ene(0,i)
479 co write (iout,*) "rescaling weights with temperature",
481 call rescale_weights(remd_t_bath(i))
482 call sum_energy(remd_ene(0,iex),.false.)
483 ene_iex_i=remd_ene(0,iex)
484 cd write (iout,*) "ene_iex_i",remd_ene(0,iex)
485 call sum_energy(remd_ene(0,i),.false.)
486 cd write (iout,*) "ene_i_i",remd_ene(0,i)
487 c write (iout,*) "rescaling weights with temperature",
489 if (real(ene_i_i).ne.real(remd_ene(0,i))) then
490 write (iout,*) "ERROR: inconsistent energies:",i,
491 & ene_i_i,remd_ene(0,i)
493 call rescale_weights(remd_t_bath(iex))
494 call sum_energy(remd_ene(0,i),.false.)
495 c write (iout,*) "ene_i_iex",remd_ene(0,i)
496 ene_i_iex=remd_ene(0,i)
497 call sum_energy(remd_ene(0,iex),.false.)
498 if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
499 write (iout,*) "ERROR: inconsistent energies:",iex,
500 & ene_iex_iex,remd_ene(0,iex)
502 c write (iout,*) "ene_iex_iex",remd_ene(0,iex)
503 c write (iout,*) "i",i," iex",iex
504 co write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
505 co & " ene_i_iex",ene_i_iex,
506 co & " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
507 delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
508 & (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
510 c write(iout,*) 'delta',delta
511 c delta=(remd_t_bath(i)-remd_t_bath(iex))*
512 c & (remd_ene(i)-remd_ene(iex))/Rb/
513 c & (remd_t_bath(i)*remd_t_bath(iex))
515 if (delta .gt. 50.0d0) then
521 else if (delta.lt.-50.0d0) then
530 iremd_tot(i2rep(i-1))=iremd_tot(i2rep(i-1))+1
531 xxx=ran_number(0.0d0,1.0d0)
532 co write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
533 if (delta .gt. xxx) then
535 remd_t_bath(i)=remd_t_bath(iex)
537 remd_ene(0,i)=ene_i_iex
538 remd_ene(0,iex)=ene_iex_i
544 ehighi(i)=ehighi(iex)
551 nupa(k,i)=nupa(k,iex)
554 ndowna(k,i)=ndowna(k,iex)
558 if (ifirst(il).eq.i) ifirst(il)=iex
560 if (nupa(k,il).eq.i) then
562 elseif (nupa(k,il).eq.iex) then
567 if (ndowna(k,il).eq.i) then
569 elseif (ndowna(k,il).eq.iex) then
575 iremd_acc(i2rep(i-1))=iremd_acc(i2rep(i-1))+1
577 i2rep(i-1)=i2rep(iex-1)
580 cd write(iout,*) 'exchange',i,iex
581 cd write (iout,'(a8,100i4)') "@ ifirst",
582 cd & (ifirst(k),k=1,remd_m(1))
584 cd write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
585 cd & (nupa(k,il),k=1,nupa(0,il))
586 cd write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
587 cd & (ndowna(k,il),k=1,ndowna(0,il))
591 remd_ene(0,iex)=ene_iex_iex
592 remd_ene(0,i)=ene_i_i
598 c-------------------------------------
600 if(iremd_tot(i).ne.0)
601 & write(iout,'(a3,i4,2f12.5)') 'ACC',i,remd_t(i)
602 & ,iremd_acc(i)/(1.0*iremd_tot(i))
605 cd write (iout,'(a6,100i4)') "ifirst",
606 cd & (ifirst(i),i=1,remd_m(1))
608 cd write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
609 cd & (nupa(i,il),i=1,nupa(0,il))
610 cd write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
611 cd & (ndowna(i,il),i=1,ndowna(0,il))
615 call mpi_scatter(remd_t_bath,1,mpi_double_precision,
616 & t_bath,1,mpi_double_precision,king,
619 call mpi_scatter(elowi,1,mpi_double_precision,
620 & elow,1,mpi_double_precision,king,
622 call mpi_scatter(ehighi,1,mpi_double_precision,
623 & ehigh,1,mpi_double_precision,king,
626 call rescale_weights(t_bath)
627 co write (iout,*) "Processor",me,
628 co & " rescaling weights with temperature",t_bath
630 stdfp=dsqrt(2*Rb*t_bath/d_time)
632 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
635 cde write(iout,*) 'REMD after',me,t_bath
637 if (me.eq.king .or. .not. out1file) then
638 write(iout,*) 'REMD exchange time=',time02-time00
644 if (restart1file) then
645 if (me.eq.king .or. .not. out1file)
646 & write(iout,*) 'writing restart at the end of run'
650 if (traj1file) call write1traj
654 if (me.eq.king .or. .not. out1file) then
655 write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))')
657 & 'MD calculations setup:',t_MDsetup,
658 & 'Energy & gradient evaluation:',t_enegrad,
659 & 'Stochastic MD setup:',t_langsetup,
660 & 'Stochastic MD step setup:',t_sdsetup,
662 write (iout,'(/28(1h=),a25,27(1h=))')
663 & ' End of MD calculation '
669 implicit real*8 (a-h,o-z)
673 include 'COMMON.IOUNITS'
674 include 'COMMON.REMD'
675 include 'COMMON.SETUP'
676 include 'COMMON.CHAIN'
677 include 'COMMON.SBRIDGE'
678 include 'COMMON.INTERACT'
680 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
681 & d_restart2(3,2*maxres*maxprocs)
689 t5_restart1(4)=t_bath
690 t5_restart1(5)=Uconst
692 call mpi_gather(t5_restart1,5,mpi_real,
693 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
701 call mpi_gather(r_d,3*2*nres,mpi_real,
702 & d_restart1,3*2*nres,mpi_real,king,
711 call mpi_gather(r_d,3*2*nres,mpi_real,
712 & d_restart2,3*2*nres,mpi_real,king,
716 open(irest1,file=mremd_rst_name,status='unknown')
717 write (irest1,*) (i2rep(i),i=0,nodes-1)
718 write (irest1,*) (ifirst(i),i=1,remd_m(1))
720 write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
721 write (irest1,*) ndowna(0,il),
722 & (ndowna(i,il),i=1,ndowna(0,il))
726 write(irest1,*) (t_restart1(j,il),j=1,4)
731 write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
736 write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
746 subroutine write1traj
747 implicit real*8 (a-h,o-z)
751 include 'COMMON.IOUNITS'
752 include 'COMMON.REMD'
753 include 'COMMON.SETUP'
754 include 'COMMON.CHAIN'
755 include 'COMMON.SBRIDGE'
756 include 'COMMON.INTERACT'
760 real xcoord(3,maxres2+2),prec
761 real r_qfrag(50),r_qpair(100)
762 real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
763 real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
765 if(.not.restart1file) then
769 t5_restart1(4)=t_bath
770 t5_restart1(5)=Uconst
771 call mpi_gather(t5_restart1,5,mpi_real,
772 & t_restart1,5,mpi_real,king,CG_COMM,ierr)
782 call mpi_gather(r_qfrag,nfrag,mpi_real,
783 & p_qfrag,nfrag,mpi_real,king,
785 call mpi_gather(qpair,nfrag,mpi_real,
786 & p_qpair,nfrag,mpi_real,king,
795 call mpi_gather(r_c,3*2*nres,mpi_real,
796 & p_c,3*2*nres,mpi_real,king,
800 call xdrfopen(ixdrf,cartname, "a", iret)
802 call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
803 call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
804 call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
805 call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
806 call xdrfint(ixdrf, nss, iret)
808 call xdrfint(ixdrf, ihpb(j), iret)
809 call xdrfint(ixdrf, jhpb(j), iret)
811 call xdrfint(ixdrf, nfrag+npair, iret)
813 call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
816 call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
821 xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
826 xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
830 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
832 call xdrfclose(ixdrf, iret)
839 subroutine read1restart
840 implicit real*8 (a-h,o-z)
844 include 'COMMON.IOUNITS'
845 include 'COMMON.REMD'
846 include 'COMMON.SETUP'
847 include 'COMMON.CHAIN'
848 include 'COMMON.SBRIDGE'
849 include 'COMMON.INTERACT'
850 real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
854 open(irest2,file=mremd_rst_name,status='unknown')
855 read (irest2,*) (i2rep(i),i=0,nodes-1)
856 read (irest2,*) (ifirst(i),i=1,remd_m(1))
858 read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
859 read (irest2,*) ndowna(0,il),
860 & (ndowna(i,il),i=1,ndowna(0,il))
863 read(irest2,*) (t_restart1(j,il),j=1,4)
866 call mpi_scatter(t_restart1,5,mpi_real,
867 & t5_restart1,5,mpi_real,king,CG_COMM,ierr)
871 t_bath=t5_restart1(4)
876 read(irest2,'(3e15.5)')
877 & (d_restart1(j,i+2*nres*il),j=1,3)
881 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
882 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
892 read(irest2,'(3e15.5)')
893 & (d_restart1(j,i+2*nres*il),j=1,3)
897 call mpi_scatter(d_restart1,3*2*nres,mpi_real,
898 & r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
904 if(me.eq.king) close(irest2)