Generating Canonical AM1-BCC Charges

Problem

You want to generate general, canonical AM1-BCC charges from just one input structure. The atomic partial charges need to be appropriate for a wide spectrum of energy-based methods including minimization, posing in an active site, SZMAP, EON, or molecular dynamics.

Ingredients

  • Omega TK - conformer generation toolkit

Note

Requires OpenEye toolkits version 2017.Feb or later.

Difficulty Level

../_images/chilly4.png ../_images/chilly4.png

Solution

The script below allows us to avoid the two big problems that easily happen from AM1-BCC charging a single structure: asymmetric charges and distorted charges (due to short-range polar interactions). In this recipe, 800 conformations are generated, and then the best one is chosen automatically and used to generate well-behaved, symmetric AM1-BCC charges. The input structure is then written out with these atomic partial charges, ready for use.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def main(argv=[__name__]):
    if len(argv) != 3:
        OEThrow.Usage("%s <infile> <outfile>" % argv[0])

    ifs = oemolistream()
    if not ifs.open(argv[1]):
        OEThrow.Fatal("Unable to open %s for reading" % argv[1])

    if not OEIs3DFormat(ifs.GetFormat()):
        OEThrow.Fatal("Invalid input format: need 3D coordinates")

    ofs = oemolostream()
    if not ofs.open(argv[2]):
        OEThrow.Fatal("Unable to open %s for writing" % argv[2])

    if ofs.GetFormat() not in [OEFormat_MOL2, OEFormat_OEB]:
        OEThrow.Error("MOL2 or OEB output file is required!")

    omega = OEOmega()
    omega.SetIncludeInput(True)
    omega.SetCanonOrder(False)
    omega.SetSampleHydrogens(True)
    eWindow = 15.0
    omega.SetEnergyWindow(eWindow)
    omega.SetMaxConfs(800)
    omega.SetRMSThreshold(1.0)

    for mol in ifs.GetOEMols():
        if omega(mol):
            OEAssignCharges(mol, OEAM1BCCELF10Charges())
            conf = mol.GetConf(OEHasConfIdx(0))
            absFCharge = 0
            sumFCharge = 0
            sumPCharge = 0.0
            for atm in mol.GetAtoms():
                sumFCharge += atm.GetFormalCharge()
                absFCharge += abs(atm.GetFormalCharge())
                sumPCharge += atm.GetPartialCharge()
            OEThrow.Info("%s: %d formal charges give total charge %d ; Sum of Partial Charges %5.4f"
                         % (mol.GetTitle(), absFCharge, sumFCharge, sumPCharge))
            OEWriteMolecule(ofs, conf)
        else:
            OEThrow.Warning("Failed to generate conformation(s) for molecule %s" % mol.GetTitle())

    return 0

Download code

can_am1_bcc.py and charges.mol2 supporting data file

Usage:

prompt > python3 can_am1_bcc.py charges.mol2 cancharges.mol2

Running the above command will generate the following output:

residue 'BFS': 0 formal charges give total charge 0 ; Sum of Partial Charges -0.0000
ceftolozane: 4 formal charges give total charge 0 ; Sum of Partial Charges -0.0000

See Also in Omega TK Manual

API

See Also in Quacpac TK Manual

API