Workflow for FFPSimul

#######################################################
# Author: Anmol Kumar                                 #
# Date: 06/23/2023                                    #
# FFPsimul is part of FFParam v2.0                    #
# Copyright: University of Maryland Baltimore 2023    #
#######################################################

Free energy of solvation in FFPSimul

After I realised that it is less hassle to run heat of vaporization/sublimation and free energy of solvation simulation using shell script, I made easy to use scripts for running these simulations. (For both additive and drude).

Structure of simulation scripts for free energy of solvation

ffparam
  |-- ffparam
      |-- ffpsimul  # This is the new directory holding all scripts
          |-- fep_0_setup.sh
          |-- fep_1_run_eqpdlrc.sh
          |-- fep_2_run_bulk.sh
          |-- fep_3_run_gas.sh
          |-- fep_4_run_wham.sh
          |-- fep_5_getresult.sh

Highlights

  • Uses packmol for generating simulation box.

  • If solvent is water, it can also use preequilibrated water box, although number of solvent molecules will vary for different systems.

  • For Drude systems, the packmol generated simulation box may result in shake error during dynamics. Thus it is advised to either use preequilibrated water box option.

  • If packmol is used to obtain simulation box for Drude systems with water as solvent, the structure of the script allows a very easy preequilibration using additive force field.

Usage

$ export PYTHONPATH=/home/anmol/ffparam

$ $PYTHONPATH/ffparam/ffpsimul/fep_0_setup.sh additive=<true/false> resi=<resname> resipdbfile=<residue pdb file> resistrfile=<strfile; can be any file containing resi topology in charmm format>

Optional Parameters:
  workdir=<working directory; default = ".">
  nsolv=<number of solvent used in bulk; default = 500 >
  solv=<residue name of solvent; default = "TIP3" >
  solvpdbfile=<if residue is other than water, then provide its pdb file>
  solvstrfile=<if residue is other than water, then provide file containing topology of solvent. Dont if resistrfile has solvent topology also.>
  autogenangdihresi=<generally true>
  autogenangdihsolv=<generally false>

Note:
  autogenangdih true means angles and dihedrals will be autogenerated internally based on bond description.
  autogenangdih false means angles and dihedrals will not be autogenerated and angle and dihedral descriptions provided in topology file will be used.
  For water topology in CHARMM, autogenangdih should be false.
  If additive=false and solvent is water then residue name is changed to SWM4 inside the script.

For following scripts optional parameters are not mentioned here. User can check all options by running script without argument.

$ $PYTHONPATH/ffparam/ffpsimul/fep_1_run_bulk_eqpdlrc.sh additive=<true/false> resi=<resname> resiboxcorfile=<box crd or pdb file> resistrfile=<strfile; Dont provide if already in toppar>

  • This script first performs equilibration and writes the equilibrated coordinates.

  • User should wait before running fep_2_run_bulk.sh script till equilibrated coordinate file are written in bulk directory.

  • The script further runs production run and obtains long range energy correction and puts in lrc_energy.csv

  • For Drude systems with water as solvent, user can create and equilibrate another simulation box using additive force field.

  • Once equilibration is done for additive system, the TIP3 resname in equilibrated coordinates can be converted to SMW4 using simple file editing.

  • This updated coordinate file can then be supplied for going ahead with Drude equilibration and following steps.

$ $PYTHONPATH/ffparam/ffpsimul/fep_2_run_bulk.sh additive=<true/false> resi=<resname> resistrfile=<strfile> resiboxcorfile=<box crd file>

  • This script creates simulation inputs for bulk system against each lambda value and runs the simulations.

  • All output is placed in bulk directory.

$ $PYTHONPATH/ffparam/ffpsimul/fep_3_run_gas.sh additive=<true/false> resi=<resname> resistrfile=<strfile> resicorfile=<box crd file>

  • This script creates simulation inputs for gas phase molecule against each lambda value and runs the simulations.

  • All output is placed in gas directory.

$ $PYTHONPATH/ffparam/ffpsimul/fep_4_run_wham.sh phase=<bulk/gas>

  • This script combines all the free energies corresponding the different lambda value,

  • and performs wham postprocessing available in charmm.

  • This script should be run for both bulk and gas phase.

  • The output is placed in bulk_freenergy.csv and gas_freenergy.csv

Note

RIGIDITY ALERT: The script assumes that bulk and gas directory contains wham files as generated by fep_2_run_bulk.sh and fep_3_run_gas.sh script.

$ $PYTHONPATH/ffparam/ffpsimul/fep_5_getresult.sh bulkGfile=<bulk_freenergy.csv> gasGfile=<gas_freenergy.csv> resicharge=<residue charge> lrcEfile=<lrc_energy.csv>

  • Uses the output csv files to get the free energy of solvation.