Merge branch 'devel' of mmka.chem.univ.gda.pl:unres into devel
[unres.git] / source / unres / src_MD-M / MREMD.F.drabinka
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 cold      integer nup(0:maxprocs),ndown(0:maxprocs)
35       integer rep2i(0:maxprocs),ireqi(maxprocs)
36       integer itime_all(maxprocs)
37       integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
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 cold       if (i2rep(me).eq.nrep) then
68 cold        nup(0)=0
69 cold       else
70 cold        nup(0)=remd_m(i2rep(me)+1)
71 cold        k=rep2i(int(i2rep(me)))+1
72 cold        do i=1,nup(0)
73 cold         nup(i)=k
74 cold         k=k+1
75 cold        enddo
76 cold       endif
77
78 cd       print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
79
80 cold       if (i2rep(me).eq.1) then
81 cold        ndown(0)=0
82 cold       else
83 cold        ndown(0)=remd_m(i2rep(me)-1)
84 cold        k=rep2i(i2rep(me)-2)+1
85 cold        do i=1,ndown(0)
86 cold         ndown(i)=k
87 cold         k=k+1
88 cold        enddo
89 cold       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 IF (.not.(rest.and.file_exist)) THEN
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(int(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(int(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 (mod(itime,ntwx).eq.0.and.traj1file) then
344           if(ntwx_cache.lt.max_cache_traj) then
345             ntwx_cache=ntwx_cache+1
346
347             totT_cache(ntwx_cache)=totT
348             EK_cache(ntwx_cache)=EK
349             potE_cache(ntwx_cache)=potE
350             t_bath_cache(ntwx_cache)=t_bath
351             Uconst_cache(ntwx_cache)=Uconst
352
353             do i=1,nfrag
354              qfrag_cache(i,ntwx_cache)=qfrag(i)
355             enddo
356             do i=1,npair
357              qpair_cache(i,ntwx_cache)=qpair(i)
358             enddo
359
360             do i=1,nres*2
361              do j=1,3
362               c_cache(j,i,ntwx_cache)=c(j,i)
363              enddo
364             enddo
365           else
366            print *,itime,"processor ",me," over cache ",ntwx_cache
367           endif
368         endif
369         if ((rstcount.eq.1000.or.itime.eq.n_timestep)
370      &                         .and..not.restart1file) then
371
372            if(me.eq.king) then
373              open(irest1,file=mremd_rst_name,status='unknown')
374              write (irest1,*) "i2rep"
375              write (irest1,*) (i2rep(i),i=0,nodes-1)
376              write (irest1,*) "ifirst"
377              write (irest1,*) (ifirst(i),i=1,remd_m(1))
378              do il=1,nodes
379               write (irest1,*) "nupa",il
380               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
381               write (irest1,*) "ndowna",il
382               write (irest1,*) ndowna(0,il),
383      &                   (ndowna(i,il),i=1,ndowna(0,il))
384              enddo
385              close(irest1)
386            endif
387            open(irest2,file=rest2name,status='unknown')
388            write(irest2,*) totT,EK,potE,totE,t_bath
389            do i=1,2*nres
390             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
391            enddo
392            do i=1,2*nres
393             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
394            enddo
395           close(irest2)
396           rstcount=0
397         endif 
398
399 c REMD - exchange
400 c forced synchronization
401         if (mod(itime,i_sync_step).eq.0 .and. me.ne.king 
402      &                                .and. .not. mremdsync) then 
403             synflag=.false.
404             call mpi_iprobe(0,101,mpi_comm_world,synflag,status,ierr)
405             if (synflag) then 
406                call mpi_recv(itime_master, 1, MPI_INTEGER,
407      &                             0,101,mpi_comm_world, status, ierr)
408                call mpi_barrier(mpi_comm_world, ierr)
409                if (out1file.or.traj1file) then
410                 call mpi_gather(itime,1,mpi_integer,
411      &             itime_all,1,mpi_integer,king,
412      &             mpi_comm_world,ierr)                 
413                endif
414                if (.not.out1file)
415      &               write(iout,*) 'REMD synchro at',itime_master,itime
416                if(itime_master.ge.n_timestep) end_of_run=.true.
417 ctime               call flush(iout)
418             endif
419         endif
420
421 c REMD - exchange
422         if ((mod(itime,nstex).eq.0.and.me.eq.king
423      &                  .or.end_of_run.and.me.eq.king )
424      &       .and. .not. mremdsync ) then
425            synflag=.true.
426            time00=MPI_WTIME()
427            do i=1,nodes-1
428               call mpi_isend(itime,1,MPI_INTEGER,i,101,
429      &                                mpi_comm_world, ireqi(i), ierr)
430 cd            write(iout,*) 'REMD synchro with',i
431 cd            call flush(iout)
432            enddo
433            call mpi_waitall(nodes-1,ireqi,statusi,ierr)
434            call mpi_barrier(mpi_comm_world, ierr)
435            time01=MPI_WTIME()
436            write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
437            if (out1file.or.traj1file) then
438             call mpi_gather(itime,1,mpi_integer,
439      &             itime_all,1,mpi_integer,king,
440      &             mpi_comm_world,ierr)
441 ctime            write(iout,'(a19,8000i8)') ' REMD synchro itime',
442 ctime     &                    (itime_all(i),i=1,nodes)
443             if(traj1file) then
444              imin_itime=itime_all(1)
445              do i=2,nodes
446                if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
447              enddo
448              ii_write=(imin_itime-imin_itime_old)/ntwx
449              imin_itime_old=int(imin_itime/ntwx)*ntwx
450              write(iout,*) imin_itime,imin_itime_old,ii_write
451             endif
452            endif
453 ctime           call flush(iout)
454         endif
455         if(mremdsync .and. mod(itime,nstex).eq.0) then
456            synflag=.true.
457            if (me.eq.king .or. .not. out1file)
458      &      write(iout,*) 'REMD synchro at',itime
459            call flush(iout)
460         endif
461         if(synflag.and..not.end_of_run) then
462            time02=MPI_WTIME()
463            synflag=.false.
464
465 cd           write(iout,*) 'REMD before',me,t_bath
466
467 c           call mpi_gather(t_bath,1,mpi_double_precision,
468 c     &             remd_t_bath,1,mpi_double_precision,king,
469 c     &             mpi_comm_world,ierr)
470            potEcomp(n_ene+1)=t_bath
471            potEcomp(n_ene+2)=iset
472            call mpi_gather(potEcomp(0),n_ene+3,mpi_double_precision,
473      &             remd_ene(0,1),n_ene+3,mpi_double_precision,king,
474      &             mpi_comm_world,ierr)
475            if(lmuca) then 
476             call mpi_gather(elow,1,mpi_double_precision,
477      &             elowi,1,mpi_double_precision,king,
478      &             mpi_comm_world,ierr)
479             call mpi_gather(ehigh,1,mpi_double_precision,
480      &             ehighi,1,mpi_double_precision,king,
481      &             mpi_comm_world,ierr)
482            endif
483
484           time03=MPI_WTIME()
485           if (me.eq.king .or. .not. out1file) then
486             write(iout,*) 'REMD gather times=',time03-time01
487      &                                        ,time03-time02
488           endif
489
490           if (restart1file) call write1rst
491
492           time04=MPI_WTIME()
493           if (me.eq.king .or. .not. out1file) then
494             write(iout,*) 'REMD writing rst time=',time04-time03
495           endif
496
497           if (traj1file) call write1traj
498
499           time05=MPI_WTIME()
500           if (me.eq.king .or. .not. out1file) then
501             write(iout,*) 'REMD writing traj time=',time05-time04
502           endif
503
504
505           if (me.eq.king) then
506             do i=1,nodes
507                remd_t_bath(i)=remd_ene(n_ene+1,i)
508             enddo
509             if(lmuca) then
510 co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
511              do i=1,nodes
512                write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
513      &            elowi(i),ehighi(i)       
514              enddo
515             else
516 cd              write(iout,*) 'REMD exchange temp,ene'
517 c             do i=1,nodes
518 co              write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
519 c              write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
520 c             enddo
521             endif
522 c-------------------------------------           
523            do irr=1,remd_m(1)
524            i=ifirst(iran_num(1,remd_m(1)))
525            do ii=1,nodes-1
526
527             if(i.gt.0.and.nupa(0,i).gt.0) then
528               iex=nupa(iran_num(1,int(nupa(0,i))),i)
529              if (lmuca) then
530                call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
531              else
532 c Swap temperatures between conformations i and iex with recalculating the free energies
533 c following temperature changes.
534               ene_iex_iex=remd_ene(0,iex)
535               ene_i_i=remd_ene(0,i)
536 co              write (iout,*) "rescaling weights with temperature",
537 co     &          remd_t_bath(i)
538               call rescale_weights(remd_t_bath(i))
539               call sum_energy(remd_ene(0,iex))
540               ene_iex_i=remd_ene(0,iex)
541 cd              write (iout,*) "ene_iex_i",remd_ene(0,iex)
542               call sum_energy(remd_ene(0,i))
543 cd              write (iout,*) "ene_i_i",remd_ene(0,i)
544 c              write (iout,*) "rescaling weights with temperature",
545 c     &          remd_t_bath(iex)
546               if (real(ene_i_i).ne.real(remd_ene(0,i))) then
547                 write (iout,*) "ERROR: inconsistent energies:",i,
548      &            ene_i_i,remd_ene(0,i)
549               endif
550               call rescale_weights(remd_t_bath(iex))
551               call sum_energy(remd_ene(0,i))
552 c              write (iout,*) "ene_i_iex",remd_ene(0,i)
553               ene_i_iex=remd_ene(0,i)
554               call sum_energy(remd_ene(0,iex))
555               if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
556                 write (iout,*) "ERROR: inconsistent energies:",iex,
557      &            ene_iex_iex,remd_ene(0,iex)
558               endif
559 c              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
560 c              write (iout,*) "i",i," iex",iex
561 co              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
562 co     &           " ene_i_iex",ene_i_iex,
563 co     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
564               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
565      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
566               delta=-delta
567 c              write(iout,*) 'delta',delta
568 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
569 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
570 c     &              (remd_t_bath(i)*remd_t_bath(iex))
571              endif
572               if (delta .gt. 50.0d0) then
573                 delta=0.0d0
574               else
575 #ifdef OSF 
576                 if(isnan(delta))then
577                   delta=0.0d0
578                 else if (delta.lt.-50.0d0) then
579                   delta=dexp(50.0d0)
580                 else
581                 delta=dexp(-delta)
582               endif
583 #else
584                 delta=dexp(-delta)
585 #endif
586               endif
587               iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
588               xxx=ran_number(0.0d0,1.0d0)
589 co              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
590               if (delta .gt. xxx) then
591                 tmp=remd_t_bath(i)       
592                 remd_t_bath(i)=remd_t_bath(iex)
593                 remd_t_bath(iex)=tmp
594                 remd_ene(0,i)=ene_i_iex
595                 remd_ene(0,iex)=ene_iex_i
596                 if(lmuca) then
597                   tmp=elowi(i)
598                   elowi(i)=elowi(iex)
599                   elowi(iex)=tmp  
600                   tmp=ehighi(i)
601                   ehighi(i)=ehighi(iex)
602                   ehighi(iex)=tmp  
603                 endif
604
605
606                 do k=0,nodes
607                   itmp=nupa(k,i)
608                   nupa(k,i)=nupa(k,iex)
609                   nupa(k,iex)=itmp
610                   itmp=ndowna(k,i)
611                   ndowna(k,i)=ndowna(k,iex)
612                   ndowna(k,iex)=itmp
613                 enddo
614                 do il=1,nodes
615                  if (ifirst(il).eq.i) ifirst(il)=iex
616                  do k=1,nupa(0,il)
617                   if (nupa(k,il).eq.i) then 
618                      nupa(k,il)=iex
619                   elseif (nupa(k,il).eq.iex) then 
620                      nupa(k,il)=i
621                   endif
622                  enddo
623                  do k=1,ndowna(0,il)
624                   if (ndowna(k,il).eq.i) then 
625                      ndowna(k,il)=iex
626                   elseif (ndowna(k,il).eq.iex) then 
627                      ndowna(k,il)=i
628                   endif
629                  enddo
630                 enddo
631
632                 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
633                 itmp=i2rep(i-1)
634                 i2rep(i-1)=i2rep(iex-1)
635                 i2rep(iex-1)=itmp
636
637 cd                write(iout,*) 'exchange',i,iex
638 cd              write (iout,'(a8,100i4)') "@ ifirst",
639 cd     &                    (ifirst(k),k=1,remd_m(1))
640 cd              do il=1,nodes
641 cd               write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
642 cd     &                    (nupa(k,il),k=1,nupa(0,il))
643 cd               write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
644 cd     &                    (ndowna(k,il),k=1,ndowna(0,il))
645 cd              enddo
646
647               else
648                remd_ene(0,iex)=ene_iex_iex
649                remd_ene(0,i)=ene_i_i
650                i=iex
651               endif 
652             endif
653            enddo
654            enddo
655 c-------------------------------------
656              do i=1,nrep
657               if(iremd_tot(i).ne.0)
658      &          write(iout,'(a3,i4,2f12.5)') 'ACC',i,remd_t(i)
659      &           ,iremd_acc(i)/(1.0*iremd_tot(i))
660              enddo
661
662 cd              write (iout,'(a6,100i4)') "ifirst",
663 cd     &                    (ifirst(i),i=1,remd_m(1))
664 cd              do il=1,nodes
665 cd               write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
666 cd     &                    (nupa(i,il),i=1,nupa(0,il))
667 cd               write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
668 cd     &                    (ndowna(i,il),i=1,ndowna(0,il))
669 cd              enddo
670           endif
671
672          time06=MPI_WTIME()
673          call mpi_scatter(remd_t_bath,1,mpi_double_precision,
674      &           t_bath,1,mpi_double_precision,king,
675      &           mpi_comm_world,ierr) 
676          time07=MPI_WTIME()
677           if (me.eq.king .or. .not. out1file) then
678             write(iout,*) 'REMD scatter time=',time07-time06
679           endif
680
681          if(lmuca) then
682            call mpi_scatter(elowi,1,mpi_double_precision,
683      &           elow,1,mpi_double_precision,king,
684      &           mpi_comm_world,ierr) 
685            call mpi_scatter(ehighi,1,mpi_double_precision,
686      &           ehigh,1,mpi_double_precision,king,
687      &           mpi_comm_world,ierr) 
688          endif
689          call rescale_weights(t_bath)
690 co         write (iout,*) "Processor",me,
691 co     &    " rescaling weights with temperature",t_bath
692
693          stdfp=dsqrt(2*Rb*t_bath/d_time)
694          do i=1,ntyp
695            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
696          enddo 
697
698 cde         write(iout,*) 'REMD after',me,t_bath
699            time08=MPI_WTIME()
700            if (me.eq.king .or. .not. out1file) then
701             write(iout,*) 'REMD exchange time=',time08-time00
702             call flush(iout)
703            endif
704         endif
705       enddo
706
707       if (restart1file) then 
708           if (me.eq.king .or. .not. out1file)
709      &      write(iout,*) 'writing restart at the end of run'
710            call write1rst
711       endif
712
713       if (traj1file) call write1traj
714
715
716       t_MD=tcpu()-tt0
717       if (me.eq.king .or. .not. out1file) then
718        write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') 
719      &  '  Timing  ',
720      & 'MD calculations setup:',t_MDsetup,
721      & 'Energy & gradient evaluation:',t_enegrad,
722      & 'Stochastic MD setup:',t_langsetup,
723      & 'Stochastic MD step setup:',t_sdsetup,
724      & 'MD steps:',t_MD
725        write (iout,'(/28(1h=),a25,27(1h=))') 
726      & '  End of MD calculation  '
727       endif
728       return
729       end
730
731       subroutine write1rst_
732       implicit real*8 (a-h,o-z)
733       include 'DIMENSIONS'
734       include 'mpif.h'
735       include 'COMMON.MD'
736       include 'COMMON.IOUNITS'
737       include 'COMMON.REMD'
738       include 'COMMON.SETUP'
739       include 'COMMON.CHAIN'
740       include 'COMMON.SBRIDGE'
741       include 'COMMON.INTERACT'
742                
743       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
744      &     d_restart2(3,2*maxres*maxprocs)
745       real t5_restart1(5)
746       integer iret,itmp
747
748
749        t5_restart1(1)=totT
750        t5_restart1(2)=EK
751        t5_restart1(3)=potE
752        t5_restart1(4)=t_bath
753        t5_restart1(5)=Uconst
754        
755        call mpi_gather(t5_restart1,5,mpi_real,
756      &      t_restart1,5,mpi_real,king,mpi_comm_world,ierr)
757
758
759        do i=1,2*nres
760          do j=1,3
761            r_d(j,i)=d_t(j,i)
762          enddo
763        enddo
764        call mpi_gather(r_d,3*2*nres,mpi_real,
765      &           d_restart1,3*2*nres,mpi_real,king,
766      &           mpi_comm_world,ierr)
767
768
769        do i=1,2*nres
770          do j=1,3
771            r_d(j,i)=dc(j,i)
772          enddo
773        enddo
774        call mpi_gather(r_d,3*2*nres,mpi_real,
775      &           d_restart2,3*2*nres,mpi_real,king,
776      &           mpi_comm_world,ierr)
777
778        if(me.eq.king) then
779          open(irest1,file=mremd_rst_name,status='unknown')
780          write (irest1,*) (i2rep(i),i=0,nodes-1)
781          write (irest1,*) (ifirst(i),i=1,remd_m(1))
782          do il=1,nodes
783               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
784               write (irest1,*) ndowna(0,il),
785      &                   (ndowna(i,il),i=1,ndowna(0,il))
786          enddo
787
788          do il=1,nodes
789            write(irest1,*) (t_restart1(j,il),j=1,4)
790          enddo
791
792          do il=0,nodes-1
793            do i=1,2*nres
794              write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
795            enddo
796          enddo
797          do il=0,nodes-1
798            do i=1,2*nres
799              write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
800            enddo
801          enddo
802          close(irest1)
803        endif
804
805       return
806       end
807
808       subroutine write1rst
809       implicit real*8 (a-h,o-z)
810       include 'DIMENSIONS'
811       include 'mpif.h'
812       include 'COMMON.MD'
813       include 'COMMON.IOUNITS'
814       include 'COMMON.REMD'
815       include 'COMMON.SETUP'
816       include 'COMMON.CHAIN'
817       include 'COMMON.SBRIDGE'
818       include 'COMMON.INTERACT'
819                
820       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
821      &     d_restart2(3,2*maxres*maxprocs)
822       real t5_restart1(5)
823       integer iret,itmp
824
825
826        t5_restart1(1)=totT
827        t5_restart1(2)=EK
828        t5_restart1(3)=potE
829        t5_restart1(4)=t_bath
830        t5_restart1(5)=Uconst
831        
832        call mpi_gather(t5_restart1,5,mpi_real,
833      &      t_restart1,5,mpi_real,king,mpi_comm_world,ierr)
834
835
836        do i=1,2*nres
837          do j=1,3
838            r_d(j,i)=d_t(j,i)
839          enddo
840        enddo
841        call mpi_gather(r_d,3*2*nres,mpi_real,
842      &           d_restart1,3*2*nres,mpi_real,king,
843      &           mpi_comm_world,ierr)
844
845
846        do i=1,2*nres
847          do j=1,3
848            r_d(j,i)=dc(j,i)
849          enddo
850        enddo
851        call mpi_gather(r_d,3*2*nres,mpi_real,
852      &           d_restart2,3*2*nres,mpi_real,king,
853      &           mpi_comm_world,ierr)
854
855        if(me.eq.king) then
856 c         open(irest1,file=mremd_rst_name,status='unknown')
857          call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
858 c         write (irest1,*) (i2rep(i),i=0,nodes-1)
859          do i=0,nodes-1
860           call xdrfint(ixdrf, i2rep(i), iret)
861          enddo
862 c         write (irest1,*) (ifirst(i),i=1,remd_m(1))
863          do i=1,remd_m(1)
864           call xdrfint(ixdrf, ifirst(i), iret)
865          enddo
866          do il=1,nodes
867 c              write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
868               do i=0,nupa(0,il)
869                call xdrfint(ixdrf, nupa(i,il), iret)
870               enddo
871
872 c              write (irest1,*) ndowna(0,il),
873 c     &                   (ndowna(i,il),i=1,ndowna(0,il))
874               do i=0,ndowna(0,il)
875                call xdrfint(ixdrf, ndowna(i,il), iret)
876               enddo
877          enddo
878
879          do il=1,nodes
880 c           write(irest1,*) (t_restart1(j,il),j=1,4)
881            do j=1,4
882             call xdrffloat(ixdrf, t_restart1(j,il), iret)
883            enddo
884          enddo
885
886          do il=0,nodes-1
887            do i=1,2*nres
888 c             write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
889             do j=1,3
890              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
891             enddo
892            enddo
893          enddo
894          do il=0,nodes-1
895            do i=1,2*nres
896 c             write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
897             do j=1,3
898              call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
899             enddo
900            enddo
901          enddo
902 c         close(irest1)
903          call xdrfclose(ixdrf, iret)
904        endif
905
906       return
907       end
908
909
910       subroutine write1traj
911       implicit real*8 (a-h,o-z)
912       include 'DIMENSIONS'
913       include 'mpif.h'
914       include 'COMMON.MD'
915       include 'COMMON.IOUNITS'
916       include 'COMMON.REMD'
917       include 'COMMON.SETUP'
918       include 'COMMON.CHAIN'
919       include 'COMMON.SBRIDGE'
920       include 'COMMON.INTERACT'
921                
922       real t5_restart1(5)
923       integer iret,itmp
924       real xcoord(3,maxres2+2),prec
925       real r_qfrag(50),r_qpair(100)
926       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
927       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
928
929       call mpi_bcast(ii_write,1,mpi_integer,
930      &           king,mpi_comm_world,ierr)
931
932       print *,'traj1file',me,ii_write,ntwx_cache
933
934       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
935       do ii=1,ii_write
936        t5_restart1(1)=totT_cache(ii)
937        t5_restart1(2)=EK_cache(ii)
938        t5_restart1(3)=potE_cache(ii)
939        t5_restart1(4)=t_bath_cache(ii)
940        t5_restart1(5)=Uconst_cache(ii)
941        call mpi_gather(t5_restart1,5,mpi_real,
942      &      t_restart1,5,mpi_real,king,mpi_comm_world,ierr)
943
944           do i=1,nfrag
945            r_qfrag(i)=qfrag_cache(i,ii)
946           enddo
947           do i=1,npair
948            r_qpair(i)=qpair_cache(i,ii)
949           enddo
950
951         call mpi_gather(r_qfrag,nfrag,mpi_real,
952      &           p_qfrag,nfrag,mpi_real,king,
953      &           mpi_comm_world,ierr)
954         call mpi_gather(r_qpair,npair,mpi_real,
955      &           p_qpair,npair,mpi_real,king,
956      &           mpi_comm_world,ierr)
957
958         do i=1,nres*2
959          do j=1,3
960           r_c(j,i)=c_cache(j,i,ii)
961          enddo
962         enddo
963
964         call mpi_gather(r_c,3*2*nres,mpi_real,
965      &           p_c,3*2*nres,mpi_real,king,
966      &           mpi_comm_world,ierr)
967
968        if(me.eq.king) then
969
970          do il=1,nodes
971           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
972           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
973           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
974           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
975           call xdrfint(ixdrf, nss, iret) 
976           do j=1,nss
977            call xdrfint(ixdrf, ihpb(j), iret)
978            call xdrfint(ixdrf, jhpb(j), iret)
979           enddo
980           call xdrfint(ixdrf, nfrag+npair, iret)
981           do i=1,nfrag
982            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
983           enddo
984           do i=1,npair
985            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
986           enddo
987           prec=10000.0
988           do i=1,nres
989            do j=1,3
990             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
991            enddo
992           enddo
993           do i=nnt,nct
994            do j=1,3
995             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
996            enddo
997           enddo
998           itmp=nres+nct-nnt+1
999           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1000          enddo
1001        endif
1002       enddo
1003       if(me.eq.king) call xdrfclose(ixdrf, iret)
1004       do i=1,ntwx_cache-ii_write
1005
1006             totT_cache(i)=totT_cache(ii_write+i)
1007             EK_cache(i)=EK_cache(ii_write+i)
1008             potE_cache(i)=potE_cache(ii_write+i)
1009             t_bath_cache(i)=t_bath_cache(ii_write+i)
1010             Uconst_cache(i)=Uconst_cache(ii_write+i)
1011
1012             do ii=1,nfrag
1013              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1014             enddo
1015             do ii=1,npair
1016              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1017             enddo
1018
1019             do ii=1,nres*2
1020              do j=1,3
1021               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1022              enddo
1023             enddo
1024       enddo
1025       ntwx_cache=ntwx_cache-ii_write
1026       return
1027       end
1028
1029
1030       subroutine read1restart
1031       implicit real*8 (a-h,o-z)
1032       include 'DIMENSIONS'
1033       include 'mpif.h'
1034       include 'COMMON.MD'
1035       include 'COMMON.IOUNITS'
1036       include 'COMMON.REMD'
1037       include 'COMMON.SETUP'
1038       include 'COMMON.CHAIN'
1039       include 'COMMON.SBRIDGE'
1040       include 'COMMON.INTERACT'
1041       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1042      &                 t5_restart1(5)
1043
1044          if(me.eq.king)then
1045               open(irest2,file=mremd_rst_name,status='unknown')
1046               read(irest2,*,err=334) i
1047               write(iout,*) "Reading old rst in ASCI format"
1048               close(irest2)
1049                call read1restart_old
1050                return
1051  334          continue
1052               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1053
1054 c             read (irest2,*) (i2rep(i),i=0,nodes-1)
1055               do i=0,nodes-1
1056                call xdrfint(ixdrf, i2rep(i), iret)
1057               enddo
1058 c             read (irest2,*) (ifirst(i),i=1,remd_m(1))
1059               do i=1,remd_m(1)
1060                call xdrfint(ixdrf, ifirst(i), iret)
1061               enddo
1062              do il=1,nodes
1063 c              read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1064               call xdrfint(ixdrf, nupa(0,il), iret)
1065               do i=1,nupa(0,il)
1066                call xdrfint(ixdrf, nupa(i,il), iret)
1067               enddo
1068
1069 c              read (irest2,*) ndowna(0,il),
1070 c     &                    (ndowna(i,il),i=1,ndowna(0,il))
1071               call xdrfint(ixdrf, ndowna(0,il), iret)
1072               do i=1,ndowna(0,il)
1073                call xdrfint(ixdrf, ndowna(i,il), iret)
1074               enddo
1075              enddo
1076              do il=1,nodes
1077 c               read(irest2,*) (t_restart1(j,il),j=1,4)
1078                do j=1,4
1079                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1080                enddo
1081              enddo
1082          endif
1083          call mpi_scatter(t_restart1,5,mpi_real,
1084      &           t5_restart1,5,mpi_real,king,mpi_comm_world,ierr)
1085          totT=t5_restart1(1)              
1086          EK=t5_restart1(2)
1087          potE=t5_restart1(3)
1088          t_bath=t5_restart1(4)
1089
1090          if(me.eq.king)then
1091               do il=0,nodes-1
1092                do i=1,2*nres
1093 c                read(irest2,'(3e15.5)') 
1094 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1095             do j=1,3
1096              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1097             enddo
1098                enddo
1099               enddo
1100          endif
1101          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1102      &           r_d,3*2*nres,mpi_real,king,mpi_comm_world,ierr)
1103
1104          do i=1,2*nres
1105            do j=1,3
1106             d_t(j,i)=r_d(j,i)
1107            enddo
1108          enddo
1109          if(me.eq.king)then 
1110               do il=0,nodes-1
1111                do i=1,2*nres
1112 c                read(irest2,'(3e15.5)') 
1113 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1114             do j=1,3
1115              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1116             enddo
1117                enddo
1118               enddo
1119          endif
1120          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1121      &           r_d,3*2*nres,mpi_real,king,mpi_comm_world,ierr)
1122          do i=1,2*nres
1123            do j=1,3
1124             dc(j,i)=r_d(j,i)
1125            enddo
1126          enddo
1127         if(me.eq.king) close(irest2)
1128         return
1129         end
1130
1131       subroutine read1restart_old
1132       implicit real*8 (a-h,o-z)
1133       include 'DIMENSIONS'
1134       include 'mpif.h'
1135       include 'COMMON.MD'
1136       include 'COMMON.IOUNITS'
1137       include 'COMMON.REMD'
1138       include 'COMMON.SETUP'
1139       include 'COMMON.CHAIN'
1140       include 'COMMON.SBRIDGE'
1141       include 'COMMON.INTERACT'
1142       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1143      &                 t5_restart1(5)
1144
1145          if(me.eq.king)then
1146              open(irest2,file=mremd_rst_name,status='unknown')
1147              read (irest2,*) (i2rep(i),i=0,nodes-1)
1148              read (irest2,*) (ifirst(i),i=1,remd_m(1))
1149              do il=1,nodes
1150               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1151               read (irest2,*) ndowna(0,il),
1152      &                    (ndowna(i,il),i=1,ndowna(0,il))
1153              enddo
1154              do il=1,nodes
1155                read(irest2,*) (t_restart1(j,il),j=1,4)
1156              enddo
1157          endif
1158          call mpi_scatter(t_restart1,5,mpi_real,
1159      &           t5_restart1,5,mpi_real,king,mpi_comm_world,ierr)
1160          totT=t5_restart1(1)              
1161          EK=t5_restart1(2)
1162          potE=t5_restart1(3)
1163          t_bath=t5_restart1(4)
1164
1165          if(me.eq.king)then
1166               do il=0,nodes-1
1167                do i=1,2*nres
1168                 read(irest2,'(3e15.5)') 
1169      &                (d_restart1(j,i+2*nres*il),j=1,3)
1170                enddo
1171               enddo
1172          endif
1173          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1174      &           r_d,3*2*nres,mpi_real,king,mpi_comm_world,ierr)
1175
1176          do i=1,2*nres
1177            do j=1,3
1178             d_t(j,i)=r_d(j,i)
1179            enddo
1180          enddo
1181          if(me.eq.king)then 
1182               do il=0,nodes-1
1183                do i=1,2*nres
1184                 read(irest2,'(3e15.5)') 
1185      &                (d_restart1(j,i+2*nres*il),j=1,3)
1186                enddo
1187               enddo
1188          endif
1189          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1190      &           r_d,3*2*nres,mpi_real,king,mpi_comm_world,ierr)
1191          do i=1,2*nres
1192            do j=1,3
1193             dc(j,i)=r_d(j,i)
1194            enddo
1195          enddo
1196         if(me.eq.king) close(irest2)
1197         return
1198         end
1199