| 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='1vp.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 |
|
|---|
