Non-polarizable embedding#

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

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

where

$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|}$

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.

Hide code cell 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.setStyle({"stick": {}})
viewer.setViewStyle({"style": "outline", "color": "black", "width": 0.1})

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

viewer.show()


You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

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:

$\boldsymbol{\mu} = \langle \hat{\boldsymbol{\mu}} \rangle + \boldsymbol{\mu}^\mathrm{n}$

where

$\boldsymbol{\mu}^\mathrm{n} = \sum_{A = 1}^M Z_A \mathbf{R}_A .$
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)

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.

Hide code cell source
# overlap
overlap_drv = vlx.OverlapIntegralsDriver()
S = overlap_drv.compute(molecule, basis).to_numpy()

# electric-dipole
dipole_drv = vlx.ElectricDipoleIntegralsDriver()

dipole_mats = dipole_drv.compute(molecule, basis)

mu_x = -1.0 * dipole_mats.x_to_numpy()
mu_y = -1.0 * dipole_mats.y_to_numpy()
mu_z = -1.0 * dipole_mats.z_to_numpy()

# kinetic energy
kinetic_drv = vlx.KineticEnergyIntegralsDriver()
T = kinetic_drv.compute(molecule, basis).to_numpy()

# nuclear attraction
nucpot_drv = vlx.NuclearPotentialIntegralsDriver()
V = -1.0 * nucpot_drv.compute(molecule, basis).to_numpy()

# one-electron Hamiltonian
h = T + V

# two-electron Hamiltonian
eri_drv = vlx.ElectronRepulsionIntegralsDriver()
g = eri_drv.compute_in_memory(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:

Element

Charge

O

-0.834

H

0.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)


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])

pot_drv = vlx.NuclearPotentialIntegralsDriver()
v_es = -1.0 * pot_drv.compute(molecule, basis, mm_charges, mm_sites).to_numpy()

# 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(
"""@environment
units: angstrom
xyz:
O        1.0991810000      2.2167050000      0.2791250000   water 1
H        2.2673500000      2.2529740000      0.6013580000   water 1
H        0.9984050000      2.8080510000     -0.5840690000   water 1
O       -0.8223690000      0.2141140000     -2.3767710000   water 2
H       -1.5990410000      0.8408210000     -2.6698860000   water 2
H       -0.4683410000      0.4151750000     -1.5365150000   water 2
O        1.7377440000     -2.0168380000     -0.4861090000   water 3
H        2.0753880000     -2.2446490000     -1.3654140000   water 3
H        1.6636970000     -2.7516570000      0.0543720000   water 3
O       -1.0758700000     -0.2692330000      2.5139440000   water 4
H       -0.5822670000     -0.1887880000      1.5058600000   water 4
H       -1.1632420000      0.7082890000      2.8568430000   water 4
@end

@charges
O  -0.834  water
H   0.417  water
H   0.417  water
@end
"""
)

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.potfile = "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.