Visualizing Molecular Dipole Moment

Problem

You want to visualize the molecular dipole moment that is a good indicator of the overall polarity of a molecule. See example in Figure 1.

../_images/dipole2img.png

Figure 1. Example of visulaizing the molecular dipole moment

Ingredients

Note

Requires OpenEye toolkits version 2013.Jun or later.

Difficulty level

../_images/chilly6.png ../_images/chilly6.png ../_images/chilly6.png

Solution

The CalculateDipoleMoment function, after assigning the partial charges of the atoms by calling the OEMMFF94PartialCharges function, calculates the center of the positive and negative charges (lines 10-22) using the partial atom charges as a weight. Then its calculates:

  • the magnitude of the dipole (lines 30-38)
  • the (unit normalized) direction of the dipole (lines 40-42)
  • the center of the dipole (lines 44-46)

The dipole inner product is then calculated for each atom.. This is the vector product of the unit dipole with the vector from the dipole center to the atom center (lines 50-52). These numbers are then normalized such that the most positive number is equal to 1.0 and the most negative to -1.0 (lines 49-53). After the normalization, the number are attached to the relevant atom as generic data with the given tag (lines 60-61). For each bond a value is also calculated by averaging the values of its end atoms See more details in the Discussion section.

Note

The dipole inner product numbers represent both the distance and correlation with direction of the dipole, i.e. an atom that is in the positive direction of the dipole will be positive and its value will be bigger the further away it is from the center of the dipole. Atoms that are in a plane orthogonal to the dipole, passing through the dipole center, will have a number close to zero.

 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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
def CalculateDipoleMoment(mol, itag):

    OEMMFFAtomTypes(mol)
    OEMMFF94PartialCharges(mol)

    ncrg = 0.0
    pcrg = 0.0
    ncen = [0.0, 0.0, 0.0]
    pcen = [0.0, 0.0, 0.0]
    for atom in mol.GetAtoms():
        charge = atom.GetPartialCharge()
        x, y, z = mol.GetCoords(atom)
        if charge < 0.0:
            ncen[0] -= charge * x
            ncen[1] -= charge * y
            ncen[2] -= charge * z
            ncrg -= charge
        elif charge > 0.0:
            pcen[0] += charge * x
            pcen[1] += charge * y
            pcen[2] += charge * z
            pcrg += charge

    for idx in range(0, 3):
        ncen[idx] = ncen[idx] / ncrg
        pcen[idx] = pcen[idx] / pcrg
    if pcrg > ncrg:
        pcrg = ncrg

    dpmag = 0.0
    for i in range(0, 3):
        dpmag += (ncen[i] - pcen[i]) * (ncen[i] - pcen[i])

    if abs(dpmag) < 0.001:
        # no dipole moment
        return False

    dpmag *= 4.80324 * pcrg

    dipdir = [0.0, 0.0, 0.0]
    for idx in range(0, 3):
        dipdir[idx] = 4.80324 * pcrg * ((pcen[idx] - ncen[idx]) / dpmag)

    dipcen = [0.0, 0.0, 0.0]
    for idx in range(0, 3):
        dipcen[idx] = 0.5 * (pcen[idx] + ncen[idx])

    dpip = OEFloatArray(mol.GetMaxAtomIdx())

    for atom in mol.GetAtoms():
        x, y, z = mol.GetCoords(atom)
        dpip[atom.GetIdx()] = sum([(x - dipcen[0]) * dipdir[0],
                                   (y - dipcen[1]) * dipdir[1],
                                   (z - dipcen[2]) * dipdir[2]])

    max_dpip = max(dpip)
    min_dpip = min(dpip)

    dpip = [-i / min_dpip if i < 0.0 else i for i in dpip]
    dpip = [+i / max_dpip if i > 0.0 else i for i in dpip]

    for atom, dipole in zip(mol.GetAtoms(), dpip):
        atom.SetData(itag, dipole)

    for bond in mol.GetBonds():
        avg = (bond.GetBgn().GetData(itag) + bond.GetEnd().GetData(itag)) / 2.0
        bond.SetData(itag, avg)

    return True

The DepictMoleculeWithDipole function below shows how the atom values calculated by the CalculateDipoleMoment function are projected onto the property map. This gives you a sense of the direction of the dipole relative to the atoms. Before rendering the molecule the OEPrepareDepictionFrom3D function is called to generate a 2D layout of the molecule that is most similar to the 3D coordinates (line 3). The atom values are then projected onto the 2D diagram using the OE2DPropMap class (lines 10-12). See example in Figure 1.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
def DepictMoleculeWithDipole(image, mol, opts, stag):

    OEPrepareDepictionFrom3D(mol, True)

    opts.SetAtomColorStyle(OEAtomColorStyle_WhiteMonochrome)
    opts.SetScale(OEGetMoleculeSurfaceScale(mol, opts))

    disp = OE2DMolDisplay(mol, opts)

    propmap = OE2DPropMap(opts.GetBackgroundColor())
    propmap.SetLegendLocation(OELegendLocation_Left)
    propmap.Render(disp, stag)

    OERenderMolecule(image, disp)

Download code

dipole2img.py

Note

The dipole2img.py Python script requires an input molecule with 3D coordinates.

Usage:

prompt > python3 dipole2img.py molecule.mol dipole.png

Discussion

In the CalculateDipoleMoment function, after calculating the atom values, a value is calculated for each bond by averaging the values of its end atoms (lines 58-60). When the property map is rendered these bond values are projected at the middle of the bond resulting in a much smoother image. See the difference between the two images shown below. Image (A) is generated by projecting only atom values, while in image (B) both the atom and the bond values are projected onto the property map before applying a Gaussian function to blur out the 2D grid underneath the molecular graph.

Table 1.
( A ) only atom values ( B ) both atom and bond values
../_images/dipole2img-nobond.png ../_images/dipole2img.png

See also in OEChem TK manual

Theory

API

See also in OEDepict TK manual

Theory

API

See also in Grapheme TK manual

Theory

API