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