From bbf213c656e02b96a44d4b22f6674619bf9c46fc Mon Sep 17 00:00:00 2001 From: jrudz Date: Tue, 10 Sep 2024 14:15:04 +0200 Subject: [PATCH] added previous ff schema from older version --- .../schema_packages/force_field.py | 480 ++++++++++++++++++ 1 file changed, 480 insertions(+) create mode 100644 src/nomad_simulations/schema_packages/force_field.py diff --git a/src/nomad_simulations/schema_packages/force_field.py b/src/nomad_simulations/schema_packages/force_field.py new file mode 100644 index 00000000..b99ed980 --- /dev/null +++ b/src/nomad_simulations/schema_packages/force_field.py @@ -0,0 +1,480 @@ +import numpy as np +import pint +import re +import typing +from structlog.stdlib import BoundLogger +from typing import Optional, List, Tuple +from ase.dft.kpoints import monkhorst_pack, get_monkhorst_pack_size_and_offset + +from nomad.units import ureg +from nomad.datamodel.data import ArchiveSection +from nomad.datamodel.metainfo.annotations import ELNAnnotation +from nomad.metainfo import ( + Quantity, + SubSection, + MEnum, + Section, + Context, + JSON, +) + +from .model_method import ModelMethod, NumericalSettings +from .model_system import ModelSystem +from .utils import is_not_representative + + +class ParamEntry(ArchiveSection): + """ + Generic section defining a parameter name and value + """ + + name = Quantity( + type=str, + shape=[], + description=""" + Name of the parameter. + """, + ) + + value = Quantity( + type=str, + shape=[], + description=""" + Value of the parameter as a string. + """, + ) + + unit = Quantity( + type=str, + shape=[], + description=""" + Unit of the parameter as a string. + """, + ) + + # TODO add description quantity + + +class ForceCalculations(NumericalSettings): + """ + Section containing the parameters for force calculations according to the referenced force field + during a molecular dynamics run. + """ + + vdw_cutoff = Quantity( + type=np.float64, + shape=[], + unit="m", + description=""" + Cutoff for calculating VDW forces. + """, + ) + + coulomb_type = Quantity( + type=MEnum( + "cutoff", + "ewald", + "multilevel_summation", + "particle_mesh_ewald", + "particle_particle_particle_mesh", + "reaction_field", + ), + shape=[], + description=""" + Method used for calculating long-ranged Coulomb forces. + + Allowed values are: + + | Barostat Name | Description | + + | ---------------------- | ----------------------------------------- | + + | `""` | No thermostat | + + | `"Cutoff"` | Simple cutoff scheme. | + + | `"Ewald"` | Standard Ewald summation as described in any solid-state physics text. | + + | `"Multi-Level Summation"` | D. Hardy, J.E. Stone, and K. Schulten, + [Parallel. Comput. **35**, 164](https://doi.org/10.1016/j.parco.2008.12.005)| + + | `"Particle-Mesh-Ewald"` | T. Darden, D. York, and L. Pedersen, + [J. Chem. Phys. **98**, 10089 (1993)](https://doi.org/10.1063/1.464397) | + + | `"Particle-Particle Particle-Mesh"` | See e.g. Hockney and Eastwood, Computer Simulation Using Particles, + Adam Hilger, NY (1989). | + + | `"Reaction-Field"` | J.A. Barker and R.O. Watts, + [Mol. Phys. **26**, 789 (1973)](https://doi.org/10.1080/00268977300102101)| + """, + ) + + coulomb_cutoff = Quantity( + type=np.float64, + shape=[], + unit="m", + description=""" + Cutoff for calculating short-ranged Coulomb forces. + """, + ) + + neighbor_update_frequency = Quantity( + type=int, + shape=[], + description=""" + Number of timesteps between updating the neighbor list. + """, + ) + + neighbor_update_cutoff = Quantity( + type=np.float64, + shape=[], + unit="m", + description=""" + The distance cutoff for determining the neighbor list. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class Potential(ArchiveSection): + """ + Section containing information about an interaction potential. + """ + + type = Quantity( + type=str, + shape=[], + description=""" + Specifies the type (i.e., functional form) of the interaction potential. + """, + ) + + parameters = Quantity( + type=dict[str, np.float64], # ? Used to be typing.ANY, not sure what I want here + shape=[], + description=""" + Dictionary of label and parameters of the interaction potential. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class BondPotential(Potential): + """ + Section containing information about bond potentials. + + Suggested types are: harmonic, cubic, Morse, fene, tabulated + """ + + equilibrium_value = Quantity( + type=np.float64, + unit='m', + shape=[], + description=""" + Specifies the equilibrium bond distance. + """, + ) + + force_constant = Quantity( + type=np.float64, + shape=[], + unit='J / m **2', + description=""" + Specifies the force constant of the bond potential. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class HarmonicBond(BondPotential): + """ + Section containing information about a Harmonic bond potential: U(x) = 1/2 k (x-x_0)^2 + C, + where k is the `force_constant` and x_0 is the `equilibrium_value` of x. + C is an arbitrary constant (not stored). + """ + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + self.type = 'harmonic' + +class CubicBond(BondPotential): + """ + Section containing information about a Cubic bond potential: U(x) = 1/2 k (x-x_0)^2 + 1/3 k_c (x-x_0)^3, + where k is the (harmonic) `force_constant`, k_c is the `force_constant_cubic`, and x_0 is the `equilibrium_value` of x. + 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.type = 'cubic' + +class MorseBond(BondPotential): + """ + Section containing information about a Morse potential: U(x) = D [1 - exp(- a (x-x_0)]^2 + C, + where a = sqrt(k/2D) is the `well_steepness`, with `force constant` k. + D is the `well_depth`, and x_0 is the `equilibrium_value` of x. C is an arbitrary constant (not stored). + """ + + well_depth = Quantity( + type=np.float64, + unit='J', + shape=[], + description=""" + Specifies the depth of the potential well. + """, + ) + + well_steepness = Quantity( + type=np.float64, + unit='1/m', + shape=[], + description=""" + Specifies the steepness of the potential well. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + self.type = 'Morse' + if self.well_depth is not None and self.well_steepness is not None: + self.force_constant = 2. * self.well_depth * self.well_steepness**2 + elif self.well_depth is not None and self.force_constant is not None: + self.well_steepness = np.sqrt(self.force_constant / (2. * self.well_depth)) + +class FenePotential(BondPotential): + """ + Section containing information about a FENE potential: U(x) = -1/2 k (x_0)^2 ln[1-(x^2)/(x_0^2)] + C, + k is the `force_constant`, and x_0 is the `equilibrium_value` of x. C is an arbitrary constant (not stored). + """ + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + self.type = 'fene' + +class TabulatedBond(BondPotential): + """ + 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='m', + shape=[], + description=""" + List of bin distances. + """, + ) + + energies = Quantity( + type=np.float64, + unit='J', + shape=[], + description=""" + List of energy values associated with each bin. + """, + ) + + forces = Quantity( + type=np.float64, + unit='J/m', + shape=[], + description=""" + List of force values associated with each bin. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + self.type = 'tabulated' + if self.bins is not None and self.energies is not None: + if len(self.bins != self.energies): + logger.error('bins and energies values have different length in TabulatedBond') + if self.bins is not None and self.forces is not None: + if len(self.bins != self.forces): + logger.error('bins and forces have different length in TabulatedBond') + +class AnglePotential(Potential): + """ + Section containing information about bond potentials. + + Suggested types are: harmonic, tabulated + """ + + equilibrium_value = Quantity( + type=np.float64, + unit='degree', + shape=[], + description=""" + Specifies the equilibrium angle. + """, + ) + + force_constant = Quantity( + type=np.float64, + shape=[], + unit='J / degree**2', + description=""" + Specifies the force constant of the angle potential. + """, + ) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class HarmonicAngle(AnglePotential): + """ + Section containing information about a Harmonic angle potential: U(x) = 1/2 k (x-x_0)^2 + C, + where k is the `force_constant` and x_0 is the `equilibrium_value` of the angle x. + C is an arbitrary constant (not stored). + """ + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + self.type = 'harmonic' + +class TabulatedAngle(AnglePotential): + """ + Section containing information about a tabulated angle 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 distances. + """, + ) + + energies = Quantity( + type=np.float64, + unit='J', + shape=[], + description=""" + List of energy values associated with each bin. + """, + ) + + 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.type = 'tabulated' + if self.bins is not None and self.energies is not None: + if len(self.bins != self.energies): + logger.error('bins and energies values have different length in TabulatedBond') + if self.bins is not None and self.forces is not None: + if len(self.bins != self.forces): + logger.error('bins and forces have different length in TabulatedBond') + +class Interactions(ArchiveSection): + """ + Section containing the list of particles involved in a particular type of interaction and/or any associated parameters. + """ + + type = Quantity( + type=MEnum("bond", "angle", "dihedral", "improper dihedral"), + shape=[], + description=""" + Denotes the classification of the interaction. + """, + ) + + name = Quantity( + type=str, + shape=[], + description=""" + Specifies the name of the interaction. Can contain information on the species, + cut-offs, potential versions, etc. + """, + ) + + n_interactions = Quantity( + type=np.dtype(np.int32), + shape=[], + description=""" + Total number of interactions of this interaction type-name. + """, + ) + + n_particles = Quantity( + type=np.int32, + shape=[], + description=""" + Number of particles included in (each instance of) the interaction. + """, + ) + + particle_labels = Quantity( + type=np.dtype(str), + shape=["n_interactions", "n_atoms"], + description=""" + Labels of the particles described by the interaction. In general, the structure is a list of list of tuples. + """, + ) + + particle_indices = Quantity( + type=np.int32, + shape=["n_interactions", "n_atoms"], + description=""" + Indices of the particles in the system described by the interaction. In general, the structure is a list of list of tuples. + """, + ) + + potential = SubSection(sub_section=Potential.m_def, repeats=False) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + +class ForceField(ModelMethod): + """ + Section containing the parameters of a (classical, particle-based) force field model. + Typical `numerical_settings` are ForceCalculations and NeighborSearching. + Lists of interactions by type and, if available, corresponding parameters can be given within `interactions`. + Additionally, a published model can be referenced with `reference`. + """ + + name = Quantity( + type=str, + shape=[], + description=""" + Identifies the name of the model. + """, + ) + + reference = Quantity( + type=str, + shape=[], + description=""" + Reference to the model e.g. DOI, URL. + """, + ) + + interactions = SubSection(sub_section=Interactions.m_def, repeats=True) + + def normalize(self, archive, logger) -> None: + super().normalize(archive, logger) + + self.name = 'ForceField'