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 typing as t 

11 

12import openmm as mm 

13 

14from ..collective_variable import CollectiveVariable 

15from .custom_writer import CustomWriter 

16 

17 

18class CVWriter(CustomWriter): 

19 """ 

20 A custom writer for reporting collective variable data. 

21 

22 Parameters 

23 ---------- 

24 variable 

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

26 value 

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

28 emass 

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

30 

31 Examples 

32 -------- 

33 >>> import cvpack 

34 >>> import openmm 

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

36 >>> from math import pi 

37 >>> from openmm import app, unit 

38 >>> from sys import stdout 

39 >>> from openmmtools import testsystems 

40 >>> model = testsystems.AlanineDipeptideVacuum() 

41 >>> phi = cvpack.Torsion(6, 8, 14, 16, name="phi") 

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

43 >>> psi = cvpack.Torsion(8, 14, 16, 18, name="psi") 

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

45 >>> reporter = StateDataReporter( 

46 ... stdout, 

47 ... 100, 

48 ... writers=[ 

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

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

51 ... ], 

52 ... step=True, 

53 ... ) 

54 >>> integrator = openmm.LangevinIntegrator( 

55 ... 300 * unit.kelvin, 

56 ... 1 / unit.picosecond, 

57 ... 2 * unit.femtosecond, 

58 ... ) 

59 >>> integrator.setRandomNumberSeed(1234) 

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

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

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

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

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

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

66 100,2.7102...,0.04970...,3.1221...,0.05386... 

67 200,2.1573...,0.04481...,2.9959...,0.05664... 

68 300,2.0859...,0.04035...,-3.001...,0.04506... 

69 400,2.8061...,0.05399...,3.0792...,0.04992... 

70 500,-2.654...,0.04784...,3.1139...,0.05592... 

71 600,3.1390...,0.05137...,-3.071...,0.05063... 

72 700,2.1133...,0.04145...,3.1047...,0.04724... 

73 800,1.7348...,0.04123...,-3.004...,0.05548... 

74 900,1.6273...,0.03007...,3.1277...,0.05271... 

75 1000,1.680...,0.03749...,2.9692...,0.04450... 

76 """ 

77 

78 def __init__( # pylint: disable=dangerous-default-value 

79 self, 

80 variable: CollectiveVariable, 

81 value: bool = False, 

82 emass: bool = False, 

83 ) -> None: 

84 if not isinstance(variable, CollectiveVariable): 

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

86 if not (value or emass): 

87 raise ValueError("At least one of value or effective_mass must be True") 

88 self._cv = variable 

89 self._value = value 

90 self._emass = emass 

91 

92 def getHeaders(self) -> t.List[str]: 

93 headers = [] 

94 if self._value: 

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

96 if self._emass: 

97 headers.append( 

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

99 ) 

100 return headers 

101 

102 def getValues(self, simulation: mm.app.Simulation) -> t.List[float]: 

103 context = simulation.context 

104 values = [] 

105 if self._value: 

106 values.append(self._cv.getValue(context) / self._cv.getUnit()) 

107 if self._emass: 

108 values.append(self._cv.getEffectiveMass(context) / self._cv.getMassUnit()) 

109 return values