| 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 | double precision :: grnd, addang, eol, energy
|
|---|
| 48 |
|
|---|
| 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()
|
|---|
| 105 | write (logString, '(a,e12.5,/)') 'Energy of start configuration: ',eol
|
|---|
| 106 | write (logString, *)
|
|---|
| 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)
|
|---|
| 160 | write (logString, *) 'last energy',eol
|
|---|
| 161 | write (logString, *) 'acceptance rate:',acz
|
|---|
| 162 |
|
|---|
| 163 | write (logString, *) 'Histogram:'
|
|---|
| 164 | do i=kmin,kmax
|
|---|
| 165 | if(xhist(i).gt.0.0d0) then
|
|---|
| 166 | write (logString, *) i,xhist(i)
|
|---|
| 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(:)
|
|---|
| 218 |
|
|---|
| 219 | double precision :: eol_old, x, grnd, addang, eol, energy, dham
|
|---|
| 220 |
|
|---|
| 221 | integer :: nresi, nswm, nhy, nhx
|
|---|
| 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()
|
|---|
| 246 | open(9,file='1vp.ref')
|
|---|
| 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
|
|---|
| 268 | write (logString, *) 'Number of contacts in reference conformation:',nci
|
|---|
| 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
|
|---|
| 284 | write (logString, *) 'Last iteration, energy:',nswm,eol_old
|
|---|
| 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()
|
|---|
| 296 | write (logString, '(e12.5,/)') eol
|
|---|
| 297 | call contacts(nhy,nhx,dham)
|
|---|
| 298 | write (logString, *) 'Number of contacts in start configuration:',nhy
|
|---|
| 299 | write (logString, *) 'Number of native contacts in start configuration:', &
|
|---|
| 300 | & nhx
|
|---|
| 301 | do i=1,nresi
|
|---|
| 302 | write (logString, '(62I1)') (ijcont(i,j), j=1,nresi)
|
|---|
| 303 | end do
|
|---|
| 304 | write (logString, *)
|
|---|
| 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)
|
|---|
| 357 | write (logString, *) 'last energy',eol
|
|---|
| 358 | write (logString, *) 'acceptance rate:',acz
|
|---|
| 359 |
|
|---|
| 360 | ! WRITE DOWN (UN-REWEIGHTED) HISTOGRAM OF MULTICANONICAL SIMULATION
|
|---|
| 361 | do i=kmin,kmax
|
|---|
| 362 | if(xhist(i).gt.0.0d0) then
|
|---|
| 363 | write (logString, *) i,xhist(i)
|
|---|
| 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 |
|
|---|
