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