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