1 | ! **************************************************************
|
---|
2 | !
|
---|
3 | ! This file contains the module multicanonical with the
|
---|
4 | ! subroutines: mulcan_par, muca_weight, mulcan_sim, and
|
---|
5 | ! mulcan_weight2
|
---|
6 | !
|
---|
7 | ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
8 | ! Shura Hayryan, Chin-Ku Hu
|
---|
9 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
10 | ! Jan H. Meinke, Sandipan Mohanty
|
---|
11 | !
|
---|
12 | ! Changed to module by Jan H. Meinke
|
---|
13 | !
|
---|
14 | ! **************************************************************
|
---|
15 |
|
---|
16 | module multicanonical
|
---|
17 |
|
---|
18 | real*8, private, allocatable :: b(:), alpha(:)
|
---|
19 | real*8, private :: ebin, beta
|
---|
20 | real*8, private :: xmin, xmax
|
---|
21 |
|
---|
22 | contains
|
---|
23 | subroutine mulcan_par(nsweep, nup, temp, kmin, kmax, binWidth, l_iter)
|
---|
24 | !
|
---|
25 | ! PURPOSE: CALCULATION OF ESTIMATORS FOR MULTICANONICAL WEIGHTS
|
---|
26 | ! USING A METHOD DESCRIBED FIRST IN: Bernd Berg,
|
---|
27 | ! {\it J.~Stat.~Phys.} {\bf 82}, 331~(1996). THE METHOD
|
---|
28 | ! IS STABLE, BUT MAY NOT ALWAYS BE THE FASTEST WAY OF
|
---|
29 | ! GETTING PARAMETERS
|
---|
30 | !
|
---|
31 | ! CALLS: addang,energy,metropolis
|
---|
32 | !
|
---|
33 | include 'INCL.H'
|
---|
34 | !
|
---|
35 | ! l_iter: .true. = for iteration of multicanonical parameters
|
---|
36 | ! .false. = starts new calculation of multicanonica parameters
|
---|
37 | ! kmin,kmax: Range of multicanonical parameter
|
---|
38 | ! ebin: bin size for multicanonical parameter
|
---|
39 | ! nup: Number of sweeps between updates of multicanonical parameters
|
---|
40 | ! nsweep: Number of sweeps for simulation run
|
---|
41 | ! temp: temperature for initial canonical run
|
---|
42 | !
|
---|
43 | integer, intent(in) :: kmin, kmax, nsweep, nup
|
---|
44 | real*8, intent(in) :: binWidth, temp
|
---|
45 | logical, intent(in) :: l_iter
|
---|
46 |
|
---|
47 | integer :: i,j
|
---|
48 | integer :: iv,nsw,muold
|
---|
49 | real*8 :: bi,bi_1,vr,g0,dv
|
---|
50 | real*8 :: acz
|
---|
51 |
|
---|
52 | real*8, allocatable :: xhist(:),g1(:),ent(:)
|
---|
53 | integer, allocatable :: ihist(:)
|
---|
54 |
|
---|
55 | xmin = kmin
|
---|
56 | xmax = kmax
|
---|
57 | ebin = binWidth
|
---|
58 |
|
---|
59 | ! Allocate the arrays
|
---|
60 | allocate(b(kmin:kmax), alpha(kmin:kmax))
|
---|
61 | allocate(xhist(kmin:kmax), g1(kmin:kmax), ent(kmin:kmax) )
|
---|
62 | allocate(ihist(kmin:kmax))
|
---|
63 |
|
---|
64 | if (.not.(allocated(b).and.allocated(alpha).and.allocated(xhist) &
|
---|
65 | .and.allocated(g1).and.allocated(ent).and.allocated(ihist))) then
|
---|
66 | stop "Unable to allocate memory in mulcan_par. Exiting."
|
---|
67 | end if
|
---|
68 |
|
---|
69 | ! Files with multicanonical parameter
|
---|
70 | open(11,file='mpar_full.d')
|
---|
71 | open(12,file='muca.d')
|
---|
72 |
|
---|
73 | beta=1.0/ ( temp * 1.98773d-3 )
|
---|
74 | do j=kmin,kmax
|
---|
75 | ihist(j)= 0
|
---|
76 | end do
|
---|
77 |
|
---|
78 |
|
---|
79 | if(l_iter) then
|
---|
80 | ! READ IN FIELDS WITH MULTICANONICAL PARAMETER
|
---|
81 | Do j=kmin,kmax
|
---|
82 | read(11,*) i,xhist(i),g1(i),b(i),alpha(i),ent(i)
|
---|
83 | end do
|
---|
84 | else
|
---|
85 | do j=kmin,kmax
|
---|
86 | b(j) = 0.0d0
|
---|
87 | alpha(j) = 0.0d0
|
---|
88 | g1(j) = 0.0d0
|
---|
89 | xhist(j) = 0.0d0
|
---|
90 | ent(j) = 0.0d0
|
---|
91 | end do
|
---|
92 | end if
|
---|
93 | !
|
---|
94 | ! _________________________________ random start
|
---|
95 | do i=1,nvr
|
---|
96 | iv=idvr(i) ! provides index of non-fixed variable
|
---|
97 | dv=axvr(iv)*(grnd()-0.5)
|
---|
98 | vr=addang(pi,dv)
|
---|
99 | vlvr(iv)=vr
|
---|
100 | enddo
|
---|
101 |
|
---|
102 | eol = energy()
|
---|
103 | write (*,'(a,e12.5,/)') 'Energy of start configuration: ',eol
|
---|
104 | write(*,*)
|
---|
105 |
|
---|
106 | call outpdb(1, 'start.pdb')
|
---|
107 |
|
---|
108 |
|
---|
109 | !======================Simulation
|
---|
110 | acz = 0.0d0
|
---|
111 | ! Loop over sweeps
|
---|
112 | do nsw=1,nsweep
|
---|
113 |
|
---|
114 | call metropolis(eol,acz,muca_weight)
|
---|
115 |
|
---|
116 | muold = int(min(xmax, max(xmin, eol/ebin+sign(0.5d0,eol))))
|
---|
117 | ihist(muold) = ihist(muold) + 1
|
---|
118 | !
|
---|
119 | ! Iterate multicanonical weights every nup sweeps
|
---|
120 | if(mod(nsw,nup).eq.0) then
|
---|
121 | xhist(kmax) = xhist(kmax) + ihist(kmax)
|
---|
122 | xhist(kmax-1) = xhist(kmax-1) + ihist(kmax-1)
|
---|
123 | rewind 11
|
---|
124 | do i=kmax-2,kmin,-1
|
---|
125 | xhist(i) = xhist(i) + ihist(i)
|
---|
126 | if((ihist(i).gt.1).and.(ihist(i+1).gt.1)) then
|
---|
127 | g0 = Float(ihist(i+1))*float(ihist(i))
|
---|
128 | g0 = g0/float(ihist(i)+ihist(i+1))
|
---|
129 | g1(i) = g1(i) + g0
|
---|
130 | g0 = g0/g1(i)
|
---|
131 | bi_1 = log(float(ihist(i+1)))
|
---|
132 | bi = log(float(ihist(i)))
|
---|
133 | b(i) = b(i) + g0*(bi_1 - bi)/ebin
|
---|
134 | else
|
---|
135 | if(xhist(i).lt.2.0) then
|
---|
136 | b(i) = max(0.0d0,b(i+1))
|
---|
137 | end if
|
---|
138 | end if
|
---|
139 | end do
|
---|
140 | b(kmax) = b(kmax-2)
|
---|
141 | b(kmax-1) = b(kmax-2)
|
---|
142 | alpha(kmax) = 0
|
---|
143 | ent(kmax) = (beta+b(kmax))*float(kmax)*ebin + alpha(kmax)
|
---|
144 | write(11,'(i5,5f15.4)') kmax,xhist(kmax),g1(kmax), &
|
---|
145 | & b(kmax),alpha(kmax),ent(kmax)
|
---|
146 | ihist(kmax) = 0
|
---|
147 | do i=kmax-1,kmin,-1
|
---|
148 | ihist(i) = 0
|
---|
149 | alpha(i) = alpha(i+1) + (b(i+1) - b(i))*(i+0.5d0)*ebin
|
---|
150 | ent(i) = (beta+b(i))*float(i)*ebin + alpha(i)
|
---|
151 | write(11,'(i5,5f15.4)') i,xhist(i),g1(i),b(i),alpha(i),ent(i)
|
---|
152 | end do
|
---|
153 | end if ! parameter update
|
---|
154 | end do ! loop over sweeps
|
---|
155 |
|
---|
156 | ! Final output
|
---|
157 | acz = acz/dble(nsw*nvr)
|
---|
158 | write(*,*) 'last energy',eol
|
---|
159 | write(*,*) 'acceptance rate:',acz
|
---|
160 |
|
---|
161 | write(*,*) 'Histogram:'
|
---|
162 | do i=kmin,kmax
|
---|
163 | if(xhist(i).gt.0.0d0) then
|
---|
164 | write(*,*) i,xhist(i)
|
---|
165 | end if
|
---|
166 | end do
|
---|
167 |
|
---|
168 | ! Save multicanonical parameters to file
|
---|
169 | alpha(kmax)=0.0d0
|
---|
170 | b(kmax) = b(kmax)+beta
|
---|
171 | do i= kmax-1,kmin,-1
|
---|
172 | b(i) = b(i) + beta
|
---|
173 | alpha(i) = alpha(i+1) + (b(i+1) - b(i))*(i+0.5d0)*ebin
|
---|
174 | end do
|
---|
175 | rewind 12
|
---|
176 | do i=kmin,kmax
|
---|
177 | write(12,*) i,b(i),alpha(i)
|
---|
178 | end do
|
---|
179 |
|
---|
180 | close(11)
|
---|
181 | close(12)
|
---|
182 |
|
---|
183 | deallocate(b, alpha)
|
---|
184 | deallocate(xhist, g1, ent )
|
---|
185 | deallocate(ihist)
|
---|
186 |
|
---|
187 | return
|
---|
188 | end subroutine mulcan_par
|
---|
189 |
|
---|
190 |
|
---|
191 | subroutine mulcan_sim(nequi, nsweep, nmes, nsave, kmin, kmax, binWidth, restart)
|
---|
192 | !
|
---|
193 | ! PURPOSE: PERFORM A MULTICANONICAL SIMULATION
|
---|
194 | ! REQUIRES AS INPUT THE MULTICANONICAL PARAMETER AS CALCULATED
|
---|
195 | ! BY THE SUBROUTINE mulcan_par
|
---|
196 | !
|
---|
197 | ! CALLS: addang, contacts,energy,metropolis
|
---|
198 | !
|
---|
199 | include 'INCL.H'
|
---|
200 |
|
---|
201 | ! restart: .true. = restart of simulation
|
---|
202 | ! .false. = start of simulation with random configuration
|
---|
203 | ! kmin,kmax: Range of multisvn+ssh://nicole.nic.kfa-juelich.de/home/meinke/repositories/smmp/smmp/trunk/SMMPcanonical parameter
|
---|
204 | ! ebin: bin size for multicanonical parameter
|
---|
205 | ! nequi: Number of sweeps for equilibrisation
|
---|
206 | ! nsweep: Number of sweeps for simulation run
|
---|
207 | ! nsave: Number of sweeps after which actual configuration is saved
|
---|
208 | ! for re-starts
|
---|
209 | ! nmes: Number of sweeps between measurments
|
---|
210 | integer, intent(in) :: kmin, kmax, nsweep, nmes, nequi, nsave
|
---|
211 | real*8, intent(in) :: binWidth
|
---|
212 | logical, intent(in) :: restart
|
---|
213 |
|
---|
214 | real*8, allocatable :: xhist(:),g1(:),ent(:)
|
---|
215 | integer, allocatable :: ihist(:)
|
---|
216 |
|
---|
217 | integer :: i,j
|
---|
218 | integer :: iv,nsw,muold
|
---|
219 | real*8 :: bi,bi_1,vr,g0,dv
|
---|
220 | real*8 :: acz
|
---|
221 |
|
---|
222 | xmin = dble(kmin)
|
---|
223 | xmax = dble(kmax)
|
---|
224 | ebin = binWidth
|
---|
225 |
|
---|
226 | allocate(b(kmin:kmax), alpha(kmin:kmax))
|
---|
227 | allocate(xhist(kmin:kmax), g1(kmin:kmax), ent(kmin:kmax) )
|
---|
228 | allocate(ihist(kmin:kmax))
|
---|
229 |
|
---|
230 | if (.not.(allocated(b).and.allocated(alpha).and.allocated(xhist) &
|
---|
231 | .and.allocated(g1).and.allocated(ent).and.allocated(ihist))) then
|
---|
232 | stop "Unable to allocate memory in mulcan_sim. Exiting."
|
---|
233 | end if
|
---|
234 |
|
---|
235 |
|
---|
236 |
|
---|
237 | ! FILE with last conformation (for re-starts)
|
---|
238 | open(8,file='start.d')
|
---|
239 | ! File with contact map of reference configuration
|
---|
240 | ! FIXME: This must go. Reference structure needs to be read in main()
|
---|
241 | open(9,file='enkefa.ref')
|
---|
242 | ! File with multicanonical parameter
|
---|
243 | open(10,file='muca.d')
|
---|
244 | ! Result file: Time series of certain quantities
|
---|
245 | open(11, file='time.d')
|
---|
246 |
|
---|
247 | do j=kmin,kmax
|
---|
248 | ihist(j)= 0
|
---|
249 | xhist(j) = 0
|
---|
250 | end do
|
---|
251 |
|
---|
252 | nresi=irsml2(1)-irsml1(1) + 1 ! nresi: Number of residues
|
---|
253 | ! Read reference contact map
|
---|
254 | nci = 0
|
---|
255 | do i=1,nresi
|
---|
256 | read(9,*) (iref(i,j) , j=1,nresi)
|
---|
257 | end do
|
---|
258 | do i=1,nresi
|
---|
259 | do j=nresi,i+3,-1
|
---|
260 | if(iref(i,j).eq.1) nci = nci + 1
|
---|
261 | end do
|
---|
262 | end do
|
---|
263 | write(*,*) 'Number of contacts in reference conformation:',nci
|
---|
264 |
|
---|
265 | ! Read in fields with multicanonical parameter
|
---|
266 | Do j=kmin,kmax
|
---|
267 | read(10,*) i,b(i),alpha(i)
|
---|
268 | end do
|
---|
269 |
|
---|
270 |
|
---|
271 | if(restart) then
|
---|
272 | read(8,*) nswm, eol_old
|
---|
273 | read(8,*) (xhist(j), j=kmin,kmax)
|
---|
274 | do i=1,nvr
|
---|
275 | read(8,*) j, x
|
---|
276 | iv = idvr(j)
|
---|
277 | vlvr(iv) = x
|
---|
278 | end do
|
---|
279 | write(*,*) 'Last iteration, energy:',nswm,eol_old
|
---|
280 | else
|
---|
281 | ! _________________________________ random start
|
---|
282 | do i=1,nvr
|
---|
283 | iv=idvr(i) ! provides index of non-fixed variable
|
---|
284 | dv=axvr(iv)*(grnd()-0.5)
|
---|
285 | vr=addang(pi,dv)
|
---|
286 | vlvr(iv)=vr
|
---|
287 | enddo
|
---|
288 | end if
|
---|
289 |
|
---|
290 | eol = energy()
|
---|
291 | write (*,'(e12.5,/)') eol
|
---|
292 | call contacts(nhy,nhx,dham)
|
---|
293 | write(*,*) 'Number of contacts in start configuration:',nhy
|
---|
294 | write(*,*) 'Number of native contacts in start configuration:', &
|
---|
295 | & nhx
|
---|
296 | do i=1,nresi
|
---|
297 | write(*,'(62I1)') (ijcont(i,j), j=1,nresi)
|
---|
298 | end do
|
---|
299 | write(*,*)
|
---|
300 |
|
---|
301 |
|
---|
302 |
|
---|
303 | if(.not.restart) then
|
---|
304 | ! _________________Equilibrization by Metropolis
|
---|
305 | do nsw=1,nequi
|
---|
306 | call metropolis(eol,acz,muca_weight2)
|
---|
307 | end do
|
---|
308 | do i=kmin,kmax
|
---|
309 | ihist(i) = 0
|
---|
310 | end do
|
---|
311 | nswm = 1
|
---|
312 | end if
|
---|
313 |
|
---|
314 | !======================Simulation
|
---|
315 | acz = 0.0d0
|
---|
316 | ! Loop over sweeps
|
---|
317 | do nsw=nswm,nsweep
|
---|
318 | call metropolis(eol,acz,muca_weight2)
|
---|
319 | muold = min(kmax,max(kmin,int(eol/ebin+sign(0.5d0,eol))))
|
---|
320 | ihist(muold) = ihist(muold) + 1
|
---|
321 |
|
---|
322 | ! SAVE ACTUAL CONFORMATIONS FOR RE-STARTS:
|
---|
323 | if(mod(nsw,nsave).eq.0) then
|
---|
324 | rewind 8
|
---|
325 | write(8,*) nswm, eol
|
---|
326 | write(8,*) (xhist(j), j=kmin,kmax)
|
---|
327 | do i=1,nvr
|
---|
328 | iv = idvr(i)
|
---|
329 | write(8,*) i,vlvr(iv)
|
---|
330 | end do
|
---|
331 | end if
|
---|
332 | ! Measurements after NMES sweeps
|
---|
333 | if(mod(nsw,nmes).eq.0) then
|
---|
334 | ! Take a histogram of energy
|
---|
335 | do i=kmin,kmax
|
---|
336 | xhist(i) = xhist(i) + ihist(i)
|
---|
337 | ihist(i) = 0
|
---|
338 | end do
|
---|
339 | ! Calculate contacts in actual configuartion and compare with reference
|
---|
340 | ! configuration
|
---|
341 | ! call contacts(nhx,nhy,dham)
|
---|
342 | ! nhx : Number of contcats in actual conformation
|
---|
343 | ! nhy : Number of contacts which are identical in actual and reference
|
---|
344 | ! configuration
|
---|
345 | ! dham: Hamming distance between actual and reference configuration
|
---|
346 | !
|
---|
347 | write(11,'(i7,f12.2,2i8,f12.4)') nsw,eol,nhx,nhy,dham
|
---|
348 | end if
|
---|
349 | end do ! End of simulation
|
---|
350 |
|
---|
351 | acz = acz/dble(nsw*nvr)
|
---|
352 | write(*,*) 'last energy',eol
|
---|
353 | write(*,*) 'acceptance rate:',acz
|
---|
354 |
|
---|
355 | ! WRITE DOWN (UN-REWEIGHTED) HISTOGRAM OF MULTICANONICAL SIMULATION
|
---|
356 | do i=kmin,kmax
|
---|
357 | if(xhist(i).gt.0.0d0) then
|
---|
358 | write(*,*) i,xhist(i)
|
---|
359 | end if
|
---|
360 | end do
|
---|
361 |
|
---|
362 | close(8)
|
---|
363 | close(9)
|
---|
364 | close(10)
|
---|
365 | close(11)
|
---|
366 |
|
---|
367 | deallocate(b, alpha)
|
---|
368 | deallocate(xhist, g1, ent )
|
---|
369 | deallocate(ihist)
|
---|
370 |
|
---|
371 | return
|
---|
372 | end subroutine mulcan_sim
|
---|
373 |
|
---|
374 | ! ************************************************************
|
---|
375 | real*8 function muca_weight(x)
|
---|
376 |
|
---|
377 | real*8, intent(in) :: x
|
---|
378 |
|
---|
379 | real*8 :: beta
|
---|
380 | common /bet/beta
|
---|
381 |
|
---|
382 | integer :: muold
|
---|
383 |
|
---|
384 |
|
---|
385 | muold = int(min(xmax,max(xmin,x/ebin+sign(0.5d0,x))))
|
---|
386 | muca_weight = (beta+b(muold))*x + alpha(muold)
|
---|
387 |
|
---|
388 | return
|
---|
389 |
|
---|
390 | end function muca_weight
|
---|
391 |
|
---|
392 | ! ************************************************************
|
---|
393 | real*8 function muca_weight2(x)
|
---|
394 |
|
---|
395 | real*8, intent(in) :: x
|
---|
396 |
|
---|
397 | real*8 :: beta
|
---|
398 | common /bet/beta
|
---|
399 |
|
---|
400 | integer :: muold
|
---|
401 |
|
---|
402 | muold = int(min(xmax,max(xmin,x/ebin+sign(0.5d0,x))))
|
---|
403 | muca_weight2 = b(muold)*x + alpha(muold)
|
---|
404 |
|
---|
405 | return
|
---|
406 |
|
---|
407 | end function muca_weight2
|
---|
408 |
|
---|
409 | end module
|
---|
410 |
|
---|