From 0af786cbda6ec74b6cee3ef8c69a9c67e0399cd6 Mon Sep 17 00:00:00 2001 From: jrudz Date: Tue, 24 Sep 2024 16:41:45 +0200 Subject: [PATCH] finished tabulated potential norm and testing --- .../schema_packages/force_field.py | 66 +++--- .../schema_packages/outputs.py | 2 - tests/test_force_field.py | 189 +++++++++++------- 3 files changed, 157 insertions(+), 100 deletions(-) diff --git a/src/nomad_simulations/schema_packages/force_field.py b/src/nomad_simulations/schema_packages/force_field.py index 4fa7960a..dfb6fbcd 100644 --- a/src/nomad_simulations/schema_packages/force_field.py +++ b/src/nomad_simulations/schema_packages/force_field.py @@ -251,49 +251,61 @@ def normalize(self, archive, logger) -> None: - np.min(self.energies.to('kJ').magnitude * MOL) ) if np.all([x < tol for x in energies_diff]): - print('consistent energies') logger.warning( - f'Forces were generated from energies in {self.name},' - f' but appear to be roughly consistent, with rtol={tol}. ' + f'Tabulated forces were generated from energies in {self.name},' + f'with consistency errors less than tol={tol}. ' ) else: - # logger.warning('inconsistent energies') - print('inconsistent energies') + logger.warning( + f'Unable to derive tabulated forces from energies in {self.name},' + f'consistency errors were greater than tol={tol}.' + ) self.forces = None - # print(energies.to('kJ').magnitude * MOL) - # print( - # self.energies.to('kJ').magnitude * MOL - # - np.min(self.energies.to('kJ').magnitude * MOL) - # ) - # print(energies_diff) - except Exception: - pass + except ValueError as e: + logger.warning( + f'Unable to derive tabulated forces from energies in {self.name},' + f'Unkown error occurred in derivation: {e}' + ) 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( - self.bins, self.forces, smoothing_factor=smoothing_factor + self.bins.magnitude, + self.forces.magnitude, + smoothing_factor=smoothing_factor, ) # re-derive forces to check consistency of the energies - forces = self.compute_forces( - self.bins, self.energies, smoothing_factor=smoothing_factor + forces = ( + self.compute_forces( + self.bins.magnitude, + self.energies.magnitude, + smoothing_factor=smoothing_factor, + ) + * ureg.J + / self.bins.units ) - forces_diff = forces.magnitude * 1000.0 * MOL - ( - self.energies.to('kJ').magnitude * 1000.0 * MOL - - np.min(self.energies.to('kJ').magnitude * 1000.0 * MOL) + forces_diff = forces.to(f'kJ/{self.bins.units}').magnitude * MOL - ( + self.forces.to(f'kJ/{self.bins.units}').magnitude * MOL ) if np.all([x < tol for x in forces_diff]): - print('consistent forces') - # logger.warning( - # f'Energies were generated from forces in {self.name},' - # f'but appear to be consistent with rtol={rtol}. ' - # ) + logger.warning( + f'Tabulated energies were generated from forces in {self.name},' + f'with consistency errors less than tol={tol}. ' + ) else: + logger.warning( + f'Unable to derive tabulated energies from forces in {self.name},' + f'consistency errors were greater than tol={tol}.' + ) self.energies = None - except Exception: - pass + except ValueError as e: + logger.warning( + f'Unable to derive tabulated energies from forces in {self.name},' + f'Unkown error occurred in derivation: {e}' + ) class BondPotential(Potential): @@ -534,7 +546,7 @@ def normalize(self, archive, logger) -> None: logger.warning('Incorrect functional form set for HarmonicAngle.') -class TabulatedAngle(BondPotential, TabulatedPotential): +class TabulatedAngle(AnglePotential, 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. diff --git a/src/nomad_simulations/schema_packages/outputs.py b/src/nomad_simulations/schema_packages/outputs.py index 4b673ef9..d1f050b3 100644 --- a/src/nomad_simulations/schema_packages/outputs.py +++ b/src/nomad_simulations/schema_packages/outputs.py @@ -66,8 +66,6 @@ class Outputs(ArchiveSection): information if the output `is_derived` from another output section or directly parsed from the simulation output files. """ - # TODO add time quantities - normalizer_level = 2 model_system_ref = Quantity( diff --git a/tests/test_force_field.py b/tests/test_force_field.py index 5516ea40..1280260d 100644 --- a/tests/test_force_field.py +++ b/tests/test_force_field.py @@ -318,6 +318,28 @@ def populate_parameters(sec_potential, parameters): 7.19013405e-21, 8.38572215e-21, ], + 'forces': [ + 1.3216958010639183e-10, + 1.0784837919852664e-10, + 8.347183948070384e-11, + 5.903996095292366e-11, + 3.521528791034704e-11, + 1.0674227406164319e-11, + -1.3922171907975959e-11, + -3.857391003207359e-11, + -6.32809869661286e-11, + -8.804340271014125e-11, + -1.1218967953067358e-10, + -1.3706127725108878e-10, + -1.6198821378146133e-10, + -1.8697048912179138e-10, + -2.1200810327207885e-10, + -2.3710105623232373e-10, + -2.6156893683081195e-10, + -2.867710717674596e-10, + -3.1202854551406456e-10, + -3.373413580706269e-10, + ], 'name': 'TabulatedBond', 'type': 'bond', 'functional_form': 'tabulated', @@ -454,51 +476,73 @@ def populate_parameters(sec_potential, parameters): 'particle_indices': [[0, 1, 2], [3, 4, 5]], 'particle_labels': [['O', 'H', 'H'], ['O', 'H', 'H']], 'bins': [ - 7.600e-11, - 7.970e-11, - 8.340e-11, - 8.710e-11, - 9.070e-11, - 9.440e-11, - 9.810e-11, - 1.018e-10, - 1.055e-10, - 1.092e-10, - 1.128e-10, - 1.165e-10, - 1.202e-10, - 1.239e-10, - 1.276e-10, - 1.313e-10, - 1.349e-10, - 1.386e-10, - 1.423e-10, - 1.460e-10, + 92.99105014973262, + 94.19727699638311, + 95.40350384303362, + 96.60973068968411, + 97.81595810929241, + 99.02218495594292, + 100.22841180259343, + 101.43463864924392, + 102.64086549589443, + 103.84709234254493, + 105.05331976215322, + 106.25954660880373, + 107.46577345545423, + 108.67200030210473, + 109.87822714875524, + 111.08445399540572, + 112.29068141501403, + 113.49690826166454, + 114.70313510831502, + 115.90936195496555, ], 'energies': [ - 1.32311751e-21, - 8.81248069e-22, - 5.28549577e-22, - 2.65354139e-22, - 9.18278089e-23, - 8.30269520e-24, - 1.47787975e-23, - 1.11422170e-22, - 2.98564919e-22, - 5.76539155e-22, - 9.45178822e-22, - 1.40498208e-21, - 1.95611499e-21, - 2.59857754e-21, - 3.33286791e-21, - 4.15882003e-21, - 5.07693206e-21, - 6.08737007e-21, - 7.19013405e-21, - 8.38572215e-21, + 2.483908793147618e-21, + 1.9871270690593166e-21, + 1.5455433097527105e-21, + 1.1591575152954348e-21, + 8.279695415419479e-22, + 5.51979703171419e-22, + 3.311878298023737e-22, + 1.6559392146862525e-22, + 5.519797819553163e-23, + 0.0, + 0.0, + 5.519797819553239e-23, + 1.6559392146862525e-22, + 3.311878298023748e-22, + 5.519797031714216e-22, + 8.2796954154194335e-22, + 1.1591575152954338e-21, + 1.545543309752715e-21, + 1.98712706905931e-21, + 2.483908793147623e-21, ], - 'name': 'TabulatedBond', - 'type': 'bond', + 'forces': [ + 4.347281042004185e-22, + 3.8896725108092945e-22, + 3.432063979614403e-22, + 2.974455448419512e-22, + 2.516846920122808e-22, + 2.0592383889279167e-22, + 1.6016298577330257e-22, + 1.1440213265381345e-22, + 6.864127953432432e-23, + 2.288042641483519e-23, + -2.288042641483519e-23, + -6.864127953432432e-23, + -1.1440213265381345e-22, + -1.6016298577330257e-22, + -2.0592383889279167e-22, + -2.516846920122808e-22, + -2.974455448419512e-22, + -3.432063979614403e-22, + -3.8896725108092945e-22, + -4.347281042004185e-22, + ], + 'name': 'TabulatedAngle', + 'type': 'angle', 'functional_form': 'tabulated', } data_tabulated_angle = ( @@ -531,30 +575,30 @@ def populate_parameters(sec_potential, parameters): 2.023, ] * ureg.radian, - 'energies': [ - 1.5, - 1.20083102, - 0.93490305, - 0.70221607, - 0.50277008, - 0.3365651, - 0.20360111, - 0.10387812, - 0.03739612, - 0.00415512, - 0.00415512, - 0.03739612, - 0.10387812, - 0.20360111, - 0.3365651, - 0.50277008, - 0.70221607, - 0.93490305, - 1.20083102, - 1.5, - ] - * ureg.kJ - / MOL, + # 'energies': [ # ! These are the exact energies. We will test the auto-generation of these energies from the forces. + # 1.5, + # 1.20083102, + # 0.93490305, + # 0.70221607, + # 0.50277008, + # 0.3365651, + # 0.20360111, + # 0.10387812, + # 0.03739612, + # 0.00415512, + # 0.00415512, + # 0.03739612, + # 0.10387812, + # 0.20360111, + # 0.3365651, + # 0.50277008, + # 0.70221607, + # 0.93490305, + # 1.20083102, + # 1.5, + # ] + # * ureg.kJ + # / MOL, 'forces': [ 15.0, 13.42105263, @@ -576,7 +620,10 @@ def populate_parameters(sec_potential, parameters): -11.84210526, -13.42105263, -15.0, - ], + ] + * ureg.kJ + / MOL + / ureg.radian, }, results_tabulated_angle, ) @@ -591,8 +638,8 @@ def populate_parameters(sec_potential, parameters): # data_fene_bond, data_tabulated_bond, # data_custom_bond, - # data_harmonic_angle, - # data_tabulated_angle, + data_harmonic_angle, + data_tabulated_angle, ], ) def test_potentials( @@ -626,10 +673,10 @@ def test_potentials( sec_FF.contributions[-1].normalize(EntryArchive, logger) # BoundLogger) potential_dict = sec_FF.contributions[-1].m_to_dict() - potential_dict = { + potential_dict = { # ! The dev is required to add new results to the dictionary, this will not be caught by the test! key: value for key, value in potential_dict.items() if key in results } print(potential_dict) print(results) - assert 1 == 2 + # assert 1 == 2 assert_dict_equal(potential_dict, results)