ISISNeutronMuon / ISISNeutronMuon/MDMC
Problem with FoM_abs.py while running water-tutorial
- Dominant language
- Python
- Stars
- 4
- Forks
- 0
- Avg merge
- 2d 4h
- Merged PRs (30d)
- 5
Description
**Description of the error**
An error was obtained while running the water tutorial, which appeared while the refinement.
(Latest master branch!)
The problem is related to operation with arrays in FoM_abs.py
**Describe the expected result**
Success
**Describe the actual result**
```
Universe created with:
Dimensions [21.73, 21.73, 21.73]
Force field None
Number of atoms 0
Total wall time: 0:17:21
Total wall time: 0:09:40
LAMMPS (29 Sep 2021 - Update 3)
using 4 OpenMP thread(s) per MPI task
LAMMPS output is captured by PyLammps wrapper
LAMMPS (29 Sep 2021 - Update 3)
using 4 OpenMP thread(s) per MPI task
LAMMPS output is captured by PyLammps wrapper
Total wall time: 0:00:00
Simulation created with lammps engine and settings:
temperature 263.0
Control created with:
Minimizer MMC
FoM type ChiSquaredExpError
Number of observables 1
Number of parameters 1
Step FoM Change state epsilon
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Input In [4], in ()
86 control.observable_pairs[0].MD_obs.independent_variables = {'E':E, 'Q':Q}
88 # Run refinement
---> 89 control.refine(n_steps=1)
File /usr/local/lib/python3.9/site-packages/MDMC/control/control.py:371, in Control.refine(self, n_steps)
368 self.equilibrate()
370 verbose_manager.header(f"Step {count + 1}")
--> 371 self.step() # advance the refinement by one step
372 count += 1
373 if self.verbose == 3: # if progress bar is there, ensure data is on new line
File /usr/local/lib/python3.9/site-packages/MDMC/control/control.py:430, in Control.step(self)
427 verbose_manager.start(4, verbose=self.verbose)
429 # Generate FoM by running MD for this step and then calculate FoM
--> 430 fom = self._generate_FoM()
432 verbose_manager.step("Selecting new parameters and updating engine")
433 # Select new parameters to consider
File /usr/local/lib/python3.9/site-packages/MDMC/control/control.py:488, in Control._generate_FoM(self)
485 self._run_MD()
486 self._calculate_observables(self.simulation, self.observable_pairs)
--> 488 FoM_value = self.FoM_calculator.calculate()
490 return FoM_value
File /usr/local/lib/python3.9/site-packages/MDMC/refinement/FoM/FoM_abs.py:396, in FigureOfMerit.calculate(self)
380 """
381 Calculates the FoM value by calculating the FoM for every
382 ``ObservablePair``
(...)
392 If calculated value of Figure of Merit is negative
393 """
395 total_weight = np.sum([obs_pair.weight for obs_pair in self.obs_pairs])
--> 396 value_unreduced = np.sum([self.calculate_single_FoM(obs_pair)
397 for obs_pair in self.obs_pairs])
398 self.value = value_unreduced / total_weight
400 assert self.value >= 0.
File /usr/local/lib/python3.9/site-packages/MDMC/refinement/FoM/FoM_abs.py:396, in (.0)
380 """
381 Calculates the FoM value by calculating the FoM for every
382 ``ObservablePair``
(...)
392 If calculated value of Figure of Merit is negative
393 """
395 total_weight = np.sum([obs_pair.weight for obs_pair in self.obs_pairs])
--> 396 value_unreduced = np.sum([self.calculate_single_FoM(obs_pair)
397 for obs_pair in self.obs_pairs])
398 self.value = value_unreduced / total_weight
400 assert self.value >= 0.
File /usr/local/lib/python3.9/site-packages/MDMC/refinement/FoM/ChiSquared_experror.py:83, in ChiSquaredExpError.calculate_single_FoM(self, obs_pair)
78 obs_pair.rescale_factor = (np.sum((MD_values / exp_errors) ** 2)
79 / np.sum(MD_values * exp_values
80 / exp_errors ** 2))
82 norm_factor = self.data_norm_factor(obs_pair=obs_pair)
---> 83 value_unreduced = np.sum((obs_pair.calculate_difference()
84 / obs_pair.calculate_exp_errors()) ** 2)
85 return obs_pair.weight * value_unreduced / norm_factor
File /usr/local/lib/python3.9/site-packages/MDMC/refinement/FoM/FoM_abs.py:293, in ObservablePair.calculate_difference(self)
280 def calculate_difference(self):
281 """
282 Assumes a single dependent variable for each ``Observable``
283
(...)
290 into account.
291 """
--> 293 diff = (np.array(*self.exp_obs.dependent_variables.values())
294 * self.rescale_factor
295 - np.array(*self.MD_obs.dependent_variables.values()))
297 return diff
ValueError: operands could not be broadcast together with shapes (1,71,373) (1,33,373)
```
**Suggested fix**
Fix arrays in a way so that the operation would become doable
**Additional details**
Actual code is below:
```
"""
An example MDMC script for optimizing spce parameters for water at 263 K
Water data provided by Bertil Halle. Ref: J. Chem. Phys. 134, 144508 (2011)
"""
from MDMC.MD.interactions import Bond, BondAngle
from MDMC.control import Control
from MDMC.MD import *
import os
from tests.test_data import data
os.environ["OMP_NUM_THREADS"] = "4"
# Build universe
# Cubic universe of side:
# 18.6270199 A is 216 water molecules
# 21.731523217 is 343 water molecules
# 24.83602653 is 512 water molecules
universe = Universe(dimensions=21.73)
H1 = Atom('H')
H2 = Atom('H', position=(0., 1.63298, 0.))
O = Atom('O', position=(0., 0.81649, 0.57736))
H_coulombic = Coulombic(atoms=[H1, H2], cutoff=10.)
O_coulombic = Coulombic(atoms=O, cutoff=10.)
water_mol = Molecule(position=(0, 0, 0),
velocity=(0, 0, 0),
atoms=[H1, H2, O],
interactions=[Bond((H1, O), (H2, O), constrained=True),
BondAngle(H1, O, H2, constrained=True)],
name='water')
shake = Shake(1e-4, 100)
universe.constraint_algorithm = shake
e_solver = PPPM(accuracy=1e-5)
universe.electrostatic_solver = e_solver
universe.fill(water_mol, num_density=0.03356718472021752)
O_dispersion = Dispersion(universe, (O.atom_type, O.atom_type), cutoff=10.,
vdw_tail_correction=True)
universe.add_force_field('SPCE')
# MD Engine setup
simulation = Simulation(universe,
engine="lammps",
time_step=1.0017365675,
temperature=263.,
traj_step=1054)
# Energy Minimization and equilibration
simulation.minimize(n_steps=5000)
simulation.run(n_steps=25000, equilibration=True)
# Setup refinement
# exp_datasets is a list of dictionaries with one dictionary per experimental dataset
exp_datasets = [{'file_name':'data/263K05Awat_LAMP',
'type':'SQw',
'reader':'LAMPSQw',
'weight':1.,
'resolution':None}]
# Fit parameters is a set(?) of all unique fit parameters in the universe which can then be filtered.
for p in universe.parameters:
if p.name != 'epsilon':
p.fixed = True
fit_parameters = universe.parameters
control = Control(simulation=simulation,
exp_datasets=exp_datasets,
fit_parameters=fit_parameters,
MC_norm=20,
minimizer_type="MMC",
MD_steps=400000,
energy_resolution=13.5)
# Bertil Halle water data is non-symmetric. Consider only a subset of the
# data in this example.
# So that the MD simulation size can be minimized, the Q min is increased and
# the Q resolution is reduced.
exp_obs = control.observable_pairs[0].exp_obs
Q_slice = slice(6, len(exp_obs.Q), 2)
Q = exp_obs.Q[Q_slice]
E_range = (exp_obs.E >=0)
E = exp_obs.E[E_range]
# copy the updated E values, and Q values back to the control.observable
control.observable_pairs[0].exp_obs.independent_variables = {'E':E, 'Q':Q}
control.observable_pairs[0].MD_obs.independent_variables = {'E':E, 'Q':Q}
# Run refinement
control.refine(n_steps=1)
```
Contributor guide
No contributing guide indexed for this repository
Assessment
This issue has not been assessed yet.