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
|
---|