wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / MREMD_nosy1traj.F
1       subroutine MREMD
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'mpif.h'
5       include 'COMMON.CONTROL'
6       include 'COMMON.VAR'
7       include 'COMMON.MD'
8 #ifndef LANG0
9       include 'COMMON.LANGEVIN'
10 #else
11       include 'COMMON.LANGEVIN.lang0'
12 #endif
13       include 'COMMON.CHAIN'
14       include 'COMMON.DERIV'
15       include 'COMMON.GEO'
16       include 'COMMON.LOCAL'
17       include 'COMMON.INTERACT'
18       include 'COMMON.IOUNITS'
19       include 'COMMON.NAMES'
20       include 'COMMON.TIME1'
21       include 'COMMON.REMD'
22       include 'COMMON.SETUP'
23       include 'COMMON.MUCA'
24       include 'COMMON.HAIRPIN'
25       double precision cm(3),L(3),vcm(3)
26       double precision energia(0:n_ene)
27       double precision remd_t_bath(maxprocs)
28       double precision remd_ene(0:n_ene+1,maxprocs)
29       integer iremd_acc(maxprocs),iremd_tot(maxprocs)
30       integer ilen,rstcount
31       external ilen
32       character*50 tytul
33       common /gucio/ cm
34       integer itime
35       integer nup(0:maxprocs),ndown(0:maxprocs)
36       integer rep2i(0:maxprocs)
37       integer itime_all(maxprocs)
38       integer status(MPI_STATUS_SIZE)
39       logical synflag,end_of_run,file_exist
40
41       time00=MPI_WTIME()
42       if(me.eq.king.or..not.out1file)
43      & write  (iout,*) 'MREMD',nodes,'time before',time00-walltime
44
45       synflag=.false.
46       mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
47
48 cd      print *,'MREMD',nodes
49
50 cd      print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
51
52 cde      write (iout,*) "Start MREMD: me",me," t_bath",t_bath
53       k=0
54       rep2i(k)=-1
55       do i=1,nrep
56         iremd_acc(i)=0
57         iremd_tot(i)=0
58         do j=1,remd_m(i)
59           i2rep(k)=i
60           rep2i(i)=k
61           k=k+1
62         enddo
63       enddo
64
65 c      print *,'i2rep',me,i2rep(me)
66 c      print *,'rep2i',(rep2i(i),i=0,nrep)
67
68        if (i2rep(me).eq.nrep) then
69         nup(0)=0
70        else
71         nup(0)=remd_m(i2rep(me)+1)
72         k=rep2i(i2rep(me))+1
73         do i=1,nup(0)
74          nup(i)=k
75          k=k+1
76         enddo
77        endif
78
79 cd       print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
80
81        if (i2rep(me).eq.1) then
82         ndown(0)=0
83        else
84         ndown(0)=remd_m(i2rep(me)-1)
85         k=rep2i(i2rep(me)-2)+1
86         do i=1,ndown(0)
87          ndown(i)=k
88          k=k+1
89         enddo
90        endif
91
92 cd       print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
93
94        
95        if(rest.and.restart1file) then 
96            inquire(file=mremd_rst_name,exist=file_exist)
97            if(file_exist) call read1restart
98        endif
99
100        if(me.eq.king) then
101         if (rest.and..not.restart1file) 
102      &          inquire(file=mremd_rst_name,exist=file_exist)
103         IF (rest.and.file_exist.and..not.restart1file) THEN
104              open(irest2,file=mremd_rst_name,status='unknown')
105              read (irest2,*) 
106              read (irest2,*) (i2rep(i),i=0,nodes-1)
107              read (irest2,*) 
108              read (irest2,*) (ifirst(i),i=1,remd_m(1))
109              do il=1,nodes
110               read (irest2,*) 
111               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
112               read (irest2,*) 
113               read (irest2,*) ndowna(0,il),
114      &                    (ndowna(i,il),i=1,ndowna(0,il))
115              enddo
116              close(irest2)
117
118              write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
119              write (iout,'(a6,1000i5)') "ifirst",
120      &                    (ifirst(i),i=1,remd_m(1))
121              do il=1,nodes
122               write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
123      &                    (nupa(i,il),i=1,nupa(0,il))
124               write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
125      &                    (ndowna(i,il),i=1,ndowna(0,il))
126              enddo
127         ELSE
128          do il=1,remd_m(1)
129           ifirst(il)=il
130          enddo
131
132          do il=1,nodes
133
134           if (i2rep(il-1).eq.nrep) then
135            nupa(0,il)=0
136           else
137            nupa(0,il)=remd_m(i2rep(il-1)+1)
138            k=rep2i(i2rep(il-1))+1
139            do i=1,nupa(0,il)
140             nupa(i,il)=k+1
141             k=k+1
142            enddo
143           endif
144
145           if (i2rep(il-1).eq.1) then
146            ndowna(0,il)=0
147           else
148            ndowna(0,il)=remd_m(i2rep(il-1)-1)
149            k=rep2i(i2rep(il-1)-2)+1
150            do i=1,ndowna(0,il)
151             ndowna(i,il)=k+1
152             k=k+1
153            enddo
154           endif
155
156          enddo
157         
158 c        write (iout,'(a6,100i4)') "ifirst",
159 c     &                    (ifirst(i),i=1,remd_m(1))
160 c        do il=1,nodes
161 c         write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
162 c     &                    (nupa(i,il),i=1,nupa(0,il))
163 c         write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
164 c     &                    (ndowna(i,il),i=1,ndowna(0,il))
165 c        enddo
166         
167         ENDIF
168        endif
169 c
170 c      t_bath=retmin+(retmax-retmin)*me/(nodes-1)
171        if(.not.(rest.and.file_exist.and.restart1file)) then
172          if (me .eq. king) then 
173             t_bath=retmin
174          else 
175             t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
176          endif
177 cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
178          if (remd_tlist) t_bath=remd_t(i2rep(me))
179        endif
180 c
181         stdfp=dsqrt(2*Rb*t_bath/d_time)
182         do i=1,ntyp
183           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
184         enddo 
185
186 c      print *,'irep',me,t_bath
187        if (.not.rest) then  
188         if (me.eq.king .or. .not. out1file)
189      &   write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
190         call rescale_weights(t_bath)
191        endif
192
193
194 c------copy MD--------------
195 c  The driver for molecular dynamics subroutines
196 c------------------------------------------------
197       t_MDsetup=0.0d0
198       t_langsetup=0.0d0
199       t_MD=0.0d0
200       t_enegrad=0.0d0
201       t_sdsetup=0.0d0
202       if(me.eq.king.or..not.out1file)
203      & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
204       tt0 = tcpu()
205 c Determine the inverse of the inertia matrix.
206       call setup_MD_matrices
207 c Initialize MD
208       call init_MD
209       if (rest) then  
210        if (me.eq.king .or. .not. out1file)
211      &  write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
212        stdfp=dsqrt(2*Rb*t_bath/d_time)
213        do i=1,ntyp
214           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
215        enddo 
216        call rescale_weights(t_bath)
217       endif
218
219       t_MDsetup = tcpu()-tt0
220       rstcount=0 
221 c   Entering the MD loop       
222       tt0 = tcpu()
223       if (lang.eq.2 .or. lang.eq.3) then
224         call setup_fricmat
225         if (lang.eq.2) then
226           call sd_verlet_p_setup        
227         else
228           call sd_verlet_ciccotti_setup
229         endif
230 #ifndef   LANG0
231         do i=1,dimen
232           do j=1,dimen
233             pfric0_mat(i,j,0)=pfric_mat(i,j)
234             afric0_mat(i,j,0)=afric_mat(i,j)
235             vfric0_mat(i,j,0)=vfric_mat(i,j)
236             prand0_mat(i,j,0)=prand_mat(i,j)
237             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
238             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
239           enddo
240         enddo
241 #endif
242         flag_stoch(0)=.true.
243         do i=1,maxflag_stoch
244           flag_stoch(i)=.false.
245         enddo  
246       else if (lang.eq.1 .or. lang.eq.4) then
247         call setup_fricmat
248       endif
249       time00=MPI_WTIME()
250       if (me.eq.king .or. .not. out1file)
251      & write(iout,*) 'Setup time',time00-walltime
252       call flush(iout)
253       t_langsetup=tcpu()-tt0
254       tt0=tcpu()
255       itime=0
256       end_of_run=.false.
257       do while(.not.end_of_run)
258         itime=itime+1
259         if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
260         if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
261         rstcount=rstcount+1
262         if (lang.gt.0 .and. surfarea .and. 
263      &      mod(itime,reset_fricmat).eq.0) then
264           if (lang.eq.2 .or. lang.eq.3) then
265             call setup_fricmat
266             if (lang.eq.2) then
267               call sd_verlet_p_setup
268             else
269               call sd_verlet_ciccotti_setup
270             endif
271 #ifndef   LANG0
272             do i=1,dimen
273               do j=1,dimen
274                 pfric0_mat(i,j,0)=pfric_mat(i,j)
275                 afric0_mat(i,j,0)=afric_mat(i,j)
276                 vfric0_mat(i,j,0)=vfric_mat(i,j)
277                 prand0_mat(i,j,0)=prand_mat(i,j)
278                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
279                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
280               enddo
281             enddo
282 #endif
283             flag_stoch(0)=.true.
284             do i=1,maxflag_stoch
285               flag_stoch(i)=.false.
286             enddo   
287           else if (lang.eq.1 .or. lang.eq.4) then
288             call setup_fricmat
289           endif
290           write (iout,'(a,i10)') 
291      &      "Friction matrix reset based on surface area, itime",itime
292         endif
293         if (reset_vel .and. tbf .and. lang.eq.0 
294      &      .and. mod(itime,count_reset_vel).eq.0) then
295           call random_vel
296           if (me.eq.king .or. .not. out1file)
297      &     write(iout,'(a,f20.2)') 
298      &     "Velocities reset to random values, time",totT       
299           do i=0,2*nres
300             do j=1,3
301               d_t_old(j,i)=d_t(j,i)
302             enddo
303           enddo
304         endif
305         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
306           call inertia_tensor  
307           call vcm_vel(vcm)
308           do j=1,3
309              d_t(j,0)=d_t(j,0)-vcm(j)
310           enddo
311           call kinetic(EK)
312           kinetic_T=2.0d0/(dimen*Rb)*EK
313           scalfac=dsqrt(T_bath/kinetic_T)
314 cd          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT     
315           do i=0,2*nres
316             do j=1,3
317               d_t_old(j,i)=scalfac*d_t(j,i)
318             enddo
319           enddo
320         endif  
321         if (lang.ne.4) then
322           if (RESPA) then
323 c Time-reversible RESPA algorithm 
324 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
325             call RESPA_step(itime)
326           else
327 c Variable time step algorithm.
328             call velverlet_step(itime)
329           endif
330         else
331           call brown_step(itime)
332         endif
333         if(ntwe.ne.0) then
334           if (mod(itime,ntwe).eq.0) call statout(itime)
335         endif
336         if (mod(itime,ntwx).eq.0.and..not.traj1file) then
337           write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
338           if(mdpdb) then
339              call hairpin(.true.,nharp,iharp)
340              call secondary2(.true.)
341              call pdbout(potE,tytul,ipdb)
342           else 
343              call cartout(totT)
344           endif
345         endif
346         if ((rstcount.eq.1000.or.itime.eq.n_timestep)
347      &                         .and..not.restart1file) then
348
349            if(me.eq.king) then
350              open(irest1,file=mremd_rst_name,status='unknown')
351              write (irest1,*) "i2rep"
352              write (irest1,*) (i2rep(i),i=0,nodes-1)
353              write (irest1,*) "ifirst"
354              write (irest1,*) (ifirst(i),i=1,remd_m(1))
355              do il=1,nodes
356               write (irest1,*) "nupa",il
357               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
358               write (irest1,*) "ndowna",il
359               write (irest1,*) ndowna(0,il),
360      &                   (ndowna(i,il),i=1,ndowna(0,il))
361              enddo
362              close(irest1)
363            endif
364            open(irest2,file=rest2name,status='unknown')
365            write(irest2,*) totT,EK,potE,totE,t_bath
366            do i=1,2*nres
367             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
368            enddo
369            do i=1,2*nres
370             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
371            enddo
372           close(irest2)
373           rstcount=0
374         endif 
375
376 c REMD - exchange
377 c forced synchronization
378         if (mod(itime,100).eq.0 .and. me.ne.king 
379      &                                .and. .not. mremdsync) then 
380             synflag=.false.
381             call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
382             if (synflag) then 
383                call mpi_recv(itime_master, 1, MPI_INTEGER,
384      &                             0,101,CG_COMM, status, ierr)
385                if (.not. out1file) then
386                 write(iout,*) 'REMD synchro at',itime_master,itime
387                else
388                 call mpi_gather(itime,1,mpi_integer,
389      &             itime_all,1,mpi_integer,king,
390      &             CG_COMM,ierr)
391                endif
392                if(itime_master.ge.n_timestep) end_of_run=.true.
393                call flush(iout)
394             endif
395         endif
396
397 c REMD - exchange
398         if ((mod(itime,nstex).eq.0.and.me.eq.king
399      &                  .or.end_of_run.and.me.eq.king )
400      &       .and. .not. mremdsync ) then
401            synflag=.true.
402            time00=MPI_WTIME()
403            do i=1,nodes-1
404             call mpi_isend(itime,1,MPI_INTEGER,i,101,
405      &                                   CG_COMM, ireq, ierr)
406 cd            write(iout,*) 'REMD synchro with',i
407 cd            call flush(iout)
408            enddo
409            time02=MPI_WTIME()
410            write(iout,*) 'REMD synchro at',itime,'time=',time02-time00
411            if (out1file) then
412             call mpi_gather(itime,1,mpi_integer,
413      &             itime_all,1,mpi_integer,king,
414      &             CG_COMM,ierr)
415             write(iout,'(a18,8000i8)') 'REMD synchro itime',
416      &                    (itime_all(i),i=1,nodes)
417            endif
418            call flush(iout)
419         endif
420         if(mremdsync .and. mod(itime,nstex).eq.0) then
421            synflag=.true.
422            if (me.eq.king .or. .not. out1file)
423      &      write(iout,*) 'REMD synchro at',itime
424            call flush(iout)
425         endif
426         if(synflag.and..not.end_of_run) then
427            time00=MPI_WTIME()
428            synflag=.false.
429
430 cd           write(iout,*) 'REMD before',me,t_bath
431
432 c           call mpi_gather(t_bath,1,mpi_double_precision,
433 c     &             remd_t_bath,1,mpi_double_precision,king,
434 c     &             CG_COMM,ierr)
435            potEcomp(n_ene+1)=t_bath
436            call mpi_gather(potEcomp(0),n_ene+2,mpi_double_precision,
437      &             remd_ene(0,1),n_ene+2,mpi_double_precision,king,
438      &             CG_COMM,ierr)
439            if(lmuca) then 
440             call mpi_gather(elow,1,mpi_double_precision,
441      &             elowi,1,mpi_double_precision,king,
442      &             CG_COMM,ierr)
443             call mpi_gather(ehigh,1,mpi_double_precision,
444      &             ehighi,1,mpi_double_precision,king,
445      &             CG_COMM,ierr)
446            endif
447
448           if (restart1file) call write1rst
449           if (traj1file) call write1traj
450
451           if (me.eq.king) then
452             do i=1,nodes
453                remd_t_bath(i)=remd_ene(n_ene+1,i)
454             enddo
455             if(lmuca) then
456 co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
457              do i=1,nodes
458                write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
459      &            elowi(i),ehighi(i)       
460              enddo
461             else
462 cd              write(iout,*) 'REMD exchange temp,ene'
463 c             do i=1,nodes
464 co              write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
465 c              write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
466 c             enddo
467             endif
468 c-------------------------------------           
469            do irr=1,remd_m(1)
470            i=ifirst(iran_num(1,remd_m(1)))
471            do ii=1,nodes-1
472
473             if(i.gt.0.and.nupa(0,i).gt.0) then
474               iex=nupa(iran_num(1,int(nupa(0,i))),i)
475              if (lmuca) then
476                call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
477              else
478 c Swap temperatures between conformations i and iex with recalculating the free energies
479 c following temperature changes.
480               ene_iex_iex=remd_ene(0,iex)
481               ene_i_i=remd_ene(0,i)
482 co              write (iout,*) "rescaling weights with temperature",
483 co     &          remd_t_bath(i)
484               call rescale_weights(remd_t_bath(i))
485               call sum_energy(remd_ene(0,iex),.false.)
486               ene_iex_i=remd_ene(0,iex)
487 cd              write (iout,*) "ene_iex_i",remd_ene(0,iex)
488               call sum_energy(remd_ene(0,i),.false.)
489 cd              write (iout,*) "ene_i_i",remd_ene(0,i)
490 c              write (iout,*) "rescaling weights with temperature",
491 c     &          remd_t_bath(iex)
492               if (real(ene_i_i).ne.real(remd_ene(0,i))) then
493                 write (iout,*) "ERROR: inconsistent energies:",i,
494      &            ene_i_i,remd_ene(0,i)
495               endif
496               call rescale_weights(remd_t_bath(iex))
497               call sum_energy(remd_ene(0,i),.false.)
498 c              write (iout,*) "ene_i_iex",remd_ene(0,i)
499               ene_i_iex=remd_ene(0,i)
500               call sum_energy(remd_ene(0,iex),.false.)
501               if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
502                 write (iout,*) "ERROR: inconsistent energies:",iex,
503      &            ene_iex_iex,remd_ene(0,iex)
504               endif
505 c              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
506 c              write (iout,*) "i",i," iex",iex
507 co              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
508 co     &           " ene_i_iex",ene_i_iex,
509 co     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
510               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
511      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
512               delta=-delta
513 c              write(iout,*) 'delta',delta
514 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
515 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
516 c     &              (remd_t_bath(i)*remd_t_bath(iex))
517              endif
518               if (delta .gt. 50.0d0) then
519                 delta=0.0d0
520               else
521 #ifdef OSF 
522                 if(isnan(delta))then
523                   delta=0.0d0
524                 else if (delta.lt.-50.0d0) then
525                   delta=dexp(50.0d0)
526                 else
527                 delta=dexp(-delta)
528               endif
529 #else
530                 delta=dexp(-delta)
531 #endif
532               endif
533               iremd_tot(i2rep(i-1))=iremd_tot(i2rep(i-1))+1
534               xxx=ran_number(0.0d0,1.0d0)
535 co              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
536               if (delta .gt. xxx) then
537                 tmp=remd_t_bath(i)       
538                 remd_t_bath(i)=remd_t_bath(iex)
539                 remd_t_bath(iex)=tmp
540                 remd_ene(0,i)=ene_i_iex
541                 remd_ene(0,iex)=ene_iex_i
542                 if(lmuca) then
543                   tmp=elowi(i)
544                   elowi(i)=elowi(iex)
545                   elowi(iex)=tmp  
546                   tmp=ehighi(i)
547                   ehighi(i)=ehighi(iex)
548                   ehighi(iex)=tmp  
549                 endif
550
551
552                 do k=0,nodes
553                   itmp=nupa(k,i)
554                   nupa(k,i)=nupa(k,iex)
555                   nupa(k,iex)=itmp
556                   itmp=ndowna(k,i)
557                   ndowna(k,i)=ndowna(k,iex)
558                   ndowna(k,iex)=itmp
559                 enddo
560                 do il=1,nodes
561                  if (ifirst(il).eq.i) ifirst(il)=iex
562                  do k=1,nupa(0,il)
563                   if (nupa(k,il).eq.i) then 
564                      nupa(k,il)=iex
565                   elseif (nupa(k,il).eq.iex) then 
566                      nupa(k,il)=i
567                   endif
568                  enddo
569                  do k=1,ndowna(0,il)
570                   if (ndowna(k,il).eq.i) then 
571                      ndowna(k,il)=iex
572                   elseif (ndowna(k,il).eq.iex) then 
573                      ndowna(k,il)=i
574                   endif
575                  enddo
576                 enddo
577
578                 iremd_acc(i2rep(i-1))=iremd_acc(i2rep(i-1))+1
579                 itmp=i2rep(i-1)
580                 i2rep(i-1)=i2rep(iex-1)
581                 i2rep(iex-1)=itmp
582
583 cd                write(iout,*) 'exchange',i,iex
584 cd              write (iout,'(a8,100i4)') "@ ifirst",
585 cd     &                    (ifirst(k),k=1,remd_m(1))
586 cd              do il=1,nodes
587 cd               write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
588 cd     &                    (nupa(k,il),k=1,nupa(0,il))
589 cd               write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
590 cd     &                    (ndowna(k,il),k=1,ndowna(0,il))
591 cd              enddo
592
593               else
594                remd_ene(0,iex)=ene_iex_iex
595                remd_ene(0,i)=ene_i_i
596                i=iex
597               endif 
598             endif
599            enddo
600            enddo
601 c-------------------------------------
602              do i=1,nrep
603               if(iremd_tot(i).ne.0)
604      &          write(iout,'(a3,i4,2f12.5)') 'ACC',i,remd_t(i)
605      &           ,iremd_acc(i)/(1.0*iremd_tot(i))
606              enddo
607
608 cd              write (iout,'(a6,100i4)') "ifirst",
609 cd     &                    (ifirst(i),i=1,remd_m(1))
610 cd              do il=1,nodes
611 cd               write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
612 cd     &                    (nupa(i,il),i=1,nupa(0,il))
613 cd               write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
614 cd     &                    (ndowna(i,il),i=1,ndowna(0,il))
615 cd              enddo
616           endif
617
618          call mpi_scatter(remd_t_bath,1,mpi_double_precision,
619      &           t_bath,1,mpi_double_precision,king,
620      &           CG_COMM,ierr) 
621          if(lmuca) then
622            call mpi_scatter(elowi,1,mpi_double_precision,
623      &           elow,1,mpi_double_precision,king,
624      &           CG_COMM,ierr) 
625            call mpi_scatter(ehighi,1,mpi_double_precision,
626      &           ehigh,1,mpi_double_precision,king,
627      &           CG_COMM,ierr) 
628          endif
629          call rescale_weights(t_bath)
630 co         write (iout,*) "Processor",me,
631 co     &    " rescaling weights with temperature",t_bath
632
633          stdfp=dsqrt(2*Rb*t_bath/d_time)
634          do i=1,ntyp
635            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
636          enddo 
637
638 cde         write(iout,*) 'REMD after',me,t_bath
639            time02=MPI_WTIME()
640            if (me.eq.king .or. .not. out1file) then
641             write(iout,*) 'REMD exchange time=',time02-time00
642             call flush(iout)
643            endif
644         endif
645       enddo
646
647       if (restart1file) then 
648           if (me.eq.king .or. .not. out1file)
649      &      write(iout,*) 'writing restart at the end of run'
650            call write1rst
651       endif
652
653       if (traj1file) call write1traj
654
655
656       t_MD=tcpu()-tt0
657       if (me.eq.king .or. .not. out1file) then
658        write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') 
659      &  '  Timing  ',
660      & 'MD calculations setup:',t_MDsetup,
661      & 'Energy & gradient evaluation:',t_enegrad,
662      & 'Stochastic MD setup:',t_langsetup,
663      & 'Stochastic MD step setup:',t_sdsetup,
664      & 'MD steps:',t_MD
665        write (iout,'(/28(1h=),a25,27(1h=))') 
666      & '  End of MD calculation  '
667       endif
668       return
669       end
670
671       subroutine write1rst
672       implicit real*8 (a-h,o-z)
673       include 'DIMENSIONS'
674       include 'mpif.h'
675       include 'COMMON.MD'
676       include 'COMMON.IOUNITS'
677       include 'COMMON.REMD'
678       include 'COMMON.SETUP'
679       include 'COMMON.CHAIN'
680       include 'COMMON.SBRIDGE'
681       include 'COMMON.INTERACT'
682                
683       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
684      &     d_restart2(3,2*maxres*maxprocs)
685       real t5_restart1(5)
686       integer iret,itmp
687
688
689        t5_restart1(1)=totT
690        t5_restart1(2)=EK
691        t5_restart1(3)=potE
692        t5_restart1(4)=t_bath
693        t5_restart1(5)=Uconst
694        
695        call mpi_gather(t5_restart1,5,mpi_real,
696      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
697
698
699        do i=1,2*nres
700          do j=1,3
701            r_d(j,i)=d_t(j,i)
702          enddo
703        enddo
704        call mpi_gather(r_d,3*2*nres,mpi_real,
705      &           d_restart1,3*2*nres,mpi_real,king,
706      &           CG_COMM,ierr)
707
708
709        do i=1,2*nres
710          do j=1,3
711            r_d(j,i)=dc(j,i)
712          enddo
713        enddo
714        call mpi_gather(r_d,3*2*nres,mpi_real,
715      &           d_restart2,3*2*nres,mpi_real,king,
716      &           CG_COMM,ierr)
717
718        if(me.eq.king) then
719          open(irest1,file=mremd_rst_name,status='unknown')
720          write (irest1,*) (i2rep(i),i=0,nodes-1)
721          write (irest1,*) (ifirst(i),i=1,remd_m(1))
722          do il=1,nodes
723               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
724               write (irest1,*) ndowna(0,il),
725      &                   (ndowna(i,il),i=1,ndowna(0,il))
726          enddo
727
728          do il=1,nodes
729            write(irest1,*) (t_restart1(j,il),j=1,4)
730          enddo
731
732          do il=0,nodes-1
733            do i=1,2*nres
734              write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
735            enddo
736          enddo
737          do il=0,nodes-1
738            do i=1,2*nres
739              write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
740            enddo
741          enddo
742          close(irest1)
743        endif
744
745
746       return
747       end
748
749       subroutine write1traj
750       implicit real*8 (a-h,o-z)
751       include 'DIMENSIONS'
752       include 'mpif.h'
753       include 'COMMON.MD'
754       include 'COMMON.IOUNITS'
755       include 'COMMON.REMD'
756       include 'COMMON.SETUP'
757       include 'COMMON.CHAIN'
758       include 'COMMON.SBRIDGE'
759       include 'COMMON.INTERACT'
760                
761       real t5_restart1(5)
762       integer iret,itmp
763       real xcoord(3,maxres2+2),prec
764       real r_qfrag(50),r_qpair(100)
765       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
766       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
767
768       if(.not.restart1file) then
769        t5_restart1(1)=totT
770        t5_restart1(2)=EK
771        t5_restart1(3)=potE
772        t5_restart1(4)=t_bath
773        t5_restart1(5)=Uconst
774        call mpi_gather(t5_restart1,5,mpi_real,
775      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
776       endif
777
778           do i=1,nfrag
779            r_qfrag(i)=qfrag(i)
780           enddo
781           do i=1,npair
782            r_qpair(i)=qpair(i)
783           enddo
784
785         call mpi_gather(r_qfrag,nfrag,mpi_real,
786      &           p_qfrag,nfrag,mpi_real,king,
787      &           CG_COMM,ierr)
788         call mpi_gather(qpair,nfrag,mpi_real,
789      &           p_qpair,nfrag,mpi_real,king,
790      &           CG_COMM,ierr)
791
792         do i=1,nres*2
793          do j=1,3
794           r_c(j,i)=c(j,i)
795          enddo
796         enddo
797
798         call mpi_gather(r_c,3*2*nres,mpi_real,
799      &           p_c,3*2*nres,mpi_real,king,
800      &           CG_COMM,ierr)
801
802        if(me.eq.king) then
803          call xdrfopen(ixdrf,cartname, "a", iret)
804          do il=1,nodes
805           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
806           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
807           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
808           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
809           call xdrfint(ixdrf, nss, iret) 
810           do j=1,nss
811            call xdrfint(ixdrf, ihpb(j), iret)
812            call xdrfint(ixdrf, jhpb(j), iret)
813           enddo
814           call xdrfint(ixdrf, nfrag+npair, iret)
815           do i=1,nfrag
816            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
817           enddo
818           do i=1,npair
819            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
820           enddo
821           prec=10000.0
822           do i=1,nres
823            do j=1,3
824             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
825            enddo
826           enddo
827           do i=nnt,nct
828            do j=1,3
829             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
830            enddo
831           enddo
832           itmp=nres+nct-nnt+1
833           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
834          enddo
835          call xdrfclose(ixdrf, iret)
836        endif
837  
838       return
839       end
840
841
842       subroutine read1restart
843       implicit real*8 (a-h,o-z)
844       include 'DIMENSIONS'
845       include 'mpif.h'
846       include 'COMMON.MD'
847       include 'COMMON.IOUNITS'
848       include 'COMMON.REMD'
849       include 'COMMON.SETUP'
850       include 'COMMON.CHAIN'
851       include 'COMMON.SBRIDGE'
852       include 'COMMON.INTERACT'
853       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
854      &                 t5_restart1(5)
855
856          if(me.eq.king)then
857              open(irest2,file=mremd_rst_name,status='unknown')
858              read (irest2,*) (i2rep(i),i=0,nodes-1)
859              read (irest2,*) (ifirst(i),i=1,remd_m(1))
860              do il=1,nodes
861               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
862               read (irest2,*) ndowna(0,il),
863      &                    (ndowna(i,il),i=1,ndowna(0,il))
864              enddo
865              do il=1,nodes
866                read(irest2,*) (t_restart1(j,il),j=1,4)
867              enddo
868          endif
869          call mpi_scatter(t_restart1,5,mpi_real,
870      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
871          totT=t5_restart1(1)              
872          EK=t5_restart1(2)
873          potE=t5_restart1(3)
874          t_bath=t5_restart1(4)
875
876          if(me.eq.king)then
877               do il=0,nodes-1
878                do i=1,2*nres
879                 read(irest2,'(3e15.5)') 
880      &                (d_restart1(j,i+2*nres*il),j=1,3)
881                enddo
882               enddo
883          endif
884          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
885      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
886
887          do i=1,2*nres
888            do j=1,3
889             d_t(j,i)=r_d(j,i)
890            enddo
891          enddo
892          if(me.eq.king)then 
893               do il=0,nodes-1
894                do i=1,2*nres
895                 read(irest2,'(3e15.5)') 
896      &                (d_restart1(j,i+2*nres*il),j=1,3)
897                enddo
898               enddo
899          endif
900          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
901      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
902          do i=1,2*nres
903            do j=1,3
904             dc(j,i)=r_d(j,i)
905            enddo
906          enddo
907         if(me.eq.king) close(irest2)
908         return
909         end
910