Coverage for openxps/couplings/harmonic_coupling.py: 100%

19 statements  

« prev     ^ index     » next       coverage.py v7.11.3, created at 2025-11-13 22:08 +0000

1""" 

2Harmonic coupling. 

3 

4.. module:: openxps.couplings.harmonic_coupling 

5 :platform: Linux, MacOS, Windows 

6 :synopsis: A harmonic coupling between a dynamical variable and a collective variable 

7 

8.. classauthor:: Charlles Abreu <craabreu@gmail.com> 

9 

10""" 

11 

12import cvpack 

13from openmm import unit as mmunit 

14 

15from ..dynamical_variable import DynamicalVariable 

16from .collective_variable_coupling import CollectiveVariableCoupling 

17 

18 

19class HarmonicCoupling(CollectiveVariableCoupling): 

20 r"""A harmonic coupling between a dynamical variable and a collective variable. 

21 

22 The coupling energy is given by: 

23 

24 .. math:: 

25 

26 U = \frac{1}{2} \kappa \left(s - q({\bf r})\right)^2 

27 

28 where :math:`s` is an extended dynamical variable, :math:`q({\bf r})` is a 

29 physical collective variable, and :math:`\kappa` is a coupling constant. 

30 

31 Parameters 

32 ---------- 

33 collective_variable 

34 The collective variable used in the coupling. 

35 dynamical_variable 

36 The dynamical variable used in the coupling. 

37 force_constant 

38 The force constant for the coupling. 

39 

40 Examples 

41 -------- 

42 >>> from copy import copy 

43 >>> import cvpack 

44 >>> import openxps as xps 

45 >>> from openmm import unit 

46 >>> dvmass = 3 * unit.dalton * (unit.nanometer / unit.radian)**2 

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

48 >>> phi_s = xps.DynamicalVariable( 

49 ... "phi_s", unit.radian, dvmass, xps.CircularBounds() 

50 ... ) 

51 >>> kappa = 1000 * unit.kilojoule_per_mole / unit.radian**2 

52 >>> xps.HarmonicCoupling(phi, phi_s, kappa) 

53 HarmonicCoupling("0.5*kappa_phi_phi_s*((phi-phi_s-6.28...*floor(...))^2)") 

54 """ 

55 

56 def __init__( 

57 self, 

58 collective_variable: cvpack.CollectiveVariable, 

59 dynamical_variable: DynamicalVariable, 

60 force_constant: mmunit.Quantity, 

61 ) -> None: 

62 self._validateArguments(collective_variable, dynamical_variable, force_constant) 

63 kappa = f"kappa_{collective_variable.getName()}_{dynamical_variable.name}" 

64 function = ( 

65 f"0.5*{kappa}*({dynamical_variable.distanceTo(collective_variable)}^2)" 

66 ) 

67 super().__init__( 

68 function=function, 

69 collective_variables=[collective_variable], 

70 dynamical_variables=[dynamical_variable], 

71 **{kappa: force_constant}, 

72 ) 

73 

74 def _validateArguments(self, cv, dv, kappa): 

75 pair = f"{cv.getName()} and {dv.name}" 

76 if not cv.getUnit().is_compatible(dv.unit): 

77 raise ValueError(f"Incompatible units for {pair}.") 

78 if (dv.isPeriodic() == cv.getPeriodicBounds() is None) or ( 

79 cv.getPeriodicBounds() != dv.bounds.asQuantity() 

80 ): 

81 raise ValueError(f"Incompatible periodicity for {pair}.") 

82 if mmunit.is_quantity(kappa) and not kappa.unit.is_compatible( 

83 mmunit.kilojoule_per_mole / dv.unit**2 

84 ): 

85 raise ValueError(f"Incompatible force constant units for {pair}.") 

86 

87 

88HarmonicCoupling.registerTag("!openxps.HarmonicCoupling")