xdrf2pdb-M.F: corrected order of reading
[unres.git] / source / unres / src_MD / old_F / MREMD.F.safe
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       integer iremd_iset(maxprocs)
28       integer*2 i_index
29      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
30       double precision remd_ene(0:n_ene+4,maxprocs)
31       integer iremd_acc(maxprocs),iremd_tot(maxprocs)
32       integer iremd_acc_usa(maxprocs),iremd_tot_usa(maxprocs)
33       integer ilen,rstcount
34       external ilen
35       character*50 tytul
36       common /gucio/ cm
37       integer itime
38 cold      integer nup(0:maxprocs),ndown(0:maxprocs)
39       integer rep2i(0:maxprocs),ireqi(maxprocs)
40       integer icache_all(maxprocs)
41       integer status(MPI_STATUS_SIZE),statusi(MPI_STATUS_SIZE,maxprocs)
42       logical synflag,end_of_run,file_exist /.false./
43
44 cdeb      imin_itime_old=0
45       ntwx_cache=0
46       time00=MPI_WTIME()
47       time01=time00
48       if(me.eq.king.or..not.out1file) then
49        write  (iout,*) 'MREMD',nodes,'time before',time00-walltime
50        write (iout,*) "NREP=",nrep
51       endif
52
53       synflag=.false.
54       if (ilen(tmpdir).gt.0 .and. (me.eq.king)) then
55         call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_mremd.rst")
56       endif
57       mremd_rst_name=prefix(:ilen(prefix))//"_mremd.rst"
58
59 cd      print *,'MREMD',nodes
60 cd      print *,'mmm',me,remd_mlist,(remd_m(i),i=1,nrep)
61 cde      write (iout,*) "Start MREMD: me",me," t_bath",t_bath
62       k=0
63       rep2i(k)=-1
64       do il=1,max0(nset,1)
65        do il1=1,max0(mset(il),1)
66         do i=1,nrep
67          iremd_acc(i)=0
68          iremd_acc_usa(i)=0
69          iremd_tot(i)=0
70          do j=1,remd_m(i)
71           i2rep(k)=i
72           i2set(k)=il
73           rep2i(i)=k
74           k=k+1
75           i_index(i,j,il,il1)=k
76          enddo
77         enddo
78        enddo
79       enddo
80
81       if(me.eq.king.or..not.out1file) then
82        write(iout,*) (i2rep(i),i=0,nodes-1)
83        write(iout,*) (i2set(i),i=0,nodes-1)
84        do il=1,nset
85         do il1=1,mset(il)
86          do i=1,nrep
87           do j=1,remd_m(i)
88            write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
89           enddo
90          enddo
91         enddo
92        enddo
93       endif
94
95 c      print *,'i2rep',me,i2rep(me)
96 c      print *,'rep2i',(rep2i(i),i=0,nrep)
97
98 cold       if (i2rep(me).eq.nrep) then
99 cold        nup(0)=0
100 cold       else
101 cold        nup(0)=remd_m(i2rep(me)+1)
102 cold        k=rep2i(int(i2rep(me)))+1
103 cold        do i=1,nup(0)
104 cold         nup(i)=k
105 cold         k=k+1
106 cold        enddo
107 cold       endif
108
109 cd       print '(i4,a4,100i4)',me,' nup',(nup(i),i=0,nup(0))
110
111 cold       if (i2rep(me).eq.1) then
112 cold        ndown(0)=0
113 cold       else
114 cold        ndown(0)=remd_m(i2rep(me)-1)
115 cold        k=rep2i(i2rep(me)-2)+1
116 cold        do i=1,ndown(0)
117 cold         ndown(i)=k
118 cold         k=k+1
119 cold        enddo
120 cold       endif
121
122 cd       print '(i4,a6,100i4)',me,' ndown',(ndown(i),i=0,ndown(0))
123
124        
125        write (*,*) "Processor",me," rest",rest,"
126      &   restart1fie",restart1file
127        if(rest.and.restart1file) then 
128            if (me.eq.king)
129      &     inquire(file=mremd_rst_name,exist=file_exist)
130 cd           write (*,*) me," Before broadcast: file_exist",file_exist
131            call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
132      &          IERR)
133 cd           write (*,*) me," After broadcast: file_exist",file_exist
134            if(file_exist) then 
135              if(me.eq.king.or..not.out1file)
136      &            write  (iout,*) 'Master is reading restart1file'
137              call read1restart(i_index)
138            else
139              if(me.eq.king.or..not.out1file)
140      &            write  (iout,*) 'WARNING : no restart1file'
141            endif
142
143            if(me.eq.king.or..not.out1file) then
144               write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
145               write(iout,*) "i_index"
146               do il=1,nset
147                do il1=1,mset(il)
148                 do i=1,nrep
149                  do j=1,remd_m(i)
150                   write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
151                  enddo
152                 enddo
153                enddo
154               enddo
155            endif 
156        endif
157
158        if(me.eq.king) then
159         if (rest.and..not.restart1file) 
160      &          inquire(file=mremd_rst_name,exist=file_exist)
161         if(.not.file_exist.and.rest.and..not.restart1file) 
162      &       write(iout,*) 'WARNING : no restart file',mremd_rst_name
163         IF (rest.and.file_exist.and..not.restart1file) THEN
164              write  (iout,*) 'Master is reading restart file',
165      &                        mremd_rst_name
166              open(irest2,file=mremd_rst_name,status='unknown')
167              read (irest2,*) 
168              read (irest2,*) (i2rep(i),i=0,nodes-1)
169              read (irest2,*) 
170              read (irest2,*) (ifirst(i),i=1,remd_m(1))
171              do il=1,nodes
172               read (irest2,*) 
173               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
174               read (irest2,*) 
175               read (irest2,*) ndowna(0,il),
176      &                    (ndowna(i,il),i=1,ndowna(0,il))
177              enddo
178              if(usampl) then
179               read (irest2,*)
180               read (irest2,*) nset
181               read (irest2,*) 
182               read (irest2,*) (mset(i),i=1,nset)
183               read (irest2,*) 
184               read (irest2,*) (i2set(i),i=0,nodes-1)
185               read (irest2,*) 
186               do il=1,nset
187                do il1=1,mset(il)
188                 do i=1,nrep
189                   read(irest2,*) (i_index(i,j,il,il1),j=1,remd_m(i))
190                 enddo
191                enddo
192               enddo
193
194               write(iout,*) "i2set",(i2set(i),i=0,nodes-1)
195               write(iout,*) "i_index"
196               do il=1,nset
197                do il1=1,mset(il)
198                 do i=1,nrep
199                  do j=1,remd_m(i)
200                   write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
201                  enddo
202                 enddo
203                enddo
204               enddo
205              endif
206
207              close(irest2)
208
209              write (iout,'(a6,1000i5)') "i2rep",(i2rep(i),i=0,nodes-1)
210              write (iout,'(a6,1000i5)') "ifirst",
211      &                    (ifirst(i),i=1,remd_m(1))
212              do il=1,nodes
213               write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
214      &                    (nupa(i,il),i=1,nupa(0,il))
215               write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
216      &                    (ndowna(i,il),i=1,ndowna(0,il))
217              enddo
218         ELSE IF (.not.(rest.and.file_exist)) THEN
219          do il=1,remd_m(1)
220           ifirst(il)=il
221          enddo
222
223          do il=1,nodes
224           if (i2rep(il-1).eq.nrep) then
225            nupa(0,il)=0
226           else
227            nupa(0,il)=remd_m(i2rep(il-1)+1)
228            k=rep2i(int(i2rep(il-1)))+1
229            do i=1,nupa(0,il)
230             nupa(i,il)=k+1
231             k=k+1
232            enddo
233           endif
234           if (i2rep(il-1).eq.1) then
235            ndowna(0,il)=0
236           else
237            ndowna(0,il)=remd_m(i2rep(il-1)-1)
238            k=rep2i(i2rep(il-1)-2)+1
239            do i=1,ndowna(0,il)
240             ndowna(i,il)=k+1
241             k=k+1
242            enddo
243           endif
244          enddo
245         
246 c        write (iout,'(a6,100i4)') "ifirst",
247 c     &                    (ifirst(i),i=1,remd_m(1))
248 c        do il=1,nodes
249 c         write (iout,'(a6,i4,a1,100i4)') "nupa",il,":",
250 c     &                    (nupa(i,il),i=1,nupa(0,il))
251 c         write (iout,'(a6,i4,a1,100i4)') "ndowna",il,":",
252 c     &                    (ndowna(i,il),i=1,ndowna(0,il))
253 c        enddo
254         
255         ENDIF
256        endif
257 c
258 c      t_bath=retmin+(retmax-retmin)*me/(nodes-1)
259        if(.not.(rest.and.file_exist.and.restart1file)) then
260          if (me .eq. king) then 
261             t_bath=retmin
262          else 
263             t_bath=retmin+(retmax-retmin)*exp(float(i2rep(me)-nrep))
264          endif
265 cd       print *,'ttt',me,remd_tlist,(remd_t(i),i=1,nrep)
266          if (remd_tlist) t_bath=remd_t(int(i2rep(me)))
267
268        endif
269        if(usampl) then
270           iset=i2set(me)
271           if(me.eq.king.or..not.out1file) 
272      &     write(iout,*) me,"iset=",iset,"t_bath=",t_bath
273        endif        
274 c
275        stdfp=dsqrt(2*Rb*t_bath/d_time)
276        do i=1,ntyp
277           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
278        enddo 
279
280 c      print *,'irep',me,t_bath
281        if (.not.rest) then  
282         if (me.eq.king .or. .not. out1file)
283      &   write (iout,'(a60,f10.5)') "REMD Temperature:",t_bath
284         call rescale_weights(t_bath)
285        endif
286
287
288 c------copy MD--------------
289 c  The driver for molecular dynamics subroutines
290 c------------------------------------------------
291       t_MDsetup=0.0d0
292       t_langsetup=0.0d0
293       t_MD=0.0d0
294       t_enegrad=0.0d0
295       t_sdsetup=0.0d0
296       if(me.eq.king.or..not.out1file)
297      & write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started"
298 #ifdef MPI
299       tt0 = MPI_Wtime()
300 #else
301       tt0 = tcpu()
302 #endif
303 c Determine the inverse of the inertia matrix.
304       call setup_MD_matrices
305 c Initialize MD
306       call init_MD
307       if (rest) then  
308        if (me.eq.king .or. .not. out1file)
309      &  write (iout,'(a60,f10.5)') "REMD restart Temperature:",t_bath
310        stdfp=dsqrt(2*Rb*t_bath/d_time)
311        do i=1,ntyp
312           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
313        enddo 
314        call rescale_weights(t_bath)
315       endif
316
317 #ifdef MPI
318       t_MDsetup = MPI_Wtime()
319 #else
320       t_MDsetup = tcpu()-tt0
321 #endif
322       rstcount=0 
323 c   Entering the MD loop       
324 #ifdef MPI
325       tt0 = MPI_Wtime()
326 #else
327       tt0 = tcpu()
328 #endif
329       if (lang.eq.2 .or. lang.eq.3) then
330         call setup_fricmat
331         if (lang.eq.2) then
332           call sd_verlet_p_setup        
333         else
334           call sd_verlet_ciccotti_setup
335         endif
336 #ifndef   LANG0
337         do i=1,dimen
338           do j=1,dimen
339             pfric0_mat(i,j,0)=pfric_mat(i,j)
340             afric0_mat(i,j,0)=afric_mat(i,j)
341             vfric0_mat(i,j,0)=vfric_mat(i,j)
342             prand0_mat(i,j,0)=prand_mat(i,j)
343             vrand0_mat1(i,j,0)=vrand_mat1(i,j)
344             vrand0_mat2(i,j,0)=vrand_mat2(i,j)
345           enddo
346         enddo
347 #endif
348         flag_stoch(0)=.true.
349         do i=1,maxflag_stoch
350           flag_stoch(i)=.false.
351         enddo  
352       else if (lang.eq.1 .or. lang.eq.4) then
353         call setup_fricmat
354       endif
355       time00=MPI_WTIME()
356       if (me.eq.king .or. .not. out1file)
357      & write(iout,*) 'Setup time',time00-walltime
358       call flush(iout)
359 #ifdef MPI
360       t_langsetup=MPI_Wtime()
361       tt0=MPI_Wtime()
362 #else
363       t_langsetup=tcpu()-tt0
364       tt0=tcpu()
365 #endif
366       itime=0
367       end_of_run=.false.
368       do while(.not.end_of_run)
369         itime=itime+1
370         if(itime.eq.n_timestep.and.me.eq.king) end_of_run=.true.
371         if(mremdsync.and.itime.eq.n_timestep) end_of_run=.true.
372         rstcount=rstcount+1
373         if (lang.gt.0 .and. surfarea .and. 
374      &      mod(itime,reset_fricmat).eq.0) then
375           if (lang.eq.2 .or. lang.eq.3) then
376             call setup_fricmat
377             if (lang.eq.2) then
378               call sd_verlet_p_setup
379             else
380               call sd_verlet_ciccotti_setup
381             endif
382 #ifndef   LANG0
383             do i=1,dimen
384               do j=1,dimen
385                 pfric0_mat(i,j,0)=pfric_mat(i,j)
386                 afric0_mat(i,j,0)=afric_mat(i,j)
387                 vfric0_mat(i,j,0)=vfric_mat(i,j)
388                 prand0_mat(i,j,0)=prand_mat(i,j)
389                 vrand0_mat1(i,j,0)=vrand_mat1(i,j)
390                 vrand0_mat2(i,j,0)=vrand_mat2(i,j)
391               enddo
392             enddo
393 #endif
394             flag_stoch(0)=.true.
395             do i=1,maxflag_stoch
396               flag_stoch(i)=.false.
397             enddo   
398           else if (lang.eq.1 .or. lang.eq.4) then
399             call setup_fricmat
400           endif
401           write (iout,'(a,i10)') 
402      &      "Friction matrix reset based on surface area, itime",itime
403         endif
404         if (reset_vel .and. tbf .and. lang.eq.0 
405      &      .and. mod(itime,count_reset_vel).eq.0) then
406           call random_vel
407           if (me.eq.king .or. .not. out1file)
408      &     write(iout,'(a,f20.2)') 
409      &     "Velocities reset to random values, time",totT       
410           do i=0,2*nres
411             do j=1,3
412               d_t_old(j,i)=d_t(j,i)
413             enddo
414           enddo
415         endif
416         if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then
417           call inertia_tensor  
418           call vcm_vel(vcm)
419           do j=1,3
420              d_t(j,0)=d_t(j,0)-vcm(j)
421           enddo
422           call kinetic(EK)
423           kinetic_T=2.0d0/(dimen*Rb)*EK
424           scalfac=dsqrt(T_bath/kinetic_T)
425 cd          write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT     
426           do i=0,2*nres
427             do j=1,3
428               d_t_old(j,i)=scalfac*d_t(j,i)
429             enddo
430           enddo
431         endif  
432         if (lang.ne.4) then
433           if (RESPA) then
434 c Time-reversible RESPA algorithm 
435 c (Tuckerman et al., J. Chem. Phys., 97, 1990, 1992)
436             call RESPA_step(itime)
437           else
438 c Variable time step algorithm.
439             call velverlet_step(itime)
440           endif
441         else
442           call brown_step(itime)
443         endif
444         if(ntwe.ne.0) then
445           if (mod(itime,ntwe).eq.0) call statout(itime)
446         endif
447         if (mod(itime,ntwx).eq.0.and..not.traj1file) then
448           write (tytul,'("time",f8.2," temp",f8.1)') totT,t_bath
449           if(mdpdb) then
450              call pdbout(potE,tytul,ipdb)
451           else 
452              call cartout(totT)
453           endif
454         endif
455         if (mod(itime,ntwx).eq.0.and.traj1file) then
456           if(ntwx_cache.lt.max_cache_traj_use) then
457             ntwx_cache=ntwx_cache+1
458           else
459            if (max_cache_traj_use.ne.1)
460      &      print *,itime,"processor ",me," over cache ",ntwx_cache
461            do i=1,ntwx_cache-1
462
463             totT_cache(i)=totT_cache(i+1)
464             EK_cache(i)=EK_cache(i+1)
465             potE_cache(i)=potE_cache(i+1)
466             t_bath_cache(i)=t_bath_cache(i+1)
467             Uconst_cache(i)=Uconst_cache(i+1)
468             iset_cache(i)=iset_cache(i+1)
469
470             do ii=1,nfrag
471              qfrag_cache(ii,i)=qfrag_cache(ii,i+1)
472             enddo
473             do ii=1,npair
474              qpair_cache(ii,i)=qpair_cache(ii,i+1)
475             enddo
476             do ii=1,nfrag_back
477               utheta_cache(ii,i)=utheta_cache(ii,i+1)
478               ugamma_cache(ii,i)=ugamma_cache(ii,i+1)
479               uscdiff_cache(ii,i)=uscdiff_cache(ii,i+1)
480             enddo
481
482
483             do ii=1,nres*2
484              do j=1,3
485               c_cache(j,ii,i)=c_cache(j,ii,i+1)
486              enddo
487             enddo
488            enddo
489           endif
490
491             totT_cache(ntwx_cache)=totT
492             EK_cache(ntwx_cache)=EK
493             potE_cache(ntwx_cache)=potE
494             t_bath_cache(ntwx_cache)=t_bath
495             Uconst_cache(ntwx_cache)=Uconst
496             iset_cache(ntwx_cache)=iset
497
498             do i=1,nfrag
499              qfrag_cache(i,ntwx_cache)=qfrag(i)
500             enddo
501             do i=1,npair
502              qpair_cache(i,ntwx_cache)=qpair(i)
503             enddo
504             do i=1,nfrag_back
505               utheta_cache(i,ntwx_cache)=utheta(i)
506               ugamma_cache(i,ntwx_cache)=ugamma(i)
507               uscdiff_cache(i,ntwx_cache)=uscdiff(i)
508             enddo
509
510             do i=1,nres*2
511              do j=1,3
512               c_cache(j,i,ntwx_cache)=c(j,i)
513              enddo
514             enddo
515
516         endif
517         if ((rstcount.eq.1000.or.itime.eq.n_timestep)
518      &                         .and..not.restart1file) then
519
520            if(me.eq.king) then
521              open(irest1,file=mremd_rst_name,status='unknown')
522              write (irest1,*) "i2rep"
523              write (irest1,*) (i2rep(i),i=0,nodes-1)
524              write (irest1,*) "ifirst"
525              write (irest1,*) (ifirst(i),i=1,remd_m(1))
526              do il=1,nodes
527               write (irest1,*) "nupa",il
528               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
529               write (irest1,*) "ndowna",il
530               write (irest1,*) ndowna(0,il),
531      &                   (ndowna(i,il),i=1,ndowna(0,il))
532              enddo
533              if(usampl) then
534               write (irest1,*) "nset"
535               write (irest1,*) nset
536               write (irest1,*) "mset"
537               write (irest1,*) (mset(i),i=1,nset)
538               write (irest1,*) "i2set"
539               write (irest1,*) (i2set(i),i=0,nodes-1)
540               write (irest1,*) "i_index"
541               do il=1,nset
542                do il1=1,mset(il)
543                 do i=1,nrep
544                   write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
545                 enddo
546                enddo
547               enddo
548
549              endif
550              close(irest1)
551            endif
552            open(irest2,file=rest2name,status='unknown')
553            write(irest2,*) totT,EK,potE,totE,t_bath
554            do i=1,2*nres
555             write (irest2,'(3e15.5)') (d_t(j,i),j=1,3)
556            enddo
557            do i=1,2*nres
558             write (irest2,'(3e15.5)') (dc(j,i),j=1,3)
559            enddo
560            if(usampl) then
561              write (irest2,*) iset
562            endif
563           close(irest2)
564           rstcount=0
565         endif 
566
567 c REMD - exchange
568 c forced synchronization
569         if (mod(itime,i_sync_step).eq.0 .and. me.ne.king 
570      &                                .and. .not. mremdsync) then 
571             synflag=.false.
572             call mpi_iprobe(0,101,CG_COMM,synflag,status,ierr)
573             if (synflag) then 
574                call mpi_recv(itime_master, 1, MPI_INTEGER,
575      &                             0,101,CG_COMM, status, ierr)
576                call mpi_barrier(CG_COMM, ierr)
577 cdeb               if (out1file.or.traj1file) then
578 cdeb                call mpi_gather(itime,1,mpi_integer,
579 cdeb     &             icache_all,1,mpi_integer,king,
580 cdeb     &             CG_COMM,ierr)                 
581                if(traj1file)
582      &          call mpi_gather(ntwx_cache,1,mpi_integer,
583      &             icache_all,1,mpi_integer,king,
584      &             CG_COMM,ierr)
585                if (.not.out1file)
586      &               write(iout,*) 'REMD synchro at',itime_master,itime
587                if(itime_master.ge.n_timestep) end_of_run=.true.
588 ctime               call flush(iout)
589             endif
590         endif
591
592 c REMD - exchange
593         if ((mod(itime,nstex).eq.0.and.me.eq.king
594      &                  .or.end_of_run.and.me.eq.king )
595      &       .and. .not. mremdsync ) then
596            synflag=.true.
597            time00=MPI_WTIME()
598            do i=1,nodes-1
599               call mpi_isend(itime,1,MPI_INTEGER,i,101,
600      &                                CG_COMM, ireqi(i), ierr)
601 cd            write(iout,*) 'REMD synchro with',i
602 cd            call flush(iout)
603            enddo
604            call mpi_waitall(nodes-1,ireqi,statusi,ierr)
605            call mpi_barrier(CG_COMM, ierr)
606            time01=MPI_WTIME()
607            write(iout,*) 'REMD synchro at',itime,'time=',time01-time00
608            if (out1file.or.traj1file) then
609 cdeb            call mpi_gather(itime,1,mpi_integer,
610 cdeb     &             itime_all,1,mpi_integer,king,
611 cdeb     &             CG_COMM,ierr)
612 cdeb            write(iout,'(a19,8000i8)') ' REMD synchro itime',
613 cdeb     &                    (itime_all(i),i=1,nodes)
614             if(traj1file) then
615 cdeb             imin_itime=itime_all(1)
616 cdeb             do i=2,nodes
617 cdeb               if(itime_all(i).lt.imin_itime) imin_itime=itime_all(i)
618 cdeb             enddo
619 cdeb             ii_write=(imin_itime-imin_itime_old)/ntwx
620 cdeb             imin_itime_old=int(imin_itime/ntwx)*ntwx
621 cdeb             write(iout,*) imin_itime,imin_itime_old,ii_write
622              call mpi_gather(ntwx_cache,1,mpi_integer,
623      &             icache_all,1,mpi_integer,king,
624      &             CG_COMM,ierr)
625 c             write(iout,'(a19,8000i8)') '     ntwx_cache',
626 c     &                    (icache_all(i),i=1,nodes)
627              ii_write=icache_all(1)
628              do i=2,nodes
629                if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
630              enddo
631 c             write(iout,*) "MIN ii_write=",ii_write
632             endif
633            endif
634 ctime           call flush(iout)
635         endif
636         if(mremdsync .and. mod(itime,nstex).eq.0) then
637            synflag=.true.
638            if (me.eq.king .or. .not. out1file)
639      &      write(iout,*) 'REMD synchro at',itime
640
641             if(traj1file) then
642              call mpi_gather(ntwx_cache,1,mpi_integer,
643      &             icache_all,1,mpi_integer,king,
644      &             CG_COMM,ierr)
645              if (me.eq.king) then
646                write(iout,'(a19,8000i8)') '     ntwx_cache',
647      &                    (icache_all(i),i=1,nodes)
648                ii_write=icache_all(1)
649                do i=2,nodes
650                  if(icache_all(i).lt.ii_write) ii_write=icache_all(i)
651                enddo
652                write(iout,*) "MIN ii_write=",ii_write
653              endif
654             endif
655
656            call flush(iout)
657         endif
658         if(synflag.and..not.end_of_run) then
659            time02=MPI_WTIME()
660            synflag=.false.
661
662 cd           write(iout,*) 'REMD before',me,t_bath
663
664 c           call mpi_gather(t_bath,1,mpi_double_precision,
665 c     &             remd_t_bath,1,mpi_double_precision,king,
666 c     &             CG_COMM,ierr)
667            potEcomp(n_ene+1)=t_bath
668            if (usampl) then
669              potEcomp(n_ene+2)=iset
670              if (iset.lt.nset) then
671                i_set_temp=iset
672                iset=iset+1
673                call EconstrQ
674                potEcomp(n_ene+3)=Uconst
675                iset=i_set_temp
676              endif
677              if (iset.gt.1) then
678                i_set_temp=iset
679                iset=iset-1
680                call EconstrQ
681                potEcomp(n_ene+4)=Uconst 
682                iset=i_set_temp
683              endif
684            endif
685            call mpi_gather(potEcomp(0),n_ene+5,mpi_double_precision,
686      &             remd_ene(0,1),n_ene+5,mpi_double_precision,king,
687      &             CG_COMM,ierr)
688            if(lmuca) then 
689             call mpi_gather(elow,1,mpi_double_precision,
690      &             elowi,1,mpi_double_precision,king,
691      &             CG_COMM,ierr)
692             call mpi_gather(ehigh,1,mpi_double_precision,
693      &             ehighi,1,mpi_double_precision,king,
694      &             CG_COMM,ierr)
695            endif
696
697           time03=MPI_WTIME()
698           if (me.eq.king .or. .not. out1file) then
699             write(iout,*) 'REMD gather times=',time03-time01
700      &                                        ,time03-time02
701           endif
702
703           if (restart1file) call write1rst(i_index)
704
705           time04=MPI_WTIME()
706           if (me.eq.king .or. .not. out1file) then
707             write(iout,*) 'REMD writing rst time=',time04-time03
708           endif
709
710           if (traj1file) call write1traj
711 cd debugging
712 cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
713 cdeb     &             icache_all,1,mpi_integer,king,
714 cdeb     &             CG_COMM,ierr)
715 cdeb            write(iout,'(a19,8000i8)') '  ntwx_cache after traj1file',
716 cdeb     &                    (icache_all(i),i=1,nodes)
717 cd end
718
719
720           time05=MPI_WTIME()
721           if (me.eq.king .or. .not. out1file) then
722             write(iout,*) 'REMD writing traj time=',time05-time04
723             call flush(iout)
724           endif
725
726
727           if (me.eq.king) then
728             do i=1,nodes
729                remd_t_bath(i)=remd_ene(n_ene+1,i)
730                iremd_iset(i)=remd_ene(n_ene+2,i)
731             enddo
732             if(lmuca) then
733 co             write(iout,*) 'REMD exchange temp,ene,elow,ehigh'
734              do i=1,nodes
735                write(iout,'(i4,4f12.5)') i,remd_t_bath(i),remd_ene(0,i),
736      &            elowi(i),ehighi(i)       
737              enddo
738             else
739 cd              write(iout,*) 'REMD exchange temp,ene'
740 c             do i=1,nodes
741 co              write(iout,'(i4,2f12.5)') i,remd_t_bath(i),remd_ene(0,i)
742 c              write(iout,'(6f12.5)') (remd_ene(j,i),j=1,n_ene)
743 c             enddo
744             endif
745 c-------------------------------------           
746            IF(.not.usampl) THEN
747             write (iout,*) "Enter exchnge, remd_m",remd_m(1),
748      &        " nodes",nodes
749             call flush(iout)
750             do irr=1,remd_m(1)
751                i=ifirst(iran_num(1,remd_m(1)))
752              call flush(iout)
753
754              do ii=1,nodes-1
755
756              if(i.gt.0.and.nupa(0,i).gt.0) then
757               iex=nupa(iran_num(1,int(nupa(0,i))),i)
758               if (lmuca) then
759                call muca_delta(remd_t_bath,remd_ene,i,iex,delta)
760               else
761 c Swap temperatures between conformations i and iex with recalculating the free energies
762 c following temperature changes.
763                ene_iex_iex=remd_ene(0,iex)
764                ene_i_i=remd_ene(0,i)
765 cd              write (iout,*) "rescaling weights with temperature",
766 cd     &          remd_t_bath(i)
767 cd              call flush(iout)
768                call rescale_weights(remd_t_bath(i))
769
770 cd               write (iout,*) "0,iex",remd_t_bath(i)
771 cd               call enerprint(remd_ene(0,iex))
772
773                call sum_energy(remd_ene(0,iex),.false.)
774                ene_iex_i=remd_ene(0,iex)
775 cd               write (iout,*) "ene_iex_i",remd_ene(0,iex)
776
777 cd               write (iout,*) "0,i",remd_t_bath(i)
778 cd               call enerprint(remd_ene(0,i))
779
780                call sum_energy(remd_ene(0,i),.false.)
781 cd               write (iout,*) "ene_i_i",remd_ene(0,i)
782 cd               call flush(iout)
783 cd              write (iout,*) "rescaling weights with temperature",
784 cd     &          remd_t_bath(iex)
785                if (real(ene_i_i).ne.real(remd_ene(0,i))) then
786                 write (iout,*) "ERROR: inconsistent energies:",i,
787      &            ene_i_i,remd_ene(0,i)
788                endif
789                call rescale_weights(remd_t_bath(iex))
790
791 cd               write (iout,*) "0,i",remd_t_bath(iex)
792 cd               call enerprint(remd_ene(0,i))
793
794                call sum_energy(remd_ene(0,i),.false.)
795 cd               write (iout,*) "ene_i_iex",remd_ene(0,i)
796 cd              call flush(iout)
797                ene_i_iex=remd_ene(0,i)
798
799 cd               write (iout,*) "0,iex",remd_t_bath(iex)
800 cd               call enerprint(remd_ene(0,iex))
801
802                call sum_energy(remd_ene(0,iex),.false.)
803                if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
804                 write (iout,*) "ERROR: inconsistent energies:",iex,
805      &            ene_iex_iex,remd_ene(0,iex)
806                endif
807 cd               write (iout,*) "ene_iex_iex",remd_ene(0,iex)
808 cd               write (iout,*) "i",i," iex",iex
809 cd               write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
810 cd     &           " ene_i_iex",ene_i_iex,
811 cd     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
812 cd              call flush(iout)
813                delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
814      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
815                delta=-delta
816 cd               write(iout,*) 'delta',delta
817 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
818 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
819 c     &              (remd_t_bath(i)*remd_t_bath(iex))
820               endif
821               if (delta .gt. 50.0d0) then
822                 delta=0.0d0
823               else
824 #ifdef OSF 
825                 if(isnan(delta))then
826                   delta=0.0d0
827                 else if (delta.lt.-50.0d0) then
828                   delta=dexp(50.0d0)
829                 else
830                   delta=dexp(-delta)
831                 endif
832 #else
833                 delta=dexp(-delta)
834 #endif
835               endif
836               iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
837               xxx=ran_number(0.0d0,1.0d0)
838 cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
839 cd              call flush(iout)
840               if (delta .gt. xxx) then
841                 tmp=remd_t_bath(i)       
842                 remd_t_bath(i)=remd_t_bath(iex)
843                 remd_t_bath(iex)=tmp
844                 remd_ene(0,i)=ene_i_iex
845                 remd_ene(0,iex)=ene_iex_i
846                 if(lmuca) then
847                   tmp=elowi(i)
848                   elowi(i)=elowi(iex)
849                   elowi(iex)=tmp  
850                   tmp=ehighi(i)
851                   ehighi(i)=ehighi(iex)
852                   ehighi(iex)=tmp  
853                 endif
854
855
856                 do k=0,nodes
857                   itmp=nupa(k,i)
858                   nupa(k,i)=nupa(k,iex)
859                   nupa(k,iex)=itmp
860                   itmp=ndowna(k,i)
861                   ndowna(k,i)=ndowna(k,iex)
862                   ndowna(k,iex)=itmp
863                 enddo
864                 do il=1,nodes
865                  if (ifirst(il).eq.i) ifirst(il)=iex
866                  do k=1,nupa(0,il)
867                   if (nupa(k,il).eq.i) then 
868                      nupa(k,il)=iex
869                   elseif (nupa(k,il).eq.iex) then 
870                      nupa(k,il)=i
871                   endif
872                  enddo
873                  do k=1,ndowna(0,il)
874                   if (ndowna(k,il).eq.i) then 
875                      ndowna(k,il)=iex
876                   elseif (ndowna(k,il).eq.iex) then 
877                      ndowna(k,il)=i
878                   endif
879                  enddo
880                 enddo
881
882                 iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
883                 itmp=i2rep(i-1)
884                 i2rep(i-1)=i2rep(iex-1)
885                 i2rep(iex-1)=itmp
886
887 cd              write(iout,*) 'exchange',i,iex
888 cd              write (iout,'(a8,100i4)') "@ ifirst",
889 cd     &                    (ifirst(k),k=1,remd_m(1))
890 cd              do il=1,nodes
891 cd               write (iout,'(a8,i4,a1,100i4)') "@ nupa",il,":",
892 cd     &                    (nupa(k,il),k=1,nupa(0,il))
893 cd               write (iout,'(a8,i4,a1,100i4)') "@ ndowna",il,":",
894 cd     &                    (ndowna(k,il),k=1,ndowna(0,il))
895 cd              enddo
896                 call flush(iout) 
897
898               else
899                remd_ene(0,iex)=ene_iex_iex
900                remd_ene(0,i)=ene_i_i
901                i=iex
902               endif 
903             endif
904            enddo
905            enddo
906 cd           write (iout,*) "exchange completed"
907 cd           call flush(iout) 
908         ELSE
909           do ii=1,nodes  
910 cd            write(iout,*) "########",ii
911
912             i_temp=iran_num(1,nrep)
913             i_mult=iran_num(1,remd_m(i_temp))
914             i_iset=iran_num(1,nset)
915             i_mset=iran_num(1,mset(i_iset))
916             i=i_index(i_temp,i_mult,i_iset,i_mset)
917
918 cd            write(iout,*) "i=",i,i_temp,i_mult,i_iset,i_mset
919
920             i_dir=iran_num(1,3)
921 cd            write(iout,*) "i_dir=",i_dir
922
923             if(i_dir.eq.1 .and. remd_m(i_temp+1).gt.0 )then            
924                
925                i_temp1=i_temp+1
926                i_mult1=iran_num(1,remd_m(i_temp1))
927                i_iset1=i_iset
928                i_mset1=iran_num(1,mset(i_iset1))
929                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
930
931             elseif(i_dir.eq.2 .and. mset(i_iset+1).gt.0)then
932
933                i_temp1=i_temp
934                i_mult1=iran_num(1,remd_m(i_temp1))
935                i_iset1=i_iset+1
936                i_mset1=iran_num(1,mset(i_iset1))
937                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
938                econstr_temp_i=remd_ene(20,i)
939                econstr_temp_iex=remd_ene(20,iex)
940                remd_ene(20,i)=remd_ene(n_ene+3,i)
941                remd_ene(20,iex)=remd_ene(n_ene+4,iex)
942
943             elseif(remd_m(i_temp+1).gt.0.and.mset(i_iset+1).gt.0)then
944
945                i_temp1=i_temp+1
946                i_mult1=iran_num(1,remd_m(i_temp1))
947                i_iset1=i_iset+1
948                i_mset1=iran_num(1,mset(i_iset1))
949                iex=i_index(i_temp1,i_mult1,i_iset1,i_mset1)
950                econstr_temp_i=remd_ene(20,i)
951                econstr_temp_iex=remd_ene(20,iex)
952                remd_ene(20,i)=remd_ene(n_ene+3,i)
953                remd_ene(20,iex)=remd_ene(n_ene+4,iex)
954
955             else
956                goto 444 
957             endif
958  
959 cd            write(iout,*) "iex=",iex,i_temp1,i_mult1,i_iset1,i_mset1
960             call flush(iout)
961
962 c Swap temperatures between conformations i and iex with recalculating the free energies
963 c following temperature changes.
964               ene_iex_iex=remd_ene(0,iex)
965               ene_i_i=remd_ene(0,i)
966 co              write (iout,*) "rescaling weights with temperature",
967 co     &          remd_t_bath(i)
968               call rescale_weights(remd_t_bath(i))
969               
970               call sum_energy(remd_ene(0,iex),.false.)
971               ene_iex_i=remd_ene(0,iex)
972 cd              write (iout,*) "ene_iex_i",remd_ene(0,iex)
973 c              call sum_energy(remd_ene(0,i),.false.)
974 cd              write (iout,*) "ene_i_i",remd_ene(0,i)
975 c              write (iout,*) "rescaling weights with temperature",
976 c     &          remd_t_bath(iex)
977 c              if (real(ene_i_i).ne.real(remd_ene(0,i))) then
978 c                write (iout,*) "ERROR: inconsistent energies:",i,
979 c     &            ene_i_i,remd_ene(0,i)
980 c              endif
981               call rescale_weights(remd_t_bath(iex))
982               call sum_energy(remd_ene(0,i),.false.)
983 cd              write (iout,*) "ene_i_iex",remd_ene(0,i)
984               ene_i_iex=remd_ene(0,i)
985 c              call sum_energy(remd_ene(0,iex),.false.)
986 c              if (real(ene_iex_iex).ne.real(remd_ene(0,iex))) then
987 c                write (iout,*) "ERROR: inconsistent energies:",iex,
988 c     &            ene_iex_iex,remd_ene(0,iex)
989 c              endif
990 cd              write (iout,*) "ene_iex_iex",remd_ene(0,iex)
991 c              write (iout,*) "i",i," iex",iex
992 cd              write (iout,'(4(a,e15.5))') "ene_i_i",ene_i_i,
993 cd     &           " ene_i_iex",ene_i_iex,
994 cd     &           " ene_iex_i",ene_iex_i," ene_iex_iex",ene_iex_iex
995               delta=(ene_iex_iex-ene_i_iex)/(Rb*remd_t_bath(iex))-
996      &              (ene_iex_i-ene_i_i)/(Rb*remd_t_bath(i))
997               delta=-delta
998 cd              write(iout,*) 'delta',delta
999 c              delta=(remd_t_bath(i)-remd_t_bath(iex))*
1000 c     &              (remd_ene(i)-remd_ene(iex))/Rb/
1001 c     &              (remd_t_bath(i)*remd_t_bath(iex))
1002               if (delta .gt. 50.0d0) then
1003                 delta=0.0d0
1004               else
1005                 delta=dexp(-delta)
1006               endif
1007               if (i_dir.eq.1.or.i_dir.eq.3)
1008      &         iremd_tot(int(i2rep(i-1)))=iremd_tot(int(i2rep(i-1)))+1
1009               if (i_dir.eq.2.or.i_dir.eq.3)
1010      &          iremd_tot_usa(int(i2set(i-1)))=
1011      &                 iremd_tot_usa(int(i2set(i-1)))+1
1012               xxx=ran_number(0.0d0,1.0d0)
1013 cd              write(iout,'(2i4,a6,2f12.5)') i,iex,' delta',delta,xxx
1014               if (delta .gt. xxx) then
1015                 tmp=remd_t_bath(i)       
1016                 remd_t_bath(i)=remd_t_bath(iex)
1017                 remd_t_bath(iex)=tmp
1018
1019                 itmp=iremd_iset(i)       
1020                 iremd_iset(i)=iremd_iset(iex)
1021                 iremd_iset(iex)=itmp
1022
1023                 remd_ene(0,i)=ene_i_iex
1024                 remd_ene(0,iex)=ene_iex_i
1025
1026                 if (i_dir.eq.1.or.i_dir.eq.3) 
1027      &           iremd_acc(int(i2rep(i-1)))=iremd_acc(int(i2rep(i-1)))+1
1028
1029                 itmp=i2rep(i-1)
1030                 i2rep(i-1)=i2rep(iex-1)
1031                 i2rep(iex-1)=itmp
1032
1033                 if (i_dir.eq.2.or.i_dir.eq.3) 
1034      &           iremd_acc_usa(int(i2set(i-1)))=
1035      &                 iremd_acc_usa(int(i2set(i-1)))+1
1036
1037                 itmp=i2set(i-1)
1038                 i2set(i-1)=i2set(iex-1)
1039                 i2set(iex-1)=itmp
1040         
1041                 itmp=i_index(i_temp,i_mult,i_iset,i_mset)
1042                 i_index(i_temp,i_mult,i_iset,i_mset)=
1043      &                i_index(i_temp1,i_mult1,i_iset1,i_mset1)
1044                 i_index(i_temp1,i_mult1,i_iset1,i_mset1)=itmp
1045
1046               else
1047                remd_ene(0,iex)=ene_iex_iex
1048                remd_ene(0,i)=ene_i_i
1049                remd_ene(20,iex)=econstr_temp_iex
1050                remd_ene(20,i)=econstr_temp_i
1051               endif
1052
1053 cd      do il=1,nset
1054 cd       do il1=1,mset(il)
1055 cd        do i=1,nrep
1056 cd         do j=1,remd_m(i)
1057 cd          write(iout,*) i,j,il,il1,i_index(i,j,il,il1)
1058 cd         enddo
1059 cd        enddo
1060 cd       enddo
1061 cd      enddo
1062
1063  444      continue           
1064
1065           enddo
1066
1067
1068         ENDIF
1069
1070 c-------------------------------------
1071              write (iout,*) "NREP",nrep
1072              do i=1,nrep
1073               if(iremd_tot(i).ne.0)
1074      &          write(iout,'(a3,i4,2f12.5,i5)') 'ACC',i,remd_t(i)
1075      &           ,iremd_acc(i)/(1.0*iremd_tot(i)),iremd_tot(i)
1076              enddo
1077
1078              if(usampl) then
1079               do i=1,nset
1080                if(iremd_tot_usa(i).ne.0)
1081      &           write(iout,'(a10,i4,f12.5,i8)') 'ACC_usampl',i,
1082      &         iremd_acc_usa(i)/(1.0*iremd_tot_usa(i)),iremd_tot_usa(i)
1083               enddo
1084              endif
1085
1086              call flush(iout)
1087
1088 cd              write (iout,'(a6,100i4)') "ifirst",
1089 cd     &                    (ifirst(i),i=1,remd_m(1))
1090 cd              do il=1,nodes
1091 cd               write (iout,'(a5,i4,a1,100i4)') "nup",il,":",
1092 cd     &                    (nupa(i,il),i=1,nupa(0,il))
1093 cd               write (iout,'(a5,i4,a1,100i4)') "ndown",il,":",
1094 cd     &                    (ndowna(i,il),i=1,ndowna(0,il))
1095 cd              enddo
1096             endif
1097
1098          time06=MPI_WTIME()
1099 cd         write (iout,*) "Before scatter"
1100 cd         call flush(iout)
1101          call mpi_scatter(remd_t_bath,1,mpi_double_precision,
1102      &           t_bath,1,mpi_double_precision,king,
1103      &           CG_COMM,ierr) 
1104 cd         write (iout,*) "After scatter"
1105 cd         call flush(iout)
1106          if(usampl)
1107      &    call mpi_scatter(iremd_iset,1,mpi_integer,
1108      &           iset,1,mpi_integer,king,
1109      &           CG_COMM,ierr) 
1110
1111          time07=MPI_WTIME()
1112           if (me.eq.king .or. .not. out1file) then
1113             write(iout,*) 'REMD scatter time=',time07-time06
1114           endif
1115
1116          if(lmuca) then
1117            call mpi_scatter(elowi,1,mpi_double_precision,
1118      &           elow,1,mpi_double_precision,king,
1119      &           CG_COMM,ierr) 
1120            call mpi_scatter(ehighi,1,mpi_double_precision,
1121      &           ehigh,1,mpi_double_precision,king,
1122      &           CG_COMM,ierr) 
1123          endif
1124          call rescale_weights(t_bath)
1125 co         write (iout,*) "Processor",me,
1126 co     &    " rescaling weights with temperature",t_bath
1127
1128          stdfp=dsqrt(2*Rb*t_bath/d_time)
1129          do i=1,ntyp
1130            stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
1131          enddo 
1132
1133 cde         write(iout,*) 'REMD after',me,t_bath
1134            time08=MPI_WTIME()
1135            if (me.eq.king .or. .not. out1file) then
1136             write(iout,*) 'REMD exchange time=',time08-time00
1137             call flush(iout)
1138            endif
1139         endif
1140       enddo
1141
1142       if (restart1file) then 
1143           if (me.eq.king .or. .not. out1file)
1144      &      write(iout,*) 'writing restart at the end of run'
1145            call write1rst(i_index)
1146       endif
1147
1148       if (traj1file) call write1traj
1149 cd debugging
1150 cdeb            call mpi_gather(ntwx_cache,1,mpi_integer,
1151 cdeb     &             icache_all,1,mpi_integer,king,
1152 cdeb     &             CG_COMM,ierr)
1153 cdeb            write(iout,'(a40,8000i8)') 
1154 cdeb     &             '  ntwx_cache after traj1file at the end',
1155 cdeb     &             (icache_all(i),i=1,nodes)
1156 cd end
1157
1158
1159 #ifdef MPI
1160       t_MD=MPI_Wtime()
1161 #else
1162       t_MD=tcpu()-tt0
1163 #endif
1164       if (me.eq.king .or. .not. out1file) then
1165        write (iout,'(//35(1h=),a10,35(1h=)/10(/a40,1pe15.5))') 
1166      &  '  Timing  ',
1167      & 'MD calculations setup:',t_MDsetup,
1168      & 'Energy & gradient evaluation:',t_enegrad,
1169      & 'Stochastic MD setup:',t_langsetup,
1170      & 'Stochastic MD step setup:',t_sdsetup,
1171      & 'MD steps:',t_MD
1172        write (iout,'(/28(1h=),a25,27(1h=))') 
1173      & '  End of MD calculation  '
1174       endif
1175       return
1176       end
1177
1178 c-----------------------------------------------------------------------
1179
1180       subroutine write1rst_oldtxt
1181       implicit real*8 (a-h,o-z)
1182       include 'DIMENSIONS'
1183       include 'mpif.h'
1184       include 'COMMON.MD'
1185       include 'COMMON.IOUNITS'
1186       include 'COMMON.REMD'
1187       include 'COMMON.SETUP'
1188       include 'COMMON.CHAIN'
1189       include 'COMMON.SBRIDGE'
1190       include 'COMMON.INTERACT'
1191                
1192       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1193      &     d_restart2(3,2*maxres*maxprocs)
1194       real t5_restart1(5)
1195       integer iret,itmp
1196
1197
1198        t5_restart1(1)=totT
1199        t5_restart1(2)=EK
1200        t5_restart1(3)=potE
1201        t5_restart1(4)=t_bath
1202        t5_restart1(5)=Uconst
1203        
1204        call mpi_gather(t5_restart1,5,mpi_real,
1205      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1206
1207
1208        do i=1,2*nres
1209          do j=1,3
1210            r_d(j,i)=d_t(j,i)
1211          enddo
1212        enddo
1213        call mpi_gather(r_d,3*2*nres,mpi_real,
1214      &           d_restart1,3*2*nres,mpi_real,king,
1215      &           CG_COMM,ierr)
1216
1217
1218        do i=1,2*nres
1219          do j=1,3
1220            r_d(j,i)=dc(j,i)
1221          enddo
1222        enddo
1223        call mpi_gather(r_d,3*2*nres,mpi_real,
1224      &           d_restart2,3*2*nres,mpi_real,king,
1225      &           CG_COMM,ierr)
1226
1227        if(me.eq.king) then
1228          open(irest1,file=mremd_rst_name,status='unknown')
1229          write (irest1,*) (i2rep(i),i=0,nodes-1)
1230          write (irest1,*) (ifirst(i),i=1,remd_m(1))
1231          do il=1,nodes
1232               write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1233               write (irest1,*) ndowna(0,il),
1234      &                   (ndowna(i,il),i=1,ndowna(0,il))
1235          enddo
1236
1237          do il=1,nodes
1238            write(irest1,*) (t_restart1(j,il),j=1,4)
1239          enddo
1240
1241          do il=0,nodes-1
1242            do i=1,2*nres
1243              write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
1244            enddo
1245          enddo
1246          do il=0,nodes-1
1247            do i=1,2*nres
1248              write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
1249            enddo
1250          enddo
1251          close(irest1)
1252        endif
1253
1254       return
1255       end
1256
1257       subroutine write1rst(i_index)
1258       implicit real*8 (a-h,o-z)
1259       include 'DIMENSIONS'
1260       include 'mpif.h'
1261       include 'COMMON.MD'
1262       include 'COMMON.IOUNITS'
1263       include 'COMMON.REMD'
1264       include 'COMMON.SETUP'
1265       include 'COMMON.CHAIN'
1266       include 'COMMON.SBRIDGE'
1267       include 'COMMON.INTERACT'
1268                
1269       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1270      &     d_restart2(3,2*maxres*maxprocs)
1271       real t5_restart1(5)
1272       integer iret,itmp
1273       integer*2 i_index
1274      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1275
1276
1277        t5_restart1(1)=totT
1278        t5_restart1(2)=EK
1279        t5_restart1(3)=potE
1280        t5_restart1(4)=t_bath
1281        t5_restart1(5)=Uconst
1282        
1283        call mpi_gather(t5_restart1,5,mpi_real,
1284      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1285
1286
1287        do i=1,2*nres
1288          do j=1,3
1289            r_d(j,i)=d_t(j,i)
1290          enddo
1291        enddo
1292        call mpi_gather(r_d,3*2*nres,mpi_real,
1293      &           d_restart1,3*2*nres,mpi_real,king,
1294      &           CG_COMM,ierr)
1295
1296
1297        do i=1,2*nres
1298          do j=1,3
1299            r_d(j,i)=dc(j,i)
1300          enddo
1301        enddo
1302        call mpi_gather(r_d,3*2*nres,mpi_real,
1303      &           d_restart2,3*2*nres,mpi_real,king,
1304      &           CG_COMM,ierr)
1305
1306        if(me.eq.king) then
1307 c         open(irest1,file=mremd_rst_name,status='unknown')
1308          call xdrfopen(ixdrf,mremd_rst_name, "w", iret)
1309 c         write (irest1,*) (i2rep(i),i=0,nodes-1)
1310          do i=0,nodes-1
1311           call xdrfint(ixdrf, i2rep(i), iret)
1312          enddo
1313 c         write (irest1,*) (ifirst(i),i=1,remd_m(1))
1314          do i=1,remd_m(1)
1315           call xdrfint(ixdrf, ifirst(i), iret)
1316          enddo
1317          do il=1,nodes
1318 c              write (irest1,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1319               do i=0,nupa(0,il)
1320                call xdrfint(ixdrf, nupa(i,il), iret)
1321               enddo
1322
1323 c              write (irest1,*) ndowna(0,il),
1324 c     &                   (ndowna(i,il),i=1,ndowna(0,il))
1325               do i=0,ndowna(0,il)
1326                call xdrfint(ixdrf, ndowna(i,il), iret)
1327               enddo
1328          enddo
1329
1330          do il=1,nodes
1331 c           write(irest1,*) (t_restart1(j,il),j=1,4)
1332            do j=1,4
1333             call xdrffloat(ixdrf, t_restart1(j,il), iret)
1334            enddo
1335          enddo
1336
1337          do il=0,nodes-1
1338            do i=1,2*nres
1339 c             write(irest1,'(3e15.5)') (d_restart1(j,i+2*nres*il),j=1,3)
1340             do j=1,3
1341              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1342             enddo
1343            enddo
1344          enddo
1345          do il=0,nodes-1
1346            do i=1,2*nres
1347 c             write(irest1,'(3e15.5)') (d_restart2(j,i+2*nres*il),j=1,3)
1348             do j=1,3
1349              call xdrffloat(ixdrf, d_restart2(j,i+2*nres*il), iret)
1350             enddo
1351            enddo
1352          enddo
1353
1354
1355              if(usampl) then
1356 c              write (irest1,*) nset
1357               call xdrfint(ixdrf, nset, iret)
1358 c              write (irest1,*) (mset(i),i=1,nset)
1359               do i=1,nset
1360                 call xdrfint(ixdrf,mset(i), iret)
1361               enddo
1362 c              write (irest1,*) (i2set(i),i=0,nodes-1)
1363               do i=0,nodes-1
1364                 call xdrfint(ixdrf,i2set(i), iret)
1365               enddo
1366 c              write (irest1,*) "i_index"
1367               do il=1,nset
1368                do il1=1,mset(il)
1369                 do i=1,nrep
1370 c                  write(irest1,*) (i_index(i,j,il,il1),j=1,remd_m(i))
1371                  do j=1,remd_m(i)
1372                    itmp=i_index(i,j,il,il1)
1373                    call xdrfint(ixdrf,itmp, iret)
1374                  enddo
1375                 enddo
1376                enddo
1377               enddo
1378            
1379              endif
1380
1381 c         close(irest1)
1382          call xdrfclose(ixdrf, iret)
1383        endif
1384
1385       return
1386       end
1387
1388
1389       subroutine write1traj
1390       implicit real*8 (a-h,o-z)
1391       include 'DIMENSIONS'
1392       include 'mpif.h'
1393       include 'COMMON.MD'
1394       include 'COMMON.IOUNITS'
1395       include 'COMMON.REMD'
1396       include 'COMMON.SETUP'
1397       include 'COMMON.CHAIN'
1398       include 'COMMON.SBRIDGE'
1399       include 'COMMON.INTERACT'
1400                
1401       real t5_restart1(5)
1402       integer iret,itmp
1403       real xcoord(3,maxres2+2),prec
1404       real r_qfrag(50),r_qpair(100)
1405       real r_utheta(50),r_ugamma(100),r_uscdiff(100)
1406       real p_qfrag(50*maxprocs),p_qpair(100*maxprocs)
1407       real p_utheta(50*maxprocs),p_ugamma(100*maxprocs),
1408      &     p_uscdiff(100*maxprocs)
1409       real p_c(3,(maxres2+2)*maxprocs),r_c(3,maxres2+2)
1410
1411       call mpi_bcast(ii_write,1,mpi_integer,
1412      &           king,CG_COMM,ierr)
1413
1414 c debugging
1415       print *,'traj1file',me,ii_write,ntwx_cache
1416 c end debugging
1417
1418       if(me.eq.king) call xdrfopen(ixdrf,cartname, "a", iret)
1419       do ii=1,ii_write
1420        t5_restart1(1)=totT_cache(ii)
1421        t5_restart1(2)=EK_cache(ii)
1422        t5_restart1(3)=potE_cache(ii)
1423        t5_restart1(4)=t_bath_cache(ii)
1424        t5_restart1(5)=Uconst_cache(ii)
1425        call mpi_gather(t5_restart1,5,mpi_real,
1426      &      t_restart1,5,mpi_real,king,CG_COMM,ierr)
1427
1428        call mpi_gather(iset_cache(ii),1,mpi_integer,
1429      &      iset_restart1,1,mpi_integer,king,CG_COMM,ierr)
1430
1431           do i=1,nfrag
1432            r_qfrag(i)=qfrag_cache(i,ii)
1433           enddo
1434           do i=1,npair
1435            r_qpair(i)=qpair_cache(i,ii)
1436           enddo
1437           do i=1,nfrag_back
1438            r_utheta(i)=utheta_cache(i,ii)
1439            r_ugamma(i)=ugamma_cache(i,ii)
1440            r_uscdiff(i)=uscdiff_cache(i,ii)
1441           enddo
1442
1443         call mpi_gather(r_qfrag,nfrag,mpi_real,
1444      &           p_qfrag,nfrag,mpi_real,king,
1445      &           CG_COMM,ierr)
1446         call mpi_gather(r_qpair,npair,mpi_real,
1447      &           p_qpair,npair,mpi_real,king,
1448      &           CG_COMM,ierr)
1449         call mpi_gather(r_utheta,nfrag_back,mpi_real,
1450      &           p_utheta,nfrag_back,mpi_real,king,
1451      &           CG_COMM,ierr)
1452         call mpi_gather(r_ugamma,nfrag_back,mpi_real,
1453      &           p_ugamma,nfrag_back,mpi_real,king,
1454      &           CG_COMM,ierr)
1455         call mpi_gather(r_uscdiff,nfrag_back,mpi_real,
1456      &           p_uscdiff,nfrag_back,mpi_real,king,
1457      &           CG_COMM,ierr)
1458
1459 #ifdef DEBUG
1460         write (iout,*) "p_qfrag"
1461         do i=1,nodes
1462           write (iout,*) i,(p_qfrag((i-1)*nfrag+j),j=1,nfrag)
1463         enddo
1464         write (iout,*) "p_qpair"
1465         do i=1,nodes
1466           write (iout,*) i,(p_qpair((i-1)*npair+j),j=1,npair)
1467         enddo
1468         call flush(iout)
1469 #endif
1470         do i=1,nres*2
1471          do j=1,3
1472           r_c(j,i)=c_cache(j,i,ii)
1473          enddo
1474         enddo
1475
1476         call mpi_gather(r_c,3*2*nres,mpi_real,
1477      &           p_c,3*2*nres,mpi_real,king,
1478      &           CG_COMM,ierr)
1479
1480        if(me.eq.king) then
1481
1482          do il=1,nodes
1483           call xdrffloat(ixdrf, real(t_restart1(1,il)), iret)
1484           call xdrffloat(ixdrf, real(t_restart1(3,il)), iret)
1485           call xdrffloat(ixdrf, real(t_restart1(5,il)), iret)
1486           call xdrffloat(ixdrf, real(t_restart1(4,il)), iret)
1487           call xdrfint(ixdrf, nss, iret) 
1488           do j=1,nss
1489            call xdrfint(ixdrf, ihpb(j), iret)
1490            call xdrfint(ixdrf, jhpb(j), iret)
1491           enddo
1492           call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
1493           call xdrfint(ixdrf, iset_restart1(il), iret)
1494           do i=1,nfrag
1495            call xdrffloat(ixdrf, p_qfrag(i+(il-1)*nfrag), iret)
1496           enddo
1497           do i=1,npair
1498            call xdrffloat(ixdrf, p_qpair(i+(il-1)*npair), iret)
1499           enddo
1500           do i=1,nfrag_back
1501            call xdrffloat(ixdrf, p_utheta(i+(il-1)*nfrag_back), iret)
1502            call xdrffloat(ixdrf, p_ugamma(i+(il-1)*nfrag_back), iret)
1503            call xdrffloat(ixdrf, p_uscdiff(i+(il-1)*nfrag_back), iret)
1504           enddo
1505           prec=10000.0
1506           do i=1,nres
1507            do j=1,3
1508             xcoord(j,i)=p_c(j,i+(il-1)*nres*2)
1509            enddo
1510           enddo
1511           do i=nnt,nct
1512            do j=1,3
1513             xcoord(j,nres+i-nnt+1)=p_c(j,i+nres+(il-1)*nres*2)
1514            enddo
1515           enddo
1516           itmp=nres+nct-nnt+1
1517           call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
1518          enddo
1519        endif
1520       enddo
1521       if(me.eq.king) call xdrfclose(ixdrf, iret)
1522       do i=1,ntwx_cache-ii_write
1523
1524             totT_cache(i)=totT_cache(ii_write+i)
1525             EK_cache(i)=EK_cache(ii_write+i)
1526             potE_cache(i)=potE_cache(ii_write+i)
1527             t_bath_cache(i)=t_bath_cache(ii_write+i)
1528             Uconst_cache(i)=Uconst_cache(ii_write+i)
1529             iset_cache(i)=iset_cache(ii_write+i)
1530
1531             do ii=1,nfrag
1532              qfrag_cache(ii,i)=qfrag_cache(ii,ii_write+i)
1533             enddo
1534             do ii=1,npair
1535              qpair_cache(ii,i)=qpair_cache(ii,ii_write+i)
1536             enddo
1537             do ii=1,nfrag_back
1538               utheta_cache(ii,i)=utheta_cache(ii,ii_write+i)
1539               ugamma_cache(ii,i)=ugamma_cache(ii,ii_write+i)
1540               uscdiff_cache(ii,i)=uscdiff_cache(ii,ii_write+i)
1541             enddo
1542
1543             do ii=1,nres*2
1544              do j=1,3
1545               c_cache(j,ii,i)=c_cache(j,ii,ii_write+i)
1546              enddo
1547             enddo
1548       enddo
1549       ntwx_cache=ntwx_cache-ii_write
1550       return
1551       end
1552
1553
1554       subroutine read1restart(i_index)
1555       implicit real*8 (a-h,o-z)
1556       include 'DIMENSIONS'
1557       include 'mpif.h'
1558       include 'COMMON.MD'
1559       include 'COMMON.IOUNITS'
1560       include 'COMMON.REMD'
1561       include 'COMMON.SETUP'
1562       include 'COMMON.CHAIN'
1563       include 'COMMON.SBRIDGE'
1564       include 'COMMON.INTERACT'
1565       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1566      &                 t5_restart1(5)
1567       integer*2 i_index
1568      &            (maxprocs/4,maxprocs/20,maxprocs/200,maxprocs/200)
1569       write (*,*) "Processor",me," called read1restart"
1570
1571          if(me.eq.king)then
1572               open(irest2,file=mremd_rst_name,status='unknown')
1573               read(irest2,*,err=334) i
1574               write(iout,*) "Reading old rst in ASCI format"
1575               close(irest2)
1576                call read1restart_old
1577                return
1578  334          continue
1579               call xdrfopen(ixdrf,mremd_rst_name, "r", iret)
1580
1581 c             read (irest2,*) (i2rep(i),i=0,nodes-1)
1582               do i=0,nodes-1
1583                call xdrfint(ixdrf, i2rep(i), iret)
1584               enddo
1585 c             read (irest2,*) (ifirst(i),i=1,remd_m(1))
1586               do i=1,remd_m(1)
1587                call xdrfint(ixdrf, ifirst(i), iret)
1588               enddo
1589              do il=1,nodes
1590 c              read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1591               call xdrfint(ixdrf, nupa(0,il), iret)
1592               do i=1,nupa(0,il)
1593                call xdrfint(ixdrf, nupa(i,il), iret)
1594               enddo
1595
1596 c              read (irest2,*) ndowna(0,il),
1597 c     &                    (ndowna(i,il),i=1,ndowna(0,il))
1598               call xdrfint(ixdrf, ndowna(0,il), iret)
1599               do i=1,ndowna(0,il)
1600                call xdrfint(ixdrf, ndowna(i,il), iret)
1601               enddo
1602              enddo
1603              do il=1,nodes
1604 c               read(irest2,*) (t_restart1(j,il),j=1,4)
1605                do j=1,4
1606                 call xdrffloat(ixdrf, t_restart1(j,il), iret)
1607                enddo
1608              enddo
1609          endif
1610          call mpi_scatter(t_restart1,5,mpi_real,
1611      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1612          totT=t5_restart1(1)              
1613          EK=t5_restart1(2)
1614          potE=t5_restart1(3)
1615          t_bath=t5_restart1(4)
1616
1617          if(me.eq.king)then
1618               do il=0,nodes-1
1619                do i=1,2*nres
1620 c                read(irest2,'(3e15.5)') 
1621 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1622             do j=1,3
1623              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1624             enddo
1625                enddo
1626               enddo
1627          endif
1628          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1629      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1630
1631          do i=1,2*nres
1632            do j=1,3
1633             d_t(j,i)=r_d(j,i)
1634            enddo
1635          enddo
1636          if(me.eq.king)then 
1637               do il=0,nodes-1
1638                do i=1,2*nres
1639 c                read(irest2,'(3e15.5)') 
1640 c     &                (d_restart1(j,i+2*nres*il),j=1,3)
1641             do j=1,3
1642              call xdrffloat(ixdrf, d_restart1(j,i+2*nres*il), iret)
1643             enddo
1644                enddo
1645               enddo
1646          endif
1647          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1648      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1649          do i=1,2*nres
1650            do j=1,3
1651             dc(j,i)=r_d(j,i)
1652            enddo
1653          enddo
1654        
1655
1656            if(usampl) then
1657              if(me.eq.king)then
1658               call xdrfint(ixdrf, nset, iret)
1659               do i=1,nset
1660                 call xdrfint(ixdrf,mset(i), iret)
1661               enddo
1662               do i=0,nodes-1
1663                 call xdrfint(ixdrf,i2set(i), iret)
1664               enddo
1665               do il=1,nset
1666                do il1=1,mset(il)
1667                 do i=1,nrep
1668                  do j=1,remd_m(i)
1669                    call xdrfint(ixdrf,itmp, iret)
1670                    i_index(i,j,il,il1)=itmp
1671                  enddo
1672                 enddo
1673                enddo
1674               enddo
1675              endif
1676
1677               call mpi_scatter(i2set,1,mpi_integer,
1678      &           iset,1,mpi_integer,king,
1679      &           CG_COMM,ierr) 
1680
1681            endif
1682
1683
1684         if(me.eq.king) close(irest2)
1685         return
1686         end
1687
1688       subroutine read1restart_old
1689       implicit real*8 (a-h,o-z)
1690       include 'DIMENSIONS'
1691       include 'mpif.h'
1692       include 'COMMON.MD'
1693       include 'COMMON.IOUNITS'
1694       include 'COMMON.REMD'
1695       include 'COMMON.SETUP'
1696       include 'COMMON.CHAIN'
1697       include 'COMMON.SBRIDGE'
1698       include 'COMMON.INTERACT'
1699       real d_restart1(3,2*maxres*maxprocs),r_d(3,2*maxres),
1700      &                 t5_restart1(5)
1701
1702          if(me.eq.king)then
1703              open(irest2,file=mremd_rst_name,status='unknown')
1704              read (irest2,*) (i2rep(i),i=0,nodes-1)
1705              read (irest2,*) (ifirst(i),i=1,remd_m(1))
1706              do il=1,nodes
1707               read (irest2,*) nupa(0,il),(nupa(i,il),i=1,nupa(0,il))
1708               read (irest2,*) ndowna(0,il),
1709      &                    (ndowna(i,il),i=1,ndowna(0,il))
1710              enddo
1711              do il=1,nodes
1712                read(irest2,*) (t_restart1(j,il),j=1,4)
1713              enddo
1714          endif
1715          call mpi_scatter(t_restart1,5,mpi_real,
1716      &           t5_restart1,5,mpi_real,king,CG_COMM,ierr)
1717          totT=t5_restart1(1)              
1718          EK=t5_restart1(2)
1719          potE=t5_restart1(3)
1720          t_bath=t5_restart1(4)
1721
1722          if(me.eq.king)then
1723               do il=0,nodes-1
1724                do i=1,2*nres
1725                 read(irest2,'(3e15.5)') 
1726      &                (d_restart1(j,i+2*nres*il),j=1,3)
1727                enddo
1728               enddo
1729          endif
1730          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1731      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1732
1733          do i=1,2*nres
1734            do j=1,3
1735             d_t(j,i)=r_d(j,i)
1736            enddo
1737          enddo
1738          if(me.eq.king)then 
1739               do il=0,nodes-1
1740                do i=1,2*nres
1741                 read(irest2,'(3e15.5)') 
1742      &                (d_restart1(j,i+2*nres*il),j=1,3)
1743                enddo
1744               enddo
1745          endif
1746          call mpi_scatter(d_restart1,3*2*nres,mpi_real,
1747      &           r_d,3*2*nres,mpi_real,king,CG_COMM,ierr)
1748          do i=1,2*nres
1749            do j=1,3
1750             dc(j,i)=r_d(j,i)
1751            enddo
1752          enddo
1753         if(me.eq.king) close(irest2)
1754         return
1755         end
1756