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

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