From b02fe384151a41843d54e596b8cfc13c04373a85 Mon Sep 17 00:00:00 2001 From: robrobo Date: Sat, 9 Aug 2025 13:43:03 +0200 Subject: [PATCH] initial commit --- create_annealing1.sh | 41 + create_emin.sh | 35 + create_sim.py | 139 ++ generate_annealing.py | 97 + generate_temperatures.py | 57 + replace_params.sh | 47 + run.sh | 53 + setup_sim.sh | 147 ++ templates/4point.gro | 4003 ++++++++++++++++++++++++++++++++++++++ templates/bak.gro | 4003 ++++++++++++++++++++++++++++++++++++++ templates/single.gro | 7 + 11 files changed, 8629 insertions(+) create mode 100755 create_annealing1.sh create mode 100755 create_emin.sh create mode 100755 create_sim.py create mode 100755 generate_annealing.py create mode 100755 generate_temperatures.py create mode 100755 replace_params.sh create mode 100755 run.sh create mode 100755 setup_sim.sh create mode 100755 templates/4point.gro create mode 100755 templates/bak.gro create mode 100755 templates/single.gro diff --git a/create_annealing1.sh b/create_annealing1.sh new file mode 100755 index 0000000..ad0230b --- /dev/null +++ b/create_annealing1.sh @@ -0,0 +1,41 @@ +#!/bin/bash + +get_ini_value() { + local file=$1 + local section=$2 + local key=$3 + + awk -F '=' -v section="$section" -v key="$key" ' + $0 ~ "\\[" section "\\]" { in_section = 1; next } + in_section && $1 ~ "^" key"[[:space:]]*$" { + gsub(/[[:space:]]+/, "", $2) + print $2 + exit + } + $0 ~ /^\[/ { in_section = 0 } + ' "$file" +} + +simdir="01_an1" + + +rm -r "$simdir" 2>/dev/null +mkdir "$simdir" + +cp topology.top "$simdir/" +../replace_params.sh params.ini topology.top --output "$simdir/topology.top" +../replace_params.sh params_an1.ini mdp_parameters.mdp --output "$simdir/mdp_parameters.mdp" +../replace_params.sh params.ini "$simdir/mdp_parameters.mdp" --output "$simdir/mdp_parameters.mdp" +cp ../run.sh "$simdir/" +cp "00_em/out/out.gro" "$simdir/gro_start.gro" + + +annealing=$(../generate_annealing.py -d 50 -l 50 -u 1000 -s 100 -e 2000) + +annealing=$(../generate_annealing.py -d 30 -l 400 -u 950 -s 50 -e 2000) + +echo "" >> "$simdir/mdp_parameters.mdp" +echo "annealing = single" >> "$simdir/mdp_parameters.mdp" +echo "$annealing" >> "$simdir/mdp_parameters.mdp" + + diff --git a/create_emin.sh b/create_emin.sh new file mode 100755 index 0000000..69d886d --- /dev/null +++ b/create_emin.sh @@ -0,0 +1,35 @@ +#!/bin/bash + +get_ini_value() { + local file=$1 + local section=$2 + local key=$3 + + awk -F '=' -v section="$section" -v key="$key" ' + $0 ~ "\\[" section "\\]" { in_section = 1; next } + in_section && $1 ~ "^" key"[[:space:]]*$" { + gsub(/[[:space:]]+/, "", $2) + print $2 + exit + } + $0 ~ /^\[/ { in_section = 0 } + ' "$file" +} + + +simdir="00_em" + + +rm -r "$simdir" 2>/dev/null +mkdir "$simdir" + +cp topology.top "$simdir/" +../replace_params.sh params.ini topology.top --output "$simdir/topology.top" +../replace_params.sh params_emin.ini mdp_parameters.mdp --output "$simdir/mdp_parameters.mdp" +../replace_params.sh params.ini "$simdir/mdp_parameters.mdp" --output "$simdir/mdp_parameters.mdp" +cp ../run.sh "$simdir/" +cp "../gro_$(get_ini_value params.ini params model).gro" "$simdir/gro_start.gro" + +echo "emstep = 0.001" >> $simdir/mdp_parameters.mdp + + diff --git a/create_sim.py b/create_sim.py new file mode 100755 index 0000000..c41cbda --- /dev/null +++ b/create_sim.py @@ -0,0 +1,139 @@ +#!/usr/bin/python3.12 +import argparse +import configparser +import math +import os + +def compute_dependent(params): + model = params["model"] + + # Charge neutrality + params["charge_B"] = round(params["charge"]/2, 5) + params["charge_A"] = -2 * params["charge_B"] + + # Relative sigma and epsilon + params["sigma_A"] = params["sigma"] + params["epsilon_A"] = params["sigma"] + sigma_A = params["sigma"] + epsilon_A = params["epsilon"] + params["sigma_B"] = round(params["sigma_rel"] * sigma_A, 5) + params["epsilon_B"] = round(params["epsilon_rel"] * epsilon_A, 5) + params["d_AB"] = round(params["r_AB"] * sigma_A, 5) + params["d_AM"] = round(params["r_AM"] * sigma_A, 5) + + # Mass: either specified or computed to give τ ≈ 0.1 + 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_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"]: + params.pop(k) + if params["model"] != "3point" and params["model"] != "4point": + 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 = [ + ("sigma", "s"), + ("epsilon", "e"), + ("charge", "q"), + ("sigma_rel", "sr"), + ("epsilon_rel", "er"), + ("r_AB", "d"), + ("angle_ABA", "a"), + ("r_AM", "m"), + ] + + parts = [params["model"]] + parts.append(f"P{params['pressure']}") + + for key, prefix in keys_order: + if key in params: + val = params[key] + parts.append(f"{prefix}{val:.5f}".rstrip("0").rstrip(".")) + + parts.append(f"v{params['version']}") + return "_".join(parts) + + +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("--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) + parser.add_argument("--q", type=float, default=1.0) + + # Relative parameters + parser.add_argument("--sigma_rel", type=float, default=0.1) + parser.add_argument("--epsilon_rel", type=float, default=0.1) + + # Mass (optional, will be calculated if omitted) + parser.add_argument("--mass", type=float, default=None) + + # 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("--r_AM", type=float, default=0.4) + + # Output + parser.add_argument("--output", type=str, default="params.ini") + + args = parser.parse_args() + + params = { + "version": args.version, + "model": args.model, + "nmol": 1000, + "pressure": args.pressure, + "sigma": args.sigma, + "sigma_rel": args.sigma_rel, + "epsilon": args.epsilon, + "epsilon_rel": args.epsilon_rel, + "charge": args.q, + "mass": args.mass, + "r_AB": args.r_AB, + "angle_ABA": args.angle_ABA, + "r_AM": args.r_AM, + } + + + # Compute dependent values + params = compute_dependent(params) + params = clean_params(params) + + dir_name = format_directory_name(params) + if os.path.exists(dir_name): + print(f"Directory {dir_name} already exists. Aborting.") + return + os.makedirs(dir_name) + + # Save to file + config_path = os.path.join(dir_name, args.output) + with open(config_path, "w") as f: + config = configparser.ConfigParser() + config["params"] = params + config.write(f) + + print(f"Parameter file written to {config_path}") + +if __name__ == "__main__": + main() + diff --git a/generate_annealing.py b/generate_annealing.py new file mode 100755 index 0000000..75aee61 --- /dev/null +++ b/generate_annealing.py @@ -0,0 +1,97 @@ +#!/usr/bin/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. + Supports negative decimals for large values. + """ + if x == 0: + return 0.0 + + # Determine the required number of decimals (can be negative) + decimals = -int(math.floor(math.log10(rel_tol * abs(x)))) + + # Clamp within allowed range + decimals = max(-max_negative_decimals, min(max_decimals, decimals)) + + if decimals >= 0: + return round(x, decimals) + else: + # Round to nearest 10^(-decimals) + 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 + + temps_down = [] + T = 300 + while T > Tmin: + T *= r + if T < Tmin: + break + temps_down.append(T) + + inv_r = 1.0 / r + temps_up = [] + T = 300 + while T < Tmax: + T *= inv_r + if T > Tmax: + break + temps_up.append(T) + + all_T = list(reversed(temps_down)) + [300.0] + temps_up + + # Apply dynamic rounding + 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): + dT = (Tmax - Tmin) / 20 + all_T = [i*dT+Tmin 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] + duration = end_time - start_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") + parser.add_argument("-l", "--Tmin", type=float, default=10.0, help="Minimum temperature") + parser.add_argument("-u", "--Tmax", type=float, default=5000.0, help="Maximum temperature") + 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") + + args = parser.parse_args() + + #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: + temps = list(reversed(temps)) + + if args.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) + else: + times = generate_times_list(len(temps), args.start_time, args.end_time) + + 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() diff --git a/generate_temperatures.py b/generate_temperatures.py new file mode 100755 index 0000000..61fc6d9 --- /dev/null +++ b/generate_temperatures.py @@ -0,0 +1,57 @@ +#!/usr/bin/python3.12 +import numpy as np +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. + Supports negative decimals for large values. + """ + if x == 0: + return 0.0 + + # Determine the required number of decimals (can be negative) + decimals = -int(math.floor(math.log10(rel_tol * abs(x)))) + + # Clamp within allowed range + decimals = max(-max_negative_decimals, min(max_decimals, decimals)) + + if decimals >= 0: + return round(x, decimals) + else: + # Round to nearest 10^(-decimals) + 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 + + temps_down = [] + T = 300 + while T > Tmin: + T *= r + if T < Tmin: + break + temps_down.append(T) + + inv_r = 1.0 / r + temps_up = [] + T = 300 + while T < Tmax: + T *= inv_r + if T > Tmax: + break + temps_up.append(T) + + all_T = np.array(temps_down[::-1] + [300.0] + temps_up) + + # Apply dynamic rounding + rounded_T = [round_dynamic(t, rel_tol=rel_tol) for t in all_T] + rounded_T = np.unique(rounded_T) + + return rounded_T + +T = generate_temperature_list_dynamic(dT=30, rel_tol=5e-3) +print(len(T)) +print(T) +print(np.linspace(0, 100000, len(T), endpoint=True)) diff --git a/replace_params.sh b/replace_params.sh new file mode 100755 index 0000000..37d810c --- /dev/null +++ b/replace_params.sh @@ -0,0 +1,47 @@ +#!/usr/bin/env bash + +set -e + +if [[ $# -lt 2 ]]; then + echo "Usage: $0