| 1 | # algorithms.py
|
|---|
| 2 | #
|
|---|
| 3 | # Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
|---|
| 4 | # Jan H. Meinke, Sandipan Mohanty
|
|---|
| 5 | #
|
|---|
| 6 | import copy
|
|---|
| 7 | import random
|
|---|
| 8 | import smmp
|
|---|
| 9 | from math import *
|
|---|
| 10 | from universe import *
|
|---|
| 11 | from protein import *
|
|---|
| 12 |
|
|---|
| 13 | class CanonicalMonteCarlo(object):
|
|---|
| 14 | """An implementation of a canonical Monte Carlo run using Metropolis weights.
|
|---|
| 15 | After a Monte-Carlo sweep over all internal variables, this implementation
|
|---|
| 16 | calculates the interactions between all Protein s in the Universe
|
|---|
| 17 | """
|
|---|
| 18 |
|
|---|
| 19 | def __init__(self, myUniverse, nequi, sweeps, weight = lambda x : x * smmp.bet.beta):
|
|---|
| 20 | self.__defSweeps = sweeps
|
|---|
| 21 | self.__defEqui = nequi
|
|---|
| 22 | self.__univ = myUniverse
|
|---|
| 23 | self.__acc = 0
|
|---|
| 24 | self.__performedSweeps = 0
|
|---|
| 25 | self.__minEnergy = myUniverse.energy()
|
|---|
| 26 | self.__l = 0.1
|
|---|
| 27 | self.__minCounter = 0
|
|---|
| 28 | self.__weight = weight
|
|---|
| 29 |
|
|---|
| 30 | ## Fraction of accepted Monte Carlo moves of the total number of attempted moves.
|
|---|
| 31 | # @pre At least one move must have been performed before this function can be
|
|---|
| 32 | # called
|
|---|
| 33 | # @return fraction of accepted moves as a number between 0 and 1.
|
|---|
| 34 | def acceptanceRate(self):
|
|---|
| 35 | return float(self.__acc) / (self.__performedSweeps * smmp.mol_par.nvr)
|
|---|
| 36 | ## The lowest energy ever seen at the end of a Monte Carlo move.
|
|---|
| 37 | # @return lowest energy seen at the end of a sweep, sweep
|
|---|
| 38 | def minimumEnergy(self):
|
|---|
| 39 | return self.__minEnergy, self.__minCounter
|
|---|
| 40 | ## Number of sweeps performed so far
|
|---|
| 41 | def performedSweeps(self):
|
|---|
| 42 | return self.__performedSweeps
|
|---|
| 43 |
|
|---|
| 44 | def run(self):
|
|---|
| 45 | """A complete Monte Carlo run including initial equilibration."""
|
|---|
| 46 | self.sweeps(self.__defSweeps)
|
|---|
| 47 |
|
|---|
| 48 | def equilibrate(self):
|
|---|
| 49 | """Do nequi sweeps to eqilibrate the system."""
|
|---|
| 50 | self.sweeps(self.__defEqui)
|
|---|
| 51 |
|
|---|
| 52 | def sweeps(self, n):
|
|---|
| 53 | """Perform n sweeps."""
|
|---|
| 54 | for i in range(0, n):
|
|---|
| 55 | self.sweep()
|
|---|
| 56 |
|
|---|
| 57 | def sweep(self):
|
|---|
| 58 | """Performs a single Monte Carlo sweep."""
|
|---|
| 59 | eol = self.__univ.energy()
|
|---|
| 60 | eol, self.__acc = smmp.metropolis(eol, self.__acc, self.__weight)
|
|---|
| 61 | #eol = self.metropolis()
|
|---|
| 62 | self.__performedSweeps += 1
|
|---|
| 63 | if eol < self.__minEnergy:
|
|---|
| 64 | self.__minEnergy = eol
|
|---|
| 65 | self.__minCounter = self.__performedSweeps
|
|---|
| 66 | smmp.outpdb(0, "best.pdb")
|
|---|
| 67 |
|
|---|
| 68 | def metropolis(self):
|
|---|
| 69 | """Implementation of the Metropolis algorithm."""
|
|---|
| 70 | eol = self.__univ.energy()
|
|---|
| 71 | # Update the internal degrees of freedom first
|
|---|
| 72 | for i in range(0, smmp.mol_par.nvr):
|
|---|
| 73 | jv = smmp.var_i.idvr[i] - 1
|
|---|
| 74 | vrol = smmp.var_r.vlvr[jv]
|
|---|
| 75 | dv=smmp.var_r.axvr[jv] * (random.random() - 0.5)
|
|---|
| 76 | smmp.var_r.vlvr[jv] = (vrol + dv + pi) % (2 * pi) - pi
|
|---|
| 77 | etrial = self.__univ.energy()
|
|---|
| 78 | dW = self.__weight(etrial) - self.__weight(eol)
|
|---|
| 79 | if dW < 0:
|
|---|
| 80 | eol = etrial
|
|---|
| 81 | self.__acc += 1
|
|---|
| 82 | else:
|
|---|
| 83 | rd =random.random()
|
|---|
| 84 | dW = min(dW, 100)
|
|---|
| 85 | if exp(-dW) > rd:
|
|---|
| 86 | eol = etrial
|
|---|
| 87 | self.__acc +=1
|
|---|
| 88 | else:
|
|---|
| 89 | smmp.var_r.vlvr[jv] = vrol
|
|---|
| 90 | # Update the global degrees of freedom
|
|---|
| 91 | return eol
|
|---|
