Package csb :: Package statistics :: Package samplers :: Package mc :: Module neqsteppropagator |
[frames] | no frames]

# Source Code for Module csb.statistics.samplers.mc.neqsteppropagator

       1  """
       2  Propagator class employing stepwise trajectories as used in the NCMC
       3  algorithm (Nilmeier et al., "Nonequilibrium candidate Monte Carlo is
       4  an efficient tool for equilibrium simulation", PNAS 2011)
       5  """
       6
       7  import csb
       8
       9  import numpy
      10
      11  from abc import ABCMeta, abstractmethod
      12  from csb.statistics.samplers.mc import TrajectoryBuilder, Trajectory, augment_state, PropagationResult
      13  from csb.statistics.samplers.mc.propagators import AbstractPropagator, MDPropagator, HMCPropagator
      14  from csb.numeric import InvertibleMatrix
      15  from csb.numeric.integrators import FastLeapFrog
      16
      17 -class NonequilibriumTrajectory(Trajectory):
      18      """
      19      Trajectory holding additional information about energy difference
      20      the Jacobian.
      21
      22      @param items: sequence of trajectory states
      23      @type items: list of L{State}s
      24
      25      @param heat: heat produced during the trajectory
      26      @type heat: float
      27
      28      @param work: work expended during the trajectory
      29      @type work: float
      30
      31      @param deltaH: energy difference between initial and final states
      32      @type deltaH: float
      33
      34      @param jacobian: product of Jacobians of perturbations applied  in the
      35                       calculation of the trajectory
      36      @type jacobian: float
      37      """
      38
      39 -    def __init__(self, items, heat=0.0, work=0.0, deltaH=0.0, jacobian=1.0, stats=None):
      40
      41          super(NonequilibriumTrajectory, self).__init__(items, heat=heat, work=work)
      42
      43          self._deltaH = None
      44          self.deltaH = deltaH
      45          self._jacobian = None
      46          self.jacobian = jacobian
      47          self._stats = None
      48          self.stats = stats
      49
      50      @property
      51 -    def jacobian(self):
      52          return self._jacobian
      53      @jacobian.setter
      54 -    def jacobian(self, value):
      55          self._jacobian = value
      56
      57      @property
      58 -    def deltaH(self):
      59          return self._deltaH
      60      @deltaH.setter
      61 -    def deltaH(self, value):
      62          self._deltaH = value
      63
      64      @property
      65 -    def stats(self):
      66          return self._stats
      67      @stats.setter
      68 -    def stats(self, value):
      69          self._stats = value
      70
      71
      72 -class AbstractSystemInfo(object):
      73      """
      74      Subclasses hold all information describing a current system setup
      75      (Hamiltonian, boundaries, ...)
      76      """
      77
      78      pass
      79
      80 -class PerturbationResult(Trajectory):
      81      """
      82      Instances hold the result of a perturbation.
      83
      84      @param items: list of states defining a phase-space trajectory
      85      @type items: list of L{AbstractState}s
      86
      87      @param work: work performed on the system during perturbation
      88      @type work: float
      89
      90      @param jacobian: jacobian of the perturbation
      91      @type jacobian: float
      92
      93      @param perturbed_sys: L{AbstractSystemInfo} instance
      94                            describing the perturbed system
      95      @type perturbed_sys: L{AbstractSystemInfo}
      96      """
      97
      98 -    def __init__(self, items, perturbed_sys, work, heat=0.0, jacobian=1.0):
      99
     100          super(PerturbationResult, self).__init__(items, heat, work)
     101
     102          self._jacobian = None
     103          self.jacobian = jacobian
     104          self._perturbed_sys = None
     105          self.perturbed_sys = perturbed_sys
     106
     107      @property
     108 -    def jacobian(self):
     109          return self._jacobian
     110      @jacobian.setter
     111 -    def jacobian(self, value):
     112          self._jacobian = value
     113
     114      @property
     115 -    def perturbed_sys(self):
     116          return self._perturbed_sys
     117      @perturbed_sys.setter
     118 -    def perturbed_sys(self, value):
     119          self._perturbed_sys = value
     120
     121
     122 -class Protocol(object):
     123      """
     124      Describes a stepwise protocol as in Nilmeier et al. (2011).
     125
     126      @param steps: the steps making up the protocol
     127      @type steps: list of L{Step}s
     128      """
     129
     130 -    def __init__(self, steps):
     131
     132          self._steps = None
     133          self.steps = steps
     134
     135      @property
     136 -    def steps(self):
     137          """
     138          The steps making up the protocol
     139          """
     140          return self._steps
     141      @steps.setter
     142 -    def steps(self, value):
     143          self._steps = value
     144
     145 -class Step(object):
     146      """
     147      Defines a step in an NCMC-like stepwise protocol.
     148
     149      @param perturbation: The perturbation of the system
     150      @type perturbation: L{AbstractPerturbation}
     151
     152      @param propagation: The propagation of the perturbed system
     153      @type propagation: L{AbstractPropagation}
     154      """
     155
     156 -    def __init__(self, perturbation, propagation):
     157
     158          self._perturbation = None
     159          self.perturbation = perturbation
     160          self._propagation = None
     161          self.propagation = propagation
     162          self._perform = None
     163          self.perform = self._perform_pert_prop
     164
     165 -    def _perform_pert_prop(self, state, extra_info=None):
     166          '''
     167          First, perform the perturbation, and then the propagation.
     168          Override this in a subclass if you want to pass on extra
     169          information to the next step in the protocol or if you want
     170          to gather some statistics on what happens in the intermediate steps.
     171
     172          @param state: state to be evolved
     173          @type state: L{State}
     174          @param extra_info: possible extra information resulting
     175                             from previous steps
     176          @type extra_info: any type
     177
     178          @rtype: L{list} containing a short trajectory consisting of
     179                  the initial and the evolved state, possible extra information
     180                  which will be passed on to the next step in the protocol and
     181                  possible subclasses of L{AbstractStepStatistics} containing
     182                  information on what happend in the step.
     183          '''
     184
     185          perturbation_result = self.perturbation(state)
     186          propagation_result = self.propagation(perturbation_result.final)
     187          result_state = propagation_result.final
     188
     189          shorttraj = NonequilibriumTrajectory([state, result_state],
     190                                               heat=propagation_result.heat,
     191                                               work=perturbation_result.work,
     192                                               jacobian=perturbation_result.jacobian)
     193
     194          return shorttraj, None, None
     195
     196 -    def _perform_prop_pert(self, state, extra_info=None):
     197          '''
     198          First, perform the propagation, and then the perturbation.
     199          Override this in a subclass if you want to pass on extra
     200          information to the next step in the protocol or if you want
     201          to gather some statistics on what happens in the intermediate steps.
     202
     203          @param state: state to be evolved
     204          @type state: L{State}
     205          @param extra_info: possible extra information resulting
     206                             from previous steps
     207          @type extra_info: any type
     208
     209          @rtype: L{list} containing a short trajectory consisting of
     210                  the initial and the evolved state, possible extra information
     211                  which will be passed on to the next step in the protocol and
     212                  possible subclasses of L{AbstractStepStatistics} containing
     213                  information on what happend in the step.
     214          '''
     215
     216          propagation_result = self.propagation(state)
     217          perturbation_result = self.perturbation(propagation_result.final)
     218          result_state = perturbation_result.final
     219
     220          shorttraj = NonequilibriumTrajectory([state, result_state],
     221                                               heat=propagation_result.heat,
     222                                               work=perturbation_result.work,
     223                                               jacobian=perturbation_result.jacobian)
     224
     225          return shorttraj, None, None
     226
     227 -    def set_perturbation_first(self):
     228          """
     229          Perform first perturbation, then propagation
     230          """
     231
     232          self.perform = self._perform_pert_prop
     233
     234 -    def set_propagation_first(self):
     235          """
     236          Perform first propagation, then perturbation
     237          """
     238
     239          self.perform = self._perform_prop_pert
     240
     241      @property
     242 -    def perturbation(self):
     243          return self._perturbation
     244      @perturbation.setter
     245 -    def perturbation(self, value):
     246          self._perturbation = value
     247
     248      @property
     249 -    def propagation(self):
     250          return self._propagation
     251      @propagation.setter
     252 -    def propagation(self, value):
     253          self._propagation = value
     254
     255
     256 -class ReducedHamiltonian(object):
     257      """
     258      Describes a reduced Hamiltonian (Hamiltonian, its position gradient
     259      and the system temperature)
     260
     261      @param log_prob: log probability of the PDF under consideration, that is,
     262                       the negative potential energy of the system
     263      @type log_prob: callable
     264
     265      @param gradient: gradient of the negative log probability of the PDF under
     266                       consideration, that is, the gradient of the potential energy;
     267                       function of position array and time
     268      @type gradient: callable
     269
     270      @param temperature: system temperature
     271      @type temperature: float
     272
     273      @param mass_matrix: system mass matrix
     274      @type mass_matrix: L{InvertibleMatrix}
     275      """
     276
     277 -    def __init__(self, log_prob, gradient=None, temperature=1.0, mass_matrix=None):
     278          self._log_prob = None
     279          self.log_prob = log_prob
     280          self._gradient = None
     281          self.gradient = gradient
     282          self._temperature = None
     283          self.temperature = temperature
     284          self._mass_matrix = None
     285          self.mass_matrix = mass_matrix
     286
     287 -    def E(self, x):
     288          """
     289          Potential energy of the system, aka negative log probability
     290
     291          @param x: position vector
     292          @type x: 1D numpy array
     293          """
     294
     295          return -self.log_prob(x)
     296
     297 -    def kinetic_energy(self, p):
     298          """
     299          Kinetic energy of the system
     300
     301          @param p: system momentum vector
     302          @type p: 1D numpy array
     303          """
     304
     305          if p is not None:
     306              if self.mass_matrix is None:
     307                  return 0.5 * sum(p ** 2)
     308              else:
     309                  if self.mass_matrix.is_unity_multiple:
     310                      return 0.5 * sum(p ** 2) / self.mass_matrix[0][0]
     311                  else:
     312                      return 0.5 * numpy.dot(p, numpy.dot(self.mass_matrix.inverse, p))
     313          else:
     314              return 0.0
     315
     316 -    def rlog_prob(self, x):
     317          """
     318          Reduced log probability
     319
     320          @param x: position vector
     321          @type x: 1D numpy array
     322          """
     323
     324          return self.log_prob(x) / self.temperature
     325
     326 -    def rkinetic_energy(self, p):
     327          """
     328          Reduced kinetic energy
     329
     330          @param p: system momentum vector
     331          @type p: 1D numpy array
     332          """
     333
     334          return self.kinetic_energy(p) / self.temperature
     335
     336 -    def __call__(self, x):
     337          """
     338          Evaluates the reduced Hamiltionian at the state x
     339
     340          @param x: system state
     341          @type x: L{State}
     342          """
     343
     344          return -self.rlog_prob(x.position) + self.rkinetic_energy(x.momentum)
     345
     346      @property
     347 -    def log_prob(self):
     348          return self._log_prob
     349      @log_prob.setter
     350 -    def log_prob(self, value):
     351          self._log_prob = value
     352
     353      @property
     354 -    def gradient(self):
     355          return self._gradient
     356      @gradient.setter
     357 -    def gradient(self, value):
     358          self._gradient = value
     359
     360      @property
     361 -    def temperature(self):
     362          return self._temperature
     363      @temperature.setter
     364 -    def temperature(self, value):
     365          self._temperature = value
     366
     367      @property
     368 -    def mass_matrix(self):
     369          return self._mass_matrix
     370      @mass_matrix.setter
     371 -    def mass_matrix(self, value):
     372          self._mass_matrix = value
     373
     374
     375 -class AbstractPerturbation(object):
     376      """
     377      Describes an abstract system perturbation
     378
     379      @param sys_before: information about the system before the perturbation
     380      @type sys_before: L{AbstractSystemInfo}
     381
     382      @param sys_after: information about the system after the perturbation
     383      @type sys_after: L{AbstractSystemInfo}
     384
     385      @param param: parameters neccessary for system perturbation
     386      @type param: L{AbstractPerturbationParam}
     387
     388      @param evaluate_work: Allows to switch off the work evaluation,
     389                            which might not always be needed, in order to
     390                            save computation time.
     391      @type evaluate_work: boolean
     392      """
     393
     394      __metaclass__ = ABCMeta
     395
     396 -    def __init__(self, sys_before, sys_after, param=None, evaluate_work=True):
     397          self._sys_before = None
     398          self.sys_before = sys_before
     399          self._sys_after = None
     400          self.sys_after = sys_after
     401          self.param = param
     402          self._evaluate_work = None
     403          self.evaluate_work = evaluate_work
     404
     405      @abstractmethod
     406 -    def _run_perturbator(self, state):
     407          """
     408          Calculates the trajectory of the system while it is being perturbed.
     409
     410          @param state: The initial system state
     411          @type state: L{State}
     412
     413          @return: The trajectory of the system while it is being perturbed
     414          @rtype: L{Trajectory}
     415          """
     416
     417          pass
     418
     419      @abstractmethod
     420 -    def _calculate_work(self, traj):
     421          """
     422          Calculates the work expended during perturbation of the system.
     423
     424          @param traj: The trajectory of the system while being perturbed
     425          @type traj: L{Trajectory}
     426
     427          @return: The work expended during perturbation
     428          @rtype: float
     429          """
     430
     431          pass
     432
     433      @abstractmethod
     434 -    def _calculate_jacobian(self, traj):
     435          """
     436          Calculates the Jacobian determinant which reflects phase
     437          space compression during perturbation.
     438
     439          @param traj: The trajectory of the system while being perturbed
     440          @type traj: L{Trajectory}
     441
     442          @return: The Jacobian determinant
     443          @rtype: float
     444          """
     445
     446          pass
     447
     448 -    def _evaluate(self, state):
     449          """
     450          Performs the perturbation of the system and / or the state
     451
     452          @param state: system state
     453          @type state: L{State}
     454          """
     455
     456          traj = self._run_perturbator(state)
     457          work = self._calculate_work(traj)
     458          jacobian = self._calculate_jacobian(traj)
     459
     460          return PerturbationResult([traj.initial, traj.final], self.sys_after,
     461                                    work, jacobian=jacobian)
     462
     463 -    def __call__(self, state):
     464          """
     465          Performs the perturbation of the system and / or the state
     466
     467          @param state: system state
     468          @type state: L{State}
     469          """
     470
     471          return self._evaluate(state)
     472
     473      @property
     474 -    def sys_before(self):
     475          return self._sys_before
     476      @sys_before.setter
     477 -    def sys_before(self, value):
     478          self._sys_before = value
     479
     480      @property
     481 -    def sys_after(self):
     482          return self._sys_after
     483      @sys_after.setter
     484 -    def sys_after(self, value):
     485          self._sys_after = value
     486
     487      @property
     488 -    def param(self):
     489          return self._param
     490      @param.setter
     491 -    def param(self, value):
     492          self._param = value
     493
     494      @property
     495 -    def evaluate_work(self):
     496          return self._evaluate_work
     497      @evaluate_work.setter
     498 -    def evaluate_work(self, value):
     499          self._evaluate_work = value
     500
     501
     502 -class AbstractPropagation(object):
     503      """
     504      Describes an abstract system propagation
     505
     506      @param sys: information about the current system setup
     507      @type sys: L{AbstractSystemInfo}
     508
     509      @param param: parameters neccessary for propagating the system
     510      @type param: L{AbstractPropagationParam}
     511
     512      @param evaluate_heat: Allows to switch off the heat evaluation,
     513                            which might not always be needed, in order to
     514                            save computation time.
     515      @type evaluate_heat: boolean
     516      """
     517
     518      __metaclass__ = ABCMeta
     519
     520 -    def __init__(self, sys, param, evaluate_heat=True):
     521
     522          self._sys = None
     523          self.sys = sys
     524          self._param = None
     525          self.param = param
     526          self._evaluate_heat = None
     527          self.evaluate_heat = evaluate_heat
     528
     529      @abstractmethod
     530 -    def _propagator_factory(self):
     531          """
     532          Factory method which returns the propagator to be used for
     533          propagating the system.
     534
     535          @return: Some propagator object
     536          @rtype: L{AbstractPropagator}
     537          """
     538
     539      @abstractmethod
     540 -    def _run_propagator(self, state):
     541          """
     542          Propagates the system using the propagator instance returned
     543          by _propagator_factory().
     544
     545          @param state: Initial state
     546          @type state: L{State}
     547
     548          @return: The result of the propagation
     549          @rtype: L{PropagationResult}
     550          """
     551
     552          pass
     553
     554      @abstractmethod
     555 -    def _calculate_heat(self, traj):
     556          """
     557          Calculates the heat resulting from system propagation.
     558
     559          @param traj: The trajectory of the system during propagation
     560          @type traj: L{Trajectory}
     561
     562          @return: The heat resulting from system propagation.
     563          @rtype: float
     564          """
     565
     566          pass
     567
     568 -    def _evaluate(self, state):
     569          """
     570          Performs the propagation of a state in the specified system
     571
     572          @param state: system state
     573          @type state: L{State}
     574          """
     575
     576          traj = self._run_propagator(state)
     577          heat = self._calculate_heat(traj)
     578
     579          return PropagationResult(traj.initial, traj.final, heat=heat)
     580
     581 -    def __call__(self, state):
     582          """
     583          Performs the propagation of a state in the specified system
     584
     585          @param state: system state
     586          @type state: L{State}
     587          """
     588
     589          return self._evaluate(state)
     590
     591      @property
     592 -    def sys(self):
     593          return self._sys
     594      @sys.setter
     595 -    def sys(self, value):
     596          self._sys = value
     597
     598      @property
     599 -    def param(self):
     600          return self._param
     601      @param.setter
     602 -    def param(self, value):
     603          self._param = value
     604
     605      @property
     606 -    def evaluate_heat(self):
     607          return self._evaluate_heat
     608      @evaluate_heat.setter
     609 -    def evaluate_heat(self, value):
     610          self._evaluate_heat = value
     611
     612
     613 -class ReducedHamiltonianPerturbation(AbstractPerturbation):
     614      """
     615      System perturbation by changing the reduced Hamiltonian
     616
     617      @param sys_before: information about the system before the perturbation
     618      @type sys_before: L{AbstractSystemInfo}
     619
     620      @param sys_after: information about the system after the perturbation
     621      @type sys_after: L{AbstractSystemInfo}
     622
     623      @param evaluate_work: Allows to switch off the work evaluation,
     624                            which might not always be needed, in order to
     625                            save computation time.
     626      @type evaluate_work: boolean
     627      """
     628
     629 -    def __init__(self, sys_before, sys_after, evaluate_work=True):
     630
     631          super(ReducedHamiltonianPerturbation, self).__init__(sys_before, sys_after,
     632                                                               evaluate_work=evaluate_work)
     633
     634 -    def _calculate_work(self, traj):
     635
     636          work = 0.0
     637          if self.evaluate_work == True:
     638              work = self.sys_after.hamiltonian(traj.final) - \
     639                     self.sys_before.hamiltonian(traj.initial)
     640
     641          return work
     642
     643 -    def _calculate_jacobian(self, traj):
     644
     645          return 1.0
     646
     647 -    def _run_perturbator(self, state):
     648
     649          return Trajectory([state, state])
     650
     651
     652 -class AbstractMCPropagation(AbstractPropagation):
     653      """
     654      System propagation by some MC algorithm.
     655
     656      @param sys: information about the current system setup
     657      @type sys: L{AbstractSystemInfo}
     658
     659      @param param: parameters neccessary for propagating the system
     660      @type param: L{AbstractPropagationParam}
     661
     662      @param evaluate_heat: Allows to switch off the heat evaluation,
     663                            which might not always be needed, in order to
     664                            save computation time.
     665      @type evaluate_heat: boolean
     666      """
     667
     668      ## Not neccessary, but otherwise epydoc complains
     669 -    def __init__(self, sys, param, evaluate_heat=True):
     670
     671          super(AbstractMCPropagation, self).__init__(sys, param, evaluate_heat=True)
     672
     673 -    def _calculate_heat(self, traj):
     674
     675          heat = 0.0
     676          if self.evaluate_heat == True:
     677              heat = self.sys.hamiltonian.E(traj.final.position) - \
     678                     self.sys.hamiltonian.E(traj.initial.position)
     679
     680          return heat
     681
     682 -    def _run_propagator(self, state):
     683
     684          gen = self._propagator_factory()
     685
     686          return gen.generate(state, self.param.iterations, False)
     687
     688
     689 -class HMCPropagation(AbstractMCPropagation):
     690      """
     691      System propagation by HMC
     692
     693      @param sys: information about the current system setup
     694      @type sys: L{AbstractSystemInfo}
     695
     696      @param param: parameters neccessary for propagating the system
     697      @type param: L{HMCPropagationParam}
     698
     699      @param evaluate_heat: Allows to switch off the heat evaluation,
     700                            which might not always be needed, in order to
     701                            save computation time.
     702      @type evaluate_heat: boolean
     703      """
     704
     705 -    def __init__(self, sys, param, evaluate_heat=True):
     706
     707          super(HMCPropagation, self).__init__(sys, param, evaluate_heat)
     708
     709          if self.param.gradient is None:
     710              self.param.gradient = self.sys.hamiltonian.gradient
     711
     712 -    def _set_mass_matrix(self, state):
     713          """
     714          Sets the mass matrix in the param object.
     715
     716          @param state: The initial state which is used to determine the dimension
     717                        of the mass matrix
     718          @type state: L{State}
     719          """
     720
     721          if self.param.mass_matrix is None:
     722              d = len(state.position)
     723              self.param.mass_matrix = InvertibleMatrix(numpy.eye(d))
     724
     725 -    def _propagator_factory(self):
     726
     727          gen = HMCPropagator(self.sys.hamiltonian, self.param.gradient,
     728                              self.param.timestep, self.param.traj_length,
     729                              temperature=self.sys.hamiltonian.temperature,
     730                              integrator=self.param.integrator,
     731                              mass_matrix=self.param.mass_matrix)
     732
     733          return gen
     734
     735 -    def _evaluate(self, state):
     736
     737          self._set_mass_matrix(state)
     738
     739          return super(HMCPropagation, self)._evaluate(state)
     740
     741      @property
     742 -    def param(self):
     743          return self._param
     744      @param.setter
     745 -    def param(self, value):
     746          self._param = value
     747
     748
     749 -class AbstractMDPropagation(AbstractPropagation):
     750      """
     751      System propagation by some MD algorithm
     752
     753      @param sys: information about the current system setup
     754      @type sys: L{AbstractSystemInfo}
     755
     756      @param param: parameters neccessary for propagating the system
     757      @type param: L{MDPropagationParam}
     758
     759      @param evaluate_heat: Allows to switch off the heat evaluation,
     760                            which might not always be needed, in order to
     761                            save computation time.
     762      @type evaluate_heat: boolean
     763      """
     764
     765      __metaclass__ = ABCMeta
     766
     767      ## Not neccessary, but otherwise epydoc complains
     768 -    def __init__(self, sys, param, evaluate_heat=True):
     769
     770          super(AbstractMDPropagation, self).__init__(sys, param, evaluate_heat=True)
     771
     772 -    def _set_mass_matrix(self, state):
     773          """
     774          Sets the mass matrix in the param object.
     775
     776          @param state: The initial state which is used to determine the dimension
     777                        of the mass matrix
     778          @type state: L{State}
     779          """
     780
     781          if self.param.mass_matrix is None:
     782              d = len(state.position)
     783              self.param.mass_matrix = InvertibleMatrix(numpy.eye(d))
     784
     785 -    def _augment_state(self, state):
     786          """
     787          Augments the initial state by a momentum if none is defined.
     788
     789          @param state: Initial state
     790          @type state: L{State}
     791          """
     792
     793          if state.momentum == None:
     794              state = augment_state(state, self.sys.hamiltonian.temperature,
     795                                    self.param.mass_matrix)
     796
     797          return state
     798
     799 -    def _run_propagator(self, state):
     800
     801          gen = self._propagator_factory()
     802          state = self._augment_state(state)
     803          traj = gen.generate(state, self.param.traj_length)
     804
     805          return traj
     806
     807
     808 -class PlainMDPropagation(AbstractMDPropagation):
     809      """
     810      System propagation by plain, microcanonical MD
     811
     812      @param sys: information about the current system setup
     813      @type sys: L{AbstractSystemInfo}
     814
     815      @param param: parameters neccessary for propagating the system
     816      @type param: L{PlainMDPropagationParam}
     817
     818      @param evaluate_heat: Allows to switch off the heat evaluation,
     819                            which might not always be needed, in order to
     820                            save computation time.
     821      @type evaluate_heat: boolean
     822      """
     823
     824      ## Not neccessary, but otherwise epydoc is complaining
     825 -    def __init__(self, sys, param, evaluate_heat=True):
     826
     827          super(PlainMDPropagation, self).__init__(sys, param, evaluate_heat=evaluate_heat)
     828
     829 -    def _propagator_factory(self):
     830
     831          return MDPropagator(self.param.gradient, self.param.timestep,
     832                              mass_matrix=self.param.mass_matrix,
     833                              integrator=self.param.integrator)
     834
     835 -    def _calculate_heat(self, traj):
     836
     837          heat = 0.0
     838          if self.evaluate_heat == True:
     839              heat = self.sys.hamiltonian(traj.final) - self.sys.hamiltonian(traj.initial)
     840
     841          return heat
     842
     843 -    def _evaluate(self, state):
     844
     845          self._set_mass_matrix(state)
     846
     847          return super(PlainMDPropagation, self)._evaluate(state)
     848
     849
     850 -class AbstractPerturbationParam(object):
     851      """
     852      Subclasses hold informations required for different kinds
     853      of system perturbation
     854      """
     855
     856      pass
     857
     858
     859 -class AbstractPropagationParam(object):
     860      """
     861      Subclasses hold informations required for different kinds
     862      of system propagation
     863      """
     864
     865      pass
     866
     867
     868 -class MDParam(object):
     869      """
     870      Holds all required information for calculating a MD trajectory.
     871
     872      @param timestep: integration timestep
     873      @type timestep: float
     874
     875      @param traj_length: MD trajectory length
     876      @type traj_length: int
     877
     878      @param gradient: gradient governing the equations of motion, function of
     879                       position array and time
     880      @type gradient: callable
     881
     882      @param temperature: System temperature for drawing momenta from the
     883                          Maxwell distribution
     884      @type temperature: float
     885
     886      @param integrator: Integration scheme to be utilized
     887      @type integrator: L{AbstractIntegrator}
     888
     889      @param mass_matrix: mass matrix for kinetic energy definition
     890      @type mass_matrix: L{InvertibleMatrix}
     891      """
     892
     893 -    def __init__(self, timestep, traj_length, gradient, temperature=1.0,
     894                   integrator=FastLeapFrog, mass_matrix=None):
     895
     896          self._timestep = None
     897          self.timestep = timestep
     898          self._traj_length = None
     899          self.traj_length = traj_length
     900          self._gradient = None
     901          self.gradient = gradient
     902          self._temperature = None
     903          self.temperature = temperature
     904          self._integrator = None
     905          self.integrator = integrator
     906          self._mass_matrix = None
     907          self.mass_matrix = mass_matrix
     908
     909      @property
     910 -    def timestep(self):
     911          return self._timestep
     912      @timestep.setter
     913 -    def timestep(self, value):
     914          self._timestep = value
     915
     916      @property
     917 -    def traj_length(self):
     918          return self._traj_length
     919      @traj_length.setter
     920 -    def traj_length(self, value):
     921          self._traj_length = value
     922
     923      @property
     924 -    def gradient(self):
     925          return self._gradient
     926      @gradient.setter
     927 -    def gradient(self, value):
     928          self._gradient = value
     929
     930      @property
     931 -    def temperature(self):
     932          return self._temperature
     933      @temperature.setter
     934 -    def temperature(self, value):
     935          self._temperature = value
     936
     937      @property
     938 -    def integrator(self):
     939          return self._integrator
     940      @integrator.setter
     941 -    def integrator(self, value):
     942          self._integrator = value
     943
     944      @property
     945 -    def mass_matrix(self):
     946          return self._mass_matrix
     947      @mass_matrix.setter
     948 -    def mass_matrix(self, value):
     949          self._mass_matrix = value
     950
     951
     952 -class HMCPropagationParam(MDParam, AbstractPropagationParam):
     953      """
     954      Holds all required information for propagating a system by HMC.
     955      The system temperature is taken from the
     956      HMCPropagation.sys.hamiltonian.temperature property.
     957
     958      @param timestep: integration timestep
     959      @type timestep: float
     960
     961      @param traj_length: MD trajectory length
     962      @type traj_length: int
     963
     964      @param gradient: gradient governing the equations of motion, function of
     965                       position array and time
     966      @type gradient: callable
     967
     968      @param iterations: number of HMC iterations to be performed
     969      @type iterations: int
     970
     971      @param integrator: Integration scheme to be utilized
     972      @type integrator: l{AbstractIntegrator}
     973
     974      @param mass_matrix: mass matrix for kinetic energy definition
     975      @type mass_matrix: L{InvertibleMatrix}
     976      """
     977
     978 -    def __init__(self, timestep, traj_length, gradient, iterations=1,
     979                   integrator=FastLeapFrog, mass_matrix=None):
     980
     981          super(HMCPropagationParam, self).__init__(timestep, traj_length, gradient,
     982                                                    integrator=integrator, mass_matrix=mass_matrix)
     983
     984          self._iterations = None
     985          self.iterations = iterations
     986
     987      @property
     988 -    def iterations(self):
     989          return self._iterations
     990      @iterations.setter
     991 -    def iterations(self, value):
     992          self._iterations = value
     993
     994
     995 -class MDPropagationParam(MDParam, AbstractPropagationParam):
     996
     997      pass
     998
     999
    1000 -class PlainMDPropagationParam(MDParam, AbstractPropagationParam):
    1001      """
    1002      Holds all required information for propagating a system by MD.
    1003      The system temperature is taken from the
    1004      MDPropagation.sys.hamiltonian.temperature property.
    1005
    1006      @param timestep: integration timestep
    1007      @type timestep: float
    1008
    1009      @param traj_length: MD trajectory length
    1010      @type traj_length: int
    1011
    1012      @param gradient: gradient governing the equations of motion, function of
    1013                       position array and time
    1014      @type gradient: callable
    1015
    1016      @param integrator: Integration scheme to be utilized
    1017      @type integrator: l{AbstractIntegrator}
    1018
    1019      @param mass_matrix: mass matrix for kinetic energy definition
    1020      @type mass_matrix: L{InvertibleMatrix}
    1021      """
    1022
    1023 -    def __init__(self, timestep, traj_length, gradient,
    1024                   integrator=FastLeapFrog, mass_matrix=None):
    1025
    1026          super(PlainMDPropagationParam, self).__init__(timestep, traj_length, gradient,
    1027                                                   integrator=integrator, mass_matrix=mass_matrix)
    1028
    1029
    1030 -class HamiltonianSysInfo(AbstractSystemInfo):
    1031      """
    1032      Holds information describing a system by its Hamiltonian only.
    1033
    1034      @param hamiltonian: the Hamiltonian of the system to be described
    1035      @type hamiltonian: L{ReducedHamiltonian}
    1036      """
    1037
    1038 -    def __init__(self, hamiltonian):
    1039
    1040          self._hamiltonian = None
    1041          self.hamiltonian = hamiltonian
    1042
    1043      @property
    1044 -    def hamiltonian(self):
    1045          return self._hamiltonian
    1046      @hamiltonian.setter
    1047 -    def hamiltonian(self, value):
    1048          self._hamiltonian = value
    1049
    1050
    1051 -class AbstractStepStatistics(object):
    1052      '''
    1053      Abstract class defining a minimal interface for objects allowing to store statistics
    1054      of what happens in L{Step} instances.
    1055      '''
    1056
    1057      @abstractmethod
    1058 -    def update(self, step_index, shorttraj, stats_data):
    1059          pass
    1060
    1061
    1062 -class DummyStepStatistics(AbstractStepStatistics):
    1063
    1064 -    def update(self, step_index, shorttraj, stats_data):
    1065          pass
    1066
    1067
    1068 -class AbstractHeatWorkJacobianLogger(object):
    1069      '''
    1070      Abstract class defining the interface for objects keeping track of and accumulating
    1071      heat, work and the Jacobian during a nonequilibrium trajectory.
    1072      '''
    1073
    1074 -    def __init__(self):
    1075
    1076          self._heat = 0.0
    1077          self._work = 0.0
    1078          self._jacobian = 1.0
    1079
    1080      @abstractmethod
    1081 -    def accumulate(self, heat=0.0, work=0.0, jacobian=1.0):
    1082          '''
    1083          Adds heat and work contribution to the so far accumulated values and
    1084          "multiply-accumulates" the Jacobian to the so far "multiply-accumulated" values.
    1085          '''
    1086          pass
    1087
    1088      @property
    1089 -    def heat(self):
    1090          return self._heat
    1091
    1092      @property
    1093 -    def work(self):
    1094          return self._work
    1095
    1096      @property
    1097 -    def jacobian(self):
    1098          return self._jacobian
    1099
    1100
    1101 -class TrivialHeatWorkJacobianLogger(AbstractHeatWorkJacobianLogger):
    1102
    1103 -    def accumulate(self, heat=0.0, work=0.0, jacobian=1.0):
    1104
    1105          self._heat += heat
    1106          self._work += work
    1107          self._jacobian *= jacobian
    1108
    1109
    1110 -class NonequilibriumStepPropagator(AbstractPropagator):
    1111      """
    1112      Propagator class which propagates a system using NCMC-like
    1113      stepwise trajectories
    1114
    1115      @param protocol: stepwise protocol to be followed
    1116      @type protocol: L{Protocol}
    1117      """
    1118
    1119 -    def __init__(self, protocol):
    1120
    1121          self._protocol = None
    1122          self.protocol = protocol
    1123
    1124 -    def _calculate_deltaH(self, traj):
    1125          """
    1126          Calculate the difference of the Hamiltonian between the initial and
    1127          the final state of a NCMC trajectory.
    1128
    1129          @param traj: The NCMC trajectory between whose initial and final states
    1130                       the Hamiltonian difference should be calculated
    1131          @type traj: L{NonequilibriumTrajectory}
    1132          """
    1133
    1134          return self.protocol.steps[-1].perturbation.sys_after.hamiltonian(traj.final) - \
    1135                 self.protocol.steps[0].perturbation.sys_before.hamiltonian(traj.initial)
    1136
    1137 -    def _create_heat_work_jacobian_logger(self):
    1138          '''
    1139          Factory method for the L{AbstractHeatWorkJacobianLogger} subclass instance
    1140          which keeps track of work, heat and Jacobian contributions during the nonequilibrium
    1141          trajectory.
    1142
    1143          @rtype: instance of an L{AbstractHeatWorkJacobianLogger}-derived class
    1144          '''
    1145
    1146          return TrivialHeatWorkJacobianLogger()
    1147
    1148 -    def _create_step_stats(self):
    1149          '''
    1150          Factory method for the L{AbstractStepStatistics} subclass instance
    1151          which can be used to collect statistics of what happens during steps.
    1152
    1153          @rtype: instance of an L{AbstractStepStatistics}-derived class
    1154          '''
    1155
    1156          return DummyStepStatistics()
    1157
    1158 -    def _set_initial_extra_info(self, init_state):
    1159          '''
    1160          Provides additional information for the first step in the protocol.
    1161
    1162          @rtype: any type
    1163          '''
    1164
    1165          return None
    1166
    1167 -    def _perform_step_iteration(self, estate, hwj_logger, step_stats, builder, extra_info):
    1168          '''
    1169          Performs the iteration over all steps in the nonequilibrium protocol.
    1170
    1171          @param estate: the state which will be evolved
    1172          @type estate: L{State}
    1173          @param hwj_logger: an instance of an L{AbstractHeatWorkJacobianLogger}-derived class,
    1174                             which keeps track of work, heat and Jacobian contributions
    1175          @type hwj_logger: subclass of L{AbstractHeatWorkJacobianLogger}
    1176          @param step_stats: an instance of an L{AbstractStepStatistics}-derived class,
    1177                             which may collect some statistics of what happens during steps
    1178          @type step_stats: subclass of L{AbstractStepStatistics}
    1179          @param builder: L{TrajectoryBuilder} instance in charge of building the L{Trajectory} object
    1180          @type builder: L{TrajectoryBuilder}
    1181          @param extra_info: Dictionary containing eventual additional information, which needs to
    1182                             be passed on from one step to the following
    1183          @type extra_info: any type
    1184
    1185          @rtype: L{Trajectory}
    1186          '''
    1187
    1188          for i in range(len(self.protocol.steps)):
    1189
    1190              shorttraj, extra_info, stats_data = self.protocol.steps[i].perform(estate, extra_info)
    1191
    1192              step_stats.update(i, shorttraj, stats_data)
    1193              hwj_logger.accumulate(shorttraj.heat, shorttraj.work, shorttraj.jacobian)
    1194
    1195              estate = shorttraj.final
    1196
    1197              if i == 0:
    1198                  builder.add_initial_state(shorttraj.initial)
    1199              elif i != len(self.protocol.steps) - 1:
    1200                  builder.add_intermediate_state(estate)
    1201              else:
    1202                  builder.add_final_state(estate)
    1203
    1204          return builder.product
    1205
    1206
    1207 -    def generate(self, init_state, return_trajectory=False):
    1208
    1209          estate = init_state.clone()
    1210
    1211          hwj_logger = self._create_heat_work_jacobian_logger()
    1212          step_stats = self._create_step_stats()
    1213          builder = TrajectoryBuilder.create(full=return_trajectory)
    1214          extra_info = self._set_initial_extra_info(estate)
    1215
    1216          traj = self._perform_step_iteration(estate, hwj_logger, step_stats,
    1217                                              builder, extra_info)
    1218
    1219          reduced_deltaH = self._calculate_deltaH(traj)
    1220
    1221          if init_state.momentum is None:
    1222              for s in traj:
    1223                  s.momentum = None
    1224
    1225          result = NonequilibriumTrajectory([x for x in traj],
    1226                                            heat=hwj_logger.heat,
    1227                                            work=hwj_logger.work,
    1228                                            deltaH=reduced_deltaH,
    1229                                            jacobian=hwj_logger.jacobian,
    1230                                            stats=step_stats)
    1231
    1232          return result
    1233
    1234      @property
    1235 -    def protocol(self):
    1236          return self._protocol
    1237      @protocol.setter
    1238 -    def protocol(self, value):
    1239          self._protocol = value
    1240