isra3l / ligpargen Goto Github PK

License: MIT License

Python 98.76% Rich Text Format 1.24%

molecular-dynamics molecular-dynamics-simulation opls opls-aa cm1a alchemical-free-energy-calculations monte-carlo-simulation

ligpargen's Introduction

LigParGen v2.1

Author:   Israel Cabeza de Vaca Lopez
Email:    [email protected] // [email protected]
Place:     William L. Jorgensen Lab at Yale University // Jens Carlsson Lab at Uppsala university
Date:    2020-2021

Description: An automatic OPLS-AA parameter generator for small organic molecules using CM1A, 1.14CM1A and CM1A-LBCC charge models. LigParGen accepts any Open Babel molecular format including SMILES, PDB, MOL, MOL2, among others. Final OPLSAA parameter outputs will be written in topology/coordinate input files for BOSS, Q, Tinker, PQR, openMM, CHARMM/NAMD, Gromacs, LAMMPS, Desmond, and xplor softwares.

Note: This new version has been written from scratch but it is based on the Leela Dodda initial ligpargen python code.

WARNING: This version of ligpargen requires BOSS5.0 to work correctly. Some features like CM1A-LBCC are not available with the BOSS version provided by the Jorgensen Lab website. I strongly recommend you request a BOSS5.0 executable by sending an email directly to Prof. Jorgensen ([email protected])

New LigParGen features:

The order and the name of the atoms will remain the same in the output files.
This new version of the ligpargen includes a robust version of the alchemical transformation method to generate single and dual topologies for four different molecular mechanics softwares (BOSS, CHARMM/NAMD, Gromacs, and Tinker).
Sanity checks to detect incorrect inputs have been implemented (incompatible charge model, net molecule charge, input format, ...).
A log file to check the inputs, outputs, and intermediate processes has been created. This allows tracking the input information (charge model used, input molecule,...) and also provides useful warnings in case of error.
Automatic net charge detection in the input molecule. If the user specifies a different charge than one automatically estimated, the log file will include a warnning.
Atom XYZ positions in the molecule input remain unchanged in the output files.
Hydrogen atoms will be added automatically just for SMILES inputs. Molecules in any other input format require to have all hydrogens to provide more flexibility of the protonation states. In this way, the user can avoid parameterization problems with tautomers or stereoisomers.
Different molecule input format, net charges, charge models, and optimization steps can be used for each molecule in alchemical transformations.
Bugs fixed (Q wrong torsion parameters, alchemical molecule overlap failure,...)
Additional default input parameters have been included such as residue name, charge model,...

INSTALATION

LigParGen requires the free BOSS software to generate the OPLSAA parameters.

1 - Download and install BOSS software from the official William L. Jorgensen lab website: http://zarbi.chem.yale.edu/software.html

BOSS is compiled for linux using 32 bits libraries so it can not run in windows using WSL. Alternatively, you can use a virtual machine in windows such as virtualBox to install a linux distro (ubuntu, centos, ...) and run LigParGen.

1.1 - Set the BOSSdir enviromental variable:

bashrc

    export BOSSdir=PATH_TO_BOSS_DIRECTORY

cshrc
```
    setenv BOSSdir PATH_TO_BOSS_DIRECTORY
```
TIP: add this command line in your ~/.bashrc or ~/.cshrc file.

2 - Download and install conda (anaconda or minicoda):

wget https://repo.anaconda.com/archive/Anaconda3-2020.11-Linux-x86_64.sh

wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh

3 - Create and activate an enviroment for python3.7 (Everything was tested with 3.7):

conda create --name py37 python=3.7

conda activate py37

Optional TIP: add this line to your .bashrc/.cshrc,....

4- Install rdkit and openbabel in py37 enviroment

conda install -c rdkit rdkit
conda install -c conda-forge openbabel

5 - Download and install LigParGen

git clone https://github.com/Isra3l/ligpargen.git

pip install -e ligpargen

Optional: Check your installation by runing the tests included in the ligpargen folder.

cd ligpargen;python -m unittest test_ligpargen/test_ligpargen.py

TIP: Do not forget to activate your py37 enviroment before using LigParGen.

conda activate py37

USAGE

Input arguments list:

' -s ': SMILES input (Ex. -s 'CCCC' )
' -i ': PDB, MOL and MOL2 files + any input format supported by Open Babel (Ex. -i benzene.pdb )
' -n ': Molecule file name (Ex. -n benzene - it will produce benzene.gmx.gro, benzene.charmm.rtf,...)
' -p ': Folder for the output files (Ex. -p bnz )
' -r ': Residue name for the output files (Ex. -r BNZ)
' -c ': Molecule net charge (Ex. -c +1)
' -o ': Number of optimizations (Ex. -o 3)
' -cgen ': Charge model to be used - CM1A or CM1A-LBCC (Ex. -cgen CM1A)

For alchemical transformations:

' -sb ': SMILES input for molecule B (Ex. -sb 'CCCCO' )
' -ib ': PDB, MOL and MOL2 files + any input format supported by Open Babel (Ex. -ib phenol.pdb )
' -cb ': Molecule net charge (Ex. -cb +1)
' -ob ': Number of optimizations (Ex. -ob 3)
' -cgenb ': Charge model to be used - CM1A or CM1A-LBCC (Ex. -cgenb CM1A)

Notes:

CM1A is automatically scaled by 1.14 in neutral molecules.
CM1A-LBCC is just for neutral molecules and it is also scaled by 1.14.

Examples

Molecule template generation:

ligpargen -i phenol.pdb -n phenol -p phenol -r MOL -c 0 -o 0 -cgen CM1A
ligpargen -s 'c1ccc(cc1)O' -n phenol -p phenol -r MOL -c 0 -o 0 -cgen CM1A
ligpargen -s 'c1ccc(cc1)O'
ligpargen -s 'c1ccc(cc1)O' -n phenol -cgen CM1A-LBCC

Alchemical transformations:

ligpargen -i phenol.pdb -ib benzene.pdb -n phenolToBenzene -p phenol2bnz -r A2B -c 0 -o 0 -cgen CM1A -cb 0 -ob 1 -cgenb CM1A-LBCC
ligpargen -i phenol.pdb -ib benzene.pdb -n phenolToBenzene
ligpargen -s 'c1ccc(cc1)O' -sb 'c1ccccc1' -n phenol_benzene
ligpargen -s 'c1ccc(cc1)O' -sb 'c1ccccc1' -n phenol_benzene -o 0 -cgen CM1A -cb 0 -ob 1 -cgenb CM1A-LBCC
ligpargen -i phenol.pdb -sb 'c1ccccc1' -n phenol_benzene -o 0 -cgen CM1A -cb 0 -ob 1 -cgenb CM1A-LBCC

Default values:

-n : molecule
-p : Current working directory
-r : MOL
-c,-cb : Automatically determined from input molecule
-o,-ob : 0
-cgen,-cgenb: CM1A-LBCC for neutral molecules and CM1A for charged molecules.

For help use the -h flag:

ligpargen -h

Please do not forget to cite the following references:

LigParGen web server: an automatic OPLS-AA parameter generator for organic ligands
Leela S. Dodda Israel Cabeza de Vaca Julian Tirado-Rives William L. Jorgensen Nucleic Acids Research, Volume 45, Issue W1, 3 July 2017, Pages W331–W336
1.14*CM1A-LBCC: Localized Bond-Charge Corrected CM1A Charges for Condensed-Phase Simulations

ligpargen's People

Contributors

Stargazers

Watchers

ligpargen's Issues

Multiple errors in file reading

Hello,
I tried to run ligpargen on the phenol file found in the exemple dir.

ligpargen -i phenol_after_BOSS.pdb
Traceback (most recent call last):
  File "/Users/marcodigennaro/miniconda3/envs/ligpargen/bin/ligpargen", line 33, in <module>
    sys.exit(load_entry_point('LigPargen', 'console_scripts', 'ligpargen')())
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/ligpargen.py", line 496, in main
    molname=args.molname, workdir= args.path, debug= args.debug)
  File "<string>", line 17, in __init__
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/ligpargen.py", line 190, in __post_init__
    self.workdir, self.debug)
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/tools/boss.py", line 208, in run
    shutil.copyfile('plt.pdb', pdbFile)
  File "/Users/marcodigennaro/miniconda3/envs/ligpargen/lib/python3.7/shutil.py", line 120, in copyfile
    with open(src, 'rb') as fsrc:
FileNotFoundError: [Errno 2] No such file or directory: 'plt.pdb'

when using my file, generated from openbabel I get the error below:

COMPND    UNNAMED
AUTHOR    GENERATED BY OPEN BABEL 2.3.90
HETATM    1  C   UNL     1      -3.420   0.763   2.248  1.00  0.00           C1-
HETATM    2  F   UNL     1      -3.829  -0.542   2.440  1.00  0.00           F
HETATM    3  F   UNL     1      -2.919   1.269   3.432  1.00  0.00           F
HETATM    4  F   UNL     1      -4.497   1.527   1.841  1.00  0.00           F
HETATM    5  N   UNL     1      -0.643  -0.104   1.322  1.00  0.00           N2-
HETATM    6  C   UNL     1       2.054  -0.880   0.123  1.00  0.00           C
HETATM    7  F   UNL     1       2.678  -1.046  -1.098  1.00  0.00           F
HETATM    8  F   UNL     1       2.071  -2.077   0.814  1.00  0.00           F
HETATM    9  F   UNL     1       2.724   0.084   0.853  1.00  0.00           F
HETATM   10  S   UNL     1      -2.128   0.811   0.975  1.00  0.00           S
HETATM   11  S   UNL     1       0.337  -0.356  -0.139  1.00  0.00           S
HETATM   12  O   UNL     1      -1.660   2.243   0.672  1.00  0.00           O
HETATM   13  O   UNL     1      -2.633   0.324  -0.393  1.00  0.00           O
HETATM   14  O   UNL     1      -0.455  -1.360  -0.992  1.00  0.00           O
HETATM   15  O   UNL     1       0.237   0.938  -0.962  1.00  0.00           O
HETATM   16 Li   UNL     2      -0.577  -2.952   1.388  1.00  0.00          Li
CONECT    1    2    3    4   10
CONECT    2    1
CONECT    3    1
CONECT    4    1
CONECT    5   10   11
CONECT    6    7    8    9   11
CONECT    7    6
CONECT    8    6
CONECT    9    6
CONECT   10    1   12   13    5
CONECT   11    6   14   15    5
CONECT   12   10
CONECT   13   10
CONECT   14   11
CONECT   15   11
MASTER        0    0    0    0    0    0    0    0   16    0   16    0
END

ligpargen -i Li-TFSI.pdb
[18:18:26] Explicit valence for atom # 0 C, 5, is greater than permitted
Traceback (most recent call last):
  File "/Users/marcodigennaro/miniconda3/envs/ligpargen/bin/ligpargen", line 33, in <module>
    sys.exit(load_entry_point('LigPargen', 'console_scripts', 'ligpargen')())
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/ligpargen.py", line 496, in main
    molname=args.molname, workdir= args.path, debug= args.debug)
  File "<string>", line 17, in __init__
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/ligpargen.py", line 174, in __post_init__
    moleculeRDkit, newIndexToOriginalIndex, atomsNameOriginal, residueNameOriginal = utilities.generateRDkitMolecule(self.ifile, self.smile, self.workdir, molnameA, self.debug)
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/tools/utilities.py", line 187, in generateRDkitMolecule
    Chem.MolToPDBFile(molecule, sfile)
Boost.Python.ArgumentError: Python argument types in
    rdkit.Chem.rdmolfiles.MolToPDBFile(NoneType, str)
did not match C++ signature:
    MolToPDBFile(RDKit::ROMol mol, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char>> filename, int confId=-1, unsigned int flavor=0)

The error is clear here: Explicit valence for atom # 0 C, 5, is greater than permitted, so I modified the input file, (C1- and N2- become respectively C N )but I get another error.

Traceback (most recent call last):
  File "/Users/marcodigennaro/miniconda3/envs/ligpargen/bin/ligpargen", line 33, in <module>
    sys.exit(load_entry_point('LigPargen', 'console_scripts', 'ligpargen')())
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/ligpargen.py", line 496, in main
    molname=args.molname, workdir= args.path, debug= args.debug)
  File "<string>", line 17, in __init__
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/ligpargen.py", line 174, in __post_init__
    moleculeRDkit, newIndexToOriginalIndex, atomsNameOriginal, residueNameOriginal = utilities.generateRDkitMolecule(self.ifile, self.smile, self.workdir, molnameA, self.debug)
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/tools/utilities.py", line 189, in generateRDkitMolecule
    atomsIndexLstWithRightOrder = buildProperOrderAtomIndexLst(molecule, molname, workdir)
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/tools/utilities.py", line 543, in buildProperOrderAtomIndexLst
    addAtomsToTheList(molecule, atomsIndexLstWithRightOrder, heavyAtomsIndexLst)
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/tools/utilities.py", line 584, in addAtomsToTheList
    else: addAtomsToTheList(molecule, atomsIndexLst, heavyAtomsIndexLst)
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/tools/utilities.py", line 584, in addAtomsToTheList
    else: addAtomsToTheList(molecule, atomsIndexLst, heavyAtomsIndexLst)
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/tools/utilities.py", line 584, in addAtomsToTheList
    else: addAtomsToTheList(molecule, atomsIndexLst, heavyAtomsIndexLst)
  [Previous line repeated 988 more times]
  File "/Users/marcodigennaro/WORK/external_packages/ligpargen/ligpargen/tools/utilities.py", line 571, in addAtomsToTheList
    neighborsIndex = [x.GetIdx() for x in atom.GetNeighbors()]
RecursionError: maximum recursion depth exceeded while calling a Python object

Thank you for your support

ligpargen running problem

I tried with BOSS 4.9 and ligpargen compiled with conda python3.7, but when I run it, I get this error, below are the error output:

ligpargen -i part2.pdb -n phenol -p phenol -r MOL -c 0 -o 0 -cgen CM1A
[16:59:22] Explicit valence for atom # 25 O, 5, is greater than permitted
Traceback (most recent call last):
File "/public1/home/scb1100/miniconda3/envs/test/bin/ligpargen", line 33, in
sys.exit(load_entry_point('LigPargen', 'console_scripts', 'ligpargen')())
File "/public1/home/scb1100/ligpargen/ligpargen-main/ligpargen/ligpargen.py", line 494, in main
molname=args.molname, workdir= args.path, debug= args.debug)
File "", line 17, in init
File "/public1/home/scb1100/ligpargen/ligpargen-main/ligpargen/ligpargen.py", line 174, in post_init
moleculeRDkit, newIndexToOriginalIndex, atomsNameOriginal, residueNameOriginal = utilities.generateRDkitMolecule(self.ifile, self.smile, self.workdir, molnameA, self.debug)
File "/public1/home/scb1100/ligpargen/ligpargen-main/ligpargen/tools/utilities.py", line 165, in generateRDkitMolecule
atomsIndexLstWithRightOrder = buildProperOrderAtomIndexLst(molecule, molname, workdir)
File "/public1/home/scb1100/ligpargen/ligpargen-main/ligpargen/tools/utilities.py", line 513, in buildProperOrderAtomIndexLst
heavyAtomsIndexLst = [atom.GetIdx() for atom in molecule.GetAtoms() if atom.GetSymbol() != 'H']
AttributeError: 'NoneType' object has no attribute 'GetAtoms'

Any ideas？

Issue with charges causing nonphysical geometries

Description of the issue

When using LigParGen with the ACE-ALA-NME molecule, I am seeing large charges on the amide N-atoms and the H-atoms attached to these N-atoms. The charges on the N-atom are -0.997258/-1.077063 and for the H-atoms they are 0.504505/0.506935. These seem rather large in magnitude. When running a short gas-phase MD and minimizing the molecule, I am getting a structure that is ~ 8 kcal/mol lower in energy then the actual crystal structure conformation. The minimized structure looks nonphysical to me as the two amide H-atoms are pointing towards each other (see structure below). It seems to me that this issue is due to the large charges on the N/H atoms in the amide groups. My question is whether I am doing something wrong (see below for reproducing this) or if there might be bug somewhere. Any thoughts and insights on this are much appreciated.

As a site note, I am seeing the same results with CM1A-LBCC option.

Code/data to reproduce above issue

Generating force field with LigParGen:

ligpargen -i input_monomer0.pdb -cgen CM1A -c 0

Running short gas-phase MD and then minimize using openmm (v7.7.0):

def OPLS_LJ(system):
    from openmm import CustomNonbondedForce
    import numpy as np
        
    forces = {system.getForce(index).__class__.__name__: system.getForce(
        index) for index in range(system.getNumForces())}
    nonbonded_force = forces['NonbondedForce']
    lorentz = CustomNonbondedForce(
        '4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=sqrt(sigma1*sigma2); epsilon=sqrt(epsilon1*epsilon2)')
    lorentz.setNonbondedMethod(nonbonded_force.getNonbondedMethod())
    lorentz.addPerParticleParameter('sigma')
    lorentz.addPerParticleParameter('epsilon')
    lorentz.setCutoffDistance(nonbonded_force.getCutoffDistance())
    system.addForce(lorentz)
    LJset = {}
    for index in range(nonbonded_force.getNumParticles()):
        charge, sigma, epsilon = nonbonded_force.getParticleParameters(index)
        LJset[index] = (sigma, epsilon)
        lorentz.addParticle([sigma, epsilon])
        nonbonded_force.setParticleParameters(
            index, charge, sigma, epsilon * 0)
    for i in range(nonbonded_force.getNumExceptions()):
        (p1, p2, q, sig, eps) = nonbonded_force.getExceptionParameters(i)
        # ALL THE 12,13 and 14 interactions are EXCLUDED FROM CUSTOM NONBONDED
        # FORCE
        lorentz.addExclusion(p1, p2)
        if eps._value != 0.0:
            #print p1,p2,sig,eps
            sig14 = np.sqrt(LJset[p1][0] * LJset[p2][0])
            eps14 = np.sqrt(LJset[p1][1] * LJset[p2][1])
            nonbonded_force.setExceptionParameters(i, p1, p2, q, sig14, eps)
    return system
    
## Created by Leela S. Dodda for LigParGen Tutorials
## Date Aug 5, 2016
from openmm import app, KcalPerKJ
import openmm as mm
from openmm import unit as u
from sys import stdout, exit

def Minimize(simulation,iters=0):
    simulation.minimizeEnergy(maxIterations=iters)
    position = simulation.context.getState(getPositions=True).getPositions()
    energy = simulation.context.getState(getEnergy=True).getPotentialEnergy()
    app.PDBFile.writeFile(simulation.topology, position,
                          open('gasmin.pdb', 'w'))
    print ('Energy at Minima is %3.3f kcal/mol' % (energy._value * KcalPerKJ))
    return simulation

temperature = 298.15 * u.kelvin
pdb = app.PDBFile('input_monomer0.pdb')

modeller = app.Modeller(pdb.topology, pdb.positions)
forcefield = app.ForceField('input_monomer0.openmm.xml')

system = forcefield.createSystem(
    modeller.topology, nonbondedMethod=app.NoCutoff,  constraints=None)
system = OPLS_LJ(system)
integrator = mm.LangevinIntegrator(
    temperature, 1 / u.picosecond,  0.0005 * u.picoseconds)
simulation = app.Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)
simulation = Minimize(simulation,1000)

simulation.reporters.append(app.PDBReporter('gas_output.pdb', 1000))
simulation.reporters.append(app.StateDataReporter('data.txt', 1000, progress=True, temperature=True, potentialEnergy=True, density=True,totalSteps=10000,speed=True))
simulation.step(100000)

simulation = Minimize(simulation,1000)

openmm xml and pdb files:
input_monomer0.pdb.gz
input_monomer0.openmm.xml.gz

CM1A-LBCC does not work

I tried with BOSS 4.9 and ligpargen compiled with conda python3.7, the charge generation algorithm only works for CM1A but not CM1A-LBCC on my system (Ubuntu 18.04), below are the comparison and error output:

ligpargen -s 'c1ccc(cc1)O' -sb 'c1ccccc1' -n phenol_benzene -o 0 -cgen CM1A -cb 0 -ob 1 -cgenb CM1A
LigParGen finished succesfully!

ligpargen -s 'c1ccc(cc1)O' -sb 'c1ccccc1' -n phenol_benzene -o 0 -cgen CM1A -cb 0 -ob 1 -cgenb CM1A-LBCC
Traceback (most recent call last):
File "/data/soft/anaconda3/envs/py37/bin/ligpargen", line 33, in
sys.exit(load_entry_point('LigPargen', 'console_scripts', 'ligpargen')())
File "/data/soft/ligpargen/ligpargen/ligpargen.py", line 494, in main
molname=args.molname, workdir= args.path, debug= args.debug)
File "", line 17, in init
File "/data/soft/ligpargen/ligpargen/ligpargen.py", line 222, in post_init
moleculeB = Molecule.fromBOSS(zmatNameB, outFileB, pdbFileB, moleculeB.shiftX, moleculeB.shiftY, moleculeB.shiftZ)
File "/data/soft/ligpargen/ligpargen/topology/Molecule.py", line 141, in fromBOSS
atoms, numberOfStructuralDummyAtoms = cls._getAtoms(cls, zmatData, pdbfile)
File "/data/soft/ligpargen/ligpargen/topology/Molecule.py", line 460, in _getAtoms
atomsLines = re.search(r'BOSS(.*?)Geometry', zmatData, re.DOTALL).group().splitlines()[1:-1]
AttributeError: 'NoneType' object has no attribute 'group'

Any ideas?

BOSS in 32 bit

There are a lot of issues with using ligpargen, as I understand it is related to BOSS coming as a 32 bit compiled code.
i.e. it will NOT run on a Mac (thanks for the recent 64 bit upgrade!) or any Linux that is 64 bit.

PS - to check what is your machine use uname -a

if you run BOSS directy, you will see errors like:

./boss/autozmat -x -i boss -o pdb -p tmppar
./boss/autozmat: Exec format error. Binary file not executable.

if you run ligpargen, you will get misguiding errors like:
FileNotFoundError: [Errno 2] No such file or directory: 'plt.pdb'
or
ERROR! Open Babel is not installed on your work station.

as mentioned in issue #1

I guess it all just boils to getting BOSS for 64 bit.

@Isra3l also tagging @leelasd as another ligpargen guru :)
**Is it possible to get BOSS as 64 bit or non-compiled? **

Thank you!

Which local version of ligpargen does the website version correspond to

Error when processing CM1A-LBCC charges

env.yaml.zip
When running the following command:

ligpargen -s 'c1ccc(cc1)O' -n phenol -cgen CM1A-LBCC

I am getting the following error message:

Traceback (most recent call last):
  File "/home/tobias/progs/anaconda3/envs/xtalmd/bin/ligpargen", line 33, in <module>
    sys.exit(load_entry_point('LigPargen', 'console_scripts', 'ligpargen')())
  File "/data/tobias_data/progs/ligpargen/ligpargen/ligpargen.py", line 496, in main
    molname=args.molname, workdir= args.path, debug= args.debug)
  File "<string>", line 17, in __init__
  File "/data/tobias_data/progs/ligpargen/ligpargen/ligpargen.py", line 192, in __post_init__
    moleculeA = Molecule.fromBOSS(zmatName, outFile, pdbFile, moleculeA.shiftX, moleculeA.shiftY, moleculeA.shiftZ)
  File "/data/tobias_data/progs/ligpargen/ligpargen/topology/Molecule.py", line 141, in fromBOSS
    atoms, numberOfStructuralDummyAtoms = cls._getAtoms(cls, zmatData, pdbfile)
  File "/data/tobias_data/progs/ligpargen/ligpargen/topology/Molecule.py", line 460, in _getAtoms
    atomsLines = re.search(r'BOSS(.*?)Geometry', zmatData, re.DOTALL).group().splitlines()[1:-1]
AttributeError: 'NoneType' object has no attribute 'group'

The temporary out file in the working directory contains

Unavailable QM Choice: CM1L - Aborted

This might be a bug related to how the input is parsed to BOSS. Any help is greatly appreciated. For reference, running ligpargen -s 'c1ccc(cc1)O' -n phenol -cgen CM1A works just fine.

Attaching yaml file from my conda environment.

--Tobias

Dihedral angles of linear groups on 5-membered rings

Compounds with a linear group (N≡C– or HC≡C–) attached to a 5-membered ring cause simulations to crash (tested with Gromacs). Before the crash, the rings become distorted. Compounds with a 6-membered ring are unaffected. By comparing the parameters generated for compounds with 5- and 6-membered rings, the issue seems to be with the coefficients of proper dihedrals that terminate with the nitrogen or carbon. In the case of 5-membered rings, the coefficients are non-zero. In some compounds, the torsion angle that terminates in the N≡C– carbon also has different coefficients than those in compounds with 6-membered rings. Making modifications to these parameters seems to fix the issue. The problem can be reproduced with the following compounds:

Mol_01: o1cc(C#C)cc1
[ dihedrals ]
; PROPER DIHEDRAL ANGLES
; ai aj ak al funct c0 c1 c2 c3 c4 c5
7 6 3 4 3 9.079 0.000 -9.079 -0.000 -0.000 0.000
6 3 4 5 3 9.079 0.000 -9.079 -0.000 -0.000 0.000

Mol_02: o1cc(C#N)cc1
7 6 3 4 3 9.079 0.000 -9.079 -0.000 -0.000 0.000
6 3 4 5 3 9.079 0.000 -9.079 -0.000 -0.000 0.000

Mol_03: o1c2ccccc2c(C#N)c1
11 8 9 10 3 30.334 0.000 -30.334 -0.000 -0.000 0.000

Mol_04: o1cccc1C#Cc2ccccc2
7 6 5 1 3 9.079 0.000 -9.079 -0.000 -0.000 0.000

Failing installation

Hi, I've installed with the following:

conda create --name ligpargen python=3.7
conda activate ligpargen
[obtained BOSS from Jorgensen  http://zarbi.chem.yale.edu/software.html]
[echo "export BOSSdir=/XXX/XXX/XXX/bin/boss" >> ~/.bashrc] where XXX is my path 
conda install -c rdkit rdkit
conda install -c conda-forge openbabel

But then, running:

cd ligpargen;python -m unittest test_ligpargen/test_ligpargen.py

gives:

....................F..
======================================================================
FAIL: test_openMM_outputs (test_ligpargen.test_ligpargen.test_ligpargen)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/XXX/XXX/bin/ligpargen/test_ligpargen/test_ligpargen.py", line 571, in test_openMM_outputs
    self.assertEqual(benzene_xml, benzene_xml_golden)
AssertionError: Lists differ: ['<!-[261 chars]class="C800" element="C" mass="12.011"/>', '<T[13538 chars]ld>'] != ['<!-[261 chars]class"C800" element="C" mass="12.011"/>', '<Ty[13527 chars]ld>']

First differing element 6:
'<Type name="opls_800" class="C800" element="C" mass="12.011"/>'
'<Type name="opls_800" class"C800" element="C" mass="12.011"/>'

Diff is 16297 characters long. Set self.maxDiff to None to see it.

----------------------------------------------------------------------
Ran 23 tests in 0.083s

FAILED (failures=1)

And running:

ligpargen -s 'c1ccc(cc1)O'

gives:

Traceback (most recent call last):
  File "/XXX/.conda/envs/ligpargen/bin/ligpargen", line 33, in <module>
    sys.exit(load_entry_point('LigPargen', 'console_scripts', 'ligpargen')())
  File "/XXX/XXX/bin/ligpargen/ligpargen/ligpargen.py", line 496, in main
    molname=args.molname, workdir= args.path, debug= args.debug)
  File "<string>", line 17, in __init__
  File "/XXX/XXX/bin/ligpargen/ligpargen/ligpargen.py", line 192, in __post_init__
    moleculeA = Molecule.fromBOSS(zmatName, outFile, pdbFile, moleculeA.shiftX, moleculeA.shiftY, moleculeA.shiftZ)
  File "/XXX/XXX/bin/ligpargen/ligpargen/topology/Molecule.py", line 141, in fromBOSS
    atoms, numberOfStructuralDummyAtoms = cls._getAtoms(cls, zmatData, pdbfile)
  File "/XXX/XXX/bin/ligpargen/ligpargen/topology/Molecule.py", line 460, in _getAtoms
    atomsLines = re.search(r'BOSS(.*?)Geometry', zmatData, re.DOTALL).group().splitlines()[1:-1]
AttributeError: 'NoneType' object has no attribute 'group'

Any help appreciated :-)

Aromatic rings warped with OPLS parametrized molecules

Hi,

When runnning a minimization invoking Tinker9 minimize with an organic molecule parametrized with OPLS (Ligpargen command-line v.2.1), I notice aromatic rings warped when the molecule is zwitterionic, like the following one...

So the result is this:

I´ve checked the .key file, comparing it with other non-zwitterionic molecules, and the parametrisation looks fine. Aromatic atoms are accordingly labelled with CA atom types, torsional values ok, etc. I haven´t escrutinize the file carefully, just a superficial check. LigParGen paraetrisation was run with the following arguments:

ligpargen -i mol.pdb -r MOL -c 1 -o 0 -n mol -cgen CM1A

So I´m wondering either this is an artifact of the parametrisation

or either

an issue of the oplsa-aa force field contained in tinker. I would appreciate any ideas in regards to this.

Thank you kindly,

E77.

A error happened when running ligpargen

Dear developer,

I installed ligpargen as suggested in the README.md. And when I ran the test, it seemed ok. But when I typed the command "ligpargen -s 'c1ccc(cc1)O'" in the terminal, a error happened, which said the file 'plt.pdb' doesn't exist. Now, the error also shows up even when I do not use the smiles string as input. I am writing to ask for help to get ligpargen installed correctly. Can you give me some advice to deal with it?

Best wishes,
Hang LYU

PS:
Test

The error

BOSS 5.0 / 64 bits executable

A big thumb-up for your great work on LigParGen, I’m currently doing an experiment about Zinc related dynamic simulation using LigParGen v2.1.

May I ask for BOSS 5.0 executable for LigParGen v2.1? I have sent emails, but Prof. Jorgensen might be busy without replying
Considering that I’m doing my experiment on Ubuntu/Mac which is based on 64 bits CPU, can you also provide 64 bits executables? It seems that creating virtual environment for 32 bit executables is not so proper for MacOS

Cannot run charge algorithm CM1A-LBCC

Hi,

When testing with small molecule this works:

$ ligpargen -i ribociclib.poltype.opt.sdf -n ribociclib -p ribociclib -r 6ZZ -c 0 -o 0 -cgen CM1A
LigParGen finished succesfully!

but using CM1A-LBCC gives traceback

$ ligpargen -i ribociclib.poltype.opt.sdf -n ribociclib -p ribociclib -r 6ZZ -c 0 -o 0 -cgen CM1A-LBCC
Traceback (most recent call last):
  File "/opt/conda/envs/py37/bin/ligpargen", line 33, in <module>
    sys.exit(load_entry_point('LigPargen', 'console_scripts', 'ligpargen')())
  File "/opt/ligpargen/ligpargen/ligpargen.py", line 496, in main
    molname=args.molname, workdir= args.path, debug= args.debug)
  File "<string>", line 17, in __init__
  File "/opt/ligpargen/ligpargen/ligpargen.py", line 192, in __post_init__
    moleculeA = Molecule.fromBOSS(zmatName, outFile, pdbFile, moleculeA.shiftX, moleculeA.shiftY, moleculeA.shiftZ)
  File "/opt/ligpargen/ligpargen/topology/Molecule.py", line 141, in fromBOSS
    atoms, numberOfStructuralDummyAtoms = cls._getAtoms(cls, zmatData, pdbfile)
  File "/opt/ligpargen/ligpargen/topology/Molecule.py", line 460, in _getAtoms
    atomsLines = re.search(r'BOSS(.*?)Geometry', zmatData, re.DOTALL).group().splitlines()[1:-1]
AttributeError: 'NoneType' object has no attribute 'group'

Test example failes and I can't do successfully even the simplest example

Hello, everyone!

I installed LigParGen into a Singularity container based on condaforge/miniforge:latest base image. I followed the LigParGen installation instructions strictly. However :

Test example fails, namely python -m unittest test_ligpargen/test_ligpargen.p gives :

.............F.........
======================================================================
FAIL: test_atoms (test_ligpargen.test_ligpargen.test_ligpargen)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/root/Work/ligpargen/ligpargen/test_ligpargen/test_ligpargen.py", line 134, in test_atoms
    self.assertEqual(atoms_benzene, golden_atoms_benzene)
AssertionError: Lists differ: ['1  [536 chars]4792045  13.471357010453586  MOL  1  0  0  0  [1567 chars] H1'] != ['1  [536 chars]4792063  13.471357010453614  MOL  1  0  0  0  [1570 chars] H1']

First differing element 4:
'5  C[46 chars]4792045  13.471357010453586  MOL  1  0  0  0  [46 chars]  C6'
'5  C[46 chars]4792063  13.471357010453614  MOL  1  0  0  0  [46 chars]  C6'

Diff is 4037 characters long. Set self.maxDiff to None to see it.

----------------------------------------------------------------------
Ran 23 tests in 0.575s

FAILED (failures=1)

the command ligpargen -s 'CC' for the simplest example I can imagine, i.e. ethane, gives the output :

Traceback (most recent call last):
  File "/opt/conda/envs/py37/bin/ligpargen", line 33, in <module>
    sys.exit(load_entry_point('LigPargen', 'console_scripts', 'ligpargen')())
  File "/root/Work/ligpargen/ligpargen/ligpargen/ligpargen.py", line 496, in main
    molname=args.molname, workdir= args.path, debug= args.debug)
  File "<string>", line 17, in __init__
  File "/root/Work/ligpargen/ligpargen/ligpargen/ligpargen.py", line 190, in __post_init__
    self.workdir, self.debug)
  File "/root/Work/ligpargen/ligpargen/ligpargen/tools/boss.py", line 208, in run
    shutil.copyfile('plt.pdb', pdbFile)
  File "/opt/conda/envs/py37/lib/python3.7/shutil.py", line 120, in copyfile
    with open(src, 'rb') as fsrc:
FileNotFoundError: [Errno 2] No such file or directory: 'plt.pdb'

What could be the reason? I digged into the code and saw that, indeed, I didn't get where the file 'plt.pdb' was created that is to be copied further into pdbFile.

UPD

Having read the previous issues, I added the 32-bit support to my 64-bit machine as suggested :

dpkg --add-architecture i386
apt-get install libc6:i386 libncurses5:i386 libstdc++6:i386

then, the output of p. 2 changed to :

Traceback (most recent call last):
  File "/opt/conda/envs/py37/bin/ligpargen", line 33, in <module>
    sys.exit(load_entry_point('LigPargen', 'console_scripts', 'ligpargen')())
  File "/root/Work/ligpargen/ligpargen/ligpargen/ligpargen.py", line 496, in main
    molname=args.molname, workdir= args.path, debug= args.debug)
  File "<string>", line 17, in __init__
  File "/root/Work/ligpargen/ligpargen/ligpargen/ligpargen.py", line 192, in __post_init__
    moleculeA = Molecule.fromBOSS(zmatName, outFile, pdbFile, moleculeA.shiftX, moleculeA.shiftY, moleculeA.shiftZ)
  File "/root/Work/ligpargen/ligpargen/ligpargen/topology/Molecule.py", line 154, in fromBOSS
    atoms, numberOfStructuralDummyAtoms = cls._getAtoms(cls, zmatData, pdbfile)
  File "/root/Work/ligpargen/ligpargen/ligpargen/topology/Molecule.py", line 473, in _getAtoms
    atomsLines = re.search(r'BOSS(.*?)Geometry', zmatData, re.DOTALL).group().splitlines()[1:-1]
AttributeError: 'NoneType' object has no attribute 'group'

whereas the test example still fails. I found in the closed issues that someone faced the same problem earlier and @Isra3l assumed this was caused by the usage of BOSS v. 4.9 and suggested to try BOSS v. 5.0 instead. Indeed, I'm using BOSS v. 4.9, since I simply can't get v 5.0, because Prof. Jorgensen doesn't reply to emails. It's not only me who contacted Prof. Jorgensen, but also a few colleagues of mine. Does anybody know how to get this BOSS v. 5.0???

Unknown error for ligand

Hi,

WHen submitting the following pdb, LigParGen crashes with following error:

Traceback (most recent call last):
File "/home/e77/.conda/envs/ligpargen/bin/ligpargen", line 33, in
sys.exit(load_entry_point('LigPargen', 'console_scripts', 'ligpargen')())
File "/home/e77/ligpargen/ligpargen/ligpargen.py", line 496, in main
molname=args.molname, workdir= args.path, debug= args.debug)
File "", line 17, in init
File "/home/e77/ligpargen/ligpargen/ligpargen.py", line 192, in post_init
moleculeA = Molecule.fromBOSS(zmatName, outFile, pdbFile, moleculeA.shiftX, moleculeA.shiftY, moleculeA.shiftZ)
File "/home/e77/ligpargen/ligpargen/topology/Molecule.py", line 149, in fromBOSS
anglesVariable, anglesAdditional = cls._getAngles(zmatData, outFileData, atoms, serial2IndexAtoms)
File "/home/e77/ligpargen/ligpargen/topology/Molecule.py", line 579, in _getAngles
atomA = atoms[serial2IndexAtoms[int(angleF[0])]]
ValueError: invalid literal for int() with base 10: 'Force'

ANy ideas of the issue? Tried with server version and again is happening...

Thank you very much,

E77

lba0000025813_lig_old.txt

Failure due to `Unavailable QM Choice: CM1L - Aborted`.

Thanks for your great job on this pretty tool!

Sadly, I am blocked in generating force field files for a molecule using CM1A-LBCC charge model.
The command to reproduce error is:

ligpargen -s 'O=C(NC1CCCCC1)[C@H](C1CCCCC1)n1c2ccccc2nc1c1ccc(OC)cc1OC' -n FXR_6 -r UNK -p tmp

The failure comes from BOSS with error recorded in file tmp/out:

Unavailable QM Choice: CM1L - Aborted

However, if I uploaded the molecule to LigParGen web server, no error happened and force field files were generated.
Could you please help debug & fix it?

Thanks a lot!

OpenMM Xml atom order Improper Dihedrals

The order of the atoms in the Improper Dihedral definition of OpenMM xml files may be incorrect. AFAIK OpenMM force field xml input expects the central atom at position 1. Attached is a case that causes the problem.

Running

ligpargen -p test -i mol_0.mol -cgen CM1A -c 0

gets the following improper dihedral definition in the OpenMM xml file:

<Improper class1="C800" class2="C801" class3="N804" class4="O806" k1="0.000000" k2="43.932000" k3="0.000000" k4="0.000000" periodicity1="1" periodicity2="2" periodicity3="3" periodicity4="4" phase1="0.00" phase2="3.141592653589793" phase3="0.00" phase4="3.141592653589793"/>

However, the atom order should be

<Improper class1="C801" class2="C800" class3="N804" class4="O806" k1="0.000000" k2="43.932000" k3="0.000000" k4="0.000000" periodicity1="1" periodicity2="2" periodicity3="3" periodicity4="4" phase1="0.00" phase2="3.141592653589793" phase3="0.00" phase4="3.141592653589793" />

When parsing the force field to openmm, the four Improper dihedrals in the molecule cannot be found. It should have 35 Torsions in total. Check out:

import openmm
from openmm import app

forcefield = app.ForceField("test/mol_0.openmm.xml")
pdb        = app.PDBFile("mol_0.pdb")
system = forcefield.createSystem(pdb.topology)
tor_force = system.getForce(2)
print(
    "Number of torsions:",
    tor_force.getNumTorsions()
)

I think there is an easy fix, but I am not sure if it is going to break anything.

Attached file: mol_0.pdb
mol_0.pdb.txt

Discrepancies between LigParGen server and ligpargen 2.1 in xyz/key files

Hi,

When submitting and contrasting the two results of parametrization of the same molecule and geometry with server version against command-line one, discrepancies appear in the tandem tinker type files:
lig.txt
lig.comline.xyz.txt
lig_server.key.txt
lig_server.xyz.txt
lig.comline.key.txt

Server options: upload pdb file + CM1A option + 0 optimizations + charge (+1)
ComLine options: ligpargen -i lig.pdb -cgen CM1A -r MOL -c 1 -o 0

There are discrepancies in both xyz and jey files when comparing. Moreover, tinker analyze routine works fine for server outputs but not for comline ones, returning a error of connectivity between first atoms.

Any hints about this? Not checked otehr ligpargen 2.1 outputs as i require tinker ones.

HOpe this helps...

Thank you,

E77

Generating parameters for water

When running the following command to obtain parameters for water

ligpargen -s O -cgen CM1A

I am getting weird angle parameters and partial charges (consistent over all parameter files for all codes supported). The angle0 parameter is much too low and the partial charges do not sum to zero. This is the openmm xml file:

<ForceField>
<AtomTypes>
<Type name="opls_800" class"O800" element="O" mass="15.999"/>
<Type name="opls_801" class"H801" element="H" mass="1.008"/>
<Type name="opls_802" class"H802" element="H" mass="1.008"/>
</AtomTypes>
<Residues>
<Residue name="MOL">
<Atom name="O00" type="opls_800"/>
<Atom name="H01" type="opls_801"/>
<Atom name="H02" type="opls_802"/>
<Bond from="1" to="0"/>
<Bond from="2" to="0"/>
</Residue>
</Residues>
<HarmonicBondForce>
<Bond class1="H801" class2="O800" length="0.094500" k="462750.400000"/>
<Bond class1="H802" class2="O800" length="0.094500" k="462750.400000"/>
</HarmonicBondForce>
<HarmonicAngleForce>
<Angle class1="H802" class2="O800" class3="H801" angle="1.433264" k="523.000000"/>
</HarmonicAngleForce>
<PeriodicTorsionForce>
</PeriodicTorsionForce>
<NonbondedForce coulomb14scale="0.5" lj14scale="0.5"/>
<Atom type="opls_800" charge="-0.838675" sigma="0.315365" epsilon="0.648520"/>
<Atom type="opls_801" charge="0.419338" sigma="0.000000" epsilon="0.000000"/>
<Atom type="opls_802" charge="0.419337" sigma="0.000000" epsilon="0.000000"/>
</NonbondedForce>
</ForceField>

I am not too familiar with the OPLS-AA force field, so I could be missing something obvious, but I think the perceived parameters should at least be physically reasonable. I know there is TIP3P, but in my case rigid bulk water models won't work that well so I am interested in getting parameters for a single water molecule.

If I run the above for say ammonia or hydrogen sulfid

ligpargen -s N -cgen CM1A

ligpargen -s S -cgen CM1A

I am getting reasonable parameters.

Any help is greatly appreciated. If this is indeed a bug I can try to help fixing it, but would need some guidance on where to look at :-)

--Tobias

Net charge is not zero due to rounding

Hi,

This error is actually discovered in either GROMACs or probably as well in APBS because for neutral molecules ligpargen will actually have slightly positive or negative charges in the *gmx.itp and *.pqr files due to rounding, i.e. total net charge for molecule is + or - (0.0001 to 0.0005). Is there a way we can correct this in ligpargen?

See:
https://gromacs.bioexcel.eu/t/ligpargen-parameterized-ligand-for-opls-aa-m-introduces-a-non-integer-charger/428/5
https://gromacs.bioexcel.eu/t/total-charge-error-in-itp-files/486
https://gromacs.bioexcel.eu/t/ligpargen-toluene-itp-file/4023

input file: does not exit

As the title,though all my input files are in the same folder and the ligpargen can be used,it still run out with input file: does not exit

Here is my command

ligpargen -i TSFI.pdb -n TFSI -P TFSI -r MOL -c -1 -0 0 -cgen CM1A