parameterization of partial charge of symmetric molecules.

[xiki-tempula]

I'm interested in using ParaMol for small molecule parameterization. For molecule such as tert-butyl alcohol, the three methyl group should have the same partial charge. I wonder how should one implement this in the RESP routine? Thank you.

[JMorado]

Hi,

Say you have the following tert-butyl molecule.

You want to perform a typical 2-stage RESP procedure where in the first step you allow all charges to vary freely without symmetry constraints, and in the second step you just optimize symmetrically-equivalent charges imposing symmetry constraints. The following script could be used to perform this in a "manual" fashion (requires manual identification of the degenerate atoms):

import os
import sys

# Path to ParaMol
sys.path.insert(0, os.path.abspath('path_to_paramol'))

# ParaMol imports
from ParaMol.System.system import *
from ParaMol.MM_engines.openmm import *

# ParaMol Tasks imports
from ParaMol.Tasks.resp_fitting import *
from ParaMol.Tasks.torsions_scan import *
from ParaMol.Utils.settings import *
from ParaMol.Utils.gaussian_esp import *
from ParaMol.Utils.Symmetrizers.amber_symmetrizer import *

# --------------------------------------------------------- #
#                         Preparation                       #
# --------------------------------------------------------- #
# Create the OpenMM engine for caffeine
openmm_system = OpenMMEngine(init_openmm=True, topology_format='AMBER', top_file='lig.prmtop', crd_format='AMBER', crd_file='lig.inpcrd')

# Create Molecular System
ligand = ParaMolSystem(name="ligand", engine=openmm_system, n_atoms=14)

# Create ParaMol's force field representation and ask to optimize charges
ligand.force_field.create_force_field(opt_charges=True)

# Create ParaMol settings instance
paramol_settings = Settings()
paramol_settings.properties["include_energies"] = False
paramol_settings.properties["include_forces"] = False
paramol_settings.properties["include_esp"] = True

# Regularization
paramol_settings.properties["include_regularization"] = True
paramol_settings.properties["regularization"]["method"] = "hyperbolic"

# Weighting
paramol_settings.objective_function["weighting_method"] = "boltzmann"
paramol_settings.objective_function["weighting_temperature"] = 300.0*unit.kelvin

# --------------------------------------------------------- #
#                 Read ESP Data into ParaMol                #
# --------------------------------------------------------- #
gaussian_esp = GaussianESP()

# Log file names
log_files = ["path_to_log_files"] 

ligand.ref_coordinates, ligand.ref_esp_grid, ligand.ref_esp, ligand.ref_energies = gaussian_esp.read_log_files(log_files)
ligand.n_structures = len(ligand.ref_coordinates)

# --------------------------------------------------------- #
#                Symmetrize ParaMol ForceField              #
# --------------------------------------------------------- #
# Symmetry ParaMol ForceField so that it respects symmetries
# In this example, we are not setting any symmetry, but we still need to do this step as we want to save a .mol2 file
amber_symmetrizer = AmberSymmetrizer(top_file="lig.prmtop", xyz="lig.inpcrd")
amber_symmetrizer.get_symmetries(ligand.force_field)
amber_symmetrizer.save("ligand_initial.mol2")

# --------------------------------------------------------- #
#                 RESP Charge Fitting - Stage  1            #
#                                                           #
#          All atoms are allowed to vary their charges      #
# --------------------------------------------------------- #
paramol_settings.properties["regularization"]["scaling_factor"] = 0.01

resp_fitting = RESPFitting()

# Perform RESP fitting
systems, parameter_space, objective_function, optimizer = resp_fitting.run_task(paramol_settings, [ligand], solver="explicit", total_charge=0)

# Write ParaMol Force Field file with final parameters
ligand.force_field.write_ff_file("ligand_resp_stage1.ff")

# Update amber symmetrizer and save .mol2 file
amber_symmetrizer.update_term_types_parameters(parameter_space.optimizable_parameters)
amber_symmetrizer.save("ligand_resp_stage1.mol2")

# --------------------------------------------------------- #
#                 RESP Charge Fitting - Stage  2            #
#                                                           #
#      Constraint all degenerate atoms to equal charge      #
# --------------------------------------------------------- #
paramol_settings.properties["regularization"]["scaling_factor"] = 0.001

# Indexes of denegerate atoms
degenerate_hydrogens = [5,6,7,8,9,10,11,12,12] 
degenerate_carbons = [2,3,4]

# Set symmetry mappings
symmetries_dict = {"Q1" : degenerate_hydrogens, 
                   "Q2" : degenerate_carbons}  

for force_field_term in ligand.force_field.force_field["NonbondedForce"][0]:
    # Iterate over all atoms of a given force field term
    for at in force_field_term.atoms:
        # Iterate over dict
        for key, value in symmetries_dict.items():
            if at in value:
                force_field_term.parameters["charge"].optimize = 1
                for parameter in force_field_term.parameters.values():
                    parameter.symmetry_group = key
                break
            else:
                force_field_term.parameters["charge"].optimize = 0
                for parameter in force_field_term.parameters.values():
                    parameter.symmetry_group = "X"

amber_symmetrizer.get_symmetries(ligand.force_field)

# Perform RESP fitting
systems, parameter_space, objective_function, optimizer = resp_fitting.run_task(paramol_settings, [ligand], solver="explicit", total_charge=0)

# Write ParaMol Force Field file with final parameters
ligand.force_field.write_ff_file("ligand_resp_stage2.ff")

# Update amber symmetrizer and save .mol2 file
amber_symmetrizer.update_term_types_parameters(parameter_space.optimizable_parameters)
amber_symmetrizer.save("ligand_resp_stage2.mol2")

A more automatic solution has to necessarily involve the automatic identification of the equivalent atoms. That probably can be done with RDKit (possible using GetSubstructMatches or something similar) - honestly, I haven't put much thought into it.

But, as you can see, once equivalent atoms have been identified and the symmetry mappings dictionary has been defined, the following bit of code can be used to set the symmetries in ParaMol's internal Force Field representation:

for force_field_term in ligand.force_field.force_field["NonbondedForce"][0]:
    # Iterate over all atoms of a given force field term
    for at in force_field_term.atoms:
        # Iterate over dict
        for key, value in symmetries_dict.items():
            if at in value:
                force_field_term.parameters["charge"].optimize = 1
                for parameter in force_field_term.parameters.values():
                    parameter.symmetry_group = key
                break
            else:
                force_field_term.parameters["charge"].optimize = 0
                for parameter in force_field_term.parameters.values():
                    parameter.symmetry_group = "X"

Hope this helps.

Best, João

[JMorado]

Alternatively, it is also possible to set the symmetries (also manually) in the .ff files that ParaMol reads/writes.

See this example:

https://paramol.readthedocs.io/en/latest/Examples/Example_5/Example_5.html

And specifications of the .ff file:

https://paramol.readthedocs.io/en/latest/Files_specification.html

[xiki-tempula]

Thank you. I will give it a try.