robin/polyamorphism_optimization
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

added pre-commit hooks; fixed whitespace and capitalization dependent parameter replacement; fixed creation of annealing and emin sim directories

Browse Source
This commit is contained in:
robrobo
2026-06-26 17:38:59 +02:00
parent 3c9edcfa74
commit 3f05cc1caf
21 changed files with 224 additions and 100 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.editorconfig
+12
View File
@@ -0,0 +1,12 @@
# .editorconfig
root = true
[*]
trim_trailing_whitespace = true
# Disable for GROMACS template files
[*.mdp]
trim_trailing_whitespace = false
[*.top]
trim_trailing_whitespace = false
.gitignore
+2 -1
View File
@@ -2,6 +2,8 @@
# gitignore template for Jupyter Notebooks
# website: http://jupyter.org/
.idea
.ipynb_checkpoints
*/.ipynb_checkpoints/*
@@ -187,4 +189,3 @@ cython_debug/
# PyPI configuration file
.pypirc
.pre-commit-config.yaml
+28
View File
@@ -0,0 +1,28 @@
repos:
- repo: https://github.com/pre-commit/pre-commit-hooks
rev: v6.0.0
hooks:
- id: trailing-whitespace
exclude: '\.(mdp|top)$'
- id: end-of-file-fixer
exclude: '\.(mdp|top)$'
- id: check-yaml
- id: check-added-large-files
- repo: https://github.com/psf/black
rev: 26.5.1
hooks:
- id: black
exclude: '\.(mdp|top)$'
- repo: https://github.com/PyCQA/isort
rev: 9.0.0a3
hooks:
- id: isort
exclude: '\.(mdp|top)$'
- repo: https://github.com/PyCQA/flake8
rev: 7.3.0
hooks:
- id: flake8
exclude: '\.(mdp|top)$'
README.md
+9 -1
View File
@@ -1,3 +1,11 @@
# polyamorphism_optimization
Code snippets for creation of simulations, analysis, and optimized parameter search for liquids with polyamorphism. Goal is to find parameters for which dynamics is still fast at the critical point and the dynamic anomaly that follows can be studied more easily.
Code snippets for creation of simulations, analysis, and optimized parameter search for liquids with polyamorphism. Goal is to find parameters for which dynamics is still fast at the critical point and the dynamic anomaly that follows can be studied more easily.
#### Setting up pre-commit hooks
```
pip install pre-commit
pre-commit autoupdate
pre-commit install
pre-commit run #--all-files
```
polyamorphism_optimization/annealing.py
+25 -15
View File
@@ -1,40 +1,50 @@
#!/usr/bin/env python3
#!/usr/bin/env python3.12
import shutil
from importlib.resources import files
from pathlib import Path
from generate_annealing import generate_annealing
from replace_params import replace_params
def main(root_dir: Path):
simdir = root_dir / "01_an1"
sim_dir = root_dir / "01_an1"
if simdir.exists():
shutil.rmtree(simdir)
simdir.mkdir()
if sim_dir.exists():
shutil.rmtree(sim_dir)
sim_dir.mkdir()
# Copy topology and replace parameters
shutil.copy(root_dir / "topology.top", simdir / "topology.top")
replace_params(root_dir / "params.ini", root_dir / "topology.top", simdir / "topology.top")
replace_params(root_dir / "params_an1.ini", root_dir / "mdp_parameters.mdp", simdir / "mdp_parameters.mdp")
replace_params(root_dir / "params.ini", simdir / "mdp_parameters.mdp", simdir / "mdp_parameters.mdp")
shutil.copy(root_dir / "topology.top", sim_dir / "topology.top")
replace_params(root_dir / "params.ini", root_dir / "topology.top", sim_dir / "topology.top", upper=True, space=True)
replace_params(
files("polyamorphism_optimization.templates.params") / "annealing1.ini",
root_dir / "mdp_parameters.mdp",
sim_dir / "mdp_parameters.mdp",
upper=True,
space=True,
)
replace_params(
root_dir / "params.ini", sim_dir / "mdp_parameters.mdp", sim_dir / "mdp_parameters.mdp", upper=True, space=True
)
# Copy run.sh
shutil.copy(root_dir / "run.sh", simdir / "run.sh")
shutil.copy(root_dir / "run.sh", sim_dir / "run.sh")
# Generate annealing section
annealing = generate_annealing(d=50, l=50, u=1000, s=50, e=2000)
annealing = generate_annealing(dT=50, Tmin=50, Tmax=1000, start_time=50, end_time=2000)
# -d 20 -l 750 -u 2100 -s 50 -e 10000 -c
with (simdir / "mdp_parameters.mdp").open("a") as f:
f.write("\n")
with (sim_dir / "mdp_parameters.mdp").open("a") as f:
f.write(annealing)
print(f"Annealing directory {simdir} created and filled.")
print(f"Annealing directory {sim_dir} created and filled.")
if __name__ == "__main__":
import sys
if len(sys.argv) != 2:
print(f"Usage: {sys.argv[0]} ROOT_DIRECTORY")
sys.exit(1)
main(Path(sys.argv[1]))
polyamorphism_optimization/create_sim.py
+21 -15
View File
@@ -1,14 +1,16 @@
#!/usr/bin/python3.12
#!/usr/bin/env python3.12
import argparse
import configparser
import math
import os
import shutil
from importlib.resources import files
from pathlib import Path
def compute_dependent(params):
model = params["model"]
# Charge neutrality
params["charge_B"] = round(params["charge"]/2, 5)
params["charge_B"] = round(params["charge"] / 2, 5)
params["charge_A"] = -2 * params["charge_B"]
# Relative sigma and epsilon
@@ -25,16 +27,16 @@ def compute_dependent(params):
if params["mass"] is None:
# Let total mass m = 2 * m_A, solve: τ = σ * sqrt(m/ε) ≈ 0.1
tau_target = 0.05 * 10
m_total = (tau_target / sigma_A)**2 * epsilon_A
m_total = (tau_target / sigma_A) ** 2 * epsilon_A
m_atom = m_total / 2
params["mass"] = round(m_atom, 5)
params["d_BB"] = round(2 * params["d_AB"] * math.sin(math.radians(params["angle_ABA"]) / 2), 5)
params["dummy_ab"] = round(params["d_AM"] / (2 * math.cos(math.radians(params["angle_ABA"]) / 2)), 5)
return params
def clean_params(params):
if params["model"] != "4point":
for k in ["r_AM", "d_AM", "d_BB", "dummy_ab"]:
@@ -43,7 +45,8 @@ def clean_params(params):
for k in ["r_AB", "d_AB", "angle_ABA"]:
params.pop(k)
return params
def format_directory_name(params):
"""Generate a short directory name from active parameters."""
keys_order = [
@@ -72,10 +75,9 @@ def format_directory_name(params):
def main():
parser = argparse.ArgumentParser(description="Generate directory and parameter file for A-B molecular systems.")
parser.add_argument("--version", type=int, default=1, help="Version identifier for parameter sweep")
parser.add_argument("--model", choices=["atomistic", "3point", "4point"], default='atomistic')
parser.add_argument("--model", choices=["atomistic", "3point", "4point"], default="atomistic")
parser.add_argument("--pressure", type=float, default=1.0)
# Independent parameters
parser.add_argument("--sigma", type=float, default=0.3)
parser.add_argument("--epsilon", type=float, default=0.8)
@@ -90,7 +92,7 @@ def main():
# Geometry for bonded/virtual_site
parser.add_argument("--r_AB", type=float, default=0.33334)
parser.add_argument("--angle_ABA", type=float, default=109.5) # 104.52
parser.add_argument("--angle_ABA", type=float, default=109.5) # 104.52
parser.add_argument("--r_AM", type=float, default=0.4)
# Output
@@ -113,13 +115,14 @@ def main():
"angle_ABA": args.angle_ABA,
"r_AM": args.r_AM,
}
if params["model"] == "atomistic":
params["nmolb"] = 2 * params["nmol"]
# Compute dependent values
params = compute_dependent(params)
params = clean_params(params)
dir_name = format_directory_name(params)
dir_name = Path(format_directory_name(params))
if os.path.exists(dir_name):
print(f"Directory {dir_name} already exists. Aborting.")
return
@@ -129,11 +132,14 @@ def main():
config_path = os.path.join(dir_name, args.output)
with open(config_path, "w") as f:
config = configparser.ConfigParser()
config.optionxform = str # keep case sensitivity
config["params"] = params
config.write(f)
print(f"Parameter file written to {config_path}")
shutil.copy(files("polyamorphism_optimization.templates") / "run.sh", dir_name / "run.sh")
if __name__ == "__main__":
main()
polyamorphism_optimization/emin.py
+24 -12
View File
@@ -1,35 +1,47 @@
#!/usr/bin/env python3
#!/usr/bin/env python3.12
import shutil
from importlib.resources import files
from pathlib import Path
from replace_params import replace_params
def main(root_dir: Path):
simdir = root_dir / "00_em"
sim_dir = root_dir / "00_em"
if simdir.exists():
shutil.rmtree(simdir)
simdir.mkdir()
if sim_dir.exists():
shutil.rmtree(sim_dir)
sim_dir.mkdir()
# Copy topology and replace parameters
shutil.copy(root_dir / "topology.top", simdir / "topology.top")
replace_params(root_dir / "params.ini", root_dir / "topology.top", simdir / "topology.top")
replace_params(root_dir / "params_emin.ini", root_dir / "mdp_parameters.mdp", simdir / "mdp_parameters.mdp")
replace_params(root_dir / "params.ini", simdir / "mdp_parameters.mdp", simdir / "mdp_parameters.mdp")
shutil.copy(root_dir / "topology.top", sim_dir / "topology.top")
replace_params(root_dir / "params.ini", root_dir / "topology.top", sim_dir / "topology.top", upper=True, space=True)
replace_params(
files("polyamorphism_optimization.templates.params") / "emin.ini",
root_dir / "mdp_parameters.mdp",
sim_dir / "mdp_parameters.mdp",
upper=True,
space=True,
)
replace_params(
root_dir / "params.ini", sim_dir / "mdp_parameters.mdp", sim_dir / "mdp_parameters.mdp", upper=True, space=True
)
# Copy run.sh
shutil.copy(root_dir / "run.sh", simdir / "run.sh")
shutil.copy(root_dir / "run.sh", sim_dir / "run.sh")
# Add finer emstep
with (simdir / "mdp_parameters.mdp").open("a") as f:
with (sim_dir / "mdp_parameters.mdp").open("a") as f:
f.write("\n")
f.write("emstep = 0.001")
print(f"Energy minimization directory {simdir} created and filled.")
print(f"Energy minimization directory {sim_dir} created and filled.")
if __name__ == "__main__":
import sys
if len(sys.argv) != 2:
print(f"Usage: {sys.argv[0]} ROOT_DIRECTORY")
sys.exit(1)
polyamorphism_optimization/generate_annealing.py
+30 -15
View File
@@ -1,7 +1,8 @@
#!/usr/bin/python3.12
#!/usr/bin/env python3.12
import argparse
import math
def round_dynamic(x, rel_tol=1e-3, max_decimals=6, max_negative_decimals=3):
"""
Dynamically round x so that the relative error is within rel_tol.
@@ -23,6 +24,7 @@ def round_dynamic(x, rel_tol=1e-3, max_decimals=6, max_negative_decimals=3):
factor = 10 ** (-decimals)
return round(x / factor) * factor
def generate_temperature_list_dynamic(dT=10, Tmin=10, Tmax=5000, rel_tol=5e-3):
r = (300 - dT) / 300.0
@@ -49,16 +51,18 @@ def generate_temperature_list_dynamic(dT=10, Tmin=10, Tmax=5000, rel_tol=5e-3):
rounded_T = [round_dynamic(t, rel_tol=rel_tol) for t in all_T]
return rounded_T
def generate_temperature_list(dT=10, Tmin=10, Tmax=5000, rel_tol=5e-3):
r = (Tmax/Tmin)**(1/20)
all_T = [Tmin*r**i for i in range(21)]
r = (Tmax / Tmin) ** (1 / 20)
all_T = [Tmin * r**i for i in range(21)]
# Apply dynamic rounding
rounded_T = [round_dynamic(t, rel_tol=rel_tol) for t in all_T]
return rounded_T
def generate_times_list(npoints, start_time, end_time):
if npoints == 1:
return [0.0]
@@ -66,6 +70,7 @@ def generate_times_list(npoints, start_time, end_time):
interval = duration / (npoints - 1)
return [round(start_time + i * interval, 3) for i in range(npoints)]
def main():
parser = argparse.ArgumentParser(description="Generate annealing times and temperatures for GROMACS.")
parser.add_argument("-d", "--dT", type=float, default=20.0, help="Temperature step size")
@@ -74,25 +79,35 @@ def main():
parser.add_argument("-r", "--rel_tol", type=float, default=5e-3, help="Relative tolerance for rounding")
parser.add_argument("-s", "--start_time", type=float, default=0, help="Start time in fs")
parser.add_argument("-e", "--end_time", type=float, required=True, help="End time in fs")
parser.add_argument("-c", "--cooling", action='store_true', help="Do cooling instead of heating")
parser.add_argument("-c", "--cooling", action="store_true", help="Do cooling instead of heating")
args = parser.parse_args()
annealing = generate_annealing(args.e, args.d, args.l, args.u, args.r, args.c, args.s)
print(annealing)
#temps = generate_temperature_list_dynamic(args.dT, args.Tmin, args.Tmax, args.rel_tol)
temps = generate_temperature_list(args.dT, args.Tmin, args.Tmax, args.rel_tol)
if args.cooling:
def generate_annealing(end_time, dT=20, Tmin=10, Tmax=5000, rel_tol=5e-3, cooling=False, start_time=0):
# temps = generate_temperature_list_dynamic(args.dT, args.Tmin, args.Tmax, args.rel_tol)
temps = generate_temperature_list(dT, Tmin, Tmax, rel_tol)
if cooling:
temps = list(reversed(temps))
if args.start_time > 0:
if start_time > 0:
temps = [temps[0]] + temps # hold first temperature
times = [0.0] + generate_times_list(len(temps) - 1, args.start_time, args.end_time)
times = [0.0] + generate_times_list(len(temps) - 1, start_time, end_time)
else:
times = generate_times_list(len(temps), args.start_time, args.end_time)
times = generate_times_list(len(temps), start_time, end_time)
annealing_config = f"""
annealing = single
annealing-npoints = {len(times)}
annealing-time = {" ".join(f'{t:.3f}' for t in times)}
annealing-temp = {" ".join(f'{t:.2f}' for t in temps)}
"""
return annealing_config
print("annealing = single")
print(f"annealing-npoints = {len(times)}")
print("annealing-time = " + " ".join(f"{t:.3f}" for t in times))
print("annealing-temp = " + " ".join(f"{t:.2f}" for t in temps))
if __name__ == "__main__":
main()
polyamorphism_optimization/initialize_sim.py
+16 -15
View File
@@ -1,10 +1,9 @@
#!/usr/bin/env python3
#!/usr/bin/env python3.12
import argparse
import os
import subprocess
from pathlib import Path
from importlib.resources import files
import configparser
import os
from importlib.resources import files
from pathlib import Path
from annealing import main as annealing
from emin import main as emin
@@ -12,19 +11,19 @@ from emin import main as emin
# --- Core function ---
def initialize_simulation_dir(directory: Path):
os.chdir(directory)
params = configparser.ConfigParser()
params = params.read('params.ini')
directory = Path(directory)
template_path = files("polyamorphism_optimization.templates") / f"topology_{params['model']}.top"
params = configparser.ConfigParser()
params.read(directory / "params.ini")
template_path = files("polyamorphism_optimization.templates") / f"top_{params['params']['model']}.top"
TOPOLOGY = template_path.read_text()
template_path = files("polyamorphism_optimization.templates") / "mdp_parameters.mdp"
MDP_PARAMETERS = template_path.read_text()
# Write files
Path("mdp_parameters.mdp").write_text(MDP_PARAMETERS)
Path("topology.top").write_text(TOPOLOGY)
Path(directory / "mdp_parameters.mdp").write_text(MDP_PARAMETERS)
Path(directory / "topology.top").write_text(TOPOLOGY)
emin(directory)
annealing(directory)
@@ -33,14 +32,16 @@ def initialize_simulation_dir(directory: Path):
for f in os.listdir(directory):
print(f)
# --- CLI entry point ---
def main():
parser = argparse.ArgumentParser(description="Create simulation input files.")
parser.add_argument("-d", "--directory", default=".", help="Base directory")
args = parser.parse_args()
create_simulation_dir(Path(args.directory))
# create_simulation_dir(Path(args.directory))
initialize_simulation_dir(Path(args.directory))
if __name__ == "__main__":
main()
polyamorphism_optimization/replace_params.py
+12 -5
View File
@@ -1,10 +1,10 @@
#!/usr/bin/env python3
import sys
#!/usr/bin/env python3.12
import configparser
import sys
from pathlib import Path
def replace_params(param_file: Path, template_file: Path, output_file: Path):
def replace_params(param_file: Path, template_file: Path, output_file: Path, upper=False, space=False):
"""Replace placeholders in template_file with values from an INI file."""
config = configparser.ConfigParser()
config.optionxform = str # keep case sensitivity
@@ -14,7 +14,15 @@ def replace_params(param_file: Path, template_file: Path, output_file: Path):
params = {}
for section in config.sections():
for key, value in config.items(section):
params[key] = value
if upper:
params[key.upper()] = value
else:
params[key] = value
if space:
for key, value in dict(params.items()).items():
params[f" {key} "] = value
params.pop(key)
content = template_file.read_text()
@@ -38,4 +46,3 @@ def main():
if __name__ == "__main__":
main()
polyamorphism_optimization/templates/NPT.ini
+1
View File
@@ -1,3 +1,4 @@
[mdp]
INTEGRATOR = md
NSTEPS = 2000000
NLOG = 1000
polyamorphism_optimization/templates/NVT.ini
+1
View File
@@ -1,3 +1,4 @@
[mdp]
INTEGRATOR = md
NSTEPS = 2000000
NLOG = 1000
polyamorphism_optimization/templates/mdp_parameters.mdp
+11 -11
View File
@@ -1,11 +1,11 @@
integrator = INTEGRATOR
integrator = INTEGRATOR
dt = 0.001
nsteps = NSTEPS
nsteps = NSTEPS
nstcomm = 10
nstlog = NLOG
nstlog = NLOG
nstcalcenergy = 10
nstenergy = NENERGY
nstxout-compressed = NSTXOUT
nstenergy = NENERGY
nstxout-compressed = NSTXOUT
compressed-x-precision = 1000
energygrps = system
rlist = 1.2
@@ -18,14 +18,14 @@ rvdw = 1.2
vdw-modifier = potential-shift-verlet
fourierspacing = 0.144
tcoupl = v-rescale
tau_t = TAUT
tau_t = TAUT
gen_vel = yes
gen_temp = TEMP
ref_t = TEMP
gen_temp = TEMP
ref_t = TEMP
tc-grps = system
pcoupl = PCOUPL
tau_p = TAUP
ref_p = PRESSURE
pcoupl = PCOUPL
tau_p = TAUP
ref_p = PRESSURE
compressibility = 4.5e-5
nstpcouple = 10
nsttcouple = 10
polyamorphism_optimization/templates/params/annealing1.ini
+1
View File
@@ -1,3 +1,4 @@
[mdp]
INTEGRATOR = md
NSTEPS = 2000000
NLOG = 1000
polyamorphism_optimization/templates/params/annealing2.ini
+1
View File
@@ -1,3 +1,4 @@
[mdp]
INTEGRATOR = md
NSTEPS = 5000000
NLOG = 1000
polyamorphism_optimization/templates/params/emin.ini
+1
View File
@@ -1,3 +1,4 @@
[mdp]
INTEGRATOR = steep
NSTEPS = 10000
NLOG = 100
polyamorphism_optimization/templates/top_3point.top
+3 -3
View File
@@ -3,15 +3,15 @@
1 2 yes 1.0 1.0
[atomtypes]
;name mass charge ptype sigma epsilon
A MASS 0.000 A SIGMA_A EPSILON_A
B MASS 0.000 A SIGMA_B EPSILON_B
A MASS 0.000 A SIGMA_A EPSILON_A
B MASS 0.000 A SIGMA_B EPSILON_B
[moleculetype]
; name nrexcl
L 1
[atoms]
; nr type resnr residu atom cgnr charge
; nr type resnr residu atom cgnr charge mass
1 A 1 L A 1 CHARGE_A MASS
2 B 1 L B1 1 CHARGE_B MASS
3 B 1 L B2 1 CHARGE_B MASS
polyamorphism_optimization/templates/top_4point.top
+3 -3
View File
@@ -3,8 +3,8 @@
1 2 yes 1.0 1.0
[atomtypes]
;name mass charge ptype sigma epsilon
A MASS 0.000 A SIGMA_A EPSILON_A
B MASS 0.000 A SIGMA_B EPSILON_B
A MASS 0.000 A SIGMA_A EPSILON_A
B MASS 0.000 A SIGMA_B EPSILON_B
D 0 0.000 D 0.0 0.0
[moleculetype]
@@ -12,7 +12,7 @@ D 0 0.000 D 0.0 0.0
L 1
[atoms]
; nr type resnr residu atom cgnr charge
; nr type resnr residu atom cgnr charge mass
1 A 1 L A 1 0 MASS
2 B 1 L B1 1 CHARGE_B MASS
3 B 1 L B2 1 CHARGE_B MASS
polyamorphism_optimization/templates/top_atomistic.top
+4 -4
View File
@@ -12,7 +12,7 @@ B MASS 0.000 A SIGMA_B EPSILON_B
A 1
[atoms]
; nr type resnr residu atom cgnr charge
; nr type resnr residu atom cgnr charge mass
1 A 1 L A 1 CHARGE_A MASS
[moleculetype]
@@ -20,12 +20,12 @@ A 1
B 1
[atoms]
; nr type resnr residu atom cgnr charge
; nr type resnr residu atom cgnr charge mass
1 B 1 L B 1 CHARGE_B MASS
[system]
MODEL
[molecules]
A NMOL
B NMOLB
A NMOL
B NMOLB
pyproject.toml
+14
View File
@@ -22,3 +22,17 @@ po_init_sim = "polyamorphism_optimization.initialize_sim:main"
[tool.setuptools.package-data]
polyamorphism_optimization = ["templates/*"]
[tool.black]
line-length = 120
target-version = ["py312"]
exclude = '''
(
\.mdp
| \.top
)
'''
[tool.isort]
profile = "black" # Makes isort compatible with black's formatting
line_length = 120
extend_skip_glob = ["*.mdp", "*.top"]
setup.cfg
+5
View File
@@ -0,0 +1,5 @@
[flake8]
max-line-length = 120
exclude =
*.mdp,
*.top