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
where
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.
Show 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.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()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
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:
where
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.
Show 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)
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])
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.