Transferable Force Fields for Six Common Solvents from Liquid Density and Surface Tension Data: A Gaussian Process Optimization Framework

Authors: Montana N. Carlozo, Alexander W. Dowling, and Edward J. Maginn

Introduction

ES-FFO is a repository used to rapidly calibrate the LJ parameters of electrolyte solvent forcefields given experimental data. The key feature of this work is using machine learning tools in the form of Gaussian processes (GPs) which allow us to cheaply estimate the results of a molecular simulation given temperature state points and thermophysical property data. This workflow combines the approaches of Wang et al., 2023, J. Chem. Theory Comput. and Carlozo et al., 2025, Digital Discovery to create one repository to build GP models and optimize parameters for electrolyte solvent FFs. ES-FFO-Public is based on the "for_sub" branch of private repository ES-FFO.

NOTE: We use Signac and signac flow (<https://signac.io/>) to manage the setup and execution of each workflow. These instructions assume a working knowledge of that software. The most complete version of this work is the for_pub branch.

Citation

This work is currently in preparation and has not yet been published. Please cite the following when referring to this work:

Carlozo, M. N., Dowling, A. W., & Maginn, E. J.. Transferable Force Fields for Six Common Solvents from Liquid Density and Surface Tension Data: A Gaussian Process Optimization Framework. In preparation.

Available Data

Repository Organization

The repository is organized as follows:
ES-FFO/ is the top level directory. It contains:

1. .gitignore prevents large files from the signac workflow and plots from being tracked by git and prevents tracking of other unimportant files.
   
2. utils/ contains functions used by all pieces of this workflow including molecule experimental data and reference parameters
   
3. block_averages/ contains a code for block averaging via the methods in H. Flyvbjerg and H.G. Peterson. Error estimates on averages of correlated data. J. Chem. Phys. 91:461-466, 1989.
   
4. fffit/ contains another set of utility functions used by multiple aspects of this workflow.
   
5. Build_GPs/ contains the workflow runs and analyzes results from applying the methods modified from those in Wang et al. 2023 to create Base FFs
   
6. Opt_ES/ contains the workflow runs and analyzes results from applying the methods in Carlozo et al. 2025. to create GP-optimized FFs
   
7. Opt_ES/ Also contains the molecular simulation results and analysis of the FF(s) developed in steps 5 and 6.
   
8. hfcs-fffit.yml is the environment for running this workflow.
   
9. submit_jobs is a shell script for submitting jobs to the cluster.
   
10. create_analysis_figs.ipynb is a jupyter notebook for creating the final analysis figures including all figures from running the scripts in 11-17.
    
11. make_corr_figs_all_molec.py is a python script for creating correlation figures for all molecules. Figure 12, Figure 14, and SI Figure S12 in the main text.
    
12. make_corr_figs_one_molec.py is a python script for creating correlation figures for individual molecules. SI Figures S1-S6 and SI Figures S13-S24 in the main text.
    
13. make_GP_vs_sim_and_sens.py is a python script for creating figures comparing GP and actual predictions and generating the results for individual molecule sensitivity analyses. Figure 9, SI Figure S11, and SI Tables S4-S15 in the main text.
    
14. make_ift_val_figs.py is a python script for creating visualizations of the interfactial tension boxes for all ST simulations.
    
15. make_param_comp_figs.py is a python script for creating parameter comparison figures. Figures 10-11 and SI Figures S7-S10 in the main text.
    
16. make_pareto_comp_figs.py is a python script for visualizing differences between pareto-optimal LJ parameter sets in ST (and LD) iterations. Figure 13 in the main text.
    
17. make_svd_sensitivity.py is a python script for creating visualizations of the SVD sensitivity analysis. SI Tables S4-S15 in the main text.
    
18. Claude_Audit/ contains the files generated by Claude to audit the experimental data is esolvs.py and generate the data for Tables 4, 8, and 9 in the main text.
    Note that all workspace/ subdirectories below include one example signac workspace with some example input data and key simulation validation figures where appropriate.

utils/

This directory contains:

1. molec_class_files/esolvs.py: A .py file which loads all electrolyte solvent data into a class
2. prep_ms_data.py: Contains scripts used to prepare data for GP training and get error data

fffit/fffit/

This directory contains functions used by both workflows:

1. models.py: Functions related to building GP models
2. pareto.py: Functions for determining pareto efficient points
3. plot.py: Functions for plotting data
4. signac.py: Functions for extracting signac data
5. utils.py: Functions for scaling and unscaling data and shuffling/splitting GP data

Build_GPs/

This directory contains all data related to running the GP building workflow, including:

1. analysis/mol/ld_iters: Results for liquid density (LD) iterations for each molecule "mol"
2. analysis/mol/vle_iters: Results for surface tension (ST) iterations
3. utils/*.py: Functions for analyzing the results of LD and ST iterations
4. ld_iters/ : The signac project directory for running the LD iterations including project.py, templates/, and workspace/
5. vle_iters/ : The signac project directory for running the ST iterations including project.py, templates/, and workspace/
6. gen_lhs_samples.py: A script for generating the LHS data necessary for LD and ST iterations
7. init_LD.py: A script to generate the signac project for LD iterations
8. init_IFT.py: A script to generate the signac project for ST iterations
9. post_ld.py: A script to analyze LD iteration data
10. post_vle.py: A script to analyze ST iteration data
11. get_CPU_hrs.py : A script to calculate the total CPU hours used by the workflow

Opt_ES/

This directory contains all data related to running the optimization workflow of Carlozo et al., 2025, Digital Discovery including:

1. analysis/: Data from analyzing the results of the optimization workflow
2. utilsOpt/: Functions for running the optimization workflow
3. gemc_val_opt/: The signac project directory for running the gemc validation for GP-optimized FF including project.py, templates/, and workspace/
4. gemc_val_no_opt/: The signac project directory for running the gemc validation for the Base FF including project.py, templates/, and workspace/
5. ift_val_opt/: The signac project directory for running the interfacial tension validation for GP-optimized FF including project.py, templates/, and workspace/
6. ift_val_no_opt/: The signac project directory for running the interfacial tension validation for the Base FF including project.py, templates/, and workspace/
7. opt_at_params_new/: The signac project directory for running the optimization of FF parameters with updated GP models including project.py, templates/, and workspace/
8. FF_params/: The directory for compiling FF parameters for each molecule for the Base and GP-optimized FFs including ff.xml, species1.mcf, species1.pdb, and species1.top files for each molecule
9. compress_gemc.py: A script to compress the gemc validation data
10. compress_ift.py: A script to compress the interfacial tension validation data
11. get_CPU_hrs.py: A script to calculate the total CPU hours used by the optimization and validation workflows
12. get_GPU_hrs.py: A script to calculate the total GPU hours used by the optimization and validation workflows
13. init_gemc_val.py: A script to initialize the gemc validation workflow
14. init_ift_val.py: A script to initialize the interfacial tension validation workflow
15. init_opt_at.py: A script to initialize the optimization workflow
16. post_analysis_opt.py: A script to analyze the results of the optimization workflow
17. post_analysis_val.py: A script to analyze the results of the validation workflows. Generates Figures 3-8 of the main text
18. get_RD2_data.py: A script to extract select RD2 property data for each molecule for the Base and GP-optimized FFs. Generates data for Table 5 in the main text.
19. all_properties.csv: A csv file containing select RD2 property data for each molecule for the Base and GP-optimized FFs. Table 5 in the main text.
20. compile_ff_params: Python file to compile FF parameters for each molecule for the Base and GP-optimized FFs.

Opt_ES/analysis/

This directory contains data from analyzing the results of the optimization workflow for each molecule "mol".:

1. at_00/mol/ExpVal/opt_res/jac_approx: Includes jacobian approximations
2. at_00/mol/ExpVal/opt_res/hess_approx: Includes hessian approximations
3. at_00/mol/ExpVal/opt_res/MAPD: Includes Mean Absolute Percentage Deviation values
4. at_00/mol/ExpVal/opt_res/ms_val_opt: Includes the error data and validation results for the GP-Optimized FF
5. at_00/mol/ExpVal/opt_res/ms_val_no_opt: Includes the error data and validation results for the Base FF
6. at_00/mol/ExpVal/opt_res/prop_pred/: Includes comparisons of GP, simulated, and experimental values for the GP-optimized and Base FFs
7. at_00/mol/ExpVal/opt_res/best_per_run.csv: The best FF optimization results for each optimization run
8. at_00/mol/ExpVal/opt_res/pareto_info.csv: The pareto-optimal FF optimization results for each molecule
9. at_00/mol/ExpVal/opt_res/unique_best_set.csv: The unique LJ parameter sets identified for each molecule during optimization
10. AT-0/ms_val_no_opt/: Base FF property prediction plots for heat of vaporization and vapor pressure (h_p_vap.pdf), surface tension (surf_tens.pdf), and density (vle.pdf)
11. AT-0/ms_val_opt/: GP-optimized FF property prediction plots for heat of vaporization and vapor pressure (h_p_vap.pdf), surface tension (surf_tens.pdf), and density (vle.pdf)
12. AT-0/ms_val_opt_comp/: GP-optimized FF and Base FF property prediction plots. (h_p_vap.pdf, surf_tens.pdf, vle.pdf)
13. comp_err_data.csv: Data comparing the errors of different FFs
14. lit_err_data.csv: Data comparing the errors of the literature FFs
15. lit_ff_data.csv: Data comparing the properties of the literature FFs
16. lit_ff_data_w_no_opt.csv: Data comparing the properties of the literature FFs with the Base FFs
17. lit_ff_data_w_opt.csv: Data comparing the properties of the literature FFs with the GP-optimized FFs
18. lit_Hvap_est_w_no_opt.csv: Data comparing the heat of vaporization of the literature FFs and Base FFs
19. lit_Hvap_est_w_opt.csv: Data comparing the heat of vaporization of the literature FFs and GP-optimized FFs

Opt_ES/utilsOpt/

This directory contains all functions related to running the optimization workflow of Carlozo et al., 2025, Digital Discovery including:

1. utilsOpt/atom_type.py: Class definitions for atom typing schemes
2. utilsOpt/opt_atom_types.py: Other functions necessary for optimizing FF parameters and analyzing workflow results
3. utilsOpt/plot.py: Plotting functions for analyzing optimization results
4. utilsOpt/signac.py: Functions for extracting signac data

Claude_Audit/

This directory contains all functions related to the Claude Code Audit of thi repository including:

1. scripts/: Includes relevant scripts generated by Claude Code
2. refprop_output/: Includes ethyleneglycol_refprop_montana.txt and methanol_refprop_montana.txt, the REFPROP property data for EG and MeOH
3. figures/: Figures (baseline_vs_exp_*.png/.pdf) comparing experimental data is esolvs.py to other data sources for each molecule 4.audit_data/: All data necessary for the files in scripts/ and all data generated from the files in scripts/
4. refs.bib: References used in the main text

Claude_Audit/audit_data

This directory contains all functions related to the Claude Code Audit of thi repository including:

1. baseline_agreement_stats.py: Comparison of experimental data sources
2. digitized_litff.csv: The digitized values for the Vahid & Maginn, Jahn et al., and Gonzalez-Salgado & Vega sources
3. litexp_*.csv: The literature experimental data from various FF development sources
4. litff_*.csv: The literature FF data from various FF development sources
5. master_ff_predictions.csv: Predictions of property data for all investigated FFs
6. montana_reference_values.csv: The values from esolvs.py laid out in csv format
7. provenance_map.csv: Data from the sources of experimental data used in esolvs.py
8. provenance_table.tex: Table 4 in the main text
9. recomputed_metrics.csv: The R1 and R2 metrics used in Tables 8 and 9 and SI Tables S2 and S3
10. reference_values_all.csv: Compiled reference values used for rulers R1 and R2
11. refprop_validation.csv: Comparison of esolvs.py and REFPROP values for EG and MeOH
12. rename_map.csv: Information related to data sources
13. student_ff_errors_repo.csv: Compiled error data from ES-FFO/Opt_ES/analysis/
14. student_ff_recompute: Error data generated by the Claude Audit. The data in Tables 8 and 9 and SI Tables S2 and S3
15. tables_6_7_reconstructed.tex: Tables 8 and 9 in the main text
16. audit_table3_refdata.csv: Claude audit for Table 3 in the main text

Claude_Audit/scripts

This directory contains all functions related to the Claude Code Audit of thi repository including:

1. build_master_ff_predictions.py: Creates a consolidated list of all FF predictions
2. build_reference_values_all.py: Generated the R1 and R2 metric reference values used in Tables 8 and 9 and SI Tables S2 and S3
3. extract_student_ff_values.py: Extracts FF ms data from ES-FFO/Opt_ES/analysis/
4. generate_provenance_table.py: Generates Table 4 in the main text
5. generate_tables.py: Generates Tables 8 and 9 and SI Tables S2 and S3 in the main text
6. plot_baseline_vs_experiment.py: Generates experimental data comparison figures
7. recompute_metrics_all.py: Recomputes error metrics for predictions. Generates the data in Tables 8 and 9 and SI Tables S2 and S3 in the main text
8. recompute_student_metrics.py: Recomputes error metrics for predictions from the Base and GP-Optimized FFs with data from ES-FFO/Opt_ES/analysis/
9. audit_table3_refdata.py: Generates audit for Table 3 in the main text
10. unit_value_audit.py: Cross-checks esolvs.py's expt_* dictionaries (the manuscript's reference values) against ruler R1
11. validate_refprop.py: Generates comparison of esolvs.py and REFPROP values for EG and MeOH

Installation

To run this software, you must have access to all packages in the hfcs-fffit environment (hfcs-fffit.yml) which can be installed using the instructions in the next section.

This package has a number of requirements that can be installed in different ways. We recommend using a conda environment to manage most of the installation and dependencies. However, some items will need to be installed from source or pip.

Running the simulations will also require an installation of GROMACS. This can be installed separately (see installation instructions here <https://manual.gromacs.org/documentation/2021.2/install-guide/index.html>).

Usage

Liquid Density Optimization (LD Iterations)

To run liquid density iterations, follow the following steps:

1. Initialize Signac workflow

     conda activate hfcs-fffit
     cd Build_GPs/
     python init_LD.py
   

2. Run LD iterations
   • Note: rm -r workspace/ signac_project_document.json signac.rc will remove everything and allow you to start fresh if you mess up

   cd Build_GPs/ld_iters/
   python project.py submit -o LD --bundle=12 --parallel

3. Check status a few times throughout the process

     python project.py status
   

Note: Step 2 operation LD runs multiple operations in series. Alternatively, the following can be run one at a time: 4. Alternative submission

  python project.py run -o create_forcefield
  python project.py submit -o create_system
  python project.py run -o fix_topology
  python project.py submit -o em_sim --bundle=12 --parallel
  python project.py submit -o nvt_eq_sim --bundle=12 --parallel
  python project.py submit -o npt_eq_sim --bundle=12 --parallel
  python project.py submit -o npt_prod_sim --bundle=12 --parallel
  python project.py submit -o calculate_props --bundle=12 --parallel

5. Extract densities, run GP optimization and get samples for the next iteration in the Build_GPs/ directory
   • Note: vle_iters/ will be populated with parameter sets automatically once the termination criteria for LD iters is satisfied

     qsub -N postLD submit_job post_ld.py
   

Surface Tension Optimization (ST Iterations)

To run surface tension iterations, follow the following steps:

1. Initialize Signac workflow

     conda activate hfcs-fffit
     cd Build_GPs/
     python init_IFT.py
   

2. Run ST iterations
   • Note: rm -r workspace/ signac_project_document.json signac.rc will remove everything and allow you to start fresh if you mess up

   cd Build_GPs/vle_iters/
   python project.py submit -o IFT #--bundle=2 --parallel

3. Check status a few times throughout the process

     python project.py status
   

Note: Step 2 operation IFT runs multiple operations in series. Alternatively, the following can be run one at a time: 4. Alternative submission

  python project.py run -o create_forcefield
  python project.py submit -o create_system
  python project.py run -o fix_topology
  python project.py submit -o em_sim #--bundle=2 --parallel
  python project.py submit -o nvt_eq_sim #--bundle=2 --parallel
  python project.py submit -o npzzat_eq_sim #--bundle=2 --parallel
  python project.py submit -o npzzat_prod_sim #--bundle=2 --parallel
  python project.py submit -o init_inter_eq_sim #--bundle=2 --parallel
  python project.py submit -o inter_eq_sim #--bundle=2 --parallel
  python project.py submit -o inter_prod_sim #--bundle=2 --parallel
  python project.py submit -o calculate_props #--bundle=2 --parallel

5. Extract properties, run GP optimization and get samples for the next iteration in the Build_GPs/ directory

     qsub -N postVLE submit_job post_vle.py
   

Note — regenerating the surface-tension GPs will not reproduce the archived DMF model. Five of the six archived surface-tension GPs use a masked linear mean function, but DMF's uses a plain linear mean function; the current code masks all six. The archived DMF model is the one behind the published results. See the "Mean functions, and a known exception for DMF" section of SENSITIVITY_ANALYSIS.md for details.

Distinct Atom Type Optimization

To run generalized FF calibration, follow the following steps:

1. Initialize Signac workflow

     conda activate hfcs-fffit
     cd Opt_ES/
     python init_opt_at.py
   

2. Run LD iterations
   • Note: rm -r workspace/ signac_project_document.json signac.rc will remove everything and allow you to start fresh if you mess up

   cd Opt_ES/opt_ff/
   python project.py submit -o OptGenFF --bundle=12 --parallel

3. Check status a few times throughout the process

     python project.py status
   

Note: Step 2 operation OptGenFF multiple operations in series. Alternatively, the following can be run one at a time: 4. Alternative submission

  python project.py submit -o gen_pareto_sets --bundle=12 --parallel
  python project.py submit -o run_obj_alg --bundle=12 --parallel

5. Extract and analyze data:

     qsub -N postOptAnaly ../submit_job post_analysis_opt.py
   

GP-optimized and Base FF Validation

To validate the generalized FF from above, follow the following steps:

1. Initialize Signac workflow

     conda activate hfcs-fffit
     cd Opt_ES/
   

2. Run GEMC validation simulations (Pvap, rho_v, Hvap) - Simulatenous with Steps 5-7
   • Note: rm -r workspace/ signac_project_document.json signac.rc will remove everything and allow you to start fresh if you mess up

   cd Opt_ES/gemc_val_xxx/
   python project.py submit -o NPT_NVT --bundle=6 --parallel
   python project.py submit -o EQ --bundle=6 --parallel
   python project.py submit -o PROD --bundle=6 --parallel
   

Note: In step 2, "xxx" represents which FF is being validated. use "no_opt" for the Base FF and "opt" for the GP-optimized FF. 3. Check status a few times throughout the process

  python project.py status

Note: Step 2 operation is multiple operations in series. Alternatively, the following can be run one at a time: 4. Alternative submission

  python project.py run -o create_forcefield
  python project.py run -o calc_boxes
  python project.py submit -o NVT_liqbox --bundle=4 --parallel
  python project.py run -o extract_final_NVT_config
  python project.py submit -o NPT_liqbox --bundle=4 --parallel
  python project.py run -o extract_final_NPT_config
  python project.py submit -o run_gemc_eq --bundle=4 --parallel
  python project.py submit -o gemc_eq_restart --bundle=4 --parallel (optional)
  python project.py submit -o run_gemc_prod --bundle=4 --parallel
  python project.py submit -o check_prod_data --bundle=4 --parallel
  python project.py submit -o calculate_props --bundle=4 --parallel
  python project.py run -o del_job (optional)

5. Run validation MD interface simulations (rho_l and surface tension) - Simultaneously with Steps 2-4
   • Note: rm -r workspace/ signac_project_document.json signac.rc will remove everything and allow you to start fresh if you mess up

   cd Opt_ES/ift_val_xxx/
   python project.py submit -o IFT #--bundle=2 --parallel

Note: In step 5, "xxx" represents which FF is being validated. use "no_opt" for the Base FF and "opt" for the GP-optimized FF. 6. Check status a few times throughout the process

  python project.py status

Note: Step 2 operation IFT runs multiple operations in series. Alternatively, the following can be run one at a time: 7. Alternative submission

  python project.py run -o create_forcefield
  python project.py submit -o create_system
  python project.py run -o fix_topology
  python project.py submit -o em_sim #--bundle=2 --parallel
  python project.py submit -o nvt_eq_sim #--bundle=2 --parallel
  python project.py submit -o npzzat_eq_sim #--bundle=2 --parallel
  python project.py submit -o npzzat_prod_sim #--bundle=2 --parallel
  python project.py submit -o init_inter_eq_sim #--bundle=2 --parallel
  python project.py submit -o inter_eq_sim #--bundle=2 --parallel
  python project.py submit -o inter_prod_sim #--bundle=2 --parallel
  python project.py submit -o calculate_props #--bundle=2 --parallel

8. Extract and analyze data:

     cd Opt_ES/
     qsub -N postOptAnaly ../submit_job post_analysis_val.py
   
     cd ../
     qsub -N makeCorrAll submit_job make_corr_figs_all_molec.py
     qsub -N makeCorrInd submit_job make_corr_figs_one_molec.py
     qsub -N makeGPSimSense submit_job make_GP_vs_sim_and_sens.py
     qsub -N makeIFTValfigs submit_job make_ift_val_figs.py
     qsub -N makeParamComp submit_job make_param_comp_figs.py
     qsub -N makeParetoComp submit_job make_pareto_comp_figs.py
     qsub -N makeSVDSens submit_job make_svd_sensitivity.py
   
     cd Claude_Audit/scripts
     python extract_student_ff_values.py
     python recompute_student_metrics.py
     python plot_baseline_vs_experiment.py
     python generate_provenance_table.py
     python build_master_ff_predictions.py
     python build_reference_values_all.py
     python recompute_metrics_all.py
     python generate_tables.py
   

Note: The second set of jobs create final figures for the SI and main text. Alternatively, these can be created by running the jupyter notebook create_analysis_figs.ipynb.

Known Issues

The instructions outlined above seem to be system-dependent. In some cases, users have the following error:

ImportError: /lib64/libstdc++.so.6: version `GLIBCXX_3.4.29' not found

If you observe this, please try the following in the terminal

export LIBRARY_PATH=$CONDA_PREFIX/lib:$LIBRARY_PATH

which should fix the problem. This is not an optimal solution and is something we would like to address. We found that related projects 1, 2 have similar issues. If you are aware of a robust solution to this issue, please let us know by raising an issue or sending an email!

Credits

This research is based upon work supported by the National Science Foundation under award number EEC-2330175 for the Engineering Research Center EARTH. Computing resources were provided by the Center for Research Computing (CRC) at the University of Notre Dame. MC acknowledges support from the Graduate Assistance in Areas of National Need fellowship from the Department of Education, grant number P200A210048.

Contact

Please contact Montana Carlozo (mcarlozo@nd.edu) with any questions, suggestions, or issues.