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

Difficulty level

../_images/chilly5.png ../_images/chilly5.png

Download

Download code

can_am1_bcc.py

See also the Usage subsection.

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
46
def main(argv=[__name__]):
    if len(argv) != 3:
        oechem.OEThrow.Usage("%s <infile> <outfile>" % argv[0])

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

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

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

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

    omega = oeomega.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):
            oequacpac.OEAssignCharges(mol, oequacpac.OEAM1BCCELF10Charges())
            conf = mol.GetConf(oechem.OEHasConfIdx(0))
            absFCharge = 0
            sumFCharge = 0
            sumPCharge = 0.0
            for atm in mol.GetAtoms():
                sumFCharge += atm.GetFormalCharge()
                absFCharge += abs(atm.GetFormalCharge())
                sumPCharge += atm.GetPartialCharge()
            print("{}: {} formal charges give total charge {}"
                  "; sum of partial charges {:5.4f}".format(mol.GetTitle(), absFCharge,
                                                            sumFCharge, sumPCharge))
            oechem.OEWriteMolecule(ofs, conf)
        else:
            print("Failed to generate conformation(s) for molecule %s" % mol.GetTitle())

    return 0

Usage

Usage

can_am1_bcc.py with charges.mol2 supporting data file

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