Setting Up a Pd₃ Nanocluster in OPC Water for GROMACS MD Simulation

Overview

This post documents the complete pipeline for setting up a molecular dynamics simulation of a Pd₃ palladium nanocluster solvated in OPC water using GROMACS on the PARAM Shakti HPC cluster at IIT Kharagpur. The goal is a reproducible, publication-quality NPT equilibration followed by a production run, with every non-obvious decision explained and justified.

The pipeline covers:

1. Force-field file preparation and Packmol box packing with fftool
2. Correcting the raw GROMACS topology for 4-site OPC (MW virtual site)
3. Fixing the Lorentz-Berthelot mixing rule for the Pd–Ow cross term
4. Inserting the MW dummy site into the coordinate file
5. EM → NVT → NPT sequential equilibration on PARAM Shakti
6. Automated post-run sanity checks and PBC correction

System: 1 Pd₃ cluster + 7,205 OPC water molecules, 60 × 60 × 60 Å cubic box, 298.15 K, 1 bar — total 28,823 atoms.

---

1. Software Stack

Tool	Version	Purpose
GROMACS	2022.2 (gmx_CPUIMPI)	MD engine — real Intel MPI build, no thread-MPI
fftool	github.com/paduagroup/fftool	Z-matrix → GROMACS files
Packmol	—	Initial box packing
Python	3.x	MW insertion script
PARAM Shakti	IIT KGP	HPC cluster (Intel MPI, 40 cores/node)

gmx_CPUIMPI is a real Intel MPI build. The thread-MPI layer is absent. This means -ntmpi and -ntomp flags are invalid for this binary. Parallelism is controlled exclusively by mpirun -np N. This distinction affects every mdrun call in the pipeline.

Installing fftool on macOS

fftool ships with a python shebang on older commits; macOS only has python3:

# Clone your fork
git clone https://github.com/<your-fork>/fftool.git ~/Downloads/fftool

# Fix the shebang on all four scripts
for f in fftool ff2xml xyztool lattice; do
  sed -i '' '1s|.*|#!/usr/bin/env python3|' ~/Downloads/fftool/$f
done

chmod +x ~/Downloads/fftool/fftool

# Add to PATH (zsh)
echo 'export PATH="$HOME/Downloads/fftool:$PATH"' >> ~/.zshrc
source ~/.zshrc

# Verify
fftool -h

---

2. Input Files

Two z-matrix (.zmat) and two force-field (.ff) files are required.

Pd₃ force field (pd3.ff)

Palladium LJ parameters from the INTERFACE force field (Heinz 2008):

Parameter	Value
σ	2.5114 Å
ε	3.58509 kcal/mol (= 15.0039 kJ/mol)
Charge	0.0 e
Pd–Pd bond	2.676 Å (equilateral D₃h, DFT geometry)

OPC water force field (opc.ff)

OPC (Izadi & Onufriev 2014) is a 4-point water model. The negative charge does not sit on the oxygen — it lives on a massless virtual site (M-site) displaced 0.1594 Å off the oxygen along the H–O–H bisector.

Parameter	Value
q(Ow)	0.0 e
q(Hw)	+0.6791 e
q(MW)	−1.3582 e
r(OH)	0.8724 Å
∠HOH	103.6°
d(OM)	0.1594 Å

Why OPC? OPC reproduces bulk water density, diffusion coefficient, and dielectric constant simultaneously better than TIP4P or SPC/E, making it the preferred choice for solvation studies of metal nanoclusters.

---

3. Packing the Box with fftool + Packmol

Step 1 — Generate per-molecule coordinate files

# Back up your custom pack.inp first — fftool overwrites it
cp pack.inp pack.inp.mine
fftool 1 pd3.zmat 7205 opc.zmat --box 60.0
cp pack.inp.mine pack.inp   # restore custom packing logic

Flag note: this fork of fftool uses --box (long form). The short form -b is not recognised and will throw an unrecognized arguments error.

Step 2 — Custom pack.inp

The custom pack.inp places Pd₃ at the exact box centre and enforces a 5.5 Å exclusion sphere around it so no water is placed inside the first coordination shell:

tolerance 2.0
filetype xyz
output simbox.xyz

structure pd3_pack.xyz
  number 1
  center
  fixed 30.0000 30.0000 30.0000 0.0 0.0 0.0
end structure

structure opc_pack.xyz
  number 7205
  inside box 1.5000 1.5000 1.5000 58.5000 58.5000 58.5000
  outside sphere 30.0000 30.0000 30.0000 5.5
end structure

The 1.5 Å buffer on each face prevents water from overlapping with the periodic image of the opposite face.

Step 3 — Run Packmol

packmol < pack.inp
# On success: writes simbox.xyz and prints "Success!"
head -2 simbox.xyz
# Expected: 21618  (3 Pd + 7205 × 3 OPC sites)

Step 4 — Generate GROMACS files

fftool 1 pd3.zmat 7205 opc.zmat --box 60.0 --gmx
# Writes: field.top  run.mdp  config.pdb

At this stage field.top has two critical errors that must be fixed before grompp will produce a physically correct simulation.

---

4. Correcting the Topology (field.top)

fftool produces a raw topology with two problems:

1. Wrong mixing rule — comb-rule 3 (geometric σ) instead of comb-rule 2 (Lorentz-Berthelot) for the Pd–Ow cross term.
2. 3-site OPC — the negative charge is left on the oxygen atom instead of being moved to the massless MW virtual site.

The corrected field.top is reproduced in full below with every column and flag documented inline.

; field.top  --  Pd3 cluster solvated in OPC water  (GROMACS topology)
; GROMACS internal units: length=nm, energy=kJ/mol, mass=u, charge=e
; Ref: GROMACS Reference Manual > Topologies > Parameter files
;      https://manual.gromacs.org/current/reference-manual/topologies/parameter-files.html

[ defaults ]
; nbfunc   comb-rule   gen-pairs   fudgeLJ   fudgeQQ
; nbfunc 1  = Lennard-Jones
; comb-rule 2 = arithmetic σ_ij=(σ_i+σ_j)/2, geometric ε_ij=√(ε_i·ε_j)
;              → Lorentz-Berthelot, correct for INTERFACE FF + OPC (AMBER lineage)
; comb-rule 3 (fftool default) = geometric σ too → WRONG for this system
  1        2           yes         0.5       0.5

[ atomtypes ]
; name  at.num     mass    charge  ptype   sigma(nm)     epsilon(kJ/mol)
; ptype A = normal atom. charge here is default/fallback; per-atom [ atoms ]
; charges override this. All set to 0 — real charges defined per atom below.
; sigma and epsilon are V and W columns interpreted under comb-rule 2.
Pd1        46  106.4200   0.0000     A     2.51140e-01   2.57316e+01
Pd2        46  106.4200   0.0000     A     2.51140e-01   2.57316e+01
Pd3        46  106.4200   0.0000     A     2.51140e-01   2.57316e+01
Ow          8   15.9990   0.0000     A     3.16655e-01   8.90360e-01
Hw          1    1.0080   0.0000     A     1.00000e-01   0.00000e+00
; MW: massless virtual site. Zero LJ, zero mass. Carries charge per [ atoms ].
MW          0    0.0000   0.0000     A     0.00000e+00   0.00000e+00

[ moleculetype ]
; name    nrexcl
; nrexcl 3: exclude non-bonded interactions up to 3 bonds apart.
; func-1 constraints count as bonds for exclusion generation —
; all three Pd-Pd pairs are auto-excluded.
pd3       3

[ atoms ]
;  nr   type  resnr  residue  atom  cgnr     charge
    1   Pd1      1   pd3      Pd1     1     0.000000
    2   Pd2      1   pd3      Pd2     2     0.000000
    3   Pd3      1   pd3      Pd3     3     0.000000

[ constraints ]
;  ai   aj   func    b0(nm)
; func 1: holonomic constraint solved by LINCS. Counts as bond for exclusions.
; Three constraints on three atoms → fully rigid equilateral triangle (D3h).
; b0 = 0.26760 nm from pd3.ff DFT geometry.
    2    1     1     0.26760
    3    1     1     0.26760
    2    3     1     0.26760

[ exclusions ]
; Explicit for the closed triangle. With func-1 constraints + nrexcl=3
; these are auto-generated, but listed here for clarity and safety.
    1    2    3
    2    1    3
    3    1    2

[ moleculetype ]
; name    nrexcl
; nrexcl 2: SETTLE generates no bonds, so explicit [ exclusions ] below
; handles all intramolecular non-bonded interactions.
OPC       2

[ atoms ]
;  nr   type  resnr  residue  atom  cgnr     charge       mass
; KEY FIX: oxygen charge = 0. Full -1.3582 moved to MW (virtual site).
; Net molecular charge: 0 + 2(+0.6791) + (-1.3582) = 0. System neutral.
    1   Ow       1   OPC      O       1     0.000000   15.99900
    2   Hw       1   OPC      H1      1     0.679100    1.00800
    3   Hw       1   OPC      H2      1     0.679100    1.00800
    4   MW       1   OPC      MW      1    -1.358200    0.00000

[ settles ]
;  OW   func    doh(nm)    dhh(nm)
; SETTLE analytically rigidifies water using doh and dhh.
; dhh = 2·doh·sin(θ/2) = 2(0.08724)sin(51.8°) = 0.13713 nm
; This encodes the HOH angle — SETTLE fixes both bond lengths AND the angle.
    1    1      0.08724    0.13713

[ exclusions ]
; MANDATORY: SETTLE creates no bonds, so nothing auto-removes intramolecular
; non-bonded interactions. All four sites must be mutually excluded.
    1    2    3    4
    2    1    3    4
    3    1    2    4
    4    1    2    3

[ virtual_sites3 ]
; Construct MW from three real atoms using func 1:
; r_MW = r_O + a(r_H1 - r_O) + b(r_H2 - r_O)
; a = b = d_OM / (2·dOH·cos(θ/2)) = 0.01594/(2·0.08724·cos51.8°) = 0.147722
; Places MW 0.1594 Å from O on the bisector. mdrun rebuilds position each step.
; site  ai   aj   ak  func      a           b
    4    1    2    3    1      0.147722    0.147722

[ system ]
Pd3 cluster in OPC water

[ molecules ]
; Order and counts MUST match the coordinate file exactly.
pd3       1
OPC       7205

Why comb-rule 2 matters for Pd–Ow

Under comb-rule 2 (Lorentz-Berthelot):

\[\sigma_{ij} = \frac{\sigma_i + \sigma_j}{2}, \quad \epsilon_{ij} = \sqrt{\epsilon_i \cdot \epsilon_j}\]

Under comb-rule 3 (OPLS geometric), both σ and ε are geometric means:

\[\sigma_{ij} = \sqrt{\sigma_i \cdot \sigma_j}\]

For Pd (σ = 0.25114 nm) and Ow (σ = 0.31666 nm):

Rule	σ(Pd–Ow) nm
comb-rule 2 (LB, correct)	0.28390
comb-rule 3 (geometric)	0.28217

A 0.0017 nm shift in σ changes the Pd–Ow potential minimum depth and position, directly altering the first solvation shell structure — the primary observable in this study.

---

5. Inserting the OPC MW Virtual Site

fftool writes config.pdb with 3 sites per water (21,618 atoms). The corrected topology expects 4 sites per water (28,823 atoms). The following Python script inserts one MW placeholder atom after each water’s H2, keyed on residue name OPC so it is reusable for any solute:

#!/usr/bin/env python3
"""
add_mw.py -- insert OPC MW virtual-site into a packed PDB.
Usage: python3 add_mw.py config.pdb config_mw.pdb
"""
import sys

WATER_RES   = "OPC"
MW_NAME     = "MW"
ATOMS_PER_W = 3

inp = sys.argv[1] if len(sys.argv) > 1 else "config.pdb"
out = sys.argv[2] if len(sys.argv) > 2 else "config_mw.pdb"

def is_atom(l):   return l.startswith("ATOM") or l.startswith("HETATM")
def resname(l):   return l[17:20].strip()
def make_mw(o):   return o[:12] + (" %-3s" % MW_NAME) + o[16:]

out_lines = []; buf = []; n_water = 0; n_other = 0

def flush_water():
    global n_water
    if not buf: return
    if len(buf) != ATOMS_PER_W:
        sys.exit(f"ERROR: OPC residue had {len(buf)} atoms, expected {ATOMS_PER_W}")
    out_lines.extend(buf)
    out_lines.append(make_mw(buf[0]))   # MW placeholder at oxygen position
    buf.clear(); n_water += 1

with open(inp) as f:
    for raw in f:
        line = raw.rstrip("\r\n")
        if is_atom(line) and resname(line) == WATER_RES:
            buf.append(line)
            if len(buf) == ATOMS_PER_W: flush_water()
        else:
            flush_water()
            if is_atom(line): n_other += 1
            out_lines.append(line)
flush_water()

serial = 0
for i, l in enumerate(out_lines):
    if is_atom(l):
        serial += 1
        out_lines[i] = l[:6] + ("%5d" % (serial % 100000)) + l[11:]

with open(out, "w") as f:
    f.write("\n".join(out_lines) + "\n")

print(f"water residues     : {n_water}   (one MW added to each)")
print(f"non-water atoms    : {n_other}")
print(f"total atoms written: {serial}   (expected {n_other + n_water*4})")

Run and verify:

python3 add_mw.py config.pdb config_mw.pdb
# Expected:
# water residues     : 7205   (one MW added to each)
# non-water atoms    : 3
# total atoms written: 28823  (expected 28823)

# Spot-check first water — must be O, H1, H2, MW in that order
grep -m4 ' OPC ' config_mw.pdb

MW coordinate = oxygen position. grompp/mdrun reconstruct MW’s real position every step from the [ virtual_sites3 ] rule. The placeholder coordinate is never used in any force calculation.

---

6. Why Pd₃ Behaves as a Single Rigid Entity

The three [ constraints ] entries fully rigidify the cluster:

[ constraints ]
    2    1     1     0.26760   ; Pd2–Pd1  fixed at 0.2676 nm
    3    1     1     0.26760   ; Pd3–Pd1  fixed at 0.2676 nm
    2    3     1     0.26760   ; Pd2–Pd3  fixed at 0.2676 nm

Three distance constraints on three atoms uniquely fix the shape of the triangle. No angle term is needed — once all three sides are fixed the angle is determined implicitly. LINCS enforces these at every step.

Degree of freedom	Status
Pd–Pd bond lengths	Frozen at 0.2676 nm (LINCS)
Pd–Pd–Pd angles	Frozen at 60° (implicit from 3 fixed sides)
Internal vibration	Zero (no internal DOF remain)
Cluster translation	Free — COM moves under net solvent force
Cluster rotation	Free — torque from solvent rotates the body

The cluster diffuses and tumbles in solution exactly as a real nanocluster would. What is frozen is only its internal geometry, not its position.

---

7. Index File

The thermostat couples Pd₃ and OPC separately to prevent the “hot solute / cold solvent” artefact:

gmx_CPUIMPI make_ndx -f config_mw.pdb -o index.ndx
# At the prompt: q

Verify:

grep '^\[' index.ndx
# Expected:
# [ System ]
# [ Other ]
# [ pd3 ]
# [ OPC ]

---

8. MDP Files — Equilibration Protocol

Three sequential phases. Key design decisions are consistent across all:

Property	EM	NVT	NPT
integrator	steep	md	md
pcoupl	no	no	Parrinello-Rahman
tcoupl	N/A	V-rescale	V-rescale
gen-vel	N/A	yes	no
continuation	no	no	yes
DispCorr	No	EnerPres	EnerPres
constraints	all-bonds	all-bonds	all-bonds

8.1 Energy Minimisation (01_em.mdp)

integrator            = steep
emtol                 = 1000.0        ; convergence: Fmax < 1000 kJ/mol/nm
emstep                = 0.01          ; initial step size [nm]
nsteps                = 50000         ; max steps (early termination expected)

cutoff-scheme         = Verlet
pbc                   = xyz
coulombtype           = PME           ; required for discrete MW charge
rcoulomb              = 1.2
ewald-rtol            = 1.0e-5
pme-order             = 4
fourierspacing        = 0.12

vdwtype               = Cut-off
vdw-modifier          = Potential-shift
rvdw                  = 1.2
DispCorr              = No            ; off: packed box density not yet converged

constraints           = all-bonds     ; Pd-Pd via LINCS + SETTLE water
constraint-algorithm  = LINCS
lincs-order           = 6             ; higher order for closed Pd3 triangle
lincs-iter            = 2

DispCorr = No during EM is intentional. The long-range dispersion tail correction assumes a converged homogeneous liquid density. A freshly Packmol-packed box has voids and inhomogeneous density — applying the correction here would give a nonsense pressure value. It is switched on from NVT onward.

8.2 NVT Equilibration (02_nvt.mdp)

integrator            = md
dt                    = 0.002         ; 2 fs — safe with all-bonds constraints
nsteps                = 500000        ; 1 ns

coulombtype           = PME
rcoulomb              = 1.2
ewald-rtol            = 1.0e-5
pme-order             = 4
fourierspacing        = 0.12

vdwtype               = Cut-off
vdw-modifier          = Potential-shift
rvdw                  = 1.2
DispCorr              = EnerPres

; V-rescale correctly samples the canonical ensemble (Bussi et al. 2007).
; Separate groups prevent hot-solute/cold-solvent artefacts.
tcoupl                = V-rescale
tc-grps               = pd3  OPC
tau-t                 = 2.0  2.0
ref-t                 = 298.15  298.15

pcoupl                = no            ; fixed volume: equilibrate T before P

gen-vel               = yes           ; fresh Maxwell-Boltzmann start after EM
gen-temp              = 298.15
gen-seed              = -1

constraints           = all-bonds
constraint-algorithm  = LINCS
lincs-order           = 6
lincs-iter            = 2
continuation          = no

8.3 NPT Equilibration (03_npt.mdp)

integrator            = md
dt                    = 0.002
nsteps                = 500000        ; 1 ns

coulombtype           = PME
rcoulomb              = 1.2
ewald-rtol            = 1.0e-5
pme-order             = 4
fourierspacing        = 0.12

vdwtype               = Cut-off
vdw-modifier          = Potential-shift
rvdw                  = 1.2
DispCorr              = EnerPres

tcoupl                = V-rescale
tc-grps               = pd3  OPC
tau-t                 = 2.0  2.0
ref-t                 = 298.15  298.15

; Parrinello-Rahman produces a correct NPT ensemble.
; The GROMACS manual explicitly states Berendsen does NOT.
pcoupl                = Parrinello-Rahman
pcoupltype            = isotropic     ; cubic box stays cubic
tau-p                 = 2.0           ; [ps] >> tau-t for P-R stability
ref-p                 = 1.0           ; [bar]
compressibility       = 4.5e-5        ; [bar-1] experimental value for water 298 K

gen-vel               = no            ; velocities inherited from nvt.cpt
continuation          = yes           ; signals continuation to grompp

constraints           = all-bonds
constraint-algorithm  = LINCS
lincs-order           = 6
lincs-iter            = 2

Critical: the NPT grompp call must include -t nvt.cpt to inherit the thermally equilibrated velocities from NVT. Without this the NPT run starts from zero velocities, discarding the NVT equilibration entirely.

---

9. SLURM Submission on PARAM Shakti

All three phases run inside a single SLURM job to avoid queue re-entry delays and to ensure each phase inherits the exact checkpoint from the previous one.

Important: gmx_CPUIMPI parallelism model

gmx_CPUIMPI is compiled with real Intel MPI only — the thread-MPI layer is completely absent. This has one critical consequence:

• -ntmpi and -ntomp flags are invalid and will crash every MPI rank with the error: “Setting the number of thread-MPI ranks is only supported with thread-MPI”
• Parallelism is controlled exclusively by mpirun -np 40
• All three phases (EM, NVT, NPT) use the same simple mpirun -np 40 launch
• GROMACS domain decomposition distributes work across all 40 ranks automatically

#!/bin/bash
#SBATCH --partition=medium          # 1-10 full nodes, max 3 days
#SBATCH --job-name=pd3_equil
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=40        # MANDATORY for medium partition
#SBATCH --time=00-12:00:00
#SBATCH --output=equil_%j.out
#SBATCH --error=equil_%j.err

set -e
set -x

WORKDIR="/scratch/20ch91r15/GROMACS/single_pd3_water/equilibration"
cd "${WORKDIR}"

source /home/apps/gromacs/gromacs-2022.2/installCPUIMPI/bin/GMXRC.bash

# Source GMXRC before set -u to avoid unbound variable errors in GMXRC itself
MPIRUN="mpirun -np 40"
GMX="gmx_CPUIMPI"

# Verify all inputs exist before consuming queue allocation
for f in field.top config_mw.pdb index.ndx 01_em.mdp 02_nvt.mdp 03_npt.mdp; do
    [[ -f "$f" ]] || { echo "ERROR: missing $f"; exit 1; }
done

# ── Phase 1: Energy Minimisation ──────────────────────────────────────────
${GMX} grompp -f 01_em.mdp -c config_mw.pdb -p field.top \
              -n index.ndx -o em.tpr -maxwarn 1

${MPIRUN} ${GMX} mdrun -v -deffnm em

# ── Phase 2: NVT Equilibration (1 ns) ─────────────────────────────────────
${GMX} grompp -f 02_nvt.mdp -c em.gro -p field.top \
              -n index.ndx -o nvt.tpr -maxwarn 1

${MPIRUN} ${GMX} mdrun -v -deffnm nvt

# ── Phase 3: NPT Equilibration (1 ns) ─────────────────────────────────────
# -t nvt.cpt inherits NVT velocities — DO NOT omit this flag
${GMX} grompp -f 03_npt.mdp -c nvt.gro -t nvt.cpt \
              -p field.top -n index.ndx -o npt.tpr -maxwarn 1

${MPIRUN} ${GMX} mdrun -v -deffnm npt

Submit:

cd /scratch/20ch91r15/GROMACS/single_pd3_water/equilibration
sbatch equilibration.sh
squeue -u 20ch91r15

PARAM Shakti filesystem policy: all job I/O must be on /scratch/$USER. The /home path is for source code and configs only and is not on the fast parallel filesystem. Never run mdrun from /home.

---

10. Post-Run Sanity Checks

After the job finishes:

bash check_equilibration.sh 2>&1 | tee sanity_report.txt

Shell compatibility note: source GMXRC.bash before set -euo pipefail in the check script. GMXRC internally references variables that may be unbound on a clean login shell, which causes -u to abort immediately.

The script checks seven things and prints PASS / WARN / FAIL for each:

Check	What it verifies	Pass criterion
File existence	All 9 output files present	All found
Atom count	.gro header = 28,823	Exact match
EM convergence	Maximum force in em.log	Fmax < 1000 kJ/mol/nm
NVT temperature	Average last 50% of trajectory	298.15 ± 5 K
NPT density	Average last 50% of trajectory	997 ± 10 kg/m³
NPT pressure	Average last 50% of trajectory	1 ± 300 bar
Pd–Pd constraint	Distance over NPT trajectory	0.2676 ± 0.001 nm

PBC correction

The script applies the two-step PBC correction automatically:

# Step A: make all molecules whole (un-wrap across box faces)
gmx trjconv -f npt.xtc -s npt.tpr -n index.ndx \
            -o npt_pbc.xtc -pbc mol -ur compact

# Step B: centre Pd3 in the box
echo -e "2\n0" | gmx trjconv -f npt_pbc.xtc -s npt.tpr -n index.ndx \
                              -o npt_center.xtc -pbc mol -center

Step order is critical. Centring before making molecules whole can re-break molecules that were just fixed. Always make whole first, then centre.

npt_center.xtc is the trajectory to use for all structural analysis (RDF, MSD, coordination number, radial density profiles).

Expected convergence

Observable	Expected	Normal fluctuation
Temperature	298.15 K	± 2–3 K
Pressure	1 bar	± 100–300 bar
Density	~997 kg/m³	± 2–3 kg/m³
Pd–Pd distance	0.2676 nm	< 0.001 nm

If density is still drifting at the end of 1 ns NPT, extend by another 1 ns before starting production.

---

11. Proceeding to Production

Once sanity_report.txt shows zero FAILs:

gmx_CPUIMPI grompp \
    -f run.mdp \
    -c npt_center.gro \
    -t npt.cpt \
    -p field.top \
    -n index.ndx \
    -o production.tpr \
    -maxwarn 1

mpirun -np 40 gmx_CPUIMPI mdrun -v -deffnm production

Key differences in run.mdp versus the NPT mdp:

• nsteps = 50000000 (100 ns at 2 fs/step)
• continuation = yes, gen-vel = no
• nstxout-compressed = 5000 (every 10 ps → 10,000 frames for analysis)

---

References

1. Izadi, S. & Onufriev, A. V. (2016). Accuracy limit of rigid 3-point water models. J. Chem. Phys. 145, 074501.
2. Heinz, H. et al. (2008). Accurate Simulation of Surfaces and Interfaces of Face-Centered Cubic Metals Using 12-6 and 9-6 Lennard-Jones Potentials. J. Phys. Chem. C 112, 17281.
3. Bussi, G., Donadio, D. & Parrinello, M. (2007). Canonical sampling through velocity rescaling. J. Chem. Phys. 126, 014101.
4. Páduai, A. A. H. fftool. https://github.com/paduagroup/fftool
5. GROMACS Reference Manual 2022.2. https://manual.gromacs.org/2022.2/reference-manual/