source: partem_p.f@ 6650a56

Last change on this file since 6650a56 was 6650a56, checked in by baerbaer <baerbaer@…>, 14 years ago

Explicitly declare variables.

git-svn-id: ​svn+ssh://svn.berlios.de/svnroot/repos/smmp/trunk@32 26dc1dd8-5c4e-0410-9ffe-d298b4865968
  • Property mode set to 100644
File size: 22.6 KB
Line 
1!**************************************************************
2!
3! This file contains the subroutines: partem_p
4! USE WITH main_p, NOT WITH main!!!!!!
5!
6! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
7! Shura Hayryan, Chin-Ku
8! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
9! Jan H. Meinke, Sandipan Mohanty
10!
11! **************************************************************
12
13 subroutine partem_p(num_rep, nequi, nswp, nmes, nsave, newsta,
14 & switch, rep_id, partem_comm)
15!
16! PURPOSE: SIMULATION OF PROTEINS BY PARALLEL TEMPERING ALGORITHM
17! ON PARALLEL COMPUTERS USING MPI
18!
19! switch: Choses the starting configuration:
20! -1 - stretched configuration
21! 0 - don't change anything
22! 1 - random start configuration
23!
24! CALLS: addang,contacts,energy,hbond,helix,iendst,metropolis,
25! outvar,(rand),rgyr
26!
27 include 'INCL.H'
28 include 'INCP.H'
29 include 'mpif.h'
30
31 logical newsta
32 integer switch, partem_comm, rep_id, nsave
33! external rand
34 external can_weight
35
36! nequi: number of Monte Carlo sweeps for thermalization
37! nswp: number of Monte Carlo sweeps
38! nmes: number of Monte Carlo sweeps between measurments
39! newsta: .true. for new simulations, .false. for re-start
40 double precision temp, eavm, sph, geavm, gsph, dv, grnd, vr
41 double precision addang, dummy, eol, energy, acz, rmsv, rmsdfun
42 double precision rgy, ee, tmhb, dham, swp, wij, rd, e_final
43
44 integer ifrrm, nmes, nswp, num_rep, i, j, nresi, iold, inode
45 integer intem, iv, jold, idum1, idum2, idum3, mpi_integer
46 integer mpi_comm_world, ierr, mpi_double_precision, nsw, nequi
47 integer nml, nhel, mhel, nbet, mbet, mhb, imhb, nctot, ncnat
48 integer mpi_comm_null, k1, k, nu, no1, in, jn
49
50 dimension eavm(MAX_PROC),sph(MAX_PROC),intem(MAX_PROC),
51 & inode(MAX_PROC), geavm(MAX_PROC), gsph(MAX_PROC)
52 double precision pbe(MAX_PROC),yol(MAX_PROC),acy(MAX_PROC),
53 & acy1(MAX_PROC),acx1(MAX_PROC),
54 & rgyrp(MAX_PROC),rmsdp(MAX_PROC), eol0,acz0
55
56 double precision e_min, e_minp(MAX_PROC), e_minpt(MAX_PROC)
57 integer h_max, h_maxp(MAX_PROC)
58! Order of replica exchange
59 integer odd
60! Counter to keep random number generators in sync
61 integer randomCount
62
63! Collect partial energies. Only the root writes to disk. We have to
64! collect the information from the different replicas and provide
65! arrays to store them.
66! eyslr storage array for solvent energy
67! eyelp - " - coulomb energy
68! eyvwp - " - van-der-Waals energy
69! eyhbp - " - hydrogen bonding energy
70! eysmi - " - intermolecular interaction energy
71! eyabp - " - Abagyan correction term
72 double precision eyslr(MAX_PROC)
73 double precision eyelp(MAX_PROC),eyvwp(MAX_PROC),eyhbp(MAX_PROC),
74 & eyvrp(MAX_PROC),eysmip(MAX_PROC), eyabp(MAX_PROC)
75! Collect information about accessible surface and van-der-Waals volume
76! asap storage array for solvent accessible surface
77! vdvolp storage array for van-der-Waals volume
78 double precision asa_p(MAX_PROC), vdvolp(MAX_PROC)
79
80 integer nhelp(MAX_PROC),nbetp(MAX_PROC), mhbp(MAX_PROC),
81 & ncnatp(MAX_PROC),nctotp(MAX_PROC)
82 integer imhbp(MAX_PROC)
83 character*80 filebase, fileNameMP, tbase0,tbase1
84! frame frame number for writing configurations
85! trackID configuration that should be tracked and written out
86! dir direction in random walk
87! -1 - visited highest temperature last
88! 1 - visited lowest temperature last
89! 0 - haven't visited the boundaries yet.
90! dirp storage array for directions.
91 integer frame, trackID, dir
92 integer dirp(MAX_PROC)
93
94 frame = ifrrm
95 trackID = 1
96 odd = 1
97 write (*,*) 'Starting parallel tempering.'
98 write (*,*) 'parameters, ',switch,newsta,nmes,nswp,nmes,
99 & rep_id, num_rep, partem_comm, myrank
100 call flush(6)
101!
102!
103! File with temperatures
104 open(11,file='temperatures',status='old')
105! File with reference conformation
106 tbase0='trj_00000'
107 open(18,file=fileNameMP(tbase0,5,9,rep_id),status='unknown')
108 if (rep_id.eq.0.and.myrank.eq.0) then
109! File with time series of simulation
110 open(14,file='ts.d',status='unknown')
111! Track weights
112! open(16, file='weights.dat', status='unknown')
113 endif
114
115! READ IN TEMPERATURES
116 do i=1,num_rep
117 read(11,*) j,temp
118 pbe(j) = 1.0d0/( temp * 1.98773d-3 )
119 end do
120 close(11)
121
122! nresi: number of residues
123 nresi=irsml2(1)-irsml1(1)+1
124!
125! Initialize variables
126 do i=1,num_rep
127 acx1(i) = 0.0d0
128 acy(i) = 0.0d0
129 eavm(i) = 0.0d0
130 sph(i) = 0.0d0
131 geavm(i) =0.0d0
132 gsph(i) = 0.0d0
133 e_minp(i) = 1.0d15
134 h_maxp(i) = 0
135 dirp(i) = 0
136 end do
137 dirp(1) = 1
138 dirp(num_rep) = -1
139 e_min = 1.0d15
140 h_max = 0
141 dir = dirp(rep_id + 1)
142
143! _________________________________ Initialize Variables
144 if(newsta) then
145 iold=0
146 do i=1,num_rep
147 inode(i) = i
148 intem(i) = i
149 end do
150! _________________________________ initialize starting configuration
151 if (switch.ne.0) then
152 do i=1,nvr
153 iv=idvr(i) ! provides index of non-fixed variable
154 if (switch.gt.0) then
155 dv=axvr(i)*(grnd()-0.5)
156 vr=addang(pi,dv)
157 else
158 vr = pi
159 endif
160 vlvr(iv)=vr
161 enddo
162 endif
163 else
164 if(rep_id.eq.0.and.myrank.eq.0) then
165 open(13,file='par_R.in', status='unknown')
166 read(13,*) iold
167 do i=1,num_rep
168 read(13,*) j,inode(i),intem(i),yol(i),e_minp(i),h_maxp(i)
169 write (*,*) "par_R.in:",i,j
170 end do
171 jold=(iold/nmes)*num_rep
172 rewind 14
173 do i=1,jold
174 read(14,*) idum1,idum2,idum3,dummy, dir
175 & ,dummy, dummy, dummy, dummy, dummy
176 & ,dummy, dummy, dummy, dummy
177 & ,dummy, idum1, idum2, idum3
178 & ,idum1, idum2, idum3, e_min
179 & ,dummy, dummy
180 write (*,*) i
181 call flush(6)
182 end do
183 close(13)
184 end if
185 CALL MPI_BCAST(IOLD,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
186 CALL MPI_BCAST(INTEM,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
187 CALL MPI_BCAST(INODE,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
188 CALL MPI_BCAST(YOL,num_rep,MPI_DOUBLE_PRECISION,0,
189 & MPI_COMM_WORLD,IERR)
190 CALL MPI_BCAST(E_MINP, num_rep, MPI_DOUBLE_PRECISION, 0,
191 & MPI_COMM_WORLD, IERR)
192 CALL MPI_BCAST(h_maxp,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
193 & IERR)
194 end if
195
196 BETA = pbe(inode(rep_id+1))
197 e_min = e_minp(rep_id+1)
198 h_max = h_maxp(rep_id+1)
199 write (*,*) "E_min=",e_min," for ", rep_id + 1
200 eol=energy()
201 if(.not.newsta.and.abs(yol(rep_id + 1) - eol).gt.0.1) then
202 write(*,*) rep_id, ' Warning: yol(rep_id).ne.eol:'
203 write(*,*) rep_id, yol(rep_id + 1), eol
204 endif
205! Start of simulation
206 write (*,*) '[',rep_id, myrank, beta, partem_comm,
207 & '] Energy before equilibration:', eol
208! =====================Equilibration by canonical Metropolis
209 do nsw=1,nequi
210 call metropolis(eol,acz,can_weight)
211 end do
212 CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
213 write (*,*) '[',rep_id,'] Energy after equilibration:', eol
214 call flush(6)
215!
216!======================Multiple Markov Chains
217 acz = 0
218 do nsw=1,nswp
219!------------First ordinary Metropolis
220 call metropolis(eol,acz,can_weight)
221 iold = iold + 1
222 eol0 = eol
223 if (myrank.eq.0.and.rep_id.eq.0) then
224 write (*,*) "Finished sweep", nsw
225 call flush(6)
226 endif
227 if(mod(iold,nmes).eq.0) then
228 if ((rep_id + 1).eq.trackID.and.myrank.eq.0) then
229 frame = iold /nmes
230 filebase = "frame_00000.pdb"
231 call outpdb(0, fileNameMP(filebase, 7, 11, frame))
232 endif
233 acz0 = acz
234! Evaluate RMSD
235 nml = 1
236 rmsv = rmsdfun(nml,irsml1(nml),irsml2(nml),ixatp,xatp,yatp, &
237 & zatp,0)
238! print *,myrank,'received RMSD,energy ',rmsv,eyab,beta
239! Measure global radius of gyration
240 call rgyr(0,rgy,ee)
241 rgyp = rgy
242! Measure Helicity and Sheetness
243 call helix(nhel,mhel,nbet,mbet)
244! Measure Number of hydrogen bonds
245 mhb = 0
246 do i = 1, ntlml
247 call hbond(i,tmhb,-1)
248 mhb = mhb + 1
249 enddo
250 call interhbond(imhb)
251! Measure total number of contacts (NCTOT) and number of
252! native contacts (NCNAT)
253 call contacts(nctot,ncnat,dham)
254! Add tracking of lowest energy configuration
255 if (eol.lt.e_min) then
256! Write out configuration
257 i=rep_id+1
258 j=inode(i)
259 e_min = eol
260 filebase = "c_emin_0000.pdb"
261 call outpdb(0, fileNameMP(filebase, 8, 11, i))
262 filebase = "c_emin_0000.var"
263 call outvar(0, fileNameMP(filebase, 8, 11, i))
264 filebase = "c_emin_0000.dat"
265 open(15, file=fileNameMP(filebase, 8, 11, i),
266 & status="unknown")
267! write(15,'(i8,2i4,f6.2,2f8.2,5i8)') iold,i,j,pbe(i),
268 write(15,*) iold,j,i,beta,
269 & eol, eyab, eysl, eyel, eyvw, eyhb, eyvr, eysmi,asa,
270 & vdvol, rgy, nhel, nbet, mhb, imhb, nctot,ncnat
271 close(15)
272 endif
273! Add tracking of configuration with larges hydrogen contents.
274 if ((mhb + imhb).gt.h_max) then
275! Write out configuration
276 i = rep_id + 1
277 j = inode(i)
278 h_max = mhb + imhb
279 filebase = "c_hmax_0000.pdb"
280 call outpdb(0,fileNameMP(filebase,8,11,i))
281 filebase = "c_hmax_0000.var"
282 call outvar(0,fileNameMP(filebase,8,11,i))
283 filebase = "c_hmax_0000.dat"
284 open(15, file=fileNameMP(filebase, 8, 11, i),
285 & status="unknown")
286! write(15,'(i8,2i4,f6.2,2f8.2,5i8)') iold,i,j,pbe(i),
287 write(15,*) iold,j,i,beta,
288 & eol, eyab, eysl, eyel, eyvw, eyhb, eyvr, eysmi,asa,
289 & vdvol, rgy, nhel, nbet, mhb, imhb, nctot,ncnat
290 close(15)
291 endif
292
293!
294!--------------------Gather measurement data
295! I only use the master node of each replica for data collection. The
296! variable partem_comm provides the appropriate communicator.
297 if (partem_comm.ne.MPI_COMM_NULL) then
298 CALL MPI_GATHER(rmsv,1,MPI_DOUBLE_PRECISION,rmsdp,1,
299 & MPI_DOUBLE_PRECISION, 0,partem_comm,IERR)
300 CALL MPI_GATHER(eyab,1,MPI_DOUBLE_PRECISION,eyabp,1,
301 & MPI_DOUBLE_PRECISION, 0,partem_comm,IERR)
302 CALL MPI_GATHER(RGYP,1,MPI_DOUBLE_PRECISION,RGYRP,1,
303 & MPI_DOUBLE_PRECISION, 0,partem_comm,IERR)
304 CALL MPI_GATHER(NHEL,1,MPI_INTEGER,NHELP,1,MPI_INTEGER,
305 & 0,partem_comm,IERR)
306 CALL MPI_GATHER(NBET,1,MPI_INTEGER,NBETP,1,MPI_INTEGER,
307 & 0,partem_comm,IERR)
308 CALL MPI_GATHER(MHB,1,MPI_INTEGER,MHBP,1,MPI_INTEGER,
309 & 0,partem_comm,IERR)
310 CALL MPI_GATHER(iMHB,1,MPI_INTEGER,iMHBP,1,MPI_INTEGER,
311 & 0,partem_comm,IERR)
312 CALL MPI_GATHER(NCTOT,1,MPI_INTEGER,NCTOTP,1,MPI_INTEGER,
313 & 0,partem_comm,IERR)
314 CALL MPI_GATHER(NCNAT,1,MPI_INTEGER,NCNATP,1,MPI_INTEGER,
315 & 0,partem_comm,IERR)
316 CALL MPI_GATHER(dir,1,MPI_INTEGER,dirp,1,MPI_INTEGER,
317 & 0,partem_comm,IERR)
318 CALL MPI_GATHER(acz0,1,MPI_DOUBLE_PRECISION,acy1,1,
319 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
320 CALL MPI_GATHER(e_min,1,MPI_DOUBLE_PRECISION,e_minp,1,
321 & MPI_DOUBLE_PRECISION,0, partem_comm,IERR)
322 CALL MPI_GATHER(EOL0,1,MPI_DOUBLE_PRECISION,YOL,1,
323 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
324 CALL MPI_GATHER(eysl,1,MPI_DOUBLE_PRECISION,eyslr,1,
325 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
326 CALL MPI_GATHER(eyel,1,MPI_DOUBLE_PRECISION,eyelp,1,
327 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
328 CALL MPI_GATHER(eyvw,1,MPI_DOUBLE_PRECISION,eyvwp,1,
329 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
330 CALL MPI_GATHER(eyhb,1,MPI_DOUBLE_PRECISION,eyhbp,1,
331 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
332 CALL MPI_GATHER(eyvr,1,MPI_DOUBLE_PRECISION,eyvrp,1,
333 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
334 CALL MPI_GATHER(eysmi,1,MPI_DOUBLE_PRECISION,eysmip,1,
335 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
336 CALL MPI_GATHER(asa,1,MPI_DOUBLE_PRECISION,asa_p,1,
337 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
338 CALL MPI_GATHER(vdvol,1,MPI_DOUBLE_PRECISION,vdvolp,1,
339 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
340
341! CALL MPI_GATHER(EOL0,1,MPI_DOUBLE_PRECISION,YOL,1,MPI_DOUBLE_PRECISION,
342! & 0,MPI_COMM_WORLD,IERR)
343! CALL MPI_GATHER(E_MIN, 1, MPI_DOUBLE_PRECISION, E_MINP, MPI_DOUBLE_PRECISION,
344! & 0,MPI_COMM_WORLD, IERR)
345
346! Write trajectory
347 write (18,*) '@@@',iold,inode(rep_id+1)
348 call outvbs(0,18)
349 write (18,*) '###'
350! call flush(18)
351! Write current configuration
352 if ((mod(iold, nsave).eq.0)) then
353 filebase = "conf_0000.var"
354 call outvar(0, fileNameMP(filebase, 6, 9, rep_id+1))
355 endif
356 endif
357
358 if(rep_id.eq.0.and.myrank.eq.0) then
359 randomCount = 0
360! Update acceptance, temperature wise average of E and E^2 used to calculate
361! specific heat.
362 do i=1,num_rep
363 j=intem(i)
364 acy(i)=0.0
365! Above: contents of acy1 are added to acy(i) a few lines down.
366! acy1(intem(i)) contains information received from the node at temperature
367! i, on how many updates have been accepted in node intem(i). Since acz
368! is not reset to 0 every cycle, acy(i) must be set to 0 here. Else, there
369! will be serious double counting and the values of acceptance printed
370! will be simply wrong.
371 end do
372 do i=1, num_rep
373 j=intem(i)
374 acy(i)=acy(i)+acy1(j)
375 eavm(i)= eavm(i)+yol(j)
376 sph(i) = sph(i)+yol(j)*yol(j)
377 enddo
378
379
380! Write measurements to the time series file ts.d
381 do i=1,num_rep
382 j=intem(i)
383 write(14,*) iold,i,j,pbe(i), dirp(j),
384 & yol(j),eyslr(j), eyelp(j), eyvwp(j), eyhbp(j),
385 & eyvrp(j),eysmip(j), asa_p(j), vdvolp(j),
386 & rgyrp(j),nhelp(j),nbetp(j),mhbp(j),
387 & imhbp(j), nctotp(j),ncnatp(j), e_minp(j),
388 & eyabp(j),rmsdp(j)
389! call flush(14)
390 end do
391! Write the current parallel tempering information into par_R.in
392! timeLeft = llwrem(2) ! Time left till hard limit
393! if ((mod(iold, nsave).eq.0).or.(timeLeft.lt.minTimeLeft)
394 if ((mod(iold, nsave).eq.0))
395 & then
396 open(13,file='par_R.in', status='unknown')
397 write(13,*) iold
398 do i=1,num_rep
399 write(13,*) i,inode(i),intem(i),yol(i),e_minp(i),
400 & h_maxp(i)
401 end do
402! -------------------------- Various statistics of current run
403! swp=nswp-nequi
404 swp=nsw
405 write(13,*) 'Acceptance rate for change of chains:'
406 do k1=1,num_rep
407 temp=1.0d0/pbe(k1)/0.00198773
408 write(13,*) temp, acx1(k1)*2.0d0*nmes/swp
409! Above: it's the acceptance rate of exchange of replicas. Since a
410! replica exchange is attempted only once every nmes sweeps, the
411! rate should be normalized with (nmes/swp).
412 end do
413 write(13,*)
414 do k1=1,num_rep
415 k = intem(k1)
416 temp=1.0d0/pbe(k1)/0.00198773
417 beta = pbe(k1)
418 geavm(k1) = nmes*eavm(k1)/swp
419 gsph(k1) = (nmes*sph(k1)/swp-geavm(k1)**2)
420 & *beta*beta/nresi
421 write(13,'(a,2f9.2,i4,f12.3)')
422 & 'Temperature, Node,local acceptance rate:',
423 & beta,temp,k,acy(k1)/dble(nsw*nvr)
424! Above: Changed (nswp-nequi) in the denominator of acceptance as
425! acceptance values are initialized to 0 after equilibration cycles are
426! finished. Note also that since this is being written in the middle of
427! the simulation, it is normalized to nsw instead of nswp.
428 write(13,'(a,3f12.2)')
429 & 'Last Energy, Average Energy, Spec. Heat:',
430 & yol(k),geavm(k1),gsph(k1)
431 write(13,*)
432 end do
433 close(13)
434! Finally, flush the time series and trajectory files to ensure that we can do
435! a proper restart.
436 call flush(14)
437 call flush(18)
438 end if
439
440!--------------------Parallel Tempering update
441! Swap with right neighbor (odd, even)
442 if(odd.eq.1) then
443 nu=1
444 no1 = num_rep-1
445! Swap with left neighbor (even, odd)
446 else
447 nu = 2
448 no1 = num_rep
449 end if
450 do i=nu,no1,2
451 j=i+1
452! Periodic bc for swaps
453 if(i.eq.num_rep) j=1
454 in=intem(i)
455 jn=intem(j)
456 wij=exp(-pbe(i)*yol(jn)-pbe(j)*yol(in)
457 & +pbe(i)*yol(in)+pbe(j)*yol(jn))
458! The random number generator is getting out of sync here, because
459! the swap is only done on node 0!
460! Keep track of number of random numbers used.
461 rd=grnd()
462 randomCount = randomCount + 1
463! write (16,*) '>', iold, i,j
464! & ,pbe(i),yol(in), pbe(j), yol(jn), wij, rd
465 if(wij.ge.rd) then
466! Next line: Replica exchange only happens after equilibration,
467! which takes place outside this loop over nsw. So, I think nsw.gt.nequi
468! is irrelevant for the calculation of acceptance of replica exchanges.
469! /Sandipan
470! if(nsw.gt.nequi)
471 acx1(i) = acx1(i)+1
472 intem(i) = jn
473 intem(j) = in
474 inode(in)= j
475 inode(jn)= i
476 end if
477 end do
478! ---------------- End Loop over nodes which creates a new temperature
479! map for all nodes, at the node with rank 0.
480!
481 odd = 1 - odd
482 end if
483! End of "if (myrank.eq.0) ...". The block above includes PT update and
484! writing of observables into the time series file etc.
485
486! Below: Communicate new temperature-node map to all nodes
487 CALL MPI_BCAST(INTEM,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
488 & IERR)
489 CALL MPI_BCAST(INODE,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
490 & IERR)
491! Synchronize random number generators for replica 0
492 if (rep_id.eq.0) then
493 CALL MPI_BCAST(randomCount,1,MPI_INTEGER,0,my_mpi_comm,
494 & IERR)
495 if (myrank.ne.0) then
496! write (*,*) '[', myrank,'] Missed', randomCount,
497! & 'random numbers.'
498 do i = 1, randomCount
499 rd = grnd()
500! write (*,*) '[', myrank,'] rd=', rd
501 enddo
502 endif
503 endif
504
505 BETA=PBE(INODE(rep_id+1))
506 if (INODE(rep_id + 1).eq.1) dir = 1
507 if (INODE(rep_id + 1).eq.num_rep) dir = -1
508
509 endif
510! End of "if (mod(iold,nmes).eq.0) ..."
511 end do
512!-----------End Loop over sweeps
513!
514! OUTPUT:
515!--------------------For Re-starts:
516 nu = rep_id + 1
517 filebase = "conf_0000.var"
518 call outvar(0, fileNameMP(filebase, 6, 9, nu))
519 e_final=energy()
520 if (partem_comm.ne.MPI_COMM_NULL) then
521 write (*,*) rep_id, ' E_final', e_final
522 endif
523 eol0 = eol
524 acz0 = acz
525 if (partem_comm.ne.MPI_COMM_NULL) then
526 CALL MPI_GATHER(EOL0,1,MPI_DOUBLE_PRECISION,YOL,1,
527 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
528 CALL MPI_GATHER(acz0,1,MPI_DOUBLE_PRECISION,acy1,1,
529 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
530 endif
531
532 if(rep_id.eq.0.and.myrank.eq.0) then
533 close(14)
534 open(13,file='par_R.in', status='unknown')
535 write(13,*) iold
536 do i=1,num_rep
537 write(13,*) i,inode(i),intem(i),yol(i),e_minp(i),h_maxp(i)
538 end do
539! -------------------------- Various statistics of current run
540 swp=nswp
541 write(13,*) 'Acceptance rate for change of chains:'
542 do k1=1,num_rep
543 temp=1.0d0/pbe(k1)/0.00198773
544 write(13,*) temp, acx1(k1)*2.0d0*nmes/swp
545 end do
546 write(13,*)
547 do k1=1,num_rep
548 k = intem(k1)
549 temp=1.0d0/pbe(k1)/0.00198773
550 beta = pbe(k1)
551 geavm(k1) = nmes*eavm(k1)/swp
552 gsph(k1) = (nmes*sph(k1)/swp-geavm(k1)**2)*beta*beta/nresi
553 write(13,'(a,2f9.2,i4,f12.3)')
554 & 'Temperature, Node,local acceptance rate:',
555 & beta,temp,k,acy(k1)/dble((nswp)*nvr)
556 write(13,'(a,3f12.2)')
557 & 'Last Energy, Average Energy, Spec. Heat:',
558 & yol(k),geavm(k1),gsph(k1)
559 write(13,*)
560 end do
561 close(13)
562! close(16)
563 end if
564 close(18)
565
566! =====================
567 CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
568
569 return
570
571 end
Note: See TracBrowser for help on using the repository browser.