initial commit

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"
+
+

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
+
+

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()
+

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()

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))

replace_params.sh

@@ -0,0 +1,47 @@
+#!/usr/bin/env bash
+
+set -e
+
+if [[ $# -lt 2 ]]; then
+    echo "Usage: $0 <params.ini> <template> [--output <output_file>]"
+    exit 1
+fi
+
+PARAMS_FILE="$1"
+TEMPLATE_FILE="$2"
+OUTPUT_FILE="${TEMPLATE_FILE}.out"
+
+# Optional --output <filename>
+if [[ "$3" == "--output" && -n "$4" ]]; then
+    OUTPUT_FILE="$4"
+fi
+
+declare -A PARAMS
+
+# Load params.ini (ignore [section] headers)
+while IFS='=' read -r key value; do
+    # Trim whitespace and ignore comments/sections
+    key="${key%%\#*}"   # remove inline comments
+    value="${value%%\#*}"
+    key="$(echo "$key" | xargs)"     # trim leading/trailing whitespace
+    value="$(echo "$value" | xargs)"
+    [[ -z "$key" || "$key" =~ ^\[.*\]$ ]] && continue
+
+    key_upper=$(echo "$key" | tr '[:lower:]' '[:upper:]')
+    PARAMS["$key_upper"]="$value"
+done < "$PARAMS_FILE"
+
+# Read template into a variable
+template=$(<"$TEMPLATE_FILE")
+
+# Replace placeholders (only when surrounded by spaces)
+for key in "${!PARAMS[@]}"; do
+    val="${PARAMS[$key]}"
+    # Use sed to do space-surrounded substitution
+    template=$(echo "$template" | sed -E "s/([[:space:]])$key([[:space:]])/$val/g")
+done
+
+# Save to output file
+echo "$template" > "$OUTPUT_FILE"
+echo "Wrote output to $OUTPUT_FILE"
+

run.sh

@@ -0,0 +1,53 @@
+#!/bin/bash -l
+
+#SBATCH --partition=normal
+#SBATCH --ntasks=1
+#SBATCH --cpus-per-task=8
+#SBATCH --gres=gpu:1
+#SBATCH --exclude=linux-05,linux-08,linux-02,linux-03,linux-04,linux-07
+
+#if [ "$(gmx --version | grep "GROMACS version")" != GROMACS*2025.2* ] ; then
+#    echo "Version loaded is not 2025.2! Exiting!"
+#    exit 1
+#fi
+
+if [ -n "$SLURM_CPUS_PER_TASK" ]; then
+    NT="-nt $SLURM_CPUS_PER_TASK"
+elif [ $(nproc) -lt 9 ] ; then
+    NT="-nt $(nproc)"
+else
+    NT="-nt 8"
+fi
+
+WORKDIR="$(dirname $(readlink -f "$0"))"
+if [ "${WORKDIR%/*}" = "/var/lib/slurm/slurmd" ] ; then
+    WORKDIR="$(pwd)"
+fi
+
+mkdir -p "$WORKDIR/out"
+
+export GMX_MAXCONSTRWARN=-1;
+gmx -nobackup grompp \
+    -f $WORKDIR/mdp_parameters.mdp \
+    -o $WORKDIR/tpr_run.tpr \
+    -c $WORKDIR/gro_start.gro \
+    -r $WORKDIR/gro_start.gro \
+    -p $WORKDIR/topology.top \
+    -po $WORKDIR/out/mdp_mdout.mdp \
+    -maxwarn 6
+
+if [ $? != 0 ] ; then
+    echo "grompp failed, exiting"
+    exit 1
+fi
+
+gmx mdrun \
+    -s $WORKDIR/tpr_run.tpr \
+    -o $WORKDIR/out/trr_traj.trr \
+    -c $WORKDIR/out/out.gro \
+    -x $WORKDIR/out/xtc_traj.xtc \
+    -e $WORKDIR/out/edr_energy.edr \
+    -g $WORKDIR/out/log.log \
+    -cpo $WORKDIR/out/state.cpt \
+    -cpi $WORKDIR/out/state.cpt \
+    -cpt 1 -notunepme -v $NT

setup_sim.sh

@@ -0,0 +1,147 @@
+#!/bin/bash -l
+
+usage() { echo "Usage: $0 [-d basedirectory ]" 1>&2; exit 1; }
+
+DIRECTORY="."
+
+while getopts "d:" opt; do
+  case $opt in
+    d )
+        DIRECTORY=$OPTARG
+        ;;
+    h )
+        usage
+        ;;
+   \? )
+        usage
+        ;;
+  esac
+done
+
+cd $DIRECTORY
+
+#=====================================================================================
+cat > mdp_parameters.mdp << EOF
+integrator               =  INTEGRATOR 
+dt                       = 0.001
+nsteps                   =  NSTEPS 
+nstcomm                  = 10
+
+nstlog                   =  NLOG 
+nstcalcenergy            = 10
+nstenergy                =  NENERGY 
+nstxout-compressed       =  NSTXOUT 
+compressed-x-precision   = 1000
+energygrps               = system
+
+rlist                    = 1.2
+nstlist                  = 10
+
+coulombtype              = pme
+rcoulomb                 = 1.2
+coulomb-modifier         = potential-shift-verlet
+vdwtype                  = pme
+rvdw                     = 1.2
+vdw-modifier             = potential-shift-verlet
+fourierspacing           = 0.144
+
+tcoupl                   = v-rescale
+tau_t                    =  TAUT 
+gen_vel                  = yes
+gen_temp                 =  TEMP 
+ref_t                    =  TEMP 
+tc-grps                  = system
+
+pcoupl                   =  PCOUPL 
+tau_p                    =  TAUP 
+ref_p                    =  PRESSURE 
+compressibility          = 4.5e-5
+
+nstpcouple               = 10
+nsttcouple               = 10
+
+EOF
+
+#=====================================================================================
+cat > params_emin.ini << EOF
+INTEGRATOR = steep
+NSTEPS = 10000
+NLOG = 100
+NENERGY = 100
+NSTXOUT = 100
+TAUT = 1.0
+TEMP = 300
+PCOUPL = no
+TCOUPL = no
+TAUP = 1.0
+EOF
+
+cat > params_an1.ini << EOF
+INTEGRATOR = md
+NSTEPS = 2000000
+NLOG = 1000
+NENERGY = 100
+NSTXOUT = 1000
+TCOUPL = v-rescale
+TAUT = 0.5
+TEMP = 10.0
+PCOUPL = c-rescale
+TAUP = 1.0
+EOF
+
+#=====================================================================================
+
+cat > topology.top << EOF
+[ defaults ]
+; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
+  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 
+D     0             0.000       D   0.0           0.0
+
+[moleculetype]
+; name nrexcl
+L  1
+
+[atoms]
+; nr type resnr residu atom cgnr charge
+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 
+4     D   1     L  D    1    CHARGE_A      0.0
+
+[settles]
+;i      funct   doh     dhh
+1       1       D_AB     D_BB 
+
+
+[dummies3]
+; Dummy from                            funct       a        b
+4       1       2       3       1       DUMMY_AB      DUMMY_AB 
+
+[exclusions]
+1       2       3       4
+2       1       3       4
+3       1       2       4
+4       1       2       3
+
+[system]
+ MODEL 
+
+[molecules]
+L                 NMOL 
+EOF
+
+#=====================================================================================
+bash ../create_emin.sh
+bash ../create_annealing1.sh
+
+
+#=====================================================================================
+
+echo "finished and simulation directory $DIRECTORY created and filled"
+ls
+exit 0

templates/single.gro

@@ -0,0 +1,7 @@
+SOL tip4p/2005 t= 8254.00000 step= 4127000
+    4
+    1SOL     OW    1   1.982   1.987   2.025
+    1SOL    HW1    2   2.020   2.075   2.021
+    1SOL    HW2    3   2.006   1.947   1.941
+    1SOL     MW    4   1.990   1.993   2.013
+   4.00000   4.00000   4.00000