source: partem_p.f@ bd2278d

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

Reformatting comments and continuation marks.

Fortran 90 and higher use ! to mark comments no matter where they are in the
code. The only valid continuation marker is &.
I also added the SMMP.kdevelop.filelist to the repository to make it easier
to use kdevelop.

git-svn-id: ​svn+ssh://svn.berlios.de/svnroot/repos/smmp/trunk@12 26dc1dd8-5c4e-0410-9ffe-d298b4865968
  • Property mode set to 100644
File size: 22.5 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,rgyp,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 asap(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(inode(rep_id+1))
189 h_max = h_maxp(inode(rep_id+1))
190 write (*,*) "E_min=",e_min," for ", intem(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, j))
253 filebase = "c_emin_0000.var"
254 call outvar(0, fileNameMP(filebase, 8, 11, j))
255 filebase = "c_emin_0000.dat"
256 open(15, file=fileNameMP(filebase, 8, 11, j),
257 & status="unknown")
258! write(15,'(i8,2i4,f6.2,2f8.2,5i8)') iold,i,j,pbe(i),
259 write(15,*) iold,i,j,pbe(i),
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,j))
272 filebase = "c_hmax_0000.var"
273 call outvar(0,fileNameMP(filebase,8,11,j))
274 filebase = "c_hmax_0000.dat"
275 open(15, file=fileNameMP(filebase, 8, 11, j),
276 & status="unknown")
277! write(15,'(i8,2i4,f6.2,2f8.2,5i8)') iold,i,j,pbe(i),
278 write(15,*) iold,i,j,pbe(i),
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,asap,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 e_minpt(i)=e_minp(intem(i))
363 end do
364 do i=1, num_rep
365 e_minp(i) = e_minpt(i)
366 j=intem(i)
367 acy(i)=acy(i)+acy1(j)
368 eavm(i)= eavm(i)+yol(j)
369 sph(i) = sph(i)+yol(j)*yol(j)
370 enddo
371
372
373! Write measurements to the time series file ts.d
374 do i=1,num_rep
375 j=intem(i)
376 write(14,*) iold,i,j,pbe(i), dirp(j),
377 & yol(j),eyslr(j), eyelp(j), eyvwp(j), eyhbp(j),
378 & eyvrp(j),eysmip(j), asap(j), vdvolp(j),
379 & rgyrp(j),nhelp(j),nbetp(j),mhbp(j),
380 & imhbp(j), nctotp(j),ncnatp(j), e_minp(i),
381 & eyabp(j),rmsdp(j)
382! call flush(14)
383 end do
384! Write the current parallel tempering information into par_R.in
385! timeLeft = llwrem(2) ! Time left till hard limit
386! if ((mod(iold, nsave).eq.0).or.(timeLeft.lt.minTimeLeft)
387 if ((mod(iold, nsave).eq.0))
388 & then
389 open(13,file='par_R.in', status='unknown')
390 write(13,*) iold
391 do i=1,num_rep
392 write(13,*) i,inode(i),intem(i),yol(i),e_minp(i),
393 & h_maxp(i)
394 end do
395! -------------------------- Various statistics of current run
396! swp=nswp-nequi
397 swp=nsw
398 write(13,*) 'Acceptance rate for change of chains:'
399 do k1=1,num_rep
400 temp=1.0d0/pbe(k1)/0.00198773
401 write(13,*) temp, acx1(k1)*2.0d0*nmes/swp
402! Above: it's the acceptance rate of exchange of replicas. Since a
403! replica exchange is attempted only once every nmes sweeps, the
404! rate should be normalized with (nmes/swp).
405 end do
406 write(13,*)
407 do k1=1,num_rep
408 k = intem(k1)
409 temp=1.0d0/pbe(k1)/0.00198773
410 beta = pbe(k1)
411 geavm(k1) = nmes*eavm(k1)/swp
412 gsph(k1) = (nmes*sph(k1)/swp-geavm(k1)**2)
413 & *beta*beta/nresi
414 write(13,'(a,2f9.2,i4,f12.3)')
415 & 'Temperature, Node,local acceptance rate:',
416 & beta,temp,k,acy(k1)/dble(nsw*nvr)
417! Above: Changed (nswp-nequi) in the denominator of acceptance as
418! acceptance values are initialized to 0 after equilibration cycles are
419! finished. Note also that since this is being written in the middle of
420! the simulation, it is normalized to nsw instead of nswp.
421 write(13,'(a,3f12.2)')
422 & 'Last Energy, Average Energy, Spec. Heat:',
423 & yol(k),geavm(k1),gsph(k1)
424 write(13,*)
425 end do
426 close(13)
427! Finally, flush the time series and trajectory files to ensure that we can do
428! a proper restart.
429 call flush(14)
430 call flush(18)
431 end if
432
433!--------------------Parallel Tempering update
434! Swap with right neighbor (odd, even)
435 if(odd.eq.1) then
436 nu=1
437 no1 = num_rep-1
438! Swap with left neighbor (even, odd)
439 else
440 nu = 2
441 no1 = num_rep
442 end if
443 do i=nu,no1,2
444 j=i+1
445! Periodic bc for swaps
446 if(i.eq.num_rep) j=1
447 in=intem(i)
448 jn=intem(j)
449 wij=exp(-pbe(i)*yol(jn)-pbe(j)*yol(in)
450 & +pbe(i)*yol(in)+pbe(j)*yol(jn))
451! The random number generator is getting out of sync here, because
452! the swap is only done on node 0!
453! Keep track of number of random numbers used.
454 rd=grnd()
455 randomCount = randomCount + 1
456! write (16,*) '>', iold, i,j
457! & ,pbe(i),yol(in), pbe(j), yol(jn), wij, rd
458 if(wij.ge.rd) then
459! Next line: Replica exchange only happens after equilibration,
460! which takes place outside this loop over nsw. So, I think nsw.gt.nequi
461! is irrelevant for the calculation of acceptance of replica exchanges.
462! /Sandipan
463! if(nsw.gt.nequi)
464 acx1(i) = acx1(i)+1
465 intem(i) = jn
466 intem(j) = in
467 inode(in)= j
468 inode(jn)= i
469 end if
470 end do
471! ---------------- End Loop over nodes which creates a new temperature
472! map for all nodes, at the node with rank 0.
473!
474 odd = 1 - odd
475 end if
476! End of "if (myrank.eq.0) ...". The block above includes PT update and
477! writing of observables into the time series file etc.
478
479! Below: Communicate new temperature-node map to all nodes
480 CALL MPI_BCAST(INTEM,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
481 & IERR)
482 CALL MPI_BCAST(INODE,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
483 & IERR)
484 CALL MPI_BCAST(E_MINP,num_rep,MPI_DOUBLE_PRECISION,0,
485 & MPI_COMM_WORLD,IERR)
486 CALL MPI_BCAST(H_MAXP,num_rep,MPI_INTEGER,0,MPI_COMM_WORLD,
487 & IERR)
488! Synchronize random number generators for replica 0
489 if (rep_id.eq.0) then
490 CALL MPI_BCAST(randomCount,1,MPI_INTEGER,0,my_mpi_comm,
491 & IERR)
492 if (myrank.ne.0) then
493! write (*,*) '[', myrank,'] Missed', randomCount,
494! & 'random numbers.'
495 do i = 1, randomCount
496 rd = grnd()
497! write (*,*) '[', myrank,'] rd=', rd
498 enddo
499 endif
500 endif
501
502 BETA=PBE(INODE(rep_id+1))
503 e_min = e_minp(inode(rep_id+1))
504 h_max = h_maxp(inode(rep_id+1))
505 if (INODE(rep_id + 1).eq.1) dir = 1
506 if (INODE(rep_id + 1).eq.num_rep) dir = -1
507
508 endif
509! End of "if (mod(iold,nmes).eq.0) ..."
510 end do
511!-----------End Loop over sweeps
512!
513! OUTPUT:
514!--------------------For Re-starts:
515 nu = rep_id + 1
516 filebase = "conf_0000.var"
517 call outvar(0, fileNameMP(filebase, 6, 9, nu))
518 e_final=energy()
519 if (partem_comm.ne.MPI_COMM_NULL) then
520 write (*,*) rep_id, ' E_final', e_final
521 endif
522 eol0 = eol
523 acz0 = acz
524 if (partem_comm.ne.MPI_COMM_NULL) then
525 CALL MPI_GATHER(EOL0,1,MPI_DOUBLE_PRECISION,YOL,1,
526 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
527 CALL MPI_GATHER(acz0,1,MPI_DOUBLE_PRECISION,acy1,1,
528 & MPI_DOUBLE_PRECISION,0,partem_comm,IERR)
529 endif
530
531 if(rep_id.eq.0.and.myrank.eq.0) then
532 close(14)
533 open(13,file='par_R.in', status='unknown')
534 write(13,*) iold
535 do i=1,num_rep
536 write(13,*) i,inode(i),intem(i),yol(i),e_minp(i),h_maxp(i)
537 end do
538! -------------------------- Various statistics of current run
539 swp=nswp
540 write(13,*) 'Acceptance rate for change of chains:'
541 do k1=1,num_rep
542 temp=1.0d0/pbe(k1)/0.00198773
543 write(13,*) temp, acx1(k1)*2.0d0*nmes/swp
544 end do
545 write(13,*)
546 do k1=1,num_rep
547 k = intem(k1)
548 temp=1.0d0/pbe(k1)/0.00198773
549 beta = pbe(k1)
550 geavm(k1) = nmes*eavm(k1)/swp
551 gsph(k1) = (nmes*sph(k1)/swp-geavm(k1)**2)*beta*beta/nresi
552 write(13,'(a,2f9.2,i4,f12.3)')
553 & 'Temperature, Node,local acceptance rate:',
554 & beta,temp,k,acy(k1)/dble((nswp)*nvr)
555 write(13,'(a,3f12.2)')
556 & 'Last Energy, Average Energy, Spec. Heat:',
557 & yol(k),geavm(k1),gsph(k1)
558 write(13,*)
559 end do
560 close(13)
561! close(16)
562 end if
563 close(18)
564
565! =====================
566 CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
567
568 return
569
570 end
Note: See TracBrowser for help on using the repository browser.