Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Non-polarizable embedding

Let us consider a system with NN electrons and MM nuclei in the quantum mechanical (QM) region embedded by a molecular mechanical (MM) region that contains point charges, qsq_s, located at site positions Rs\mathbf{R}_s. In the non-polarizable embedding model, the electrostatic interaction operator between the QM and MM regions is given as

H^QM/MM=Vn,es+∑i=1Nv^es(i)\hat{H}^\mathrm{QM/MM} = V^\mathrm{n,es} + \sum_{i=1}^N \hat{v}^\mathrm{es}(i)
(1)

where

Vn,es=∑A=1M∑sZAe qs4πε0∣RA−Rs∣;v^es=∑s−eqs4πε0∣r−Rs∣V^\mathrm{n,es} = \sum_{A=1}^M \sum_s \frac{Z_A e \, q_s}{4\pi\varepsilon_0 |\mathbf{R}_A - \mathbf{R}_s|} ; \qquad % \hat{v}^\mathrm{es} = \sum_s \frac{-e q_s}{4\pi\varepsilon_0 |\mathbf{r} - \mathbf{R}_s|}
(2)

To illustrate the concept, we consider the QM/MM embedding of water in water using the Hartree–Fock method. We determine the energy contributions and dipole moment of a QM water molecule embedded by four MM water molecules in a tetrahedral configuration. The oxygen in the QM region acts as donor and acceptor of two hydrogen bonds each.

Let us first define the system coordinates and view the system.

Source
import numpy as np
import py3Dmol
import scipy
import veloxchem as vlx

h2o_xyz = """3
water
O        0.0000000000      0.0000000000      0.0000000000                 
H        0.6891400000      0.8324710000      0.0000000000                 
H        0.7224340000     -0.8726890000      0.0000000000
"""

solvent_xyz = """12
water solvent
O        1.0991810000      2.2167050000      0.2791250000                 
H        2.2673500000      2.2529740000      0.6013580000                 
H        0.9984050000      2.8080510000     -0.5840690000                 
O       -0.8223690000      0.2141140000     -2.3767710000                 
H       -1.5990410000      0.8408210000     -2.6698860000                 
H       -0.4683410000      0.4151750000     -1.5365150000                 
O        1.7377440000     -2.0168380000     -0.4861090000                 
H        2.0753880000     -2.2446490000     -1.3654140000                 
H        1.6636970000     -2.7516570000      0.0543720000                 
O       -1.0758700000     -0.2692330000      2.5139440000                 
H       -0.5822670000     -0.1887880000      1.5058600000                 
H       -1.1632420000      0.7082890000      2.8568430000 
"""

donors_xyz = """6
water solvent                                  
O       -0.8223690000      0.2141140000     -2.3767710000                 
H       -1.5990410000      0.8408210000     -2.6698860000                 
H       -0.4683410000      0.4151750000     -1.5365150000
O       -1.0758700000     -0.2692330000      2.5139440000                 
H       -0.5822670000     -0.1887880000      1.5058600000                 
H       -1.1632420000      0.7082890000      2.8568430000
"""

acceptors_xyz = """6
water solvent
O        1.0991810000      2.2167050000      0.2791250000                 
H        2.2673500000      2.2529740000      0.6013580000                 
H        0.9984050000      2.8080510000     -0.5840690000
O        1.7377440000     -2.0168380000     -0.4861090000                 
H        2.0753880000     -2.2446490000     -1.3654140000                 
H        1.6636970000     -2.7516570000      0.0543720000              
"""
viewer = py3Dmol.view(width=300, height=300, linked=False)

viewer.addModel(h2o_xyz)
viewer.addModel(solvent_xyz)

viewer.setStyle({"stick": {}})
viewer.setViewStyle({"style": "outline", "color": "black", "width": 0.1})

# rotate for a better initial view
viewer.rotate(45, "x")

viewer.show()
Loading...

As a probe of the changes in the electron density of the solute water in the QM region as due to the embedding, we will study the molecular dipole moment. The dipole moment is the sum of its electronic and nuclear parts, where the electronic part is the expectation value of a one-electron operator and the nuclear part is obtained as a sum over point charges:

μ=⟨μ^⟩+μn\boldsymbol{\mu} = \langle \hat{\boldsymbol{\mu}} \rangle + \boldsymbol{\mu}^\mathrm{n}
(3)

where

μn=∑A=1MZARA.\boldsymbol{\mu}^\mathrm{n} = \sum_{A = 1}^M Z_A \mathbf{R}_A .
(4)
element_charges = {"H": 1.0, "O": 8.0}

We implement a function to determine the dipole moment given the one-particle density matrix in the atomic orbital (AO) basis.

def dipmom(D):
    mu_e = np.zeros(3)
    mu_n = np.zeros(3)

    # electronic part
    mu_e[0] = np.einsum("ab, ab", D, mu_x)
    mu_e[1] = np.einsum("ab, ab", D, mu_y)
    mu_e[2] = np.einsum("ab, ab", D, mu_z)

    # nuclear part
    for A, molecule_atom_label in enumerate(molecule.get_labels()):
        R_A = np.array(molecule.get_atom_coordinates(A))
        Z_A = element_charges[molecule_atom_label]

        mu_n += Z_A * R_A

    # total dipole moment
    mu = mu_e + mu_n

    return mu

Gas phase water

For reference, we determine the properties of an isolated water molecule.

molecule = vlx.Molecule.read_xyz_string(h2o_xyz)
basis = vlx.MolecularBasis.read(molecule, "cc-pvdz")

nocc = molecule.number_of_alpha_electrons()
V_nuc = molecule.nuclear_repulsion_energy()

SCF calculation

We first retrieve needed one- and two-electron integrals in the AO basis.

Source
# overlap
S = vlx.compute_overlap_integrals(molecule, basis)

# electric-dipole
dipole_mats = vlx.compute_electric_dipole_integrals(molecule, basis)

mu_x = dipole_mats[0]
mu_y = dipole_mats[1]
mu_z = dipole_mats[2]

# kinetic energy
T = vlx.compute_kinetic_energy_integrals(molecule, basis)

# nuclear attraction
V = vlx.compute_nuclear_potential_integrals(molecule, basis)

# one-electron Hamiltonian
h = T + V

# two-electron Hamiltonian
fock_drv = vlx.FockDriver()
g = fock_drv.compute_eri(molecule, basis)

The SCF optimization is performed with the following solver, see the Hartree-Fock section for details.

def scf_solver(h, V_nuc, C):
    max_iter = 50
    conv_thresh = 1e-4

    print("iter      SCF energy    Error norm")

    for iter in range(max_iter):
        D_alpha = np.einsum("ik,jk->ij", C[:, :nocc], C[:, :nocc])

        J = np.einsum("ijkl,kl->ij", g, D_alpha)
        K = np.einsum("ilkj,kl->ij", g, D_alpha)
        F = h + 2 * J - K

        E = np.einsum("ij,ij->", h + F, D_alpha) + V_nuc

        # compute convergence metric
        F_MO = np.einsum("ki,kl,lj->ij", C, F, C)
        e_vec = np.reshape(F_MO[:nocc, nocc:], -1)
        error = np.linalg.norm(e_vec)
        print(f"{iter:>2d}  {E:16.8f}  {error:10.2e}")

        if error < conv_thresh:
            print("SCF iterations converged!")
            break

        epsilon, C = scipy.linalg.eigh(F, S)

    return E, C

We perform an SCF optimization to get the Hartree–Fock wave function of an isolated water molecule.

epsilon, C = scipy.linalg.eigh(h, S)  # initial guess
E_HF, C_HF = scf_solver(h, V_nuc, C)

Dipole moment

From the molecular orbital (MO) coefficients, we determine the density matrix and compute the dipole moment.

D = 2 * np.einsum("ik,jk->ij", C_HF[:, :nocc], C_HF[:, :nocc])
mu_HF = dipmom(D)

print(f"Vacuum dipole moment: {np.linalg.norm(mu_HF) : 8.6f}")
Vacuum dipole moment:  0.878333

Liquid water

The TIP3P model provides a classical description of liquid water in terms of a three-site model with the following partial charges:

ElementCharge
O-0.834
H0.417
tip3p_charges = {"H": 0.417, "O": -0.834}

Let us calculate the dipole moment for a water molecule in the TIP3P model.

mu_tip3p = np.zeros(3)

for A, molecule_atom_label in enumerate(molecule.get_labels()):
    R_A = np.array(molecule.get_atom_coordinates(A))
    q_A = tip3p_charges[molecule_atom_label]

    mu_tip3p += q_A * R_A

print(f"Dipole moment: {np.linalg.norm(mu_tip3p) : 8.6f}")
Dipole moment:  1.112794

This dipole moment is substantially larger than the gas phase result.

In order to determine a QM/MM value for the dipole moment of water, we adopt the TIP3P parameters for the description of the four water molecules in the environment (or solvent) that consists of two donors and two acceptors.

solvent = vlx.Molecule.read_xyz_string(solvent_xyz)
donors = vlx.Molecule.read_xyz_string(donors_xyz)
acceptors = vlx.Molecule.read_xyz_string(acceptors_xyz)

Embedding operator

We implement a function to construct the embedding operators given a Molecule object that contains the embedding atoms.

def get_embedding_operator(embedding):
    # electronic part
    mm_sites = []
    mm_charges = []
    for atom_idx, atom_label in enumerate(embedding.get_labels()):
        mm_sites.append(embedding.get_atom_coordinates(atom_idx))
        mm_charges.append(tip3p_charges[atom_label])

    v_es = vlx.compute_nuclear_potential_integrals(molecule, basis, mm_charges, mm_sites)

    # nuclear part
    V_n_es = 0.0
    for A, molecule_atom_label in enumerate(molecule.get_labels()):
        R_A = np.array(molecule.get_atom_coordinates(A))
        Z_A = element_charges[molecule_atom_label]

        for s, solvent_atom_label in enumerate(embedding.get_labels()):
            R_s = np.array(embedding.get_atom_coordinates(s))
            q_s = mm_charges[s]

            V_n_es += Z_A * q_s / np.linalg.norm(R_A - R_s)

    return V_n_es, v_es

Solvation energies

Using the initial guess for gas phase water, the SCF optimizations for three embedded systems are performed corresponding to (i) all four water molecules in the solvent, (ii) the two donors, and (iii) the two acceptors.

V_n_s, v_s = get_embedding_operator(solvent)
E_s, C_s = scf_solver(h + v_s, V_nuc + V_n_s, C_HF)

V_n_d, v_d = get_embedding_operator(donors)
E_d, C_d = scf_solver(h + v_d, V_nuc + V_n_d, C_HF)

V_n_a, v_a = get_embedding_operator(acceptors)
E_a, C_a = scf_solver(h + v_a, V_nuc + V_n_a, C_HF)

The solvation energy of the fully coordinated pentamer is calculated as

D = 2 * np.einsum("ai,bi->ab", C_s[:, :nocc], C_s[:, :nocc])

v_s_aver = np.einsum("ab, ab", D, v_s)

print(f"Nuclear QM/MM energy   : {V_n_s : 16.8f}")
print(f"Electronic QM/MM energy: {v_s_aver : 16.8f}")
print(41 * "-")
print(f"Total QM/MM energy     : {V_n_s + v_s_aver : 16.8f}")
Nuclear QM/MM energy   :      -0.08545956
Electronic QM/MM energy:       0.00987780
-----------------------------------------
Total QM/MM energy     :      -0.07558176

As can be seen, the pure electronic contributions destabilize the system, but the nuclear contributions stabilize it even more. As such, solvation is seen to be beneficial, as a result of the nuclear contributions.

Let us next determine these energies for all three coordination complexes.

print("Structure   Total energy      Nuclear        Electronic   Diff. to vacuum  ")
print(73 * "-")
print(f"   vacuum:  {E_HF : 12.8f}")
print(
    f" solution:  {E_s : 12.8f}  {V_n_s: 14.8f}  {E_s - V_n_s: 14.8f}  {E_s - V_n_s - E_HF: 13.8f}"
)
print(
    f"   donors:  {E_d : 12.8f}  {V_n_d: 14.8f}  {E_d - V_n_d: 14.8f}  {E_d - V_n_d - E_HF: 13.8f}"
)
print(
    f"acceptors:  {E_a : 12.8f}  {V_n_a: 14.8f}  {E_a - V_n_a: 14.8f}  {E_a - V_n_a - E_HF: 13.8f}"
)
Structure   Total energy      Nuclear        Electronic   Diff. to vacuum  
-------------------------------------------------------------------------
   vacuum:  -75.98692530
 solution:  -76.05504275     -0.08545956    -75.96958319     0.01734211
   donors:  -76.01370114      0.65101990    -76.66472104    -0.67779574
acceptors:  -76.02542388     -0.73647946    -75.28894441     0.69798088

As can be seen, the purely electronic contribution stabilizes the system for the donors, but not for the acceptors and the fully coordinated system. This can be understood by noting that the donors add positive point charges (hydrogens) close to the QM oxygen atom, which provides a beneficial potential for the electronic structure. Adding two acceptors instead puts two negative (oxygen) charges close to the electron density, and thus destabilizes the system.

Dipole moments

Next, we determine the dipole moments of the solute in the different configurations.

D = 2 * np.einsum("ai,bi->ab", C_s[:, :nocc], C_s[:, :nocc])
mu_s = dipmom(D)

D = 2 * np.einsum("ai,bi->ab", C_d[:, :nocc], C_d[:, :nocc])
mu_d = dipmom(D)

D = 2 * np.einsum("ai,bi->ab", C_a[:, :nocc], C_a[:, :nocc])
mu_a = dipmom(D)

print(f"   vacuum: {np.linalg.norm(mu_HF) : 8.6f}")
print(f" solution: {np.linalg.norm(mu_s) : 8.6f}")
print(f"   donors: {np.linalg.norm(mu_d) : 8.6f}")
print(f"acceptors: {np.linalg.norm(mu_a) : 8.6f}")
   vacuum:  0.878333
 solution:  1.147602
   donors:  1.006242
acceptors:  1.031627

The dipole moment is seen to substantially increase due to solvation. The donors pull and acceptors push the electron density of the solute in the direction from the hydrogens toward the oxygen and thus increase the polarization of the QM system.

We note that the QM/MM dipole moment of the fully solvated water molecule is in good agreement with that of the TIP3P model.

Reference calculation

Let us verify our results with a reference calculation using VeloxChem.

# create a veloxchem *.pot file
with open("embedding.pot", "w") as pot_file:
    pot_file.write("""12
xyz
O        1.0991810000      2.2167050000      0.2791250000    -0.834
H        2.2673500000      2.2529740000      0.6013580000     0.417
H        0.9984050000      2.8080510000     -0.5840690000     0.417
O       -0.8223690000      0.2141140000     -2.3767710000    -0.834
H       -1.5990410000      0.8408210000     -2.6698860000     0.417
H       -0.4683410000      0.4151750000     -1.5365150000     0.417
O        1.7377440000     -2.0168380000     -0.4861090000    -0.834
H        2.0753880000     -2.2446490000     -1.3654140000     0.417
H        1.6636970000     -2.7516570000      0.0543720000     0.417
O       -1.0758700000     -0.2692330000      2.5139440000    -0.834
H       -0.5822670000     -0.1887880000      1.5058600000     0.417
H       -1.1632420000      0.7082890000      2.8568430000     0.417
""")
scf_drv = vlx.ScfRestrictedDriver()
scf_drv.point_charges = "embedding.pot"
scf_results = scf_drv.compute(molecule, basis)
print(f"Total energy: {scf_results['scf_energy']:12.8f}")
Total energy: -76.05504276

We note that this energy is in perfect agreement with our previous result.

Environment
Embedding potentials
Environment
Polarizable embedding