source: partem_p.f@ c021909

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

Changed asap to asa_p to avoid shadowing

The global variable asap was shadowed in partem_p by the area used to collect
results from different MPI tasks. This change takes care of it.
TODO: Change all variable names that are used to collect data from the different
MPI tasks so they end in _p.

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