Non-polarizable embedding#

To illustrate this concept, here we consider the QM/MM version using the Hartree-Fock method. The total energy of the molecular system in the QM/MM method is partitioned into three contributions:

\[ E = E_\mathrm{MM} + E_\mathrm{MM/QM} + E_\mathrm{QM} , \]

where the first term is the energy of the MM region, the second term is the interaction energy of the QM and MM regions, and the last term is the QM region energy. In this section we use a non-polarizable force field to describe the MM region, and the second term can be expanded as:

\[ E_\mathrm{MM/QM} = \sum_i q_i^\mathrm{perm} (\phi_{i}^\mathrm{ele}+ \phi_{i}^{nuc} ), \]

where \(\{ q_i^\mathrm{perm}\}\) is the set of permanent charges in MM region, and the \(\phi_{i}^\mathrm{ele}\) and \(\phi_{i}^\mathrm{nuc}\) are potential components are generated by electrons and nuclei in MM region.

This is now to be illustrated by considering the energy contributions and dipole moment of one QM water molecule, as considered individually, when coordinated with four MM water molecules, or coordinated with two donors or two acceptors:

import numpy as np
import scipy
import veloxchem as vlx
import py3Dmol
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(viewergrid=(2, 2), width=600, height=500, linked=False)
viewer.setViewStyle({"style": "outline", "color": "black", "width": 0.1})
viewer.addModel(h2o_xyz, viewer=(0, 0))
viewer.addModel(h2o_xyz, viewer=(0, 1))
viewer.addModel(h2o_xyz, viewer=(1, 0))
viewer.addModel(h2o_xyz, viewer=(1, 1))
viewer.addModel(solvent_xyz, viewer=(0, 1))
viewer.addModel(donors_xyz, viewer=(1, 0))
viewer.addModel(acceptors_xyz, viewer=(1, 1))
viewer.setStyle({"stick": {}})
# 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

Dipole moments#

The molecular dipole moment is the sum of its electronic and nuclear parts. 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 . \]

QM#

Setting up a calculation for a single molecule:

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

norb = vlx.MolecularBasis.get_dimensions_of_basis(basis, molecule)
nocc = molecule.number_of_alpha_electrons()
V_nuc = molecule.nuclear_repulsion_energy()

print("Number of contracted basis functions:", norb)
print("Number of doubly occupied molecular orbitals:", nocc)
print(f"Nuclear repulsion energy (in a.u.): {V_nuc : 14.12f}")
* Info * Reading basis set from file: /home/thomas/Notebook/anaconda/envs/echem/lib/python3.9/site-packages/veloxchem/basis/CC-PVDZ
                                                                                                                          
                                              Molecular Basis (Atomic Basis)                                              
                                             ================================                                             
                                                                                                                          
                                  Basis: CC-PVDZ                                                                          
                                                                                                                          
                                  Atom Contracted GTOs          Primitive GTOs                                            
                                                                                                                          
                                   O   (3S,2P,1D)               (17S,4P,1D)                                               
                                   H   (2S,1P)                  (4S,1P)                                                   
                                                                                                                          
                                  Contracted Basis Functions : 24                                                         
                                  Primitive Basis Functions  : 48                                                         
                                                                                                                          
Number of contracted basis functions: 24
Number of doubly occupied molecular orbitals: 5
Nuclear repulsion energy (in a.u.):  7.964303024713

SCF calculation#

In preparation for the embedded systems, we here set up and use a tailored SCF solver:

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 = np.einsum("ik,jk->ij", C[:, :nocc], C[:, :nocc])

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

        E = np.einsum("ij,ij->", h + F, D) + 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

Calculating the overlap matrix and one- and two-electron Hamiltonians:

# overlap matrix
overlap_drv = vlx.OverlapIntegralsDriver()
S = overlap_drv.compute(molecule, basis).to_numpy()

# one-electron Hamiltonian
kinetic_drv = vlx.KineticEnergyIntegralsDriver()
T = kinetic_drv.compute(molecule, basis).to_numpy()

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

h = T + V

# two-electron Hamiltonian
eri_drv = vlx.ElectronRepulsionIntegralsDriver()
g = eri_drv.compute_in_memory(molecule, basis)

The SCF is calculated with:

epsilon, C = scipy.linalg.eigh(h, S)  # initial guess
E_HF, C_HF = scf_solver(h, V_nuc, C)
iter      SCF energy    Error norm
 0      -68.97278671    1.97e+00
 1      -69.54016846    1.83e+00
 2      -72.78571482    1.77e+00
 3      -72.80872293    1.48e+00
 4      -74.05300746    1.49e+00
 5      -74.77710497    1.14e+00
 6      -75.34750222    9.29e-01
 7      -75.66180142    6.38e-01
 8      -75.83364767    4.64e-01
 9      -75.91583567    3.07e-01
10      -75.95476888    2.12e-01
11      -75.97240552    1.41e-01
12      -75.98042144    9.52e-02
13      -75.98400915    6.33e-02
14      -75.98562135    4.25e-02
15      -75.98634173    2.84e-02
16      -75.98666440    1.90e-02
17      -75.98680860    1.27e-02
18      -75.98687312    8.50e-03
19      -75.98690196    5.68e-03
20      -75.98691487    3.80e-03
21      -75.98692063    2.54e-03
22      -75.98692321    1.70e-03
23      -75.98692437    1.14e-03
24      -75.98692488    7.60e-04
25      -75.98692511    5.08e-04
26      -75.98692522    3.40e-04
27      -75.98692526    2.27e-04
28      -75.98692528    1.52e-04
29      -75.98692529    1.02e-04
30      -75.98692530    6.79e-05
SCF iterations converged!

Dipole moment#

With this, the dipole moment can be calculated as:

def dipmom(D):

    mu_e = np.zeros(3)
    mu_n = np.zeros(3)

    # electronic part
    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()

    mu_e[0] = 2 * np.einsum("ab, ab", D, mu_x)
    mu_e[1] = 2 * np.einsum("ab, ab", D, mu_y)
    mu_e[2] = 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

    mu = mu_e + mu_n

    return mu

Using the element charges, we obtain

element_charges = {"H": 1.0, "O": 8.0}

D = 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

QM/MM#

Embedding potentials#

The hybrid QM/MM approach for multiscale modeling of complex systems was recognized by the Nobel Prize in 2013. In this model, the Hamiltonian takes the form

\[ \hat{H} = \hat{H}^\mathrm{QM} + \hat{H}^\mathrm{MM} + \hat{H}^\mathrm{QM/MM} \]

The electrostatic interaction operator between the QM and MM regions takes the form

\[ \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|} \]

We will use this model to study water in the condensed phase with hydrogen bonds. The TIP3P model provides a classical description of liquid water in terms of three-site model with partial charges.

Element

Charge

O

-0.834

H

0.417

Setting up these charges and molecule objects:

tip3p_charges = {"H": 0.417, "O": -0.834}

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

The embedding operator is constructed as:

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

SCF optimizations#

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)
iter      SCF energy    Error norm
 0      -76.04721904    6.37e-02
 1      -76.05404597    3.00e-02
 2      -76.05480572    1.40e-02
 3      -76.05496833    9.79e-03
 4      -76.05501549    5.76e-03
 5      -76.05503195    3.89e-03
 6      -76.05503832    2.44e-03
 7      -76.05504091    1.61e-03
 8      -76.05504198    1.04e-03
 9      -76.05504243    6.77e-04
10      -76.05504262    4.37e-04
11      -76.05504270    2.85e-04
12      -76.05504273    1.85e-04
13      -76.05504274    1.20e-04
14      -76.05504275    7.79e-05
SCF iterations converged!
iter      SCF energy    Error norm
 0      -76.01164058    3.40e-02
 1      -76.01344618    1.41e-02
 2      -76.01364388    6.42e-03
 3      -76.01368402    4.49e-03
 4      -76.01369499    2.65e-03
 5      -76.01369869    1.83e-03
 6      -76.01370012    1.16e-03
 7      -76.01370071    7.85e-04
 8      -76.01370096    5.11e-04
 9      -76.01370106    3.42e-04
10      -76.01370111    2.25e-04
11      -76.01370113    1.49e-04
12      -76.01370114    9.87e-05
SCF iterations converged!
iter      SCF energy    Error norm
 0      -76.02250376    3.98e-02
 1      -76.02508400    1.79e-02
 2      -76.02534197    8.27e-03
 3      -76.02539710    5.90e-03
 4      -76.02541372    3.52e-03
 5      -76.02541974    2.40e-03
 6      -76.02542215    1.53e-03
 7      -76.02542314    1.02e-03
 8      -76.02542357    6.58e-04
 9      -76.02542375    4.34e-04
10      -76.02542382    2.84e-04
11      -76.02542386    1.86e-04
12      -76.02542387    1.22e-04
13      -76.02542388    8.00e-05
SCF iterations converged!

Solvation energies#

D = np.einsum("ik,jk->ij", C_s[:, :nocc], C_s[:, :nocc])

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

print(f"Nuclear QM/MM energy   : {V_n_s : 16.12f}")
print(f"Electronic QM/MM energy: {v_s_aver : 16.12f}")
print(41 * "-")
print(f"Total QM/MM energy     : {V_n_s + v_s_aver : 16.12f}")
Nuclear QM/MM energy   :  -0.085459561697
Electronic QM/MM energy:   0.009877803071
-----------------------------------------
Total QM/MM energy     :  -0.075581758627
print(f"   vacuum: {E_HF : 12.8f}")
print(f" solution: {E_s : 12.8f} {E_s - V_n_s: 14.8f} {E_s - V_n_s - E_HF: 13.8f}")
print(f"   donors: {E_d : 12.8f} {E_d - V_n_d: 14.8f} {E_d - V_n_d - E_HF: 13.8f}")
print(f"acceptors: {E_a : 12.8f} {E_a - V_n_a: 14.8f} {E_a - V_n_a - E_HF: 13.8f}")
   vacuum: -75.98692530
 solution: -76.05504275   -75.96958319    0.01734211
   donors: -76.01370114   -76.66472104   -0.67779574
acceptors: -76.02542388   -75.28894441    0.69798088

Dipole moments#

D = np.einsum("ik,jk->ij", C_s[:, :nocc], C_s[:, :nocc])
mu_s = dipmom(D)

D = np.einsum("ik,jk->ij", C_d[:, :nocc], C_d[:, :nocc])
mu_d = dipmom(D)

D = np.einsum("ik,jk->ij", 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

Compare with the TIP3P dipole moment.

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 momemnt: {np.linalg.norm(mu_tip3p) : 8.6f}")
Dipole momemnt:  1.112794