martimunicoy/peleffy

The peleffy (PELE Force Field Yielder) is a Python package that builds PELE-compatible force field templates.

About
Status
Installation

PELE Force Field Yielder

The peleffy (PELE Force Field Yielder) is a Python package that builds PELE-compatible force field templates. The current supported force fields are:

• Any force field from the Open Force Field toolkit.
• OPLS2005.
• A combination of them.

Documentation

The documentation for the peleffy package is available at GitHub Pages.

Main features

Small molecules parameterization

It can take a molecular structure from a PDB file or a SMILES tag and parameterize it with any of the supported force fields. The resulting parameters are stored in a dictionary-like object that can be easily manipulated.

from peleffy.topology import Molecule
from peleffy.forcefield import OpenForceField

molecule = Molecule('path_to_pdb_file.pdb')
openff = OpenForceField('openff_unconstrained-1.3.0.offxml')

parameters = openff.parameterize(molecule)

Template generation

A peleffy's molecular representation can be employed, along with its parameters, to build a Topology file. A Topology file is a wrapper of topological elements that can be written as an Impact template file, compatible with PELE.

from peleffy.topology import Topology
from peleffy.template import Impact

topology = Topology(molecule, parameters)

impact_template = Impact(topology)
impact_template.to_file('ligand_parameters.txt')

A template with the OBC parameters for the implicit solvent of PELE can also be generated.

from peleffy.solvent import OBC2

obc2_solvent = OBC2(topology)
obc2_solvent.to_file('ligand_obc.txt')

It can also generate the rotamer library file for the side chain exploration of PELE.

from peleffy.topology import RotamerLibrary

rotamer_library = RotamerLibrary(molecule)
rotamer_library.to_file('ligand_rotamers.txt')

User-friendly CLI

All the commands above can be run with a single CLI directive.

python -m peleffy.main path_to_pdb_file.pdb -f 'openff_unconstrained-1.3.0.offxml' --with_solvent

Contributors

For a full list of contributors, see the GitHub Contributors page.

peleffy's Issues

Default molecule tags do not match

Description

The default tag for a molecule is UNK. However, when writing an untagged molecule to a PDB file, it is written as UNL. This mismatch causes some problems when using these PDB files with parameter templates created for a UNK ligand.

Solution

Consider using UNL as the default molecule tag or always apply the molecule tag from the offpele representation to the RDKit molecule object to obtain a correct PDB file from it.

Ignore terminal rotatable bonds

Allow the user to ignore terminal rotatable bonds when generating the rotamer library of a ligand. By default, PlopRotTemp ignores the terminal rotatable bond of each branch, so, if we want to obtain the same results as PlopRotTemp we need to add this feature to offpele.

Proper dihedrals with idivf different from 1

Description

We were expecting that all OpenFF's proper dihedrals had idivf term equal to 1.0. However, some proper dihedrals like t84 and t86 are not satisfying this condition.

Solution

The condition to force idivf term in offpele to equal 1 will be removed.

Use the PELE's default DataLocal path for OPLS

Description

When an Impact template with OPLS parameters is built and the output files are saved following the DataLocal hierarchy, we should take into account that the output path changes. So, there are two possible paths and offpele should choose between them depending on the source of the parameters that are written into the Impact template:

• OPLS template: DataLocal/Templates/OPLS2005/HeteroAtoms/
• OFF template: DataLocal/Templates/OFF/Parsley/HeteroAtoms/

Solution

Add two possible output paths and let the Impact template builder choose where to save it depending on the source of the parameters.

Wrong molecule tag is written to PDB files

Description

Molecule tags are used as residue names in the PDB writer. However, the current implementation does not update the RDKit's default residue name after a custom tag is defined in the main Molecule class.

Solution

Modify the implementation to assign the residue name in the PDB according to the Molecule's tag.

Extend support to openff-1.3.0 and test it

Description

OpenFF 1.3.0 has been released. It should be included to our toolkit. Besides, it would be interesting to have a suitable test to check that every new force field addition works as expected.

Specify dependency versions

Description

To avoid unexpected incompatibilities or bugs in the future, we should specify versions for the main dependencies of offpele.

Is RDKit correctly assigning connectivities from ligand PDBs?

We should double check some of the "known tricky" cases below:

• Tertiary nitrogens like NH4
• Carboxylic acids
• Carboxylates

Unexpected dihedral parameters from Open Force Field Toolkit

Description

We noticed that some of the parameters assigned to torsions by the Open Force Field Toolkit are unexpected. The current implementation of offpele expects the following values for torsional parameters:

• phase constant: either 0.0 or 180.0 degrees
• periodicity: one of (1, 2, 3, 4, 6)

See issue #33 for more details.

While the requirement in the periodicity value is trivial, the phase constant does really need to be either 0.0 or 180.0 degrees. That is because of PELE's current dihedral equation (explained here), which is:

Since, there is not a phase constant term, we need to use the prefactor which can equal 1 or -1. The sign determines whether the shape of the torsion behaves as in Open Force Field with a phase constant equal to 0.0 or 180.0 degrees. Thus, in case that the Open Force Field Toolkit works with phase values other than 0.0 and 180.0 degrees we will need to change the equation of dihedrals in PELE.

Integrate method for the OPLS OBC parameters

Description

PELE has a script to obtain the OBC parameters given a small molecule. This script is quite old though and presents many compatibility problems. Maybe we can incorporate this script into the offpele package to automatize the parameterization workflow. We could also look at its compatibility issues to try to improve it.

Solution

This would be a new Solvent class, that given a Molecule representation, parameterized with OPLS2005 would use atom types to assign the OBC parameters. Then, it should write the template file in the corresponding format file expected by PELE.

Compatibility with OpenFF allow_undefined_stereo parameter

Description

OpenFF allow the initialization of molecules with undefined stereochemistry data (although the default behaviour is always complain about it). In some cases, (mainly during benchmarking) we might happen to skip this checking.

Solution

Allow the user to skip the stereochemistry verification through a parameter in the initialization of the Molecule class.

Missing charge module in conda package

After installing the package using conda, the charge module is missing. Full traceback running the main module as shown in the documentation:

Traceback (most recent call last):
  File "/home/jgilaber/anaconda3/envs/offpele_env/lib/python3.7/runpy.py", line 193, in _run_module_as_main
    "__main__", mod_spec)
  File "/home/jgilaber/anaconda3/envs/offpele_env/lib/python3.7/runpy.py", line 85, in _run_code
    exec(code, run_globals)
  File "/home/jgilaber/anaconda3/envs/offpele_env/lib/python3.7/site-packages/offpele/main.py", line 277, in <module>
    main()
  File "/home/jgilaber/anaconda3/envs/offpele_env/lib/python3.7/site-packages/offpele/main.py", line 273, in main
    args.output, args.with_solvent, args.as_datalocal)
  File "/home/jgilaber/anaconda3/envs/offpele_env/lib/python3.7/site-packages/offpele/main.py", line 209, in run_offpele
    from offpele.topology import Molecule
  File "/home/jgilaber/anaconda3/envs/offpele_env/lib/python3.7/site-packages/offpele/topology/__init__.py", line 2, in <module>
    from .molecule import Atom, Molecule
  File "/home/jgilaber/anaconda3/envs/offpele_env/lib/python3.7/site-packages/offpele/topology/molecule.py", line 13, in <module>
    from offpele.charge import (Am1bccCalculator, GasteigerCalculator,
ModuleNotFoundError: No module named 'offpele.charge'

Version is 0.3.1, I checked in the installation folder and the charge subfolder was not installed

Combine OPLS parameters with OFF partial charges

Description

We can combine bonding parameters from the Open Force Field Toolkit with OPLS2005 parameters. When OPLS is used to parameterize nonbonding parameters, charges are automatically assigned from this source. However, a nice add-on would be to allow the user to assign partial charges from one of the OFF methods although the rest of the nonbonding parameters are assigned with OPLS.

Check out input PDB files prior parameterizing

Description

For instance, offpele would fail if the input PDB that contains the ligand to parameterize has atoms with unconsistent residue names. The exception raised in this case is totally misleading.

Solution

Implement a safety check for PDB files in order to properly handle exceptions related with its format prior running the parameterization.

Allow the assignment of external partial charges

Description

PlopRotTemp allows the assignment of partial charges coming from external suites such as Jaguar. We need to implement this feature in offpele too.

Solution

A new Charge method that would read an external file containing the partial charges to assign to the Molecule representation.

Construct offpele's molecule from SMILES

We can construct molecules from SMILES tags. This method would ease the automatization workflows that employ offpele. For this purpose, we need three main features:

• Build RDKit's molecule from smiles using the method Chem.MolFromSmiles.
• New function to automatically assign atom names to the molecule (as we no longer will feed from a PDB file).
• A function to write the new molecule to a PDB file for PELE (PELE will still need a PDB file with the correct conformation and atom names assigned).

Implement a conda install

Many dependencies of offpele can be easily installed through conda. However, some of them lack of an pip install method. For this reason, we are sort of forced to use conda to install offpele.

Generate Sphinx documentation

We need to documentate the package with Sphinx.

Wrong default output name for Impact template and rotamer library

Description

Impact templates and rotamer libraries are saved with the molecule's name instead of the tag.

Solution

Replace name by tag.

At least one rotatable bond must be inside the rotamer library core

Description

Right now, the core only holds a set of atoms, centered in the molecular graph, with no rotatable bonds. A single rotatable bond should be inside the core to increase the performance of PELE when sampling the rotamer library.

Solution

Change the implementation of the core assignment in the MolecularGraph class.

Test offpele with complex non-standard residues

Description

Let's see how OpenFF works when parameterizing a set of cofactors / prosthetic groups

Solution

@SergiR1996 could provide several examples of compounds to test. It would also be nice to rely on QM calculations for these compounds to use as reference when checking the resulting forcefield parameters.

Workflow to automatically parameterize protein residues

Design a workflow to generate automatically the parameters of standard protein residues with Open Force Field toolkit.

Change prints to logs

In order to be able to use the package via API and easily change the stdout level let's switch from prints to logging. https://docs.python.org/3/library/logging.html

Consider replacing the use of label_molecules() by create_openmm_system() from OFF's ForceField

Description

create_openmm_system runs the parameter assignment for OpenMM and it might be a safer option than label_molecules(). For example, the latter does not care about aspects such as exclusion lists (see issue #24).

Solution

Check the differences between both methods before performing any substitution. In case, it happens to be a better method to obtain the parameters for PELE, replace it in: toolkits.py

Allow the combination of PELE’s VDGBNP terms with OFF terms

Benchmarks with OBC solvent parameters exhibited disparate energies when predicting hydration free energies. For this reason, we will try to combine OFF parameters with the VDGBNP terms (which displayed a more reasonable correlation).

However, the VDGBNP implementation requires OPLS2005 atom types to be assigned to small molecules. Therefore, we will need to add as a dependency Schrodinger’s ffld_server, which is the tool that is used for this purpose.

Handle exclusions in PELE 1-4 list from dihedral definitions

Update about PELE

We realized that PELE uses the proper dihedral definitions in the Impact template to generate the 1-4 list. For this purpose, it uses first and fourth atoms from the dihedral definition.

Dihedral: [ 1,  2,  6,  8 ] --> new 1-4 pair: 1-8

However, proper definitions are not unique. Thus, to avoid duplicates in the 1-4 list it only adds the corresponding pair when a new proper dihedral (with at least 1 difference in any of the atoms that are involved in the dihedral) is found.

Dihedral: [ 1,  2,  6,  8 ] --> new 1-4 pair: 1-8
Dihedral: [ 1,  2,  6,  8 ] --> X (since it has been already included before)
Dihedral: [ 1,  2,  6,  9 ] --> new 1-4 pair: 1-9

Finally, we must exclude from the 1-4 list those pairs that may actually belong to 1-2 or 1-3 relations. This filtering is not internally done in PELE and has to be done when generating the list of proper dihedrals for the Impact template. So, by looping over bonds angles and previous dihedrals we need to detect these cases. PELE will exclude from the 1-4 list dihedrals whose second or third atom indexes are defined in the Impact template with a minus sign in front.

Dihedral: [ 1,  2,  6,  8 ] --> new 1-4 pair: 1-8
Dihedral: [ 1,  2,  6,  8 ] --> X (since it has been already included before)
Dihedral: [ 1,  2,  6,  9 ] --> new 1-4 pair: 1-9
Dihedral: [ 1,  3, -5,  8 ] --> X (we need to add the minus sign in front of the second or third indexes to avoid adding twice pair 1-8 in the 1-4 list)
Dihedral: [ 1,  2, -6,  7 ] --> X (if we imagine that 1-7 pair belongs to a bond or an angle, then, we would need to add the minus sign in front of the second or third indexes to avoid adding this pair that actually is not a 1-4 relation)

Originally posted by @martimunicoy in #33 (comment)

Are OFF torsions translated correctly to PELE?

Description

The equation that OFF is using is explained here, which is of the form:

On the other side, the dihedrals equation in PELE is explained here, which is:

The function that is responsible for translating OFF torsions to PELE is this one. In this way, it expects the following values:

• periodicity: any integer from 1 to 6
• phase: either 0 or 180 degrees
• k: any quantity in kcal/mol
• idivf: any integer

Then it assigns, PELE's parameters as:

• periodicity: it equals the analogous OFF term
• k: it equals the analogous OFF terms divided by idivf
• prefactor: it is set to 1 when OFF phase equals 180 degrees and it is set to -1 when OFF phase equals 0 degrees.

This OFF-PELE translation is tested at test_parameters.py. However, we should deepen its validation.

Multiple partial charges calculations

It is possible to use Antechamber with different calculation methods to extract partial charges from a structure. We should investigate if there are other methods that can be used to generate partial charges faster than the current method (am1bcc).

Add individual tests for the RDKit toolkit

Description

Methods in the RDKit toolkit are not individually tested. We should test them

Impact class needs to be tested

Description

Currently, Impact class is not tested. We should design some tests to check it out.

Improve documentation for version 3.0

Description

At least, the usage of CLI of main.py should be documented. Besides, some docstrings can be updated to the API documentation.

Correction of bad geometries

We were wondering if peleffy can correct ligands with bad geometries but correct connectivity. This feature would be very useful for FragPELE simulations. I provide here a PDB example with an hydrogen atom in a bad position.

growing_result.txt

Molecule from SMILES is unnamed

Description

When an offpele Molecule is initialized from a SMILES tag, no name is set to it. This behavior can cause some problems later when generating the template files.

Solution

We should assign a default name in Molecule objects.

Allow ligand core to be in an edge

Description

The new feature of PELE, called side chain perturbation, will need a custom template with a core defined in the region covalently attached to the protein backbone.

Solution

offpele should support the creation of a ligand template with a custom core containing a specific set of atoms (which would be in a molecule's edge, not in the center).

Test Coverage

When unit-test are implemented let's try to setup test coverage to check all lines are validated. Also, as github is linked with codecov (codecov.io) every time we PR or merge a branch we will be able to see whether we improve coverage or not, and therefore, merge or not consequently.

Wrong output name for the solvent parameter file

handle_output_paths() function in main.py assigns an invalid output name to the solvent parameter file. PELE looks for a ligandParams.txt file, not a LIG_solv.json.

Github Actions (CI)

Let's implement github actions workflows so we know whether we broke any code or have any missing dependency!

• Deployment pypi
• Deployment conda
• Source code deployment
• Testing

Warn the user about the Schrodinger's environment variable

The Schrodinger Toolkit introduced in #29 requires an environment variable to be set with the path to Schrodinger's installation path. There is no warning about this requirement and offpele simply crashes when no such variable is found. We should warn the user about this prerequisite.

Assign Molecule's connectivity with an RDKit template

Description

If we want to load PDB files with no connectivity data (such as those coming from PELE simulations) we need a method to assign the missing connectivity.

Solution

We can use an RDKit molecule representation with the correct connectivity assigned as the connectivity template. Then, the RDKit's method AllChem.AssignBondOrdersFromTemplate can be employed to assign the connectivity to a third molecular representation using the previous connectivity template.

Describe all modules through docstrings

We need to describe all modules, classes and functions of this first release through docstrings.

Allow a different RDKit conformer to be set to offpele's molecule

Description

Some benchmarks need to apply conformational variations to an offpele's molecule instance. A method to set a new RDKit's conformer would be useful in these cases.

Change repository name

This package is no longer exclusively employing OpenFF force fields (as it can also employ OPLS2005). Thus, its name should be more general.

Some proposals:

• pelep: PELE Parameterizer
• pelepar: PELE PARameterizer
• peleup: PELE Universal Parameterizer
• param2pele: PARAMeters to PELE
• peleffy: PELE Force Field Yielder
• pelet: PELE Templatizer

Load parameters from an Impact file

Description

A nice add-on to the package would be a method to get a BaseParameterWrapper from an Impact file.

Solution

Implement a static method in the BaseParameterWrapper class to load the parameters from an Impact file.

Missing dummy propers in OPLS

Description

The OPLS implementation skips any proper that has a null force constant. However, we discovered that PELE still needs the definitions of these propers (dummy propers because they have this null force constant) for the consider them in the 1-4 list.

Solution

Prevent offpele to skipping these dummy propers, despite having no effect in the proper energetic term.

Insert unit tests

We need to unit test as many modules as possible. Right now, only a single test for dihedral equations has been develop.

Check compatibility between 1-4 scaling factors between OpenFF and PELE

PELE uses its built-in implementation of the exclusion lists and 1-4 scaling factors. Therefore, there is no need to previously exclude covalently linked atoms nor apply scaling factors on 1-4 interactions from OpenFF.

Regarding the implementation of scaling factors in OpenFF: [1] [2]

Test main.py function

Description

main.py function contains many untested lines. We should increase the coverage on this function

Full compatibility with OPLS2005 force field

Description

Currently, we are partially supporting OPLS2005 parameters. We should allow the standalone use of OPLS2005 parameters.

Solution

Define a new method to assign the OPLS2005 parameters to the ligand employing the already implemented SchrodingerToolkitWrapper.

OPLS2005 would not keep PDB atom names

Description

When using the OPLS2005 force field, PDB atom names are not kept but they ones from ffld_server are used.

Solution

Instead of assigning PDB names from the ffld_server, we must use the get_atom_names method from the RDKitToolkitWrapper.

Torsions terms -1

In the different torsions terms, the last element of the column which controls the sign of the sin/cos term is always -1 while in OPLS2005 variates between +-1. I haven't looked at the complete trigonometric function to check whether is ok or not but let's have a look when we have a moment!