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