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