Embedded many-body expansions

The embedded many-body expansion (EMBE) program allows the calculation of interacting fragment energies, such as terms in a many-body expansion, in the presence of an approximate model environment created by a self-consistent charge-density-embedding scheme. The implementation allows for embedding in both finite systems (e.g. a water cluster or clathrate) and three-dimensionally periodic systems (e.g. a molecular crystal).

The environment density is constructed from Hartree-Fock (or DFT) calculations on the isolated fragment molecules; Gaussian basis functions are then fitted to the fragment densities (using standard density-fitting techniques) to create an initial guess at the environment charge density. This process can then be iterated until a converged embedding density is obtained.

A complete embedded many-body expansion calculation proceeds in the following steps:

(1) Initialisation and determination of fragments in input geometry.
(2) Self-consistent embedding environment generation, via iterated calculations on each individual fragment in the input geometry.
(3) Many-body expansion of fragment interactions (up to second-order) in the presence of the converged embedding environment.

The embedded MBE method as implemented in Molpro was originally described in
P.J. Bygrave, N.L. Allan, and F.R. Manby, J. Chem. Phys, 137, (2012)
and was generalised to involve environment basis functions of higher angular momentum in
C.R. Taylor, P.J. Bygrave, N.L. Allan, and F.R. Manby, in preparation.
Any publications resulting from the use of the approximate Gaussian embedding scheme must acknowledge the aforementioned work.

The exchange-repulsion contribution to the embedding potential is modelled using the overlap-based model of Wheatley and Price, see
R.J. Wheatley and S.L. Price, Mol. Phys. 69, 507 (1990).

An EMBE calculation consists of two principal program calls in Molpro; the embed,method=DENSITY program and the mbe program. The former is responsible for computing the embedding density, while the latter performs the many-body expansion.

The embed,method=DENSITY program additionally calls the ewald subprogram for calculating electrostatic interactions when periodic boundary conditions (PBC) are employed using the pbc program.

The basic structure of an EMBE input file is:

Basis specification, including embedding basis

Geometry

{ PBC, [lattice vectors]} if running a periodic calculation

{ embed,method=DENSITY, [options], [ewald, (ewald options)]} Ewald directive and options only if running periodic calculation

{ mbe,crystal, [options] }

MBEENG many-body fragment energy calculation procedure

DFCPROC procedure for generating density-fitted charge embedding environment

The reader is directed to the examples for full displays of the input structure.

The basis set to be used to fit the environment density must be specified in input as a basis block named DFCSET (case insensitive), for “density-fitted charges set“.

Only one fitting set may be specified, and must have this name, but is otherwise of flexible construction (e.g. different sets on different atoms or manually-specified primitives are acceptable)

The use of the Weigend-type JKFIT density-fitting sets is recommended, matched to the basis set used to compute the charge density.

The angular momentum of functions available in the basis may be limited via the standard basis input syntax, e.g.

basis={ set,dfcset default,avdz/jkfit }

produces a fitting set with the full flexibility of the aVDZ fitting basis on all atoms present, while

basis={ set,dfcset default,(p)avdz/jkfit }

produces a fitting set limited to angular momentum up to $p$ (i.e. including both $s$- and $p$-type functions).

As of June 2015, the use of higher than $d$-type functions in periodic embedding is discouraged due to a lack of Ewald screening routines for functions beyond $d$-type. Finite system embedding, however, should present no problems when higher functions are employed.

The overall system geometry is specified as usual in Molpro input using a geometry block. When dealing with periodic boundary conditions (PBC), the atomic positions can be specified in either absolute (Cartesian) coordinates or fractional coordinates as the user desires (but see discussion of space group specification below). Cubic (’CUBIC’ and parallelipiped (’PARAP’) unit cells are supported.

Atomic coordinates for the geometry block may also be stored in a separate XYZ file as per usual Molpro input.

Symmetry must be disabled (symmetry=nosym) in input for the EMBE program to run correctly.

The lattice vectors for a periodic system are specified as described in the PBC program documentation. As an example, for input geometry specified in Cartesian coordinates,

{PBC,LATT_TYPE=’PARAP’; vec, x_a , y_a, z_a vec, x_b , y_b, z_b vec, x_c , y_c, z_c }

would generate a parallelipiped (’PARAP’) unit cell with lattice vectors $\mathbf{a},\mathbf{b},\mathbf{c}$ each defined by $x$, $y$, and $z$ components $x_a, y_a, z_a$, etc (in Ångstroms).

The choice of Cartesian or fractional coordinates in the geometry affects how the PBC definitions must be given; if fractional coordinates are chosen, then the space group expansion routines must be invoked to generate the unit cell from the asymmetric unit cell provided. If the space group is unknown and/or the full primitive unit cell geometry is provided in fractional coordinates, then the space group may simply be set to the least symmetric possible (sg=1) to avoid applying any symmetry operations.

For example, if a full set of fractional atomic coordinates are specified for all atoms in the primitive unit cell, the equivalent unit cell to the above Cartesian case would be defined by:

{PBC,sg=1,LATT_TYPE=’PARAP’; vec, x_a , y_a, z_a vec, x_b , y_b, z_b vec, x_c , y_c, z_c }

i.e. the only necessary change is the addition of the sg=1 flag to invoke the routine to convert the fractional coordinates specified by the user. Expansion of an asymmetric unit cell input is not fully implemented (space groups 1-142 and 195-230 are available); however, any primitive unit cell may be specified in fractional coordinates through the use of sg=1.

The specification of the particular method to be used to calculate the embedding density and the embedded fragments’ energies is done through the standard Molpro procedure block format. One block, DFCPROC, specifies the method that will calculate the fragment densities to construct the embedding environment (the “density-fitted charges” procedure). The other, MBEENG, is used for performing the fully-embedded interaction energy calculations between fragments in the MBE.

A typical set-up for the DFCPROC would use, for example, density-fitted Hartree-Fock calculations to obtain the fragment densities that will be used to construct the environment:

PROC DFCPROC {df-hf;start,atden} ENDPROC

to be followed by an embedded MBE calculation in which fragments are treated with e.g. MP2:

PROC MBEENG {df-hf;start,atden} {df-mp2} ENDPROC

The use of F12 methods as an MBEENG energy procedure is temporarily disabled due to issues with the supporting basis set infrastructure. Restoration of this functionality is a priority in future updates of the EMBE code.

EMBED,method=DENSITY options:

ksab=[value] Defines the exchange-overlap parameter $k_{S_{AB}}$ which determines the strength of exchange repulsion as a function of density overlap.

EWALD options:

rcut=[value] Cut-off for the real-space part of the Ewald summation (in Å).

kcut=[value] Cut-off for the Fourier- (or $k$-)space part of the Ewald summation (in number of reciprocal unit cells).

Deprecated EWALD options:

These are included for reference and comparison with previous (undocumented) functionality of the EMBE code, but are no longer relevant.

gamma=[value] The value of the screening exponent for the Ewald screening density basis functions. Default value of zero, i.e. for a given atom, the widest Gaussian (smallest exponent) is chosen (see below). [The use of specified exponents in the screening density as described by Bygrave et al. is deprecated, and if a value of $\gamma$ is specified, it is simply used as an approximate target for the functions selected to form the screening density from the defined basis.]

eps=[value] Value of the dielectric constant $\varepsilon$ for the dipolar unit cell (surface term) correction to the Ewald sum. [The surface term correction is not currently implemented for the generalised embedding procedure.]

The MBE program must be invoked using {mbe,crystal,... }, even if performing a calculation in a finite embedding environment (the crystal directive initiates the embedded MBE, of which the finite form is a special case).

maxlev=[0, 1, or 2] Specifies the level (order) of many-body expansion to be performed – no MBE, first-order (monomer terms only), or second-order (monomer and dimer terms). Orders higher than 2 are not supported in EMBE. Order 0 means that no calculation will be performed after the embedding environment is generated – this is mainly for testing purposes.

cutoff=[value] The distance cut-off (in Å) for inclusion of dimer interactions. Any interactions between dimers with all atoms more than cutoff apart will be excluded.

savedb=[on or off] Specifies whether to save files (*.idb and *.edb) containing a database of all interaction energies (including monomer terms) after the job has finished. This file can then be read in using the usedb option if a calculation is restarted with a larger dimer cut-off, avoiding recalculation of the nearer interactions. Please note that the files do not store any information about the embedding environment – if a calculation is restarted with parameters that affect the embedding (such as a change of basis or change of Ewald sum cut-offs), reading this file will result in energies that are inconsistent with the new environment!

usedb=[on or off] Specifies whether to use a saved set of database files if they are present. If no files are present matching the job name with the *.idb and *.edb suffixes, this option is ignored and a new, fresh calculation is started. Note the above warning against using this option when restarting a job with different embedding parameters to the previous run!

dfciter=[value] The number of iterations to perform in the self-consistent charge embedding procedure. As a rough guide, 10 iterations typically gives convergence of the interaction of the fragments with the environment below 1 microhartree ($<$ a few J/mol).

binres=[value] At the end of an MBE calculation, a table of the total energy as a function of cut-off distance is printed; this option controls the interval (in Å) between rows of the table (e.g. binres=0.5 would print a table with increments of 0.5 Å between rows).

There are a few special options for the MBE program controlling the presence and type of counterpoise correction used in the dimer calculations:

cp=[on or off] Specifies whether to perform any counterpoise correction; cp=on requests counterpoise correction, of the type determined by voidcp.

voidcp=[on or off] If this option is off, the standard counterpoise correction is performed (i.e. dummy basis functions on the inactive fragment). Setting voidcp=on uses the counterpoise scheme of Kamiya et al. (M. Kamiya, S. Hirata, and M. Valiev, J. Chem. Phys., 128, 2008), termed the void counterpoise correction, which modifies the standard counterpoise correction to account for the influence of the embedding environment on the dummy functions on an inactive fragment. If cp=off, then voidcp is ignored.

The MBE program prints a table of energy versus radial cut-off as described in the Options section, followed by a single value corresponding to the total energy at the maximum cut-off specified in input (in the form MBE/[basis] energy $=$ [value]). These tables are printed for each of the reference (HF or KS-DFT) energy ENERGR, the correlation energy (from MP2, CCSD(T), etc.) EMP2,ECCSD etc., and the embedding energy EMBEDENERGY, broken down into its one- and two-body contributions.

Note that when performing an MP2 calculation, the related but distinct SCS-MP2 correlation energy is printed automatically also. This is to match standard Molpro functionality, as when an MP2 calculation is done on any syste, the SCS-MP2 energy is printed as well. These correlation energies are mutually-exclusive, i.e. they should not be summed together when calculating the total energy at a given level of theory.

The final total energy of an e.g. MP2 calculation is the sum of each of the tabulated terms (the reference and correlation energy at both the 1-body and 2-body levels). However, note that the embedding energy as printed in the table double-counts the interaction of the orbitals with the embedding environment. Thus, the correct energy expression for the result of an EMBE calculation is $$E_{\text{total}} = E_{\text{HF}}+E_{\text{MP2}} - \frac{1}{2}E_{\text{embed}}.$$ Note that the “final” single energy value printed at the very end of the calculation already has this correction for double-counting applied (i.e. the above sum only needs to be performed when calculating the energy at a particular cut-off distance using the values from the summary tables).

In this example, a finite system embedded two-body expansion is performed on the water hexamer in the prism structure, using DF-HF for both the embedding environment generation and the fragment interaction energies. The embedding environment consists of $s$- and $p$-functions and the calculations are not counterpoise-corrected.

examples/embe_water_hexamer.inp

skipped ! bug5219
gdirect
BASIS=avdz
basis={
set,dfcset
default,(p)avdz/jkfit
}

symmetry,nosym
noorient
geomtyp=xyz
! E(RHF/6-31G(d,p))=-456.2206685  # water-6 prism from Cambridge Cluster Database
geometry={
O                   1.526924   -0.407825    1.467087
H                   0.615293   -0.694811    1.569578
H                   2.001971   -0.712098    2.222674
O                   0.831994    2.016713   -0.073306
H                   1.169971    1.534473    0.670375
H                   1.181010    1.553886   -0.823001
O                   1.709464   -0.492590   -1.370853
H                   2.421392   -0.900609   -1.835542
H                   1.839304   -0.657788   -0.441680
O                  -1.891336    1.277844   -0.075746
H                  -1.013035    1.654737   -0.078267
H                  -2.500751    1.992871   -0.154652
O                  -1.172193   -1.306705   -1.416807
H                  -0.263214   -1.117130   -1.610277
H                  -1.591494   -0.462719   -1.314422
O                  -1.135583   -1.129560    1.451997
H                  -1.200738   -1.536124    0.593247
H                  -1.613872   -0.317697    1.342986
}

{embed,ksab=7.5d0}

{mbe,crystal,cp=off,cutoff=20.0,maxlev=2,usedb=off,binres=0.5,dfciter=10}

PROC MBEENG
 {df-hf;start,atden}
ENDPROC

PROC DFCPROC
 {df-hf;start,atden}
ENDPROC

The anisotropy of the embedding environment’s representation of individual atoms can be controlled by altering the maximum angular momentum available to the embedding basis DFCSET. In this example, a periodic Hartree-Fock MBE calculation is performed on crystalline hydrogen fluoride, using an embedding density which is constructed only using $s$-functions (i.e. spherical atomic densities only). The dimer calculations are not counterpoise-corrected. This input also demonstrates the use of the space-group expansion routine to provide an asymmetric unit cell structure in fractional coordinates.

examples/embe_hfcrystal.inp

memory,200,m
gdirect
BASIS=avdz
basis={
set,dfcset
default,(s)avdz/jkfit
}

symmetry,nosym
noorient
geomtyp=xyz

! Hydrogen fluoride crystal geometry from Atoji and Lipscomb, Acta. Cryst., 7, 1954.
a=4.394930691
b=3.411518525
c=5.268740272
z1=0.635051640/c
x2=1.904707822/a
z2=1.124859274/c

geometry={
 f  0.25  0.25  z1
 h  x2    0.25  z2
}

{PBC,sg=62
VEC,$a, 0 ,0;
VEC, 0,$b, 0;
VEC, 0, 0,$c}

{embed,method=DENSITY,ksab=4.0,
 ewald,rcut=7.5,kcut=8,gamma=0.00,eps=-1}

{mbe,crystal,cp=off,cutoff=5,maxlev=2,usedb=off,uselog=on,savedb=off,binres=0.5,dfciter=8}

PROC MBEENG
 {df-hf;start,atden}
ENDPROC

PROC DFCPROC
 {df-hf;start,atden}
ENDPROC

Increasing the angular momentum available to the fitting set gives, in principle, a more physically-realistic description of the surrounding environment. This example demonstrates an MP2 MBE calculation on the cubic phase of ice XI, using an embedding environment constructed from $s$-,$p$-, and $d$-functions. In this calculation, the void counterpoise-correction scheme is used.

examples/embe_ice.inp

memory,200,m
gdirect
BASIS=avdz
basis={
set,dfcset
default,(d)avdz/jkfit
}

symmetry,nosym
noorient
geometry={
O     0.000000000000000     0.000000000000000     0.000000000000000
H     0.000000000000000     0.780449156624030     0.590499143827740
H     0.000000000000000    -0.780449156624030     0.590499143827740
O    -2.240586682500000    -2.240586682500000    -3.242747568500000
H    -2.240586682500000    -1.460137525875970    -2.652248424672259
H    -2.240586682500000    -3.021035839124030    -2.652248424672259
O     0.000000000000000    -2.240586682500000     1.621373784250000
H    -0.780449156624030    -2.240586682500000     2.211872928077740
H     0.780449156624030    -2.240586682500000     2.211872928077740
O    -2.240586682500000     0.000000000000000    -1.621373784250000
H    -3.021035839124030     0.000000000000000    -1.030874640422259
H    -1.460137525875970     0.000000000000000    -1.030874640422259
}

{PBC,LATT_TYPE='PARAP';
vec,    4.481173365000000,     0.000000000000000,     0.000000000000000
vec,    0.000000000000000,     4.481173365000000,     0.000000000000000
vec,    0.000000000000000,     0.000000000000000,     6.485495137000000
}

! Change infile for each job (current file name)
{embed,method=DENSITY,dfcmethod=restart,ksab=7.7d0,
 ewald,rcut=8,kcut=8,gamma=0.00,eps=-1}

{mbe,crystal,cp=on,cutoff=10.0,maxlev=2,usedb=off,uselog=on,binres=0.5,dfciter=10,voidcp=on}

PROC MBEENG
 {df-hf;start,atden}
 {df-mp2}
ENDPROC

PROC DFCPROC
 {df-hf;start,atden}
ENDPROC