Coverage for cvpack/reporting/cv_writer.py: 93%

1""" 

2.. class:: CVWriter 

3 :platform: Linux, MacOS, Windows 

4 :synopsis: A custom writer for reporting collective variable data 

5 

6.. moduleauthor:: Charlles Abreu <craabreu@gmail.com> 

7 

8""" 

9 

10import openmm as mm 

11 

12from ..collective_variable import CollectiveVariable 

13from .custom_writer import CustomWriter 

14 

15 

16class CVWriter(CustomWriter): 

17 """ 

18 A custom writer for reporting collective variable data. 

19 

20 Parameters 

21 ---------- 

22 variable 

23 The collective variable whose values and/or effective masses will be reported. 

24 value 

25 If `True`, report the values of the collective variable. 

26 emass 

27 If `True`, report the effective masses of the collective variable. 

28 allowReinitialization 

29 Whether to allow context reinitialization when the collective variable shares 

30 its force group with other forces (see :meth:`CollectiveVariable.getValue`). 

31 

32 Examples 

33 -------- 

34 >>> import cvpack 

35 >>> import openmm 

36 >>> from cvpack.reporting import StateDataReporter, CVWriter 

37 >>> from math import pi 

38 >>> from openmm import app, unit 

39 >>> from sys import stdout 

40 >>> from openmmtools import testsystems 

41 >>> model = testsystems.AlanineDipeptideVacuum() 

42 >>> phi = cvpack.Torsion(4, 6, 8, 14, name="phi") 

43 >>> phi.addToSystem(model.system) 

44 >>> psi = cvpack.Torsion(6, 8, 14, 16, name="psi") 

45 >>> psi.addToSystem(model.system) 

46 >>> reporter = StateDataReporter( 

47 ... stdout, 

48 ... 100, 

49 ... writers=[ 

50 ... CVWriter(phi, value=True, emass=True), 

51 ... CVWriter(psi, value=True, emass=True), 

52 ... ], 

53 ... step=True, 

54 ... ) 

55 >>> integrator = openmm.LangevinIntegrator( 

56 ... 300 * unit.kelvin, 

57 ... 1 / unit.picosecond, 

58 ... 2 * unit.femtosecond, 

59 ... ) 

60 >>> integrator.setRandomNumberSeed(1234) 

61 >>> # Integrate only the physical forces, so phi and psi do not bias the dynamics 

62 >>> integrator.setIntegrationForceGroups({0}) 

63 >>> simulation = app.Simulation(model.topology, model.system, integrator) 

64 >>> simulation.context.setPositions(model.positions) 

65 >>> simulation.context.setVelocitiesToTemperature(300 * unit.kelvin, 5678) 

66 >>> simulation.reporters.append(reporter) 

67 >>> simulation.step(1000) # doctest: +SKIP 

68 #"Step","phi (rad)",...,"emass[psi] (nm**2 Da/(rad**2))" 

69 100,-2.4309...,0.0404...,2.7368...,0.0497... 

70 200,-2.0470...,0.0386...,2.2298...,0.0456... 

71 300,-2.4721...,0.0424...,2.2528...,0.0422... 

72 400,-2.4032...,0.0406...,-3.1341...,0.0546... 

73 500,-2.8165...,0.0460...,-2.5845...,0.0458... 

74 600,-2.6256...,0.0401...,3.1382...,0.0509... 

75 700,-2.4156...,0.0463...,2.1037...,0.0406... 

76 800,-2.3871...,0.0372...,1.8683...,0.0408... 

77 900,-1.8631...,0.0322...,2.0393...,0.0374... 

78 1000,-1.0804...,0.0240...,1.8363...,0.0354... 

79 """ 

80 

81 def __init__( 

82 self, 

83 variable: CollectiveVariable, 

84 value: bool = False, 

85 emass: bool = False, 

86 allowReinitialization: bool = False, 

87 ) -> None: 

88 if not isinstance(variable, CollectiveVariable): 

89 raise TypeError("variable must be a CollectiveVariable object") 

90 if not (value or emass): 

91 raise ValueError("At least one of value or emass must be True") 

92 self._cv = variable 

93 self._value = value 

94 self._emass = emass 

95 self._allow_reinitialization = allowReinitialization 

96 

97 def getHeaders(self) -> list[str]: 

98 headers = [] 

99 if self._value: 

100 headers.append(f"{self._cv.getName()} ({self._cv.getUnit().get_symbol()})") 

101 if self._emass: 

102 headers.append( 

103 f"emass[{self._cv.getName()}] ({self._cv.getMassUnit().get_symbol()})" 

104 ) 

105 return headers 

106 

107 def getValues(self, simulation: mm.app.Simulation) -> list[float]: 

108 context = simulation.context 

109 values = [] 

110 if self._value: 

111 value = self._cv.getValue(context, self._allow_reinitialization) 

112 values.append(value / self._cv.getUnit()) 

113 if self._emass: 

114 mass = self._cv.getEffectiveMass(context, self._allow_reinitialization) 

115 values.append(mass / self._cv.getMassUnit()) 

116 return values