diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 3c2e91f..0b6eedc 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,6 +1,13 @@ repos: - - repo: https://github.com/psf/black - rev: 23.7.0 +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.4.0 + hooks: + - id: trailing-whitespace + - id: end-of-file-fixer + - id: check-yaml + - id: check-added-large-files +- repo: https://github.com/psf/black + rev: 23.9.1 hooks: - id: black-jupyter # It is recommended to specify the latest version of Python @@ -8,3 +15,9 @@ repos: # pre-commit's default_language_version, see # https://pre-commit.com/#top_level-default_language_version language_version: python3.11 +- repo: https://github.com/astral-sh/ruff-pre-commit + # Ruff version. + rev: v0.0.292 + hooks: + - id: ruff + types_or: [python, pyi, jupyter] diff --git a/examples/lennard.py b/examples/lennard.py index 222a5ca..be4d512 100644 --- a/examples/lennard.py +++ b/examples/lennard.py @@ -23,9 +23,7 @@ def create_system(): for xyz in pos: particles.add_particle(pos=xyz, mass=1.0, ptype=0, name="Ar") - generate_maxwell_velocities( - particles, rgen=default_rng(), temperature=0.8 - ) + generate_maxwell_velocities(particles, rgen=default_rng(), temperature=0.8) lennard = LennardJonesCut(dim=3, shift=True, mixing="geometric") parameters = { diff --git a/examples/movie-doublewell-pair.py b/examples/movie-doublewell-pair.py index ba40769..84bab0c 100644 --- a/examples/movie-doublewell-pair.py +++ b/examples/movie-doublewell-pair.py @@ -13,7 +13,6 @@ from turtlemd.tools import generate_lattice sns.set_context("talk", font_scale=0.8) -# matplotlib.rcParams.update({'font.size': 10}) def create_system() -> System: diff --git a/pyproject.toml b/pyproject.toml index bf51ae2..a0710a6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "turtlemd" -version = "2023.3.1" +version = "2024.1" description = "A slow molecular dynamics library for testing" authors = [ "Daniel Zhang ", @@ -13,6 +13,8 @@ readme = "README.md" [tool.poetry.dependencies] python = ">=3.10,<3.13" numpy = "^1.24.2" +toml = "^0.10.2" +pint = "^0.23" [tool.poetry.group.dev.dependencies] ruff = "^0.0.261" @@ -24,6 +26,7 @@ pytest-randomly = "^3.12.0" jupyter-black = "^0.3.4" pre-commit = "^3.3.3" flake8 = "^6.0.0" +types-toml = "^0.10.8.7" [tool.ruff] line-length = 79 @@ -39,6 +42,9 @@ warn_unused_ignores = true show_error_codes = true plugins = "numpy.typing.mypy_plugin" +[tool.poetry.scripts] +turtlemd = "turtlemd.bin:main" + [build-system] requires = ["poetry-core"] build-backend = "poetry.core.masonry.api" diff --git a/tests/inout/box1.toml b/tests/inout/box1.toml new file mode 100644 index 0000000..a3bab67 --- /dev/null +++ b/tests/inout/box1.toml @@ -0,0 +1,4 @@ +[box] +low = [10, 11, 12] +high = [23, 24, 25] +periodic = [true, true, true] diff --git a/tests/inout/box2.toml b/tests/inout/box2.toml new file mode 100644 index 0000000..8a1d982 --- /dev/null +++ b/tests/inout/box2.toml @@ -0,0 +1,4 @@ +[box] +low = [10, 11, 12] +high = [23, 24, 25] +periodic = true diff --git a/tests/inout/box3.toml b/tests/inout/box3.toml new file mode 100644 index 0000000..bc9157c --- /dev/null +++ b/tests/inout/box3.toml @@ -0,0 +1,6 @@ +[box] +low = [10, 11, 12] +high = [23, 24, 25] +periodic = [true, true, true] +alpha = 75 +beta = 60 diff --git a/tests/inout/config/start.xyz b/tests/inout/config/start.xyz new file mode 100644 index 0000000..bd1be80 --- /dev/null +++ b/tests/inout/config/start.xyz @@ -0,0 +1,4 @@ +2 +# Step 2 box: 1.0 2.0 3.0 +A 1.0 1.0 1.0 1.0 1.0 1.0 +B 2.0 2.0 2.0 1.0 1.0 1.0 diff --git a/tests/inout/config/system_one_down.toml b/tests/inout/config/system_one_down.toml new file mode 100644 index 0000000..8b3ab13 --- /dev/null +++ b/tests/inout/config/system_one_down.toml @@ -0,0 +1,7 @@ +[box] +periodic = [true, true, true] +low = [0.0, 0.0, 0.0] +high = [1.0, 1.0, 1.0] + +[particles] +file = "../start2.xyz" diff --git a/tests/inout/integrators/integrator.toml b/tests/inout/integrators/integrator.toml new file mode 100644 index 0000000..e1b8d83 --- /dev/null +++ b/tests/inout/integrators/integrator.toml @@ -0,0 +1,3 @@ +[integrator] +class = "this-is-integrator-hello" +timestep = 1234.5678 diff --git a/tests/inout/integrators/langevin1.toml b/tests/inout/integrators/langevin1.toml new file mode 100644 index 0000000..c43e94b --- /dev/null +++ b/tests/inout/integrators/langevin1.toml @@ -0,0 +1,6 @@ +[integrator] +class = "LangevinOverdamped" +timestep = 1 +gamma = 2 +beta = 3 +seed = 123 diff --git a/tests/inout/integrators/langevin2.toml b/tests/inout/integrators/langevin2.toml new file mode 100644 index 0000000..3444555 --- /dev/null +++ b/tests/inout/integrators/langevin2.toml @@ -0,0 +1,6 @@ +[integrator] +class = "LangevinInertia" +timestep = 1 +gamma = 2 +beta = 3 +seed = 123 diff --git a/tests/inout/integrators/velocityverlet.toml b/tests/inout/integrators/velocityverlet.toml new file mode 100644 index 0000000..c1a7f09 --- /dev/null +++ b/tests/inout/integrators/velocityverlet.toml @@ -0,0 +1,3 @@ +[integrator] +class = "VelocityVerlet" +timestep = 1234.5678 diff --git a/tests/inout/integrators/verlet.toml b/tests/inout/integrators/verlet.toml new file mode 100644 index 0000000..0c2ead0 --- /dev/null +++ b/tests/inout/integrators/verlet.toml @@ -0,0 +1,3 @@ +[integrator] +class = "Verlet" +timestep = 1234.5678 diff --git a/tests/inout/md.toml b/tests/inout/md.toml new file mode 100644 index 0000000..bacf46a --- /dev/null +++ b/tests/inout/md.toml @@ -0,0 +1,40 @@ +# Molecular Dynamics Configuration + +# Basic MD settings +[md] +total_steps = 200 # Total simulation steps +temperature = 298.15 # Kelvin + +# Integrator settings +[integrator] +class = "VelocityVerlet" +timestep = 0.005 + +# Initial structure file: +[particles] +file = "config/start.xyz" + +[particles.masses] +Ar = 1.0 + +[particles.types] +Ar = 1 + +# Periodic boundary conditions: +[box] +periodic = true +high = [10.0, 10.0, 10.0] + +# Potential/Forcefield settings +[potential] +type = "Lennard-Jones" # Could also be "AMBER", "CHARMM", etc. + + [potential.parameters] + epsilon = 1.0 # LJ epsilon parameter (for Lennard-Jones potential) + sigma = 3.40 # LJ sigma parameter (for Lennard-Jones potential) + cut = 2.5 # Cut-off for the potential + +# Output settings: +[output] +trajectory_file = "traj.xyz" +log_file = "md.log" diff --git a/tests/inout/nesting.toml b/tests/inout/nesting.toml new file mode 100644 index 0000000..96e96e1 --- /dev/null +++ b/tests/inout/nesting.toml @@ -0,0 +1,16 @@ +[box] +periodic = [true, true, true] +low = [0.0, 0.0, 0.0] +high = [1.0, 1.0, 1.0] + +[system] +file = "system" + +[particles] +file = "particles" + +[system.particles] +file = "system.particles" + +[x.y.z.w] +file = "x.y.z.w" diff --git a/tests/inout/particles/masses.txt b/tests/inout/particles/masses.txt new file mode 100644 index 0000000..2d4a60e --- /dev/null +++ b/tests/inout/particles/masses.txt @@ -0,0 +1,2 @@ +123. +456. diff --git a/tests/inout/particles/particles1.toml b/tests/inout/particles/particles1.toml new file mode 100644 index 0000000..ccea99d --- /dev/null +++ b/tests/inout/particles/particles1.toml @@ -0,0 +1,15 @@ +[box] +periodic = [true, true, true] +low = [0.0, 0.0, 0.0] +high = [1.0, 1.0, 1.0] + +[particles] +file = "../config/start.xyz" + +[particles.masses] +A = 2.0 +B = 10.0 + +[particles.types] +A = 9 +B = 55 diff --git a/tests/inout/particles/particles2.toml b/tests/inout/particles/particles2.toml new file mode 100644 index 0000000..bfbd7e1 --- /dev/null +++ b/tests/inout/particles/particles2.toml @@ -0,0 +1,9 @@ +[box] +periodic = [true, true, true] +low = [0.0, 0.0, 0.0] +high = [1.0, 1.0, 1.0] + +[particles] +file = "../config/start.xyz" +mass_list = [10.0, 40.0] +type_list = [99, 101] diff --git a/tests/inout/particles/particles3.toml b/tests/inout/particles/particles3.toml new file mode 100644 index 0000000..97a0a9b --- /dev/null +++ b/tests/inout/particles/particles3.toml @@ -0,0 +1,9 @@ +[box] +periodic = [true, true, true] +low = [0.0, 0.0, 0.0] +high = [1.0, 1.0, 1.0] + +[particles] +file = "../config/start.xyz" +mass_file = "masses.txt" +type_file = "types.txt" diff --git a/tests/inout/particles/types.txt b/tests/inout/particles/types.txt new file mode 100644 index 0000000..d36fe97 --- /dev/null +++ b/tests/inout/particles/types.txt @@ -0,0 +1,2 @@ +33. +77. diff --git a/tests/inout/start.xyz b/tests/inout/start.xyz new file mode 100644 index 0000000..cc5967a --- /dev/null +++ b/tests/inout/start.xyz @@ -0,0 +1,12 @@ +2 +# Step 2 box: 1.0 2.0 3.0 +A 1.0 1.0 1.0 1.0 1.0 1.0 +B 1.0 1.0 1.0 1.0 1.0 1.0 +2 +# Step 3 box: 1.0 2.0 3.0 +A 4.0 4.0 4.0 4.0 4.0 4.0 +B 4.0 4.0 4.0 4.0 4.0 4.0 +2 +# Step 4 box: 1.0 2.0 3.0 +A 8.0 8.0 8.0 8.0 8.0 8.0 +B 8.0 8.0 8.0 8.0 8.0 8.0 diff --git a/tests/inout/start2.xyz b/tests/inout/start2.xyz new file mode 100644 index 0000000..b939764 --- /dev/null +++ b/tests/inout/start2.xyz @@ -0,0 +1,4 @@ +2 +# Step 4 box: 1.0 2.0 3.0 +A 8.0 8.0 8.0 8.0 8.0 8.0 +B 8.0 8.0 8.0 8.0 8.0 8.0 diff --git a/tests/inout/system.toml b/tests/inout/system.toml new file mode 100644 index 0000000..5364ede --- /dev/null +++ b/tests/inout/system.toml @@ -0,0 +1,7 @@ +[box] +periodic = [true, true, true] +low = [0.0, 0.0, 0.0] +high = [1.0, 1.0, 1.0] + +[particles] +file = "config/start.xyz" diff --git a/tests/inout/system_missing_xyz.toml b/tests/inout/system_missing_xyz.toml new file mode 100644 index 0000000..ea6be26 --- /dev/null +++ b/tests/inout/system_missing_xyz.toml @@ -0,0 +1,7 @@ +[box] +periodic = [true, true, true] +low = [0.0, 0.0, 0.0] +high = [1.0, 1.0, 1.0] + +[particles] +file = "config/start_not_present.xyz" diff --git a/tests/inout/system_units1.toml b/tests/inout/system_units1.toml new file mode 100644 index 0000000..4b551d8 --- /dev/null +++ b/tests/inout/system_units1.toml @@ -0,0 +1,10 @@ +[system] +units = "real" + +[box] +periodic = [true, true, true] +low = [0.0, 0.0, 0.0] +high = [1.0, 1.0, 1.0] + +[particles] +file = "config/start.xyz" diff --git a/tests/inout/system_units2.toml b/tests/inout/system_units2.toml new file mode 100644 index 0000000..d23a1d0 --- /dev/null +++ b/tests/inout/system_units2.toml @@ -0,0 +1,10 @@ +[system] +units = "units-that-does-not-exist" + +[box] +periodic = [true, true, true] +low = [0.0, 0.0, 0.0] +high = [1.0, 1.0, 1.0] + +[particles] +file = "config/start.xyz" diff --git a/tests/inout/test_common.py b/tests/inout/test_common.py new file mode 100644 index 0000000..50f140f --- /dev/null +++ b/tests/inout/test_common.py @@ -0,0 +1,149 @@ +"""Test the common methods for classes.""" +import inspect +import logging + +import pytest + +from turtlemd.inout.common import ( + get_args_and_kwargs, + get_parameter_type, + initiate_instance, + inspect_function, +) + +# Define some functions that we can inspect to get +# the types of arguments. + + +def foo1(arg1, arg2, *arg3, arg4=10, arg5=101, **arg6): + pass + + +def foo2(a, b=10, *args, c, **kwargs): + pass + + +def foo3(): + pass + + +def test_get_parameter_type(): + """Test that we can get the type of parameters.""" + params = inspect.signature(foo1).parameters.values() + param_types = [get_parameter_type(i) for i in params] + correct = ["args", "args", "varargs", "kwargs", "kwargs", "keywords"] + assert param_types == correct + + +def test_inspect_function(): + """Test that we can inspect functions to get parameters.""" + params = inspect_function(foo2) + assert params == { + "args": ["a"], + "kwargs": ["b", "c"], + "varargs": ["args"], + "keywords": ["kwargs"], + } + params = inspect_function(foo3) + assert params == {} + + +class ClassNoArg: + def __init__(self): + self.a = 123 + + +class ClassArg: + def __init__(self, a, b): + self.a = a + self.b = b + + +class ClassArgKwarg: + def __init__(self, a, b=5): + self.a = a + self.b = b + + +class ClassKwarg: + def __init__(self, a=10, b=5, c=101): + self.a = a + self.b = b + self.c = c + + +class ClassVarArgs: + def __init__(self, a, *args): + self.a = a + self.args = args + + +class ClassKeyword: + def __init__(self, a, **kwargs): + self.a = a + self.kwargs = kwargs + + +def test_get_args_and_kwargs(caplog): + """Test that we can get the args and kwargs of constructors.""" + # No arguments: + settings = {} + args, kwargs = get_args_and_kwargs(ClassNoArg, settings) + assert not args + assert not kwargs + # Only positional arguments: + settings = {"a": 1, "b": 2} + args, kwargs = get_args_and_kwargs(ClassArg, settings) + assert args == [1, 2] + assert not kwargs + # Only keyword: + settings = {"a": 1, "c": 2} + args, kwargs = get_args_and_kwargs(ClassKwarg, settings) + assert not args + assert kwargs == {"a": 1, "c": 2} + # Positional and keyword: + settings = {"a": 1, "b": 2} + args, kwargs = get_args_and_kwargs(ClassArgKwarg, settings) + assert args == [1] + assert kwargs == {"b": 2} + # Varargs: + with pytest.raises(ValueError): + get_args_and_kwargs(ClassVarArgs, settings) + # Keywords: + with pytest.raises(ValueError): + get_args_and_kwargs(ClassKeyword, settings) + # Missing arguments: + with pytest.raises(ValueError): + settings = {"a": 10} + get_args_and_kwargs(ClassArg, settings) + # Extra arguments: + settings = {"a": 10, "b": 10, "__extra__": 101, "__stop__": 102} + with caplog.at_level(logging.WARNING): + get_args_and_kwargs(ClassArg, settings) + assert "Ignored extra argument" in caplog.text + assert "__extra__" in caplog.text + assert "__stop__" in caplog.text + + +def test_initiate_instance(): + """Test that we can initiate classes.""" + settings = {} + instance = initiate_instance(ClassNoArg, settings) + assert isinstance(instance, ClassNoArg) + assert instance.a == 123 + settings = {"a": -1, "b": -2} + instance = initiate_instance(ClassArg, settings) + assert isinstance(instance, ClassArg) + assert instance.a == -1 + assert instance.b == -2 + settings = {"a": -1} + instance = initiate_instance(ClassArgKwarg, settings) + assert isinstance(instance, ClassArgKwarg) + assert instance.a == -1 + assert instance.b == 5 + settings = {"a": -1, "b": -5} + instance = initiate_instance(ClassArgKwarg, settings) + assert instance.b == -5 + settings = {"b": -5} + with pytest.raises(ValueError): + initiate_instance(ClassArgKwarg, settings) diff --git a/tests/inout/test_settings.py b/tests/inout/test_settings.py new file mode 100644 index 0000000..4e9ca43 --- /dev/null +++ b/tests/inout/test_settings.py @@ -0,0 +1,165 @@ +"""Test that we can read and interpret the input file.""" +import logging +import pathlib +import shutil + +import numpy as np +import pytest + +from turtlemd.inout.settings import ( + create_box_from_settings, + create_integrator_from_settings, + create_particles_from_settings, + create_system_from_settings, + read_settings_file, + search_for_setting, +) +from turtlemd.integrators import ( + LangevinInertia, + LangevinOverdamped, + VelocityVerlet, + Verlet, +) + +HERE = pathlib.Path(__file__).resolve().parent + + +def test_search_for_setting(): + """Test that we can search for settings.""" + settings = read_settings_file(HERE / "nesting.toml") + assert settings["system"]["particles"]["file"] == "system.particles" + assert settings["system"]["file"] == "system" + assert settings["x"]["y"]["z"]["w"]["file"] == "x.y.z.w" + match = search_for_setting(settings, "file") + for item in match: + item["file"] = "updated" + assert settings["system"]["particles"]["file"] == "updated" + assert settings["system"]["file"] == "updated" + assert settings["x"]["y"]["z"]["w"]["file"] == "updated" + + +def test_create_box(): + """Test that we can create a box from settings.""" + settings = read_settings_file(HERE / "box1.toml") + box = create_box_from_settings(settings) + assert pytest.approx(box.low) == np.array([10.0, 11.0, 12.0]) + assert pytest.approx(box.high) == np.array([23.0, 24.0, 25.0]) + assert len(box.periodic) == 3 + assert all(box.periodic) + # Read a similar file, but periodic is given as only true + # (that is, not a list). + settings = read_settings_file(HERE / "box2.toml") + box = create_box_from_settings(settings) + assert len(box.periodic) == 3 + assert all(box.periodic) + # Read a triclinic box: + settings = read_settings_file(HERE / "box3.toml") + boxt = create_box_from_settings(settings) + assert pytest.approx(boxt.alpha) == 75 + assert pytest.approx(boxt.beta) == 60 + assert pytest.approx(boxt.gamma) == 90 + + +def test_create_integrator(caplog: pytest.LogCaptureFixture): + """Test the creation of integrators from settings.""" + settings = read_settings_file(HERE / "integrators" / "verlet.toml") + integ = create_integrator_from_settings(settings) + assert integ is not None + assert integ.timestep == 1234.5678 + assert isinstance(integ, Verlet) + + settings = read_settings_file(HERE / "integrators" / "verlet.toml") + settings["integrator"].pop("class") + with pytest.raises(ValueError): + with caplog.at_level(logging.ERROR): + create_integrator_from_settings(settings) + assert 'No "class" given for integrator' in caplog.text + + settings = read_settings_file(HERE / "integrators" / "integrator.toml") + with pytest.raises(ValueError): + with caplog.at_level(logging.ERROR): + create_integrator_from_settings(settings) + assert "Could not create unknown class" in caplog.text + + # Test that we can create all integrators: + classes = (Verlet, VelocityVerlet, LangevinOverdamped, LangevinInertia) + files = ( + "verlet.toml", + "velocityverlet.toml", + "langevin1.toml", + "langevin2.toml", + ) + + for klass, filei in zip(classes, files): + settings = read_settings_file(HERE / "integrators" / filei) + integ = create_integrator_from_settings(settings) + assert integ is not None + assert isinstance(integ, klass) + + +def test_create_system(tmp_path: pathlib.PosixPath): + """Test that we can create systems.""" + settings_file = HERE / "system.toml" + settings = read_settings_file(settings_file) + system = create_system_from_settings(settings) + correct = np.array([[1.0, 1.0, 1.0], [2.0, 2.0, 2.0]]) + assert pytest.approx(system.particles.pos) == correct + # Test that we can read relative paths: + settings_file = HERE / "config" / "system_one_down.toml" + settings = read_settings_file(settings_file) + system = create_system_from_settings(settings) + correct2 = np.full((2, 3), 8.0) + assert pytest.approx(system.particles.pos) == correct2 + # Test a absolute path: + new_file = (tmp_path / "start_absolute.xyz").resolve() + shutil.copy(HERE / "config" / "start.xyz", new_file) + settings["particles"] = {"file": new_file} + system = create_system_from_settings(settings) + assert pytest.approx(system.particles.pos) == correct + # Test a missing file: + with pytest.raises(FileNotFoundError): + new_file = (tmp_path / "missing_file.xyz").resolve() + settings["particles"] = {"file": new_file} + create_system_from_settings(settings) + # Test with units: + settings_file = HERE / "system_units1.toml" + settings = read_settings_file(settings_file) + system = create_system_from_settings(settings) + assert system.units.name == "real" + # Test with units that does not exist: + settings_file = HERE / "system_units2.toml" + settings = read_settings_file(settings_file) + with pytest.raises(KeyError): + system = create_system_from_settings(settings) + + +def help_with_particles(settings_file): + """Help with creating the particles.""" + settings = read_settings_file(settings_file) + box = create_box_from_settings(settings) + return create_particles_from_settings(settings, dim=box.dim) + + +def test_create_particles(): + """Test that we can create particles.""" + particles = help_with_particles(HERE / "system.toml") + correct = np.array([[1.0, 1.0, 1.0], [2.0, 2.0, 2.0]]) + assert pytest.approx(particles.pos) == correct + + particles = help_with_particles(HERE / "particles" / "particles1.toml") + assert pytest.approx(particles.mass[0][0]) == 2.0 + assert pytest.approx(particles.mass[1][0]) == 10.0 + assert particles.ptype[0] == 9 + assert particles.ptype[1] == 55 + + particles = help_with_particles(HERE / "particles" / "particles2.toml") + assert pytest.approx(particles.mass[0][0]) == 10.0 + assert pytest.approx(particles.mass[1][0]) == 40.0 + assert particles.ptype[0] == 99 + assert particles.ptype[1] == 101 + + particles = help_with_particles(HERE / "particles" / "particles3.toml") + assert pytest.approx(particles.mass[0][0]) == 123.0 + assert pytest.approx(particles.mass[1][0]) == 456.0 + assert particles.ptype[0] == 33 + assert particles.ptype[1] == 77 diff --git a/tests/inout/test_xyz.py b/tests/inout/test_xyz.py new file mode 100644 index 0000000..432ba2b --- /dev/null +++ b/tests/inout/test_xyz.py @@ -0,0 +1,173 @@ +"""Test that we can read xyz-files.""" +import logging +import pathlib + +import numpy as np +import pytest + +from turtlemd.inout.xyz import ( + configuration_from_xyz_file, + pad_to_nd, + read_xyz_file, + system_to_xyz, +) +from turtlemd.system import Box, Particles, System +from turtlemd.tools import generate_lattice + +HERE = pathlib.Path(__file__).resolve().parent +XYZDIR = HERE / "xyz" + +CORRECT_XYZ = np.array( + [ + [0.0, 0.0, 0.0], + [0.5, 0.5, 0.5], + [0.5, 0.5, 0.0], + [0.5, 0.0, 0.5], + [0.0, 0.5, 0.5], + ] +) + + +def test_read_single_config(): + """Test that we can read a single config from a xyz-file.""" + xyz_file = XYZDIR / "config.xyz" + snapshot = next(read_xyz_file(xyz_file)) + assert snapshot.natoms == 5 + assert snapshot.atoms == ["Ba", "Hf", "O", "O", "O"] + assert ( + snapshot.comment + == "Example from: https://en.wikipedia.org/wiki/XYZ_file_format" + ) + assert pytest.approx(snapshot.xyz) == CORRECT_XYZ + + +def test_read_trajectory(): + """Test that we can read a xyz-trajectory.""" + xyz_file = XYZDIR / "traj.xyz" + for i, snapshot in enumerate(read_xyz_file(xyz_file)): + assert snapshot.natoms == 3 + xyz = np.full_like(snapshot.xyz, 500 - i) + assert pytest.approx(snapshot.xyz) == xyz + assert snapshot.xyz.shape == (3, 6) + assert snapshot.comment.startswith(f"# Step: {500-i} Box:") + assert snapshot.atoms == ["A", "B", "C"] + + +def test_malformed_xyz(caplog): + xyz_file = XYZDIR / "error.xyz" + with pytest.raises(ValueError): + with caplog.at_level(logging.ERROR): + next(read_xyz_file(xyz_file)) + assert "Could not read the number of atoms" in caplog.text + + +def test_configuration_from_xyz(): + """Test that we can read configurations from a given xyz-file.""" + xyz_file = XYZDIR / "config.xyz" + atoms, pos, _, _ = configuration_from_xyz_file(xyz_file, dim=3) + assert atoms == ["Ba", "Hf", "O", "O", "O"] + assert pytest.approx(pos) == CORRECT_XYZ + # Test what happens if we use a 2D system: + _, pos, vel, _ = configuration_from_xyz_file(xyz_file, dim=2) + assert pos.shape == (5, 2) + assert pytest.approx(vel) == np.zeros((5, 2)) + xyz_file = XYZDIR / "config2D.xyz" + _, pos, _, _ = configuration_from_xyz_file(xyz_file, dim=2) + assert pos.shape == (5, 2) + # Test what happens if we have more columns: + xyz_file = XYZDIR / "traj.xyz" + _, pos, vel, _ = configuration_from_xyz_file(xyz_file, dim=3) + assert pytest.approx(pos) == np.full((3, 3), 500) + assert pytest.approx(vel) == np.full((3, 3), 500) + _, pos, vel, force = configuration_from_xyz_file(xyz_file, dim=2) + assert pytest.approx(pos) == np.full((3, 2), 500) + assert pytest.approx(vel) == np.full((3, 2), 500) + assert pytest.approx(force) == np.full((3, 2), 500) + + +def test_pad_to_nd(): + """Test that we can pad 1D and 2D coordinates.""" + vec1 = np.array( + [ + 1, + ] + ) + vec2 = np.array([1, 2]) + vec3 = np.array([1, 2, 3]) + vec4 = np.array([1, 2, 3, 4]) + pad_to_nd(vec1, dim=3) + assert pytest.approx(np.array([1, 0, 0])) == pad_to_nd(vec1, dim=3) + assert pytest.approx(np.array([1, 2, 0])) == pad_to_nd(vec2, dim=3) + assert pytest.approx(vec3) == pad_to_nd(vec3, dim=3) + assert pytest.approx(vec4) == pad_to_nd(vec4, dim=3) + + +def create_test_system(lattice="fcc"): + """Create a test system.""" + if lattice == "fcc": + xyz, size = generate_lattice("fcc", [3, 3, 3], density=0.9) + elif lattice == "sq2": + xyz, size = generate_lattice(lattice, [3, 3], density=0.9) + else: + xyz = [ + 1.0, + ] + size = [None, None] + box = Box(low=size[:, 0], high=size[:, 1]) + particles = Particles(dim=box.dim) + for pos in xyz: + particles.add_particle(pos=pos, mass=1.0, name="Ar", ptype=0) + return System(box=box, particles=particles) + + +def test_system_to_xyz(tmp_path: pathlib.PosixPath): + """Test that we can create a XYZ file from a system.""" + system = create_test_system() + xyz_file = (tmp_path / "system.xyz").resolve() + system_to_xyz(system, xyz_file) + _, pos, _, _ = configuration_from_xyz_file(xyz_file) + assert pytest.approx(pos) == system.particles.pos + positions = [system.particles.pos.copy()] + # Test that we can append to a file to create two frames: + system.particles.pos += 1.234 + positions.append(system.particles.pos.copy()) + system_to_xyz(system, xyz_file, filemode="a", title="Second frame") + for i, snapshot in enumerate(read_xyz_file(xyz_file)): + assert pytest.approx(snapshot.xyz) == positions[i] + if i == 0: + assert snapshot.comment.startswith( + "# TurtleMD system. Box: 4.9324" + ) + elif i == 1: + assert snapshot.comment.strip() == "Second frame" + + +def test_system_to_xyz_2dim(tmp_path: pathlib.PosixPath): + """Test that we can write XYZ-files when the system is 2D or 1D.""" + # For 2D: + system = create_test_system(lattice="sq2") + xyz_file = (tmp_path / "systemsq2.xyz").resolve() + system_to_xyz(system, xyz_file) + _, pos, _, _ = configuration_from_xyz_file(xyz_file) + assert pytest.approx(pos[:, 0]) == system.particles.pos[:, 0] + assert pytest.approx(pos[:, 1]) == system.particles.pos[:, 1] + assert pos.shape == (18, 3) + assert system.particles.pos.shape == (18, 2) + assert pytest.approx(pos[:, 2]) == np.zeros(18) + + # For 1D: + box = Box(periodic=[False]) + particles = Particles(dim=box.dim) + for i in range(3): + particles.add_particle( + pos=[i], + mass=1.0, + name="Ar", + ptype=0, + ) + system_1d = System(box=box, particles=particles) + xyz_file = (tmp_path / "system1D.xyz").resolve() + system_to_xyz(system_1d, xyz_file) + _, pos, _, _ = configuration_from_xyz_file(xyz_file) + pos1 = pos[:, 0] + assert pytest.approx(pos1) == system_1d.particles.pos.flatten() diff --git a/tests/tools/config.xyz b/tests/inout/xyz/config.xyz similarity index 100% rename from tests/tools/config.xyz rename to tests/inout/xyz/config.xyz diff --git a/tests/inout/xyz/config2D.xyz b/tests/inout/xyz/config2D.xyz new file mode 100644 index 0000000..a0cb4e2 --- /dev/null +++ b/tests/inout/xyz/config2D.xyz @@ -0,0 +1,7 @@ +5 +Example from: https://en.wikipedia.org/wiki/XYZ_file_format +Ba 0.0 0.0 +Hf 0.5 0.5 +O 0.5 0.5 +O 0.5 0.0 +O 0.0 0.5 diff --git a/tests/tools/error.xyz b/tests/inout/xyz/error.xyz similarity index 100% rename from tests/tools/error.xyz rename to tests/inout/xyz/error.xyz diff --git a/tests/tools/traj.xyz b/tests/inout/xyz/traj.xyz similarity index 86% rename from tests/tools/traj.xyz rename to tests/inout/xyz/traj.xyz index 69d14ab..e5bdb2b 100644 --- a/tests/tools/traj.xyz +++ b/tests/inout/xyz/traj.xyz @@ -1,20 +1,20 @@ 3 -# Step: 500 Box: 0.0000 1.0000 2.0000 +# Step: 500 Box: 0.0000 1.0000 2.0000 A 500.000000000 500.000000000 500.000000000 500.000000000 500.000000000 500.000000000 B 500.000000000 500.000000000 500.000000000 500.000000000 500.000000000 500.000000000 C 500.000000000 500.000000000 500.000000000 500.000000000 500.000000000 500.000000000 3 -# Step: 499 Box: 0.0000 1.0000 2.0000 +# Step: 499 Box: 0.0000 1.0000 2.0000 A 499.000000000 499.000000000 499.000000000 499.000000000 499.000000000 499.000000000 B 499.000000000 499.000000000 499.000000000 499.000000000 499.000000000 499.000000000 C 499.000000000 499.000000000 499.000000000 499.000000000 499.000000000 499.000000000 3 -# Step: 498 Box: 0.0000 1.0000 2.0000 +# Step: 498 Box: 0.0000 1.0000 2.0000 A 498.000000000 498.000000000 498.000000000 498.000000000 498.000000000 498.000000000 B 498.000000000 498.000000000 498.000000000 498.000000000 498.000000000 498.000000000 C 498.000000000 498.000000000 498.000000000 498.000000000 498.000000000 498.000000000 3 -# Step: 497 Box: 0.0000 1.0000 2.0000 +# Step: 497 Box: 0.0000 1.0000 2.0000 A 497.000000000 497.000000000 497.000000000 497.000000000 497.000000000 497.000000000 B 497.000000000 497.000000000 497.000000000 497.000000000 497.000000000 497.000000000 C 497.000000000 497.000000000 497.000000000 497.000000000 497.000000000 497.000000000 diff --git a/tests/potentials/test_lennardjones.py b/tests/potentials/test_lennardjones.py index 576154b..26ac5ea 100644 --- a/tests/potentials/test_lennardjones.py +++ b/tests/potentials/test_lennardjones.py @@ -109,7 +109,7 @@ def test_mix_sixthpower(): def test_unkown_mixing(): - """Test that we fail for an unknow mixing rule.""" + """Test that we fail for an unknown mixing rule.""" with pytest.raises(ValueError): mix_parameters(1, 2, 3, 4, 5, 6, mixing="?") diff --git a/tests/potentials/test_wellwca.py b/tests/potentials/test_wellwca.py index be18cfd..279e4db 100644 --- a/tests/potentials/test_wellwca.py +++ b/tests/potentials/test_wellwca.py @@ -85,7 +85,7 @@ def test_activate(): def test_potential(): - """Test that we can calulate the potential energy correctly.""" + """Test that we calculate the potential energy correctly.""" system = create_system() pot = DoubleWellPair(types=(0, 0)) pot.set_parameters(PARAMETERS) diff --git a/tests/simulation/test_mdsimulation.py b/tests/simulation/test_mdsimulation.py index 86d55e8..9219d56 100644 --- a/tests/simulation/test_mdsimulation.py +++ b/tests/simulation/test_mdsimulation.py @@ -127,7 +127,7 @@ def test_md_simulation(): traj = np.load(HERE / "md-traj.npy") # Accumulate energies: for i, systemi in enumerate(simulation.run()): - therm = systemi.thermo(boltzmann=1) + therm = systemi.thermo() assert ( pytest.approx(therm["temperature"]) == CORRECT_ENERGIES["temperature"][i] diff --git a/tests/system/test_box.py b/tests/system/test_box.py index 3ab24b8..d23fb0e 100644 --- a/tests/system/test_box.py +++ b/tests/system/test_box.py @@ -112,17 +112,17 @@ def test_volume(): assert vol == float("inf") -def test_print(capfd): +def test_print(capsys): """Test that we can print box information""" box = Box(high=[5]) print(box) - captured = capfd.readouterr() + captured = capsys.readouterr() assert "Hello, this is box" in captured.out assert f"and my matrix is:\n{box.box_matrix}" in captured.out box = TriclinicBox(high=[10.0, 10.0], alpha=None, beta=None, gamma=45.0) print(box) - captured = capfd.readouterr() + captured = capsys.readouterr() assert "Hello, this is triclinic box and my matrix" in captured.out diff --git a/tests/system/test_system.py b/tests/system/test_system.py index 00bc509..3c5384d 100644 --- a/tests/system/test_system.py +++ b/tests/system/test_system.py @@ -8,7 +8,7 @@ from turtlemd.system.system import System -def test_system_setup(): +def test_system_setup(capsys): """Test that we can initiate a system.""" box = Box() particles = Particles() @@ -19,6 +19,10 @@ def test_system_setup(): assert system.potentials[0] is potentials[0] assert len(system.potentials) == 1 + print(system) + captured = capsys.readouterr() + assert "TurtleMD system with" in captured.out + def test_potential(): """Test that we can evaluate the potential via the system.""" diff --git a/tests/test_bin.py b/tests/test_bin.py new file mode 100644 index 0000000..572f7c8 --- /dev/null +++ b/tests/test_bin.py @@ -0,0 +1,46 @@ +import logging +import pathlib + +import pytest + +from turtlemd.bin import main +from turtlemd.version import __version__ + +HERE = pathlib.Path(__file__).resolve().parent + + +def test_main_info(monkeypatch, capsys): + """Test that we can run turtlemd with "-v".""" + + for arg in ("-v", "--version"): + test_args = ["turtlemd", arg] + + monkeypatch.setattr("sys.argv", test_args) + + with pytest.raises(SystemExit): + main() + + captured = capsys.readouterr() + assert __version__ in captured.out + + +def test_main_input(monkeypatch, caplog, capsys): + """Test that we can read an input file.""" + for arg in ("-i", "--input_file"): + test_args = ["turtlemd", arg] + monkeypatch.setattr("sys.argv", test_args) + + with pytest.raises(SystemExit): + main() + + captured = capsys.readouterr() + assert "expected one argument" in captured.err + + tomlfile = HERE / "inout" / "md.toml" + + test_args = ["turtlemd", arg, str(tomlfile)] + monkeypatch.setattr("sys.argv", test_args) + with caplog.at_level(logging.INFO): + main() + assert "Reading settings from file" in caplog.text + assert "Created system" in caplog.text diff --git a/tests/test_integrators.py b/tests/test_integrators.py index d051108..ad71c88 100644 --- a/tests/test_integrators.py +++ b/tests/test_integrators.py @@ -5,6 +5,7 @@ import turtlemd.integrators from turtlemd.integrators import ( + INTEGRATORS, LangevinInertia, LangevinOverdamped, LangevinParameter, @@ -483,11 +484,19 @@ def test_verlet(): def test_langevin_brownian(): """Test the overdamped Langevin integrator.""" system = create_test_system() + # Test that we create a default random generator: + integrator = LangevinOverdamped( + timestep=0.002, + gamma=0.3, + beta=1.0, + ) + assert integrator.rgen is not None + integrator = LangevinOverdamped( timestep=0.002, gamma=0.3, - rgen=FakeRandomGenerator(seed=1), beta=1.0, + rgen=FakeRandomGenerator(seed=1), ) for i in range(51): assert ( @@ -560,3 +569,11 @@ def test_langevin_integration(monkeypatch): assert pytest.approx(pos[0][0]) == TISMOL_POS_LANG[i] assert pytest.approx(vel[0][0]) == TISMOL_VEL_LANG[i] integrator(system) + + +def test_integrator_registry(): + """Test that we can access the integrator registry.""" + assert "verlet" in INTEGRATORS + assert issubclass(INTEGRATORS["verlet"], Verlet) + vverlet = INTEGRATORS.get("velocityverlet", None) + assert issubclass(vverlet, VelocityVerlet) diff --git a/tests/test_units.py b/tests/test_units.py new file mode 100644 index 0000000..881a76a --- /dev/null +++ b/tests/test_units.py @@ -0,0 +1,23 @@ +import pytest +from pint import Quantity +from scipy.constants import Avogadro as AVOGADRO # type: ignore +from scipy.constants import Boltzmann as BOLTZMANN + +from turtlemd.units import UnitSystem + + +def test_units(capsys): + """Test that we can create create a unit system with correct time.""" + unit = UnitSystem( + name="gromacs", + length=Quantity(1.0, "nm"), + mass=Quantity(1.0 / AVOGADRO, "g"), + energy=Quantity(1.0 / AVOGADRO, "kJ"), + boltzmann=(Quantity(BOLTZMANN, "J/K").to("kJ/K") * AVOGADRO).magnitude, + input_time_unit="ps", + ) + assert pytest.approx(unit.time_factor) == 1.0 + assert pytest.approx(unit.boltzmann) == 0.00831446261815324 + print(unit) + captured = capsys.readouterr() + assert "gromacs" in captured.out diff --git a/tests/tools/config2D.xyz b/tests/tools/config2D.xyz deleted file mode 100644 index d552047..0000000 --- a/tests/tools/config2D.xyz +++ /dev/null @@ -1,7 +0,0 @@ -5 -Example from: https://en.wikipedia.org/wiki/XYZ_file_format -Ba 0.0 0.0 -Hf 0.5 0.5 -O 0.5 0.5 -O 0.5 0.0 -O 0.0 0.5 diff --git a/tests/tools/test_tools.py b/tests/tools/test_tools.py index 85d0379..06282a1 100644 --- a/tests/tools/test_tools.py +++ b/tests/tools/test_tools.py @@ -86,7 +86,7 @@ def test_lattice_generation(): def test_lattice_special(): """Test all the 'special' cases for lattice generation.""" - # Test that we fail for a unknow lattice. + # Test that we fail for a unknown lattice. with pytest.raises(ValueError): generate_lattice(lattice="this is not a lattice!") # Test that we set lcon to just one and repeat to ones, if diff --git a/tests/tools/test_xyz.py b/tests/tools/test_xyz.py deleted file mode 100644 index b4e4514..0000000 --- a/tests/tools/test_xyz.py +++ /dev/null @@ -1,89 +0,0 @@ -"""Test that we can read xyz-files.""" -import logging -import pathlib - -import numpy as np -import pytest - -from turtlemd.tools.xyz import particles_from_xyz_file, read_xyz_file - -HERE = pathlib.Path(__file__).resolve().parent - -CORRECT_XYZ = np.array( - [ - [0.0, 0.0, 0.0], - [0.5, 0.5, 0.5], - [0.5, 0.5, 0.0], - [0.5, 0.0, 0.5], - [0.0, 0.5, 0.5], - ] -) - - -def test_read_single_config(): - """Test that we can read a single config from a xyz-file.""" - xyz_file = HERE / "config.xyz" - snapshot = next(read_xyz_file(xyz_file)) - assert snapshot.natoms == 5 - assert snapshot.atoms == ["Ba", "Hf", "O", "O", "O"] - assert ( - snapshot.comment - == "Example from: https://en.wikipedia.org/wiki/XYZ_file_format" - ) - assert pytest.approx(snapshot.xyz) == CORRECT_XYZ - - -def test_read_trajectory(): - """Test that we can read a xyz-trajectory.""" - xyz_file = HERE / "traj.xyz" - for i, snapshot in enumerate(read_xyz_file(xyz_file)): - assert snapshot.natoms == 3 - xyz = np.full_like(snapshot.xyz, 500 - i) - assert pytest.approx(snapshot.xyz) == xyz - assert snapshot.xyz.shape == (3, 6) - assert snapshot.comment.startswith(f"# Step: {500-i} Box:") - assert snapshot.atoms == ["A", "B", "C"] - - -def test_malformed_xyz(caplog): - xyz_file = HERE / "error.xyz" - with pytest.raises(ValueError): - with caplog.at_level(logging.ERROR): - next(read_xyz_file(xyz_file)) - assert "Could not read the number of atoms" in caplog.text - - -def test_particles_from_xyz(): - """Test that we can create particles from a given xyz-file.""" - xyz_file = HERE / "config.xyz" - # Set up some masses: - masses = { - "O": 16.0, - "Hf": 178.49, - "Ba": 137.33, - } - particles = particles_from_xyz_file(xyz_file, dim=3, masses=masses) - assert list(particles.name) == ["Ba", "Hf", "O", "O", "O"] - assert pytest.approx(particles.pos) == CORRECT_XYZ - mass_table = np.array([137.33, 178.49, 16.0, 16.0, 16.0]).reshape(5, 1) - assert pytest.approx(particles.mass) == mass_table - assert particles.ptype[0] != particles.ptype[1] - assert particles.ptype[0] != particles.ptype[2] - assert particles.ptype[2] == particles.ptype[3] - assert particles.ptype[2] == particles.ptype[4] - # Test what happens if we use a 2D system: - particles = particles_from_xyz_file(xyz_file, dim=2) - assert particles.pos.shape == (5, 2) - assert pytest.approx(particles.vel) == np.zeros((5, 2)) - xyz_file = HERE / "config2D.xyz" - particles = particles_from_xyz_file(xyz_file, dim=2) - assert particles.pos.shape == (5, 2) - # Test what happens if we have more columns: - xyz_file = HERE / "traj.xyz" - particles = particles_from_xyz_file(xyz_file, dim=3) - assert pytest.approx(particles.pos) == np.full((3, 3), 500) - assert pytest.approx(particles.vel) == np.full((3, 3), 500) - particles = particles_from_xyz_file(xyz_file, dim=2) - assert pytest.approx(particles.pos) == np.full((3, 2), 500) - assert pytest.approx(particles.vel) == np.full((3, 2), 500) - assert pytest.approx(particles.force) == np.full((3, 2), 500) diff --git a/turtlemd/bin.py b/turtlemd/bin.py new file mode 100644 index 0000000..7c7a6bb --- /dev/null +++ b/turtlemd/bin.py @@ -0,0 +1,53 @@ +"""Command line execution of TurtleMD.""" +import argparse +import logging + +from turtlemd.inout.settings import ( + create_integrator_from_settings, + create_system_from_settings, + read_settings_file, +) +from turtlemd.simulation import MDSimulation +from turtlemd.version import __version__ + +LOGGER = logging.getLogger(__name__) +LOGGER.addHandler(logging.NullHandler()) + + +def main(): + """Set up the parser and read the input file.""" + parser = argparse.ArgumentParser( + # prog="TurtleMD", + description="Run TurtleMD from input files.", + ) + parser.add_argument( + "-i", + "--input_file", + help="path to the input TOML file.", + type=str, + required=True, + ) + parser.add_argument( + "-v", + "--version", + action="version", + version=f"%(prog)s {__version__}", + help="show program version and exit.", + ) + args = parser.parse_args() + + LOGGER.info(f"Reading settings from file {args.input_file}") + settings = read_settings_file(args.input_file) + system = create_system_from_settings(settings) + LOGGER.info(f"Created system: {system}") + integrator = create_integrator_from_settings(settings) + LOGGER.info(f"Created integrator {integrator}") + + simulation = MDSimulation( + system=system, + integrator=integrator, + steps=settings["md"]["steps"], + start=0, + ) + + return settings, simulation diff --git a/turtlemd/inout/__init__.py b/turtlemd/inout/__init__.py new file mode 100644 index 0000000..cbdb58f --- /dev/null +++ b/turtlemd/inout/__init__.py @@ -0,0 +1,4 @@ +"""Make inout into a package.""" +from .xyz import configuration_from_xyz_file, read_xyz_file + +__all__ = ["configuration_from_xyz_file", "read_xyz_file"] diff --git a/turtlemd/inout/common.py b/turtlemd/inout/common.py new file mode 100644 index 0000000..c6da684 --- /dev/null +++ b/turtlemd/inout/common.py @@ -0,0 +1,195 @@ +"""Methods for creating objects from settings.""" +from __future__ import annotations + +import inspect +import logging +from collections import defaultdict +from inspect import Parameter +from typing import TYPE_CHECKING, Any + +if TYPE_CHECKING: # pragma: no cover + from collections.abc import Callable + + +LOGGER = logging.getLogger(__name__) +LOGGER.addHandler(logging.NullHandler()) + + +def generic_factory( + settings: dict[str, Any], + registry: dict[str, type], + name: str = "generic", +) -> Any: + """Create an instance of a class based on the given settings. + + Args: + settings: Specifies arguments and keyword arguments for object + initiation. Must contain a "class" key to specify which + class should be instantiated. + registry: A mapping from class names to classes. + name: Labels the type of classes we create. Mainly + used for error messages. Defaults to "generic". + + Returns: + An instance of the specified class if successful. Returns + None otherwise. + + Raises: + ValueError: If the "class" key is missing in settings or if the + specified class is not found in the registry. + """ + klass_name = settings.get("class", "").lower() + + if not klass_name: + msg = f'No "class" given for {name}! Could not create object!' + LOGGER.error(msg) + raise ValueError(msg) + + klass = registry.get(klass_name, None) + + if klass is None: + msg = f'Could not create unknown {name} class "{settings["class"]}"!' + LOGGER.error(msg) + raise ValueError(msg) + settings.pop("class", None) + return initiate_instance(klass, settings) + + +def get_parameter_type(arg: Parameter) -> str: + """Determine the type of a method's parameter. + + This method categorizes the given parameter based on its kind + (e.g., positional, keyword-only, etc.). + + Args: + arg: The parameter whose type will be determined. + + Returns: + The type of the parameter. Can be one of the following: + - "args": Positional parameters or positional-or-keyword + parameters without default values. + - "kwargs": Keyword-only parameters or positional-or-keyword + parameters with default values. + - "varargs": For variable positional parameters (i.e., *args). + - "keywords": For variable keyword parameters (i.e., **kwargs). + """ + if arg.kind == Parameter.POSITIONAL_OR_KEYWORD: + return "args" if arg.default is Parameter.empty else "kwargs" + kind_mapping = { + Parameter.POSITIONAL_ONLY: "args", + Parameter.VAR_POSITIONAL: "varargs", + Parameter.KEYWORD_ONLY: "kwargs", + Parameter.VAR_KEYWORD: "keywords", + } + return kind_mapping[arg.kind] + + +def get_args_and_kwargs( + klass: Any, settings: dict[str, Any] +) -> tuple[list[Any], dict[str, Any]]: + """Extract arguments for a class's constructor from given settings. + + This method identifies the required arguments for the constructor + of `klass` and retrieves their values from the `settings` dictionary. + + Args: + klass: The class for which the constructor arguments need to be + extracted. + settings: A dictionary containing potential positional and keyword + arguments for `klass.__init__()`. + + Returns: + A tuple containing: + - Positional arguments extracted from `settings`. + - Keyword arguments extracted from `settings`. + + Raises: + ValueError: If a required positional argument for `klass` is not + present in `settings` or if the method needs *args or **kwargs. + """ + info = inspect_function(klass.__init__) + used_args, positional_args, keyword_args = set(), [], {} + missing_args = [] + + # We disallow the usage of *args when creating instances from settings: + if info.get("varargs"): + raise ValueError( + "*args are not supported when initiating classes from settings!" + ) + # We disallow the usage of **kwargs when creating instances from settings: + if info.get("keywords"): + raise ValueError( + "*kwargs are not supported when initiating classes from settings!" + ) + + for arg in info.get("args", []): + if arg == "self": + continue + if arg in settings: + positional_args.append(settings[arg]) + used_args.add(arg) + else: + missing_args.append(arg) + + if missing_args: + raise ValueError( + f"Required arguments {', '.join(missing_args)} for " + f"{klass} were not found!" + ) + + for kwarg in info.get("kwargs", []): + if kwarg != "self" and kwarg in settings: + keyword_args[kwarg] = settings[kwarg] + used_args.add(kwarg) + + for key in settings: + if key not in used_args: + LOGGER.warning( + 'Ignored extra argument "%s" when initiating %s', key, klass + ) + return positional_args, keyword_args + + +def inspect_function(method: Callable) -> dict[str, list[Any]]: + """Return arguments/kwargs of a given function. + + This function categorizes the parameters of the given function or method + based on their types (e.g., positional, keyword-only, etc.). + + Args: + method: The function or method to inspect. + + Returns: + A dictionary categorizing parameters: + - "args": Positional or positional-or-keyword parameters + without default values. + - "kwargs": Keyword-only or positional-or-keyword parameters + with default values. + - "varargs": Variable positional parameters (*args). + - "keywords": Variable keyword parameters (**kwargs). + """ + arguments = inspect.signature(method) + parameters_by_type = defaultdict(list) + for name, param in arguments.parameters.items(): + parameters_by_type[get_parameter_type(param)].append(name) + return parameters_by_type + + +def initiate_instance(klass: type[Any], settings: dict[str, Any]) -> Any: + """Initialize a class with arguments from settings. + + Args: + klass: The class to initialize. + settings): A dictionary containing positional and keyword arguments + (if any) for `klass.__init__()`. + + Returns: + An instance of the given class, initialized with the arguments + extracted from `settings`. + """ + args, kwargs = get_args_and_kwargs(klass, settings) + msg = 'Initiated "%s" from "%s"' + msg += "\n-arguments: %s" + msg += "\n-kwargs: %s" + LOGGER.debug(msg, klass.__name__, klass.__module__, args, kwargs) + return klass(*args, **kwargs) diff --git a/turtlemd/inout/default.toml b/turtlemd/inout/default.toml new file mode 100644 index 0000000..03c85d0 --- /dev/null +++ b/turtlemd/inout/default.toml @@ -0,0 +1,25 @@ +[run] +directory = "." + +[md] +steps = 1000 + +[integrator] +class = "VelocityVerlet" +timestep = 0.005 + +[system] +units = "reduced" + +[box] +periodic = [true, true, true] +low = [0.0, 0.0, 0.0] +high = [1.0, 1.0, 1.0] + +[particles] +file = "initial.xyz" + +[particles.velocity] +temperature = 1.0 +momentum = true +seed = 1 diff --git a/turtlemd/inout/settings.py b/turtlemd/inout/settings.py new file mode 100644 index 0000000..fa1b1f0 --- /dev/null +++ b/turtlemd/inout/settings.py @@ -0,0 +1,211 @@ +"""Define the input file.""" +from __future__ import annotations + +import logging +import pathlib +from typing import TYPE_CHECKING, Any + +import numpy as np +import toml + +from turtlemd.inout.common import generic_factory +from turtlemd.inout.xyz import configuration_from_xyz_file +from turtlemd.integrators import INTEGRATORS, MDIntegrator +from turtlemd.random import create_random_generator +from turtlemd.system import Box, Particles, System +from turtlemd.system.box import TriclinicBox +from turtlemd.system.particles import generate_maxwell_velocities + +if TYPE_CHECKING: # pragma: no cover + pass + +DEFAULT = pathlib.Path(__file__).resolve().parent / "default.toml" + +LOGGER = logging.getLogger(__name__) +LOGGER.addHandler(logging.NullHandler()) + + +def read_settings_file(settings_file: pathlib.Path | str) -> dict[str, Any]: + """Read settings from the given file.""" + settings_file = pathlib.Path(settings_file) + default = toml.load(DEFAULT) + settings = toml.load(settings_file) + settings = deep_update(default, settings) + settings["run"]["directory"] = settings_file.resolve().parent + return settings + + +def deep_update( + target: dict[str, Any], source: dict[str, Any] +) -> dict[str, Any]: + """Update the target dictionary with settings from the source.""" + for key, val in source.items(): + if isinstance(val, dict): + target[key] = deep_update(target.get(key, {}), val) + else: + target[key] = val + return target + + +def search_for_setting( + settings: dict[str, Any], target: str +) -> list[dict[str, Any]]: + """Search for dictionary items by a given key.""" + stack = [settings] + found = [] + while stack: + current = stack.pop() + for key, val in current.items(): + if isinstance(val, dict): + stack.append(val) + elif key == target: + found.append(current) + return found + + +def create_box_from_settings(settings: dict[str, Any]) -> Box | TriclinicBox: + """Create a simulation box from settings.""" + low = settings["box"]["low"] + periodic = settings["box"]["periodic"] + + try: + len(periodic) + except TypeError: + periodic = [periodic] * len(low) + settings["box"]["periodic"] = periodic + + if any(angle in settings["box"] for angle in ("alpha", "beta", "gamma")): + return TriclinicBox(**settings["box"]) + return Box(**settings["box"]) + + +def create_integrator_from_settings(settings: dict[str, Any]) -> MDIntegrator: + """Create an integrator from settings.""" + return generic_factory( + settings["integrator"], INTEGRATORS, name="integrator" + ) + + +def look_for_file(settings: dict[str, Any], filename: str) -> pathlib.Path: + """Find the file from a give input string.""" + file_path = pathlib.Path(filename) + + if file_path.is_absolute() and file_path.is_file(): + return file_path + else: + base_dir = settings["run"]["directory"] + file_path = (base_dir / file_path).resolve() + if not file_path.is_file(): + msg = "File %s not found." + LOGGER.critical(msg, file_path) + raise FileNotFoundError(msg, file_path) + return file_path + + +def get_particle_data_from_settings( + settings: dict[str, Any], + dict_key: str, + list_key: str, + file_key: str, + dtype: type = float, +) -> tuple[dict[str, int], np.ndarray | None]: + """Get masses or types from the settings.""" + data_dict = settings["particles"].get(dict_key, {}) + + data_list = settings["particles"].get(list_key) + if data_list: + data_list = np.array(data_list, dtype=dtype) + + if file_key in settings["particles"]: + filename = look_for_file(settings, settings["particles"][file_key]) + data_list = np.loadtxt(filename, dtype=dtype) + + return data_dict, data_list + + +def create_particles_from_settings( + settings: dict[str, Any], dim: int = 3 +) -> Particles: + """Create particles from settings. + + Args: + settings: The settings to create particles from. + dim: The dimensionality of the system. + """ + particles = Particles(dim=dim) + + xyz_file = look_for_file(settings, settings["particles"]["file"]) + msg = "Loading initial coordinates from file: %s" + LOGGER.info(msg, xyz_file) + atoms, pos, vel, force = configuration_from_xyz_file(xyz_file, dim=dim) + + mass_dict, mass_list = get_particle_data_from_settings( + settings, + dict_key="masses", + list_key="mass_list", + file_key="mass_file", + dtype=float, + ) + type_dict, type_list = get_particle_data_from_settings( + settings, + dict_key="types", + list_key="type_list", + file_key="type_file", + dtype=float, + ) + + if not type_dict: + type_dict = {atom: idx for idx, atom in enumerate(set(atoms))} + + for i, (atomi, posi, veli, forcei) in enumerate( + zip(atoms, pos, vel, force) + ): + if mass_list is not None: + massi = mass_list[i] + else: + massi = mass_dict.get(atomi, 1.0) + if type_list is not None: + typei = type_list[i] + else: + typei = type_dict.get(atomi, -1) + particles.add_particle( + pos=posi, + vel=veli, + force=forcei, + mass=massi, + name=atomi, + ptype=typei, + ) + return particles + + +def create_velocities(settings: dict[str, Any], system: System) -> None: + """Create velocities for the particles in a system.""" + vel_settings = settings.get("particles", {}).get("velocity") + if vel_settings is None: + return + + rgen = create_random_generator(seed=vel_settings.get("seed")) + + generate_maxwell_velocities( + system.particles, + rgen, + temperature=vel_settings["temperature"], + boltzmann=system.units.boltzmann, + dof=system.dof(), + momentum=vel_settings["momentum"], + ) + + +def create_system_from_settings(settings: dict[str, Any]) -> System: + """Create a system from the given settings.""" + # Set up the box: + box = create_box_from_settings(settings) + particles = create_particles_from_settings(settings, dim=box.dim) + system = System( + box=box, + particles=particles, + units=settings["system"]["units"], + ) + create_velocities(settings, system) + return system diff --git a/turtlemd/tools/xyz.py b/turtlemd/inout/xyz.py similarity index 51% rename from turtlemd/tools/xyz.py rename to turtlemd/inout/xyz.py index 2fd68df..bc88617 100644 --- a/turtlemd/tools/xyz.py +++ b/turtlemd/inout/xyz.py @@ -1,17 +1,28 @@ """Methods to read particles from xyz-files.""" + +from __future__ import annotations + import logging import pathlib from collections.abc import Iterator from dataclasses import dataclass, field +from typing import TYPE_CHECKING import numpy as np -from turtlemd.system.particles import Particles +if TYPE_CHECKING: # pragma: no cover + from turtlemd.system import System LOGGER = logging.getLogger(__name__) LOGGER.addHandler(logging.NullHandler()) +def pad_to_nd(x: np.ndarray, dim: int = 3) -> np.ndarray: + """Pad 1D and 2D vectors to 3D.""" + length = max(0, dim - len(x)) + return np.pad(x, (0, length), mode="constant") + + @dataclass class Snapshot: """Store coordinates and atoms for a snapshot.""" @@ -24,11 +35,10 @@ class Snapshot: def read_xyz_file(filename: str | pathlib.Path) -> Iterator[Snapshot]: - """Read configurations from a xyz-file""" - + """Read configurations from a xyz-file.""" lines_to_read = 0 snapshot = None - with open(filename) as fileh: + with open(filename, encoding="utf-8") as fileh: for lines in fileh: if lines_to_read == 0: # This is a new frame @@ -61,23 +71,19 @@ def read_xyz_file(filename: str | pathlib.Path) -> Iterator[Snapshot]: yield snapshot -def particles_from_xyz_file( +def configuration_from_xyz_file( filename: str | pathlib.Path, dim: int = 3, - masses: dict[str, float] | None = None, -) -> Particles: - """Create particles from a given xyz-file. +) -> tuple[list[str], np.ndarray, np.ndarray, np.ndarray]: + """Get atoms, positions, velocities, and forces from a xyz-file. Args: - dim: The number of dimensions to consider. - masses: dict[str, float] + filename: The file to read the configuration from. + dim: The number of dimensions. """ - if masses is None: - masses = {} snapshot = next(read_xyz_file(filename)) - particles = Particles(dim=dim) - # We will just assign particle types from the atom name: - ptypes = {atom: idx for idx, atom in enumerate(set(snapshot.atoms))} + atoms, positions, velocities, forces = [], [], [], [] + for atom, xyz in zip(snapshot.atoms, snapshot.xyz): pos = xyz[:dim] vel = None @@ -85,16 +91,43 @@ def particles_from_xyz_file( if len(xyz) > dim: vel = xyz[dim : dim * 2] if len(vel) != len(pos): - vel = None + vel = np.zeros_like(pos) force = xyz[dim * 2 : dim * 3] if len(force) != len(pos): - force = None - particles.add_particle( - pos=pos, - vel=vel, - force=force, - mass=masses.get(atom, 1.0), - name=atom, - ptype=ptypes.get(atom, -1), - ) - return particles + force = np.zeros_like(pos) + atoms.append(atom) + positions.append(pos) + velocities.append(vel) + forces.append(force) + return atoms, np.array(positions), np.array(velocities), np.array(forces) + + +def system_to_xyz( + system: System, + filename: str | pathlib.Path, + filemode: str = "w", + title: str | None = None, +) -> None: + """Write the system configuration to a xyz-file. + + Args: + system: The system to get the particles from. + filename: The path to the file to write. + filemode: If "w" this method will overwrite if the `filename` + already exists. Otherwise it will append to it. + title: A string to use as the title for the frame. + """ + if filemode != "w": + filemode = "a" + if title is None: + box = " ".join([f"{i}" for i in system.box.box_matrix.flatten()]) + txt_title = f"# TurtleMD system. Box: {box}\n" + else: + txt_title = f"{title}\n" + with open(filename, filemode, encoding="utf-8") as output_xyz: + output_xyz.write(f"{system.particles.npart}\n") + output_xyz.write(txt_title) + for part in system.particles: + name = f"{part['name']:5s}" + pos = " ".join([f"{i:15.9f}" for i in pad_to_nd(part["pos"])]) + output_xyz.write(f"{name} {pos}\n") diff --git a/turtlemd/integrators.py b/turtlemd/integrators.py index b72f4ca..219af1a 100644 --- a/turtlemd/integrators.py +++ b/turtlemd/integrators.py @@ -1,17 +1,36 @@ """Definition of time integrators.""" +from __future__ import annotations + import logging from abc import ABC, abstractmethod from dataclasses import dataclass, field import numpy as np -from numpy.random import Generator, default_rng +from numpy.random import Generator +from turtlemd.random import create_random_generator from turtlemd.system.system import System LOGGER = logging.getLogger(__name__) LOGGER.addHandler(logging.NullHandler()) +INTEGRATORS = {} + + +def register_integrator(name=None): + """Register integrators in the global registry.""" + + def decorator(cls): + nonlocal name + if name is None: + name = cls.__name__ + INTEGRATORS[name.lower()] = cls + return cls + + return decorator + + class MDIntegrator(ABC): """Base class for MD integrators.""" @@ -30,9 +49,11 @@ def integration_step(self, system: System): """Integrate the system a single time step.""" def __call__(self, system: System): + """Perform the integration step.""" return self.integration_step(system) +@register_integrator() class Verlet(MDIntegrator): """The Verlet integrator.""" @@ -41,6 +62,7 @@ class Verlet(MDIntegrator): previous_pos: np.ndarray | None # Positions at the previous step. def __init__(self, timestep: float): + """Set up the Verlet integrator.""" super().__init__( timestep=timestep, description="Verlet integrator", @@ -51,6 +73,7 @@ def __init__(self, timestep: float): self.previous_pos = None def integration_step(self, system: System): + """Do one step with the Verlet integrator.""" particles = system.particles if self.previous_pos is None: self.previous_pos = particles.pos - particles.vel * self.timestep @@ -68,12 +91,14 @@ def integration_step(self, system: System): system.potential_and_force() +@register_integrator() class VelocityVerlet(MDIntegrator): """The Velocity Verlet integrator.""" half_timestep: float # Half of the time step def __init__(self, timestep: float): + """Set ut the Velocity Verlet integrator.""" super().__init__( timestep=timestep, description="Velocity Verlet integrator", @@ -82,6 +107,7 @@ def __init__(self, timestep: float): self.half_timestep = self.timestep * 0.5 def integration_step(self, system: System): + """Do one integration step with the Velocity Verlet integrator.""" particles = system.particles imass = particles.imass # Update velocity @@ -94,6 +120,7 @@ def integration_step(self, system: System): particles.vel += self.half_timestep * particles.force * imass +@register_integrator() class LangevinOverdamped(MDIntegrator): """Overdamped version of the Langevin integrator. @@ -120,8 +147,14 @@ class LangevinOverdamped(MDIntegrator): _initiate: bool # If True, we still need to set some parameters def __init__( - self, timestep: float, gamma: float, rgen: Generator, beta: float + self, + timestep: float, + gamma: float, + beta: float, + rgen: Generator | None = None, + seed: int | None = None, ): + """Set up the overdamped Langevin integrator.""" super().__init__( timestep=timestep, description="Langevin overdamped integrator", @@ -129,7 +162,10 @@ def __init__( ) self.gamma = gamma self.sigma = 0.0 - self.rgen = rgen + if rgen is None: + self.rgen = create_random_generator(seed=seed) + else: + self.rgen = rgen self.bddt = 0.0 self.beta = beta self._initiate = True @@ -177,6 +213,7 @@ class LangevinParameter: cho: list[np.ndarray] = field(default_factory=list) +@register_integrator() class LangevinInertia(MDIntegrator): """The `Langevin`_ integrator. @@ -204,16 +241,17 @@ def __init__( gamma: float, beta: float, rgen: Generator | None = None, - seed: int = 0, + seed: int | None = None, ): + """Set up the Langevin integrator.""" super().__init__( timestep=timestep, - description="Langevin overdamped integrator", + description="Langevin integrator", dynamics="stochastic", ) self.gamma = gamma if rgen is None: - self.rgen = default_rng(seed=seed) + self.rgen = create_random_generator(seed=seed) else: self.rgen = rgen self.beta = beta @@ -284,7 +322,7 @@ def integration_step(self, system: System): def draw_random_numbers( self, system: System ) -> tuple[np.ndarray, np.ndarray]: - """This method draws random numbers for the integration step.""" + """Draws random numbers for the integration step.""" particles = system.particles dim = particles.dim pos_rand = np.zeros_like(particles.pos) diff --git a/turtlemd/potentials/__init__.py b/turtlemd/potentials/__init__.py index ab3547f..b1cf397 100644 --- a/turtlemd/potentials/__init__.py +++ b/turtlemd/potentials/__init__.py @@ -1,3 +1,4 @@ +"""Import and make potentials available.""" from .lennardjones import LennardJonesCut from .well import DoubleWell, DoubleWellPair diff --git a/turtlemd/potentials/potential.py b/turtlemd/potentials/potential.py index bef810f..3570163 100644 --- a/turtlemd/potentials/potential.py +++ b/turtlemd/potentials/potential.py @@ -53,7 +53,7 @@ def potential_and_force( ) -> tuple[float, np.ndarray, np.ndarray]: """Evaluate potential & force. - It may be more efficient to do both togehter. If not changed + It may be more efficient to do both together. If not changed in the subclasses, it will just use the other functions. """ vpot = self.potential(system) diff --git a/turtlemd/random.py b/turtlemd/random.py new file mode 100644 index 0000000..089f173 --- /dev/null +++ b/turtlemd/random.py @@ -0,0 +1,13 @@ +"""Make random number generators available. + +This module is just present to make sure that all random +number generators are created by the same method. +""" +from __future__ import annotations + +from numpy.random import Generator, default_rng + + +def create_random_generator(seed: int | None = None) -> Generator: + """Create a random generator.""" + return default_rng(seed=seed) diff --git a/turtlemd/simulation/__init__.py b/turtlemd/simulation/__init__.py index 62e55e4..1b53c2f 100644 --- a/turtlemd/simulation/__init__.py +++ b/turtlemd/simulation/__init__.py @@ -1,3 +1,4 @@ +"""Make the main simulation class(es) available.""" from .mdsimulation import MDSimulation __all__ = ["MDSimulation"] diff --git a/turtlemd/system/particles.py b/turtlemd/system/particles.py index 46458b4..4dae6c1 100644 --- a/turtlemd/system/particles.py +++ b/turtlemd/system/particles.py @@ -1,6 +1,6 @@ """Define a list of particles. -Tha particle list collects the mass, positions, velocities, +The particle list collects the mass, positions, velocities, forces, types, etc. for a collection of particles. """ from collections.abc import Iterator @@ -56,7 +56,16 @@ def add_particle( name: str = "?", ptype: int = 1, ): - """Add a single particle to the particle list.""" + """Add a single particle to the particle list. + + Args: + pos: Particle's position. + vel: Particle's velocity. Defaults to zero. + force: Force on the particle. Defaults to zero. + mass: Particle's mass. Defaults to 1.0. + name: Name of the particle. Defaults to "?". + ptype: Type of the particle. Defaults to 1. + """ if vel is None: vel = np.zeros_like(pos) if force is None: @@ -79,10 +88,10 @@ def add_particle( self.npart += 1 def __iter__(self) -> Iterator[dict]: - """Yield the properties of the particles. + """Iterate over the properties of the particles. Yields: - out: The information in `self.pos`, `self.vel`, ... etc. + The information in `self.pos`, `self.vel`, ... etc. """ for i, pos in enumerate(self.pos): part = { @@ -97,7 +106,7 @@ def __iter__(self) -> Iterator[dict]: yield part def __len__(self) -> int: - """Just give the number of particles.""" + """Return the number of particles.""" return self.npart def __getitem__(self, key: SupportsIndex | tuple[SupportsIndex, ...]): @@ -121,10 +130,11 @@ def pairs(self) -> Iterator[tuple[int, int, int, int]]: """Iterate over all pairs of particles. Yields: - out[0]: The index for the first particle in the pair. - out[1]: The index for the second particle in the pair. - out[2]: The particle type of the first particle. - out[3]: The particle type of the second particle. + A tuple with + - the index for the first particle in the pair, + - the index for the second particle in the pair, + - the particle type of the first particle, + - the particle type of the second particle. """ for i, itype in enumerate(self.ptype[:-1]): @@ -133,7 +143,7 @@ def pairs(self) -> Iterator[tuple[int, int, int, int]]: def linear_momentum(particles: Particles) -> np.ndarray: - """Return linear momentum of the particles.""" + """Calculate the linear momentum of the particles.""" return np.sum(particles.vel * particles.mass, axis=0) @@ -171,13 +181,14 @@ def kinetic_energy(particles: Particles) -> tuple[np.ndarray, float]: def kinetic_temperature( particles: Particles, boltzmann: float, - dof: list[float] | None = None, + dof: np.ndarray | list[float] | None = None, kin_tensor: np.ndarray | None = None, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Calculate the kinetic temperature of a collection of particles. Args: - boltzmann: This is the Boltzmann factor/constant in correct units. + particles: The particles to calculate the kinetic temperature of. + boltzmann: This is the Boltzmann constant in simulation units. dof: The degrees of freedom to subtract. Its shape should be equal to the number of dimensions. kin_tensor: The kinetic energy tensor. If the kinetic energy @@ -187,7 +198,6 @@ def kinetic_temperature( out[0]: The temperature averaged over all dimensions. out[1]: The temperature for each spatial dimension. out[2]: The kinetic energy tensor. - """ ndof = particles.npart * np.ones(particles.vel[0].shape) @@ -215,7 +225,6 @@ def pressure_tensor( Returns: out[0]: The symmetric pressure tensor. out[1] : The scalar pressure. - """ if kin_tensor is None: kin_tensor, _ = kinetic_energy(particles) @@ -229,7 +238,7 @@ def generate_maxwell_velocities( rgen: Generator, temperature: float = 1.0, boltzmann: float = 1.0, - dof: list[float] | None = None, + dof: np.ndarray | list[float] | None = None, momentum: bool = True, ): """Generate velocities from a Maxwell distribution. @@ -250,13 +259,11 @@ def generate_maxwell_velocities( particles: The particles we will set the velocity for. rgen: The random number generator used for drawing velocities. temperature: The desired temperature. - boltzmann: This is the Boltzmann factor/constant in - correct units. + boltzmann: This is the Boltzmann constant in simulation units. dof: The degrees of freedom to subtract. Its shape should be equal to the number of dimensions. momentum: If False, we will not zero the linear momentum. """ - vel = np.sqrt(particles.imass) * rgen.normal( loc=0.0, scale=1.0, size=particles.vel.shape ) diff --git a/turtlemd/system/system.py b/turtlemd/system/system.py index fc9ee7e..18ba317 100644 --- a/turtlemd/system/system.py +++ b/turtlemd/system/system.py @@ -5,15 +5,17 @@ import numpy as np -if TYPE_CHECKING: # pragma: no cover - from turtlemd.system.box import Box - from turtlemd.system.particles import Particles - from turtlemd.system.particles import ( kinetic_energy, kinetic_temperature, pressure_tensor, ) +from turtlemd.units import UNIT_SYSTEMS + +if TYPE_CHECKING: # pragma: no cover + from turtlemd.system.box import Box, TriclinicBox + from turtlemd.system.particles import Particles + from turtlemd.units import UnitSystem class Thermo(TypedDict): @@ -28,19 +30,21 @@ class Thermo(TypedDict): class System: """A system the MD is run on. - The system bridges some other objects together, like the particles - and the box. + A system consist of a simulation box, the particles, the + potential functions, and definition of units. """ - box: Box # The simulation box. + box: Box | TriclinicBox # The simulation box. particles: Particles # The particles in the system. potentials: list[Any] # The force field. + units: UnitSystem # Conversion factors and Boltzmann's constant. def __init__( self, - box: Box, + box: Box | TriclinicBox, particles: Particles, + units: str = "reduced", potentials: list[Any] | None = None, ): """Initialize a new system. @@ -54,6 +58,7 @@ def __init__( self.potentials = [] if potentials is not None: self.potentials = list(potentials) + self.units = UNIT_SYSTEMS[units] def potential_and_force(self): """Evaluate the potential energy and the force.""" @@ -76,7 +81,7 @@ def potential_and_force(self): return v_pot, force, virial def potential(self): - """ "Evaluate the potential energy.""" + """Evaluate the potential energy.""" v_pot = sum(pot.potential(self) for pot in self.potentials) self.particles.v_pot = v_pot return v_pot @@ -98,7 +103,7 @@ def force(self): self.particles.virial = virial return force, virial - def thermo(self, boltzmann: float = 1.0) -> Thermo: + def thermo(self) -> Thermo: """Evaluate simple thermodynamic properties for the system.""" thermo: Thermo = { "ekin": None, @@ -124,7 +129,19 @@ def thermo(self, boltzmann: float = 1.0) -> Thermo: ) thermo["pressure"] = pressure - dof = getattr(self.box, "dof", None) - temp, _, _ = kinetic_temperature(particles, boltzmann, dof=dof) + dof = self.dof() + temp, _, _ = kinetic_temperature( + particles, self.units.boltzmann, dof=dof + ) thermo["temperature"] = float(temp) return thermo + + def dof(self) -> np.ndarray | None: + """Extract the degrees of freedom of the system.""" + dof = getattr(self.box, "dof", None) + return dof + + def __str__(self): + """Write some info about the system.""" + msg = f"TurtleMD system with {len(self.particles)} particles." + return msg diff --git a/turtlemd/tools/__init__.py b/turtlemd/tools/__init__.py index 2893ecd..f13fef7 100644 --- a/turtlemd/tools/__init__.py +++ b/turtlemd/tools/__init__.py @@ -1,4 +1,6 @@ +"""Make tools available.""" from .tools import generate_lattice -from .xyz import particles_from_xyz_file, read_xyz_file -__all__ = ["generate_lattice", "particles_from_xyz_file", "read_xyz_file"] +__all__ = [ + "generate_lattice", +] diff --git a/turtlemd/units.py b/turtlemd/units.py new file mode 100644 index 0000000..622d32a --- /dev/null +++ b/turtlemd/units.py @@ -0,0 +1,111 @@ +"""Define unit conversions for TurtleMD.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np +from pint import UnitRegistry +from scipy.constants import Avogadro as AVOGADRO # type: ignore +from scipy.constants import Boltzmann as BOLTZMANN + +if TYPE_CHECKING: # pragma: no cover + from pint.facets.plain.quantity import PlainQuantity + +Quantity = UnitRegistry().Quantity + + +def time_unit( + length: PlainQuantity, mass: PlainQuantity, energy: PlainQuantity +) -> PlainQuantity: + """Calculate the time unit given the length, mass, and energy units. + + A system of units will typically define the length (L), mass (M), + and time (T) units. The energy unit (E) is then given as E = ML²/T². + With MD, it is often convenient to define the length, mass, and + energy units. The time unit is then given by: T = (ML²/E)^(1/2). + This method calculates the time unit as stated above and returns + the answer in fs. + + Args: + length: The length unit. + mass: The mass unit. + energy: The energy unit. + + Returns: + The time unit in fs. + """ + time_s = np.sqrt(length.to("m") ** 2 * mass.to("kg") / energy.to("J")) + return Quantity(time_s.magnitude, "s").to("fs") + + +class UnitSystem: + """Define a simple unit system.""" + + def __init__( + self, + name: str, + length: PlainQuantity, + mass: PlainQuantity, + energy: PlainQuantity, + boltzmann: float, + input_time_unit: str = "internal", + ): + """Set up the unit system. + + In particular, the time unit is calculated here. + """ + self.name = name + self.base_units = { + "length": length, + "mass": mass, + "energy": energy, + } + self.time_fs = time_unit(length, mass, energy) + self.boltzmann = boltzmann + if input_time_unit == "internal": + self.time_factor = 1.0 + else: + self.time_factor = 1.0 / self.time_fs.to(input_time_unit).magnitude + + def __str__(self): + """Return information on values and units as text.""" + msg = [ + f'# Unit system "{self.name}"', + f"\t* Time: {self.time_fs}", + f"\t* Time factor: {self.time_factor}", + f"\t* Boltzmann: {self.boltzmann}", + ] + for key, val in self.base_units.items(): + msg.append(f"\t* {key.capitalize()}: {val}") + return "\n".join(msg) + + +UNIT_SYSTEMS: dict[str, UnitSystem] = {} + +UNIT_SYSTEMS["reduced"] = UnitSystem( + name="reduced", + length=Quantity(1.0, "Å"), + mass=Quantity(1.0 / AVOGADRO, "g"), + energy=Quantity(BOLTZMANN, "J"), + boltzmann=1.0, + input_time_unit="internal", +) + +UNIT_SYSTEMS["lj"] = UnitSystem( + name="lj", + length=Quantity(3.405, "Å"), + mass=Quantity(39.948 / AVOGADRO, "g"), + energy=Quantity(BOLTZMANN, "J") * 119.8, + boltzmann=1.0, + input_time_unit="internal", +) + +UNIT_SYSTEMS["real"] = UnitSystem( + name="real", + length=Quantity(1.0, "Å"), + mass=Quantity(1.0 / AVOGADRO, "g"), + energy=Quantity(1.0 / AVOGADRO, "kcal"), + boltzmann=(Quantity(BOLTZMANN, "J/K").to("kcal/K") * AVOGADRO).magnitude, + input_time_unit="fs", +)