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/chilly.png ../_images/chilly.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.

 1def main(argv=[__name__]):
 2    if len(argv) != 3:
 3        oechem.OEThrow.Usage("%s <infile> <outfile>" % argv[0])
 4
 5    ifs = oechem.oemolistream()
 6    if not ifs.open(argv[1]):
 7        oechem.OEThrow.Fatal("Unable to open %s for reading" % argv[1])
 8
 9    if not oechem.OEIs3DFormat(ifs.GetFormat()):
10        oechem.OEThrow.Fatal("Invalid input format: need 3D coordinates")
11
12    ofs = oechem.oemolostream()
13    if not ofs.open(argv[2]):
14        oechem.OEThrow.Fatal("Unable to open %s for writing" % argv[2])
15
16    if ofs.GetFormat() not in [oechem.OEFormat_MOL2, oechem.OEFormat_OEB]:
17        oechem.OEThrow.Error("MOL2 or OEB output file is required!")
18
19    omegaOpts = oeomega.OEOmegaOptions()
20    omegaOpts.SetIncludeInput(True)
21    omegaOpts.SetCanonOrder(False)
22    omegaOpts.SetSampleHydrogens(True)
23    eWindow = 15.0
24    omegaOpts.SetEnergyWindow(eWindow)
25    omegaOpts.SetMaxConfs(800)
26    omegaOpts.SetRMSThreshold(1.0)
27
28    omega = oeomega.OEOmega(omegaOpts)
29    
30    for mol in ifs.GetOEMols():
31        if omega(mol):
32            oequacpac.OEAssignCharges(mol, oequacpac.OEAM1BCCELF10Charges())
33            conf = mol.GetConf(oechem.OEHasConfIdx(0))
34            absFCharge = 0
35            sumFCharge = 0
36            sumPCharge = 0.0
37            for atm in mol.GetAtoms():
38                sumFCharge += atm.GetFormalCharge()
39                absFCharge += abs(atm.GetFormalCharge())
40                sumPCharge += atm.GetPartialCharge()
41            print("{}: {} formal charges give total charge {}"
42                  "; sum of partial charges {:5.4f}".format(mol.GetTitle(), absFCharge,
43                                                            sumFCharge, sumPCharge))
44            oechem.OEWriteMolecule(ofs, conf)
45        else:
46            print("Failed to generate conformation(s) for molecule %s" % mol.GetTitle())
47
48    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