From 8c0399eafe13dcb580124dbab26e42c95fa9b4ea Mon Sep 17 00:00:00 2001 From: jrudz Date: Wed, 25 Sep 2024 15:55:54 +0200 Subject: [PATCH] more angle progress --- .../schema_packages/force_field.py | 153 +++++++++++++----- tests/test_force_field.py | 146 ++++++++++++++++- 2 files changed, 260 insertions(+), 39 deletions(-) diff --git a/src/nomad_simulations/schema_packages/force_field.py b/src/nomad_simulations/schema_packages/force_field.py index b09b7050..d638cbcf 100644 --- a/src/nomad_simulations/schema_packages/force_field.py +++ b/src/nomad_simulations/schema_packages/force_field.py @@ -268,7 +268,6 @@ def normalize(self, archive, logger) -> None: ) if self.forces is not None and self.energies is None: - print('in gen energies') try: # generated energies from forces numerically using spline self.energies = self.compute_energies( @@ -308,8 +307,38 @@ def normalize(self, archive, logger) -> None: ) -class PolynomialPotential(Potential): +class PolynomialForceConstant(ParameterEntry): + """ + Section defining a force constant for a potential of polynomial form. """ + + name = Quantity( + type=str, + shape=[], + description=""" + Name of the force constant. + """, + ) + + exponent = Quantity( + type=np.int32, + shape=[], + description=""" + Exponent for this term in the polynomial. + """, + ) + + value = Quantity( + type=np.float64, + shape=[], + description=""" + Value of the force constant. + """, + ) + + +class PolynomialPotential(Potential): + r""" Abstract class for potentials with polynomial form: $V(x) = [\left k_1 (x - x_0) + k_2 (x - x_0)^2 + x_3 (x - x_0)^3 + \dots + C$, where $\{x_1, x_2, x_3 \dots}$ are the `force_constants` for each term in the polynomial @@ -318,7 +347,7 @@ class PolynomialPotential(Potential): """ force_constants = SubSection( - sub_section=ParameterEntry.m_def, + sub_section=PolynomialForceConstant.m_def, repeats=True, description=""" List of force constants value and corresponding unit for polynomial potentials. @@ -329,13 +358,17 @@ def normalize(self, archive, logger) -> None: super().normalize(archive, logger) self.name = 'PolynomialPotential' + if not self.functional_form: + self.functional_form = 'polynomial' + elif self.functional_form != 'polynomial': + logger.warning('Incorrect functional form set for PolynomialPotential.') class BondPotential(Potential): """ Section containing information about bond potentials. - Suggested types are: harmonic, cubic, Morse, fene, tabulated + Suggested types are: harmonic, cubic, polynomial, Morse, fene, tabulated """ equilibrium_value = Quantity( @@ -419,47 +452,19 @@ def normalize(self, archive, logger) -> None: logger.warning('Incorrect functional form set for CubicBond.') -# TODO Add Fourth Power potential from gromacs, might want to make it a more general quartic polynomial, even though it's not as general -# class QuarticBond(BondPotential): -# r""" -# Section containing information about a Cubic bond potential: -# $V(r) = \frac{1}{2} k_r (r - r_0)^2 + \frac{1}{3} k_c (r - r_0)^3 + C$, -# where $k_r$ is the (harmonic) `force_constant`, $k_c$ is the `force_constant_cubic`, -# and $r_0$ is the `equilibrium_value` of $r$. -# C is an arbitrary constant (not stored). -# """ - -# force_constant_cubic = Quantity( -# type=np.float64, -# shape=[], -# unit='J / m**3', -# description=""" -# Specifies the cubic force constant of the bond potential. -# """, -# ) - -# def normalize(self, archive, logger) -> None: -# super().normalize(archive, logger) - -# self.name = 'CubicBond' -# if not self.functional_form: -# self.functional_form = 'cubic' -# elif self.functional_form != 'cubic': -# logger.warning('Incorrect functional form set for CubicBond.') - - class PolynomialBond(PolynomialPotential, BondPotential): """ Section containing information about a polynomial bond potential: """ def __init__(self): + super().__init__() docstring = PolynomialPotential.__doc__ pattern = r'\$V\(x\)(.*?)(\(not stored\)\.)' match = re.search(pattern, docstring, re.DOTALL) extracted_text = '' if match: - extracted_text = match.group(1).strip() + extracted_text = match.group().strip() # .group(1).strip() self.__doc__ = rf"""{self.__doc__} {extracted_text}. Here the dependent variable of the potential, $x$, corresponds to the bond distance.""" @@ -572,7 +577,7 @@ class AnglePotential(Potential): """ Section containing information about angle potentials. - Suggested types are: harmonic, cosine, fourier_series, urey_bradley, tabulated + Suggested types are: harmonic, cosine, restricted_cosinse, fourier_series, urey_bradley, polynomial, tabulated """ equilibrium_value = Quantity( @@ -745,12 +750,13 @@ class PolynomialAngle(PolynomialPotential, AnglePotential): """ def __init__(self): + super().__init__() docstring = PolynomialPotential.__doc__ pattern = r'\$V\(x\)(.*?)(\(not stored\)\.)' match = re.search(pattern, docstring, re.DOTALL) extracted_text = '' if match: - extracted_text = match.group(1).strip() + extracted_text = match.group().strip() # .group(1).strip() self.__doc__ = rf"""{self.__doc__} {extracted_text}. Here the dependent variable of the potential, $x$, corresponds to the angle between three particles.""" @@ -790,6 +796,81 @@ def normalize(self, archive, logger) -> None: self.name = 'TabulatedAngle' +class DihedralPotential(Potential): + """ + Section containing information about dihedral potentials. + + Suggested types are: fourier_series, tabulated + + # ? Something about angle convention? + """ + + equilibrium_value = Quantity( + type=np.float64, + unit='degree', + shape=[], + description=""" + Specifies the equilibrium dihedral angle. + """, + ) + + force_constant = Quantity( + type=np.float64, + shape=[], + unit='J / degree**2', + description=""" + Specifies the force constant of the dihedral angle potential. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + self.name = 'DihedralPotential' + if not self.type: + self.type = 'dihedral' + elif self.type != 'dihedral': + logger.warning('Incorrect type set for DihedralPotential.') + + if self.n_particles: + if self.n_particles != 4: + logger.warning( + 'Incorrect number of particles set for DihedralPotential.' + ) + else: + self.n_particles = 4 + + +class TabulatedDihedral(DihedralPotential, TabulatedPotential): + """ + Section containing information about a tabulated bond potential. The value of the potential and/or force + is stored for a set of corresponding bin distances. + """ + + bins = Quantity( + type=np.float64, + unit='degree', + shape=[], + description=""" + List of bin dihedral angles. + """, + ) + + forces = Quantity( + type=np.float64, + unit='J/degree', + shape=[], + description=""" + List of force values associated with each bin. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + self.name = 'TabulatedDihedral' + + class ForceField(ModelMethod): """ Section containing the parameters of a (classical, particle-based) force field model. diff --git a/tests/test_force_field.py b/tests/test_force_field.py index f0021d2b..951132aa 100644 --- a/tests/test_force_field.py +++ b/tests/test_force_field.py @@ -14,6 +14,8 @@ BondPotential, HarmonicBond, CubicBond, + PolynomialForceConstant, + PolynomialBond, MorseBond, FeneBond, TabulatedBond, @@ -23,6 +25,7 @@ RestrictedCosineAngle, FourierSeriesAngle, UreyBradleyAngle, + PolynomialAngle, TabulatedAngle, ) from nomad_simulations.schema_packages.numerical_settings import ForceCalculations @@ -217,7 +220,56 @@ def populate_parameters(sec_potential, parameters): results_cubic_bond, ) -# TODO add polynomial bond +# polynomial +results_polynomial_bond = { + 'n_interactions': 3, + 'n_particles': 2, + 'particle_indices': [[0, 1], [2, 3], [4, 5]], + 'particle_labels': [['O', 'H'], ['O', 'H'], ['O', 'H']], + 'equilibrium_value': 9.6e-11, + 'force_constants': [ + { + 'name': 'k_2', + 'value': 8.302695202135821e-21, + 'unit': 'kilojoule / nanometer ** 2', + 'exponent': 2, + }, + { + 'name': 'k_4', + 'value': 3.4871319848970445e-21, + 'unit': 'kilojoule / nanometer ** 4', + 'exponent': 4, + }, + ], + 'name': 'PolynomialBond', + 'type': 'bond', + 'functional_form': 'polynomial', +} +data_polynomial_bond = ( + PolynomialBond, + n_interactions, + n_particles, + particle_labels, + particle_indices, + { + 'equilibrium_value': 0.96 * ureg.angstrom, + 'force_constants': [ + PolynomialForceConstant( + name='k_2', + exponent=2, + value=5000.0 / MOL, + unit=str(ureg.kJ / ureg.nanometer**2), + ), + PolynomialForceConstant( + name='k_4', + exponent=4, + value=2100.0 / MOL, + unit=str(ureg.kJ / ureg.nanometer**4), + ), + ], + }, + results_polynomial_bond, +) # morse results_morse_bond = { @@ -500,8 +552,30 @@ def populate_parameters(sec_potential, parameters): results_cosine_angle, ) -# TODO add RestrictedCosineAngle -# TODO add PolynomialAngle +# restricted cosine +results_restrictedcosine_angle = { + 'n_interactions': 2, + 'n_particles': 3, + 'particle_indices': [[0, 1, 2], [3, 4, 5]], + 'particle_labels': [['O', 'H', 'H'], ['O', 'H', 'H']], + 'equilibrium_value': 104.45020605234907, + 'force_constant': 8.70026082874034e-23, + 'name': 'RestrictedCosineAngle', + 'type': 'angle', + 'functional_form': 'restricted_cosine', +} +data_restrictedcosine_angle = ( + RestrictedCosineAngle, + n_interactions, + n_particles, + particle_labels, + particle_indices, + { + 'equilibrium_value': 1.823 * ureg.radian, + 'force_constant': 172 * ureg.kJ / MOL / ureg.radian**2, + }, + results_restrictedcosine_angle, +) # fourier_series results_fourier_angle = { @@ -557,6 +631,69 @@ def populate_parameters(sec_potential, parameters): results_ureybradley_angle, ) +# polynomial +results_polynomial_angle = { + 'n_interactions': 2, + 'n_particles': 3, + 'particle_indices': [[0, 1, 2], [3, 4, 5]], + 'particle_labels': [['O', 'H', 'H'], ['O', 'H', 'H']], + 'equilibrium_value': 104.45020605234907, + 'force_constants': [ + { + 'name': 'k_2', + 'value': 1.245404280320373e-22, + 'unit': 'kilojoule / radian ** 2', + 'exponent': 2, + }, + { + 'name': 'k_3', + 'value': -3.3210780808543285e-24, + 'unit': 'kilojoule / radian ** 4', + 'exponent': 3, + }, + { + 'name': 'k_4', + 'value': 7.472425681922239e-24, + 'unit': 'kilojoule / radian ** 4', + 'exponent': 4, + }, + ], + 'name': 'PolynomialAngle', + 'type': 'angle', + 'functional_form': 'polynomial', +} +data_polynomial_angle = ( + PolynomialAngle, + n_interactions, + n_particles, + particle_labels, + particle_indices, + { + 'equilibrium_value': 1.823 * ureg.radian, + 'force_constants': [ + PolynomialForceConstant( + name='k_2', + exponent=2, + value=75 / MOL, + unit=str(ureg.kJ / ureg.radian**2), + ), + PolynomialForceConstant( + name='k_3', + exponent=3, + value=-2.0 / MOL, + unit=str(ureg.kJ / ureg.radian**4), + ), + PolynomialForceConstant( + name='k_4', + exponent=4, + value=4.5 / MOL, + unit=str(ureg.kJ / ureg.radian**4), + ), + ], + }, + results_polynomial_angle, +) + # tabulated results_tabulated_angle = { 'n_interactions': 2, @@ -698,14 +835,17 @@ def populate_parameters(sec_potential, parameters): [ data_harmonic_bond, data_cubic_bond, + data_polynomial_bond, data_morse_bond, data_fene_bond, data_tabulated_bond, data_custom_bond, data_harmonic_angle, data_cosine_angle, + data_restrictedcosine_angle, data_fourier_angle, data_ureybradley_angle, + data_polynomial_angle, data_tabulated_angle, ], )