This repository contains information and instructions on how to install and run relative alchemical transofrmations with position restraints for MLIPs code with the MACE-OFF23 potential. We need a combination of several packages. This repository only contains required information, input data, sample and analysis scripts. The actual code and used packages are stored in two other (forked) repositories:
We forked the openmm-ml and the mace github repositories and made some adjustments. You need both packages to get the entire set-up running. You will need:
https://github.com/cbc-univie/openmm-ml-anchor-4D and https://github.com/cbc-univie/mace
git clone git@github.com:cbc-univie/anchor-4D.git
git clone git@github.com:cbc-univie/openmm-ml-anchor-4D.git
git clone git@github.com:cbc-univie/mace.git
mamba create -n anchor-4D python=3.12
mamba activate anchor-4D
mamba install pytorch=2.5.1 pytorch-gpu openmm-torch cudatoolkit nnpops -c conda-forge
cd openmm-ml-anchor-4D/
pip install .
cd ../mace/
pip install .
The data folder contains pdb files for all equilibrated small solute boxes we used for ASFE calculations. The data/make_overlapping_boxes/ contains a script that constructs initial coordinates for RSFE computation from the equilibrated boxes. Details on how to use this are stored in data/make_overlapping_boxes/docu.txt. To run this script you need to install mdtraj
mamba install mdtraj -c conda-forge
The script folder contains both the sampling scripts (slurm submit and python script) and the analysis scripts (slurm submit and python scripts.) The submit script writes and runs submit scripts for all lambda states. To run this, go to the anchor-4D/scripts/ folder and execute:
bash generate_samples.sh
Note that in the generate_samples.sh file you need to specify the anchor atoms for both solutes (we recommend using the anchor atoms that were used to contruct initial coordinates, see data/make_overlapping_boxes/sepcs.yaml). Also, currently, only the linear_to_cutoff shifting scheme is implemented. So you should not change this option.
If you want to run a single job for one lambda state, you can run (e.g. for ethane to water):
python sample_states.py --lamb 0.0 --pdb ../data/ethane_water/input/ethane_water_overlapping_box.pdb --anchor_sol1 C1 --anchor_sol2 O --shifting linear_to_cutoff --model_size small --default_dtype float32 --timestep 0.001 --tag v1
Note that this only works if you have first constructed initial coordinates in data/make_overlapping_boxes/ .
In the generate_samples.sh script you can define your lambda-schedule, the shifting style, model size, precision and the time step for the simulation.
Before running the analysis, make sure the required dependencies are installed:
mamba install -c conda-forge pymbar mdtraj seaborn tqdm
The free energy calculation can be performed using either the pairwise BAR or MBAR method (in this application, MBAR can sometimes produce incorrect results and/or fail to converge; in such cases you can apply the mbar_filtered method, see calculate_asfe.py for more details)
Free energy calculations can be performed using the provided run_analysis.sh script.
Make sure to set the correct base directory path inside the script.
Example usage:
sbatch run_analysis.sh --system_name methane_ethane --version v1 --method bar --run 1
sbatch run_analysis.sh --system_name methane_ethane --version v1 --method mbar --run 1
sbatch run_analysis.sh --system_name methane_ethane --version v1 --method mbar_filtered --run 1 --start_index 6000 --lambda_range 0 4
The MBAR analysis produces an overlap matrix and a weight matrix, as well as a .pkl file containing the u_kn matrix and N_k array, which were used for the free energy estimation. The free energy result is printed to a .txt file. The BAR analysis saves .pkl files for all lambda pairs containing the u_kn matrix and N_k array. The pairwiese free energy differences are printed to a .txt file and need to be summed up to get the final free energy estimate.