esol.py: Compute solvation energies with the E-sol Workflow
The script esol.py is a workflow for calculating the energy of solvation with QM methods [324]
The workflow contains the following stages:
- A conformer search is performed with Macromodel. The phases that are considered is specified with the conf_search_phases keyword, and the number of conformers generated for each phase is specified with the n_conformers_conf_search keyword.
- State penalty is calculated with Epik. If use_input_state_penalties is set, structures that already have the
r_epik_State_Penaltyproperty defined and are neutral with no formal charges skip the Epik stage and use the precomputed state penalty. If the structures have reference pKas listed under thes_j_ref_pkasproperty in themaefile, they are used to correct the Epik state penalty. Note: If you set conf_sets_provided=True, there can be no clashings_j_ref_pkasvalues within conformer sets. It is ok if not all of the conformers in the set have thes_j_ref_pkasproperty, but any that do must have the exact same value. - Gas phase QM geometry optimizations are performed with B3LYP/LACVP* on each conformer.
- PBF water single point energy calculations are performed with M06-2X/LACVP** on the lowest energy geometry optimization structures. The number of conformers for which single point energy calculations are performed are specified with n_conformers_to_score.
- Using the energies from Step 4, the hydration energy dE (kcal/mol) is calculated as Lowest solution energy - Lowest gas energy.
- The energy of solvation (E-sol) is calculated as Hydration Energy - State Penalty.
- The output file
esol_results_{jobname}.csvcontains relevant properties computed for each input molecule.{jobname}-out.maecontains the lowest energy gas and solution phase structures, with the respective structures titled as{mol_name}_gas_phase_conformerand{mol_name}_solution_phase_conformer.{jobname}.outcontains a summary of the various stages of the jobs submitted by the workflow.
By default, the geometry optimization and single point stages use the QM methods described above. You can instead run these stages with MLFF or xTB methods, specified through the geopt_functional, geopt_solvation, sp_functional, and sp_solvation keywords. The protocol keyword provides a convenient way to preset these and other keywords for common E-sol configurations, such as using MLFF for some or all stages, rather than setting each keyword individually. Any keyword that you set explicitly overrides the protocol settings.
Using the esol.py script
The esol.py script uses an input file that contains keywords to specify the parameters of the calculation. The command syntax is:
jaguar run esol.py inputfile.in
See Table 1 for a full list of the available esol.py keywords and their descriptions.
The input file has &JaguarKeywords_QMStageGeopt and &JaguarKeywords_QMStageSPE sections which allow the user to apply Jaguar keywords that do not have a corresponding inputfile keyword to the calculations in steps 3 and 4. The legal keywords in this section are the same as those found in the &gen section of a typical Jaguar job. You can also use the &JaguarKeywords section to specify keywords for both stages. If keywords are listed in both the &JaguarKeywords and a stage specific section (e.g., &JaguarKeywords_QMStageGeopt), the stage specific keyword value takes precedence. You cannot set the basis/functional/solvation method in these sections since they are already controlled by E-sol input file options. These sections should be used with care since the QM jobs can easily crash with incompatible Jaguar keywords (for instance ecd=1 and itddft=0 will fail since ecd jobs require tddft to be run).
Here is an example .in file. In the example, for each molecule in molecules.mae, the calculation get up to 80 conformers out of the conformer search (40 for gas, 40 for water) and up to 80 geometry optimization jobs are launched. However, only 8 conformers should go on to the PBF single point energy stage.
input=molecules.mae n_conformers_conf_search=40 n_conformers_to_score=8 conf_search_phases = ['gas','water'] &JaguarKeywords iuhf=2 & &JaguarKeywords_QMStageGeopt maxitg=10 & &JaguarKeywords_QMStageSPE maxit=200 &
| Keyword | Type | Allowed Values | Default Values | Description |
|---|---|---|---|---|
| conf_energy_window | real number | Any | 5.0 | Window to filter high-energy structures at end of QRNN/MLFF conformational filtering stage. Should be a positive number. |
| conf_filter | string | Any | NONE | Use MLFF geopts to filter conformers obtained from Macromodel conformer search. Multiplies the Macromodel energy threshold (5 kcal/mol) and "n_conformers_conf_search" values by "conf_filter_factor" and then uses MLFF energies to filter the conformers with respect to "conf_energy_window". A maximum of "n_conformers_to_score" conformers will pass the filter. If molecules are unsupported by QRNN/MLFF, fall back to xTB conformer filtering protocol. QRNN/MLFF require a forcefield license. Available options are: ['QRNN', 'MLFF', 'XTB', 'NONE']. |
| conf_filter_factor | real number | Any | 1.5 | Factor by which to loosen energy threshold and n_conformers_conf_search for Macromodel conformer search job if QRNN/MLFF conformational filtering stage has been requested. Should be a positive number. |
| conf_search_phases | string | ['gas', 'water'] |
[gas, water] | Phase(s) to run conformer search in. Only used when `conf_search_program=macromodel`. |
| conf_search_program | string | macromodel prime |
macromodel | Program to use for conformer search. Options are "macromodel" or "prime". "prime" search is only applicable for systems containing macrocycles (rings with 10 or more atoms). |
| conf_sets_provided | boolean | Any | False | If set, groups the molecules in "input" into sets of conformers and skips conformer search. |
| geopt_basis | string | Any | LACVP* | Basis to use during QM geopt stage |
| geopt_functional | string | Any | B3LYP | Functional to use during QM geopt stage |
| geopt_solvation | string | NONE PBF PCM SMD QRNN_ALPB MLFF_ALPB XTB_ALPB |
NONE | Solvation method to use during QM geopt stage |
| input | string | Any | None | Path (absolute or relative) to a file (.mae/.maegz/.sdf/.sdfgz) containing molecular structures. |
| mlff_auto_fallback | string | NONE QM XTB |
XTB | Controls fallback behavior for molecules unsupported by QRNN/MLFF. Only applies if QRNN/MLFF is requested for a QM stage of E-Sol. `NONE`: Unsupported molecules are written to a separate .mae file. `QM`: Unsupported molecules fall back to default QM settings. `XTB`: Unsupported molecules fall back to xTB (default). |
| n_conformers_conf_search | integer | Any | 12 | Maximum number of conformers to keep from each conformer search. |
| n_conformers_to_score | integer | Any | 12 | Maximum number of conformers to score with QM. |
| protocol | string | default mlff_filter mlff_geopt mlff_energy full_exhaustive_dft |
default | Select a protocol for the E-Sol calculation. `default`: Default settings. The `mlff_<stage>` protocols use MLFF up through the named stage and QM for any later stages: `mlff_filter`: MLFF conformer filtering, then QM geopt and single-point. `mlff_geopt`: MLFF conformer filtering and geopt, then QM single-point. `mlff_energy`: MLFF for all stages. `full_exhaustive_dft`: QM for all stages with exhaustive conformer search. Individual keywords override the protocol settings when both are specified. |
| sp_basis | string | Any | LACVP** | Basis to use during QM single-point stage |
| sp_functional | string | Any | M06-2X | Functional to use during QM single-point stage |
| sp_solvation | string | PBF SM8 SMD PCM QRNN_CPCMX MLFF_CPCMX XTB_CPCMX |
PBF | Solvation method to use during QM single-point stage |
| use_input_state_penalties | boolean | Any | False | Skip the Epik run and use the existing Epik state penalty for totally neutral structures (no formal charges) with the `r_epik_State_Penalty` property. |
| prepare_input_structures | boolean | Any | True | Preprocess input structures, converting coordinates from 2D to 3D and adding any missing hydrogens. |