source: EXAMPLES/partem_p.f

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