Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
73 changes: 39 additions & 34 deletions pack_mm/cli/packmm.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
from typer import Exit, Option, Typer
from typer_config import use_config

from pack_mm.core.core import pack_molecules
from pack_mm.core.core import pack_molecules, setup_file_logger


class InsertionMethod(str, Enum):
Expand Down Expand Up @@ -133,41 +133,45 @@ def packmm(
),
geometry: bool = Option(True, help="Perform geometry optimization at the end."),
out_path: str = Option(".", help="path to save various outputs."),
log: str = Option("pack-mm.log", help="file to save logs."),
):
"""Pack molecules into a system based on the specified parameters."""
print("Script called with following input")
print(f"{system=}")
print(f"{nmols=}")
print(f"{molecule=}")
print(f"{ntries=}")
print(f"{seed=}")
print(f"where={where.value}")
print(f"{centre=}")
print(f"{radius=}")
print(f"{height=}")
print(f"{a=}")
print(f"{b=}")
print(f"{c=}")
print(f"{cell_a=}")
print(f"{cell_b=}")
print(f"{cell_c=}")
print(f"{arch=}")
print(f"{model=}")
print(f"{dispersion=}")
print(f"{device=}")
print(f"{temperature=}")
print(f"{fmax=}")
print(f"{threshold=}")
print(f"{geometry=}")
print(f"{out_path=}")
print(f"{every=}")
print(f"insert_strategy={insert_strategy.value}")
print(f"relax_strategy={relax_strategy.value}")
print(f"{md_steps=}")
print(f"{md_timestep=}")
print(f"{md_temperature=}")
logger = setup_file_logger(log_file=log)
logger.info("Script called with following input")
logger.info(f"{system=}")
logger.info(f"{nmols=}")
logger.info(f"{molecule=}")
logger.info(f"{ntries=}")
logger.info(f"{seed=}")
logger.info(f"where={where.value}")
logger.info(f"{centre=}")
logger.info(f"{radius=}")
logger.info(f"{height=}")
logger.info(f"{a=}")
logger.info(f"{b=}")
logger.info(f"{c=}")
logger.info(f"{cell_a=}")
logger.info(f"{cell_b=}")
logger.info(f"{cell_c=}")
logger.info(f"{arch=}")
logger.info(f"{model=}")
logger.info(f"{dispersion=}")
logger.info(f"{device=}")
logger.info(f"{temperature=}")
logger.info(f"{fmax=}")
logger.info(f"{threshold=}")
logger.info(f"{geometry=}")
logger.info(f"{out_path=}")
logger.info(f"{every=}")
logger.info(f"insert_strategy={insert_strategy.value}")
logger.info(f"relax_strategy={relax_strategy.value}")
logger.info(f"{md_steps=}")
logger.info(f"{md_timestep=}")
logger.info(f"{md_temperature=}")
logger.info(f"log_file={log}")
print(f"Output is logger to {log}")
if nmols == -1:
print("nothing to do, no molecule to insert")
logger.info("nothing to do, no molecule to insert")
raise Exit(0)

center = centre
Expand All @@ -176,7 +180,7 @@ def packmm(
lc = [x < 0.0 for x in center]
if len(center) != 3 or any(lc):
err = "Invalid centre 3 coordinates expected!"
print(f"{err}")
logger.info(f"{err}")
raise Exception("Invalid centre 3 coordinates expected!")

pack_molecules(
Expand Down Expand Up @@ -209,4 +213,5 @@ def packmm(
md_steps=md_steps,
md_timestep=md_timestep,
md_temperature=md_temperature,
logger=logger,
)
101 changes: 74 additions & 27 deletions pack_mm/core/core.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@

from __future__ import annotations

from copy import copy
import logging
from pathlib import Path

from ase import Atoms
Expand All @@ -18,6 +20,29 @@
from numpy import cos, exp, pi, random, sin, sqrt


def setup_file_logger(log_file="pack-mm.log", log_level=logging.INFO):
"""
Set up a basic file logger.

Args:
log_file (str): The name of the file to write logs to.
log_level (int): The minimum level of messages to log
(e.g., logging.DEBUG, logging.INFO).
"""
logger = logging.getLogger(__name__)
logger.setLevel(log_level)

file_handler = logging.FileHandler(log_file)

formatter = logging.Formatter(
"%(asctime)s - %(name)s - %(levelname)s - %(message)s"
)
file_handler.setFormatter(formatter)
logger.addHandler(file_handler)

return logger


def random_point_in_sphere(c: (float, float, float), r: float) -> (float, float, float):
"""
Generate a random point inside a sphere of radius r, centered at c.
Expand Down Expand Up @@ -131,11 +156,14 @@ def random_point_in_cylinder(
return (x, y, z)


def validate_value(label: str, x: float | int) -> None:
def validate_value(label: str, x: float | int, logger: logging.logger = None) -> None:
"""Validate input value, and raise an exception."""
if x is not None and x < 0.0:
err = f"Invalid {label}, needs to be positive"
print(err)
if logger:
logger.info(err)
else:
print(err)
raise Exception(err)


Expand Down Expand Up @@ -216,6 +244,7 @@ def pack_molecules(
md_steps: int = 10,
md_timestep: float = 1.0,
md_temperature: float = 100.0,
logger: logging.logger = None,
) -> tuple(float, Atoms):
"""
Pack molecules into a system based on the specified parameters.
Expand Down Expand Up @@ -259,22 +288,22 @@ def pack_molecules(

"""
kbt = temperature * kB
validate_value("temperature", temperature)
validate_value("radius", radius)
validate_value("height", height)
validate_value("fmax", fmax)
validate_value("seed", seed)
validate_value("box a", a)
validate_value("box b", b)
validate_value("box c", c)
validate_value("ntries", ntries)
validate_value("cell box cell a", cell_a)
validate_value("cell box cell b", cell_b)
validate_value("cell box cell c", cell_c)
validate_value("nmols", nmols)
validate_value("MD steps", md_steps)
validate_value("MD timestep", md_timestep)
validate_value("MD temperature", md_temperature)
validate_value("temperature", temperature, logger)
validate_value("radius", radius, logger)
validate_value("height", height, logger)
validate_value("fmax", fmax, logger)
validate_value("seed", seed, logger)
validate_value("box a", a, logger)
validate_value("box b", b, logger)
validate_value("box c", c, logger)
validate_value("ntries", ntries, logger)
validate_value("cell box cell a", cell_a, logger)
validate_value("cell box cell b", cell_b, logger)
validate_value("cell box cell c", cell_c, logger)
validate_value("nmols", nmols, logger)
validate_value("MD steps", md_steps, logger)
validate_value("MD timestep", md_timestep, logger)
validate_value("MD temperature", md_temperature, logger)

set_random_seed(seed)

Expand All @@ -289,9 +318,13 @@ def pack_molecules(
sysname = Path(system).stem + "+"

# Print summary
print(f"Inserting {nmols} {molecule} molecules in {sysname}.")
print(f"Using {arch} model {model} on {device}.")
print(f"Insert in {where}.")
summary = f"""Inserting {nmols} {molecule} molecules in {sysname}.
Using {arch} model {model} on {device}.
Insert in {where}."""
if logger:
logger.info(summary)
else:
print(summary)

cell = sys.cell.lengths()

Expand All @@ -305,12 +338,12 @@ def pack_molecules(
device=device,
dispersion=dispersion,
)
sys.calc = calc
sys.calc = copy(calc)

e = sys.get_potential_energy() if len(sys) > 0 else 0.0

mol = load_molecule(molecule)
mol.calc = calc
mol.calc = copy(calc)
emol = mol.get_potential_energy()

csys = sys.copy()
Expand Down Expand Up @@ -351,13 +384,17 @@ def pack_molecules(
relax_strategy=relax_strategy,
)

tsys.calc = calc
tsys.calc = copy(calc)
en = tsys.get_potential_energy()
de = en - e

acc = exp(-de / kbt)
u = random.random()
print(f"Old energy={e}, new energy={en}, {de=}, {acc=}, random={u}")
message_ene = f"Old energy={e}, new energy={en}, {de=}, {acc=}, random={u}"
if logger:
logger.info(message_ene)
else:
print(message_ene)

if abs(de / emol) > threshold and u <= acc:
accept = True
Expand All @@ -366,12 +403,22 @@ def pack_molecules(
csys = tsys.copy()
e = en
i += 1
print(f"Inserted particle {i}")
message_insert = f"Inserted particle {i}"
if logger:
logger.info(message_insert)
else:
print(message_insert)
write(Path(out_path) / f"{sysname}{i}{Path(molecule).stem}.cif", csys)
else:
# Things are bad, maybe geomatry optimisation saves us
# once you hit here is bad, this can keep looping
print(f"Failed to insert particle {i + 1} after {ntries} tries")
message_fail = f"""Failed to insert particle {i + 1} after {ntries} tries.
Trying alternative methods, {relax_strategy}, to save the day.
"""
if logger:
logger.info(message_fail)
else:
print(message_fail)
csys = save_the_day(
csys,
device,
Expand Down
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[project]
name = "pack-mm"
version = "0.1.12"
version = "0.1.14"
description = "packing materials and molecules in boxes using for machine learnt interatomic potentials"
authors = [
{ name = "Alin M. Elena" },
Expand Down
4 changes: 2 additions & 2 deletions tests/test_advanced.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@

runner = CliRunner()

err = 1.0e-8
err = 1.0e-5


def test_packmm_hmc(tmp_path):
Expand Down Expand Up @@ -84,4 +84,4 @@ def test_packmm_every(tmp_path):
assert (tmp_path / "system+1H2.cif").exists()
assert (tmp_path / "system+2H2.cif").exists()
f = read(tmp_path / "system+2H2.cif")
assert f[0].position == pytest.approx([1.24701529, 0.024216, 0.11632905], abs=err)
assert f[0].position == pytest.approx([1.247015, 0.024216, 0.116329], abs=err)
Loading