diff --git a/environment-dev.yml b/environment-dev.yml index 7a0dd6a..1ffc9ee 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -3,7 +3,7 @@ channels: - conda-forge dependencies: - cuda-version=11.8 -- freud >=2.13 +- freud =2.13.2 - gsd >=3.0 - hoomd >=4.0 - python >=3.9 diff --git a/environment.yml b/environment.yml index b060d7d..1bd48ea 100644 --- a/environment.yml +++ b/environment.yml @@ -3,7 +3,7 @@ channels: - conda-forge dependencies: - cuda-version=11.8 -- freud >=2.13 +- freud =2.13.2 - gsd >=3.0 - hoomd >=4.0 - python >=3.9 diff --git a/msibi/forces.py b/msibi/forces.py index bdb2719..9bfb0a8 100644 --- a/msibi/forces.py +++ b/msibi/forces.py @@ -84,12 +84,6 @@ def __init__( self._tail_correction_history = [] self._learned_potential_history = [] - if optimize and nbins is None: - raise ValueError( - "If this force is set to be optimized, the nbins " - "must be set as a non-zero value" - ) - def __repr__(self): return ( f"Type: {self.__class__}; " @@ -131,7 +125,7 @@ def smoothing_window(self) -> int: @smoothing_window.setter def smoothing_window(self, value: int): - if not isinstance(value, int): + if not isinstance(value, int) or value <= 0: raise ValueError("The smoothing window must be an integer.") self._smoothing_window = value for state in self._states: @@ -144,7 +138,7 @@ def smoothing_order(self) -> int: @smoothing_order.setter def smoothing_order(self, value: int): - if not isinstance(value, int): + if not isinstance(value, int) or value <= 0: raise ValueError("The smoothing order must be an integer.") self._smoothing_order = value for state in self._states: @@ -157,7 +151,7 @@ def nbins(self) -> int: @nbins.setter def nbins(self, value: int): - if not isinstance(value, int): + if not isinstance(value, int) or value <= 0: raise ValueError("nbins must be an integer.") self._nbins = value for state in self._states: diff --git a/msibi/tests/test_forces.py b/msibi/tests/test_forces.py index e1d2c7d..b3d0424 100644 --- a/msibi/tests/test_forces.py +++ b/msibi/tests/test_forces.py @@ -31,6 +31,7 @@ def test_potential_setter(self, bond): initial_pot = np.copy(bond.potential) bond.potential = bond.potential * 2 assert np.allclose(bond.potential, initial_pot * 2) + assert bond.format == "table" def test_smooth_potential(self, bond): bond.set_quadratic( @@ -49,6 +50,37 @@ def test_smooth_potential(self, bond): for i, j in zip(bond.potential, noisy_pot): assert i != j + def test_set_from_file(self): + bond = Bond(type1="A", type2="B", optimize=True, nbins=60) + bond.set_quadratic(x_min=0.0, x_max=3.0, x0=1, k2=200, k3=0, k4=0) + bond.save_potential("test.csv") + bond2 = Bond(type1="A", type2="B", optimize=True, nbins=60) + bond2.set_from_file("test.csv") + assert np.allclose(bond.potential, bond2.potential) + os.remove("test.csv") + + def test_fit_scores(self, msibi, stateX, stateY): + msibi.gsd_period = 10 + bond = Bond(type1="A", type2="B", optimize=True, nbins=60) + bond.set_quadratic(x_min=0.0, x_max=3.0, x0=1, k2=200, k3=0, k4=0) + angle = Angle(type1="A", type2="B", type3="A", optimize=False) + angle.set_harmonic(k=500, t0=2) + angle2 = Angle(type1="B", type2="A", type3="B", optimize=False) + angle2.set_harmonic(k=500, t0=2) + msibi.add_state(stateX) + msibi.add_state(stateY) + msibi.add_force(bond) + msibi.add_force(angle) + msibi.add_force(angle2) + init_bond_pot = np.copy(bond.potential) + msibi.run_optimization(n_steps=500, n_iterations=1) + assert len(bond._states[stateY]["f_fit"]) == 1 + assert len(bond._states[stateX]["f_fit"]) == 1 + bond.plot_fit_scores(state=stateY) + bond.plot_fit_scores(state=stateX) + with pytest.raises(RuntimeError): + angle.plot_fit_scores(state=stateY) + def test_smoothing_window(self, bond): bond.smoothing_window = 5 assert bond.smoothing_window == 5 @@ -61,6 +93,48 @@ def test_nbins(self, bond): bond.nbins = 60 assert bond.nbins == 60 + def test_set_potential_error(self): + bond = Bond(type1="A", type2="B", optimize=False) + bond.set_harmonic(k=500, r0=2) + with pytest.raises(ValueError): + bond.potential = np.array([1, 2, 3]) + + def test_plot_target_dist(self, stateX): + angle = Angle(type1="A", type2="B", type3="A", optimize=True, nbins=60) + angle.set_quadratic(x0=2, k4=0, k3=0, k2=100, x_min=0, x_max=np.pi) + angle._add_state(stateX) + angle.plot_target_distribution(state=stateX) + + def test_static_warnings(self): + bond = Bond(type1="A", type2="B", optimize=False) + bond.set_harmonic(k=500, r0=2) + with pytest.warns(UserWarning): + bond.potential + bond.force + + def test_bad_smoothing_args(self, bond): + with pytest.raises(ValueError): + bond.smoothing_window = 0 + with pytest.raises(ValueError): + bond.smoothing_window = 4.5 + with pytest.raises(ValueError): + bond.smoothing_order = 0 + with pytest.raises(ValueError): + bond.smoothing_order = 2.2 + with pytest.raises(ValueError): + bond.nbins = 20.5 + with pytest.raises(ValueError): + bond.nbins = 0 + angle = Angle(type1="A", type2="B", type3="A", optimize=False) + angle.set_harmonic(k=500, t0=2) + with pytest.raises(RuntimeError): + angle.smooth_potential() + + def test_save_static_force(self): + angle = Angle(type1="A", type2="B", type3="A", optimize=False) + angle.set_harmonic(k=500, t0=2) + with pytest.raises(RuntimeError): + angle.save_potential("test.csv") class TestBond(BaseTest): def test_bond_name(self, bond): diff --git a/msibi/tests/test_optimize.py b/msibi/tests/test_optimize.py index 6d47069..5f2db61 100644 --- a/msibi/tests/test_optimize.py +++ b/msibi/tests/test_optimize.py @@ -49,6 +49,11 @@ def test_run(self, msibi, stateX, stateY): assert msibi.n_iterations == 1 ff = msibi._build_force_objects() assert len(ff) == 1 + assert len(bond.distribution_history(state=stateX)) == 1 + assert len(bond.potential_history) == 2 + assert len(bond._head_correction_history) == 1 + assert len(bond._tail_correction_history) == 1 + assert len(bond._learned_potential_history) == 1 def test_run_with_static_force(self, msibi, stateX, stateY): msibi.gsd_period = 10 @@ -68,4 +73,55 @@ def test_run_with_static_force(self, msibi, stateX, stateY): assert not np.array_equal(bond.potential, init_bond_pot) assert msibi.n_iterations == 1 ff = msibi._build_force_objects() - assert len(ff) == 2 \ No newline at end of file + assert len(ff) == 2 + + def test_run_with_all_forces(self, msibi, stateX, stateY): + msibi.add_state(stateX) + msibi.add_state(stateY) + + bond = Bond(type1="A", type2="B", optimize=False, nbins=60) + bond.set_harmonic(r0=1.1, k=100) + msibi.add_force(bond) + + angle = Angle(type1="A", type2="B", type3="A", optimize=False) + angle.set_harmonic(t0=1.9, k=100) + msibi.add_force(angle) + angle2 = Angle(type1="B", type2="A", type3="B", optimize=False) + angle2.set_harmonic(t0=2.3, k=100) + msibi.add_force(angle2) + + pair = Pair(type1="A", type2="B", r_cut=3.0, nbins=60, optimize=True, exclude_bonded=True) + pair.set_lj(sigma=1.5, epsilon=1, r_cut=3.0, r_min=0.1) + msibi.add_force(pair) + pair2 = Pair(type1="A", type2="A", r_cut=3.0, nbins=60, optimize=True, exclude_bonded=True) + pair2.set_lj(sigma=2, epsilon=2, r_cut=3.0, r_min=0.1) + msibi.add_force(pair2) + pair3 = Pair(type1="B", type2="B", r_cut=3.0, nbins=60, optimize=True, exclude_bonded=True) + pair3.set_lj(sigma=1.5, epsilon=1, r_cut=3.0, r_min=0.1) + msibi.add_force(pair3) + dihedral = Dihedral(type1="A", type2="B", type3="A", type4="B", optimize=False) + dihedral.set_harmonic(k=100, phi0=0, d=-1, n=1) + msibi.add_force(dihedral) + + msibi.run_optimization(n_steps=500, n_iterations=0) + + def test_raise_errors(self, msibi, stateX, stateY): + with pytest.raises(RuntimeError): + msibi.pickle_forces(file_path="test.pkl") + + bond = Bond(type1="A", type2="B", optimize=True, nbins=60) + angle = Angle(type1="A", type2="B", type3="A", optimize=True, nbins=60) + with pytest.raises(RuntimeError): + msibi.add_force(bond) + msibi.add_force(angle) + + with pytest.raises(ValueError): + msibi = MSIBI( + nlist=hoomd.md.nlist.Cell, + integrator_method=hoomd.md.methods.DisplacementCapped, + method_kwargs=dict(), + thermostat=hoomd.md.methods.thermostats.MTTK, + thermostat_kwargs=dict(tau=0.01), + dt=0.003, + gsd_period=int(1e3), + ) \ No newline at end of file diff --git a/msibi/utils/exceptions.py b/msibi/utils/exceptions.py deleted file mode 100644 index f27cf2a..0000000 --- a/msibi/utils/exceptions.py +++ /dev/null @@ -1,10 +0,0 @@ -SUPPORTED_ENGINES = ["hoomd"] - -class UnsupportedEngine(Exception): - def __init__(self, engine): - message = ( - 'Unsupported engine: "{0}". Supported engines are: {1}'.format( - engine, ", ".join(SUPPORTED_ENGINES) - ) - ) - super(UnsupportedEngine, self).__init__(message) diff --git a/msibi/utils/hoomd_run_template.py b/msibi/utils/hoomd_run_template.py deleted file mode 100755 index 9527985..0000000 --- a/msibi/utils/hoomd_run_template.py +++ /dev/null @@ -1,53 +0,0 @@ -HOOMD3_HEADER = """ -import hoomd - -device = hoomd.device.auto_select() -sim = hoomd.Simulation(device=device) -sim.create_state_from_snapshot({0}) -nl = {1}() -integrator = hoomd.md.Integrator(dt={2}) -method_kwargs = {3} -method = {method}( - kT={4}, - filter=hoomd.filter.All(), - **method_kwargs -) -integrator.methods.append(method) - -sim.operations.integrator = integrator -sim.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT={4}) -""" - -HOOMD3_TEMPLATE = """ - -""" - -HOOMD2_HEADER = """ -import hoomd -import hoomd.md -from hoomd.init import read_gsd - -hoomd.context.initialize("") -system = read_gsd("{0}", frame=-1, time_step=0) - -nl = {1}() -nl.reset_exclusions(exclusions={2}) -""" - -HOOMD_TEMPLATE = """ -_all = hoomd.group.all() -hoomd.md.integrate.mode_standard({dt}) -integrator_kwargs = {integrator_kwargs} -integrator = {integrator}(group=_all, **integrator_kwargs) - - -hoomd.dump.gsd( - filename="query.gsd", - group=_all, - period={gsd_period}, - overwrite=True, - dynamic=["momentum"] - ) - -hoomd.run({n_steps}) -"""