Table of Contents

Projection-based WF-in-DFT embedding

Embedding methods allow a system to be divided into two smaller subsystems, each of which can be treated using a different level of theory. For example, in wavefunction-in-DFT (WF-in-DFT) embedding, a WF-level (e.g. HF, MP2, CCSD(T), CASSCF, etc.) calculation is performed on one subsystem, while a DFT-level calculation is performed on the other subsystem. The interactions between the two subsystems are calculated at the DFT level. The primary advantage of WF-in-DFT embedding is that it facilitates the application of an accurate, systematically improvable WF method to regions where such accuracy is required, while a more efficient DFT method is applied to the remainder of the system. The overall strategy for projection-based WF-in-DFT embedding is described as follows:

(1) A DFT calculation is performed on the full system to obtain a set of canonical molecular orbitals (MOs).

(2) The canonical orbitals are localized using a localization method such as Pipek-Mezey or intrinsic bond orbitals (IBOs).

(3) The occupied localized MOs (LMOs) are partitioned into two subsystems, labeled the active subsystem (A) and the frozen subsystem (B).

(4) The interaction potential between the two subsystems is calculated at the DFT level.

(5) A WF-level calculation is performed on the LMOs in the active subsystem, embedded in the DFT-level interaction potential produced by the electrons in the frozen subsystem (B).

Note: WF-in-DFT embedding is not implemented with symmetry.

It is possible to replace the WF-level calculation on subsystem A with a DFT-level calculation, which corresponds to DFT-in-DFT embedding. Since the interaction potential between the subsystems is calculated at the DFT level, the result of a DFT-in-DFT embedding calculation is numerically identical to the results of a DFT calculation on the full system. Similarly, it is possible to replace the DFT-level calculation with a HF-level calculation, which corresponds to WF-in-HF embedding.

Molpro implements the numerically exact projection-based WF-in-DFT embedding method for open and closed shell systems developed in the following papers:

F. R. Manby, M. Stella, J. D. Goodpaster, and T. F. Miller III, J. Chem. Theory Comput. 8, 2564 (2012)

T. A. Barnes, J. D. Goodpaster, F. R. Manby, and T. F. Miller III, J. Chem. Phys. 139, 024103 (2013)

J. D. Goodpaster, T. A. Barnes, F. R. Manby, and T. F. Miller III, J. Chem. Phys., 140, 18A507 (2014)

S. J. Bennie, M. Stella, T. F. Miller III and F. R. Manby, J. Chem. Phys. 143, 024105 (2015)

S. J. R. Lee, F. Ding, F. R. Manby, T. F. Miller III, arXiv e-prints, (2019), arXiv:1903.05830

S. J. R. Lee, M. Welborn, F. R. Manby, and T. F. Miller III, Acc. Chem. Res., acs.accounts.8b00672 (2019).

All publications resulting from the use of this method must acknowledge the above.

Getting Started

The WF-in-DFT embedding program is called by the EMBED directive. When using this directive, the user must specify which MOs are associated with the active subsystem. This can be most directly accomplished through the ORBS option:

The KS calculation is performed on the full system. The EMBED directive indicates the beginning of an embedding calculation in which subsystem A corresponds to the single MO orbital1. The HF and CCSD(T) calculations are performed on subsystem A, embedded in a subsystem interaction potential calculated at the KS level.

Alternatively to setting the ORBS array, the user can instead specify the ATOMS option:

CASSCF calculations can be used as the WF method, but some care must be taken for the specification of the active space. In particular, the OCC, CLOSED, and FROZEN cards must be specified, and they must correspond to the number of OCC, CLOSED, and FROZEN orbitals within the active subsystem:

proc multi_proc

endproc

{ks}

{[atom1,atom2,...]}

Options

EMBED Options

The following options may be specified on the EMBED command line:

Density fitting

Density fitting may be enabled by specifying DF-EMBED in place of the EMBED command. It is strongly recommended that both EMBED and the mean-field calculation that preceded it be run with density fitting.

{df-embed,atoms=[atom1,atom2,...]}

Further density fitting options may be supplied to the CFIT directive (see section density fitting).

Saving orbitals and densities

The following options may be used to save densities and orbitals for later analysis. If AO truncation is employed (truncate), these densities and orbitals represent the truncated basis.

Improved efficiency using basis set truncation

By default, the calculation on the active subsystem (subsystem A) is performed using all basis functions employed by the calculation on the full system. The computational cost of a high-level WF calculation on even a small active subsystem can thus be prohibitively high if the full system is sufficiently large. This cost can be greatly reduced by performing the calculation on the active subsystem using only the basis functions that are most important for representation of the MOs associated with the active subsystem. The density threshold method is implemented to facilitate truncation of the active subsystem basis set, and this method is described in the subsections below.

Density threshold method

The density threshold truncation method requires just one parameter to be adjusted for all types of molecular system and functions that are important in the long range are automatically kept. At its low truncation limit it produces the full basis embedding answer or at high truncation it’s limit is a result constructed from just the functions on the active embedding region. Application of this method involves slightly different input options:

When using this method, the following should be taken into consideration:
(1) The accuracy of this truncation method is sensitive to the size of the basis set. For this reason, it is recommended to use at least a triple-zeta basis set.
(2) It is recommended to use Knizia’s IBO localization method, which reduces the orbital tails.
(3) For accurate HF-in-HF and DFT-in-DFT results, it is recommended to use a value for truncate of less than 0.0001. This is especially important when embedding across covalent bonds.
(4) The geometry must be in Cartesian coordinates.

All publications resulting from the use of this method must acknowledge the following:

S. J. Bennie, M. Stella, T. F. Miller III and F. R. Manby, J. Chem. Phys. 143, 024105 (2015).

Gradients for Projection-based WF-in-DFT Embedding

Analytical nuclear gradients have been implemented for projection-based wavefunction-in-DFT (WF-in-DFT) embedding with and without density threshold truncation. The current available methods that can be used for the WF method on the active subsystem (subsystem A) are CCSD(T), CCSD, MP2, and HF. The current available methods that can be used for the low-level SCF method are LDA, GGAs, hybrid GGAs and HF. ECP gradients should also work for WF-in-DFT both with and without density threshold truncation.

Important Notes:

  1. Only closed-shell gradients are supported right now.
  2. Embedding gradients currently only supports Pipek-Mezey localized orbitals such that the core and valence are localized together. Therefore, when calling Pipek-Mezey localizaton the core option must be specified to be 0 so core and valence orbitals are localized together.
  3. Embedding gradients do not include the projector correction so the option HF_COR must be set to 0.
  4. If density threshold truncation (density threshold method) is used a lowproc must be specified for the type-in-type correction.
  5. Density fitting embedding gradients are not supported right now.
  6. Please note the following embedding gradient methods have not yet been implemented: Multireference-WF-in-DFT, and DFT-in-DFT embedding.

Example input files for gradients are located under the examples section (examples)

All publications resulting from the use of this method must acknowledge the analytical WF-in-DFT nuclear gradients paper.

S. J. R. Lee, F. Ding, F. R. Manby, T. F. Miller III, arXiv e-prints, (2019), arXiv:1903.05830

Embedding Gradient Options

Embedding Gradient Best Practices

Examples

Projection-Based Embedding

This example performs a HF-in-HF and CCSD(T)-in-HF calculation on methanol. The HF-in-HF energy is numerically identical to the HF energy of the full system. By replacing the initial HF calculation on the full system with a DFT calculation, it is possible to perform HF-in-DFT and CCSD(T)-in-DFT calculations.

examples/embed_proj_methanol.inp
!methanol
memory,100,m
symmetry,nosym ! Embedding is not implemented with symmetry

r1=2.648205499 ! Bohr
r2=1.780518852 ! Bohr
r3=2.064847045 ! Bohr
a1=110.61344
d1=0
d2=119.9573

{gthresh,grid=1e-12,orbital=1e-8,coeff=1e-7}

geometry= {C;
           O,C,r1;
           H1,O,r2,C,a1;
           H2,C,r3,O,a1,H1,d1;
           H3,C,r3,O,a1,H1,d2;
           H4,C,r3,O,a1,H1,-d2}

basis,def2-svp

! Step 1: Perform a HF calculation on the full system.
{df-hf}

! Step 2: Perform orbital localization.
! Localized MOs are saved to record 3100.2
{ibba;save,3100.2}

! Step 3: Define the procedure that will be used for the embedded calculation.
!         Here, we define a standard CCSD(T) calculation.
proc embedded_ccsdt
    ! Step 3a: Call HF. In this procedure, only subsystem A will be calculated.
    ! This corresponds to a HF-in-HF calculation.
    {df-hf}

    ! Step 3b: Call CCSD(T). This will perform a CCSD(T)-in-HF calculation
    {df-ccsd(t)}
endproc

! Step 4: Perform the embedding calculation. This requires:
! Specification of embedded MOs
! Location of localized orbitals
! The procedure to run on the embedded subsystem A.
{df-embed,orbs=[1,6,7,8,9],highproc=embedded_ccsdt      ! The active subsystem corresponds to the indicated MOs
orbital,3100.2,local}             ! Use the localized MOs from record 3100.2

Embedded Multiconfigurational Calculations

This example performs an MCSCF-in-DFT calculation on butene.

examples/embed_mrci.inp
!butene
symmetry,nosym ! Embedding is not implemented with symmetry

{gthresh,grid=1e-12,orbital=1e-8,coeff=1e-7}

geometry={
C1,,          0.835432384 ,   0.049538050 ,  -1.305177530
C2,,         -0.187431703 ,  -1.125601622 ,   1.070665830
C3,,         -0.143380420 ,   0.687719249 ,   3.341106035
H1,,         -2.118667450 ,  -1.789600571 ,   0.735591607
H2,,          0.942690000 ,  -2.802115823 ,   1.537268881
H3,,         -0.860542916 ,  -0.240216354 ,   5.042940132
H4,,         -1.309266799 ,   2.354788168 ,   2.978399494
H5,,          1.778149207 ,   1.339614495 ,   3.742333311
C4,,         -0.440429873 ,   0.344384131 ,  -3.450328416
H6,,          2.782544497 ,   0.715697009 ,  -1.206788342
H7,,         -0.434433198 ,  -1.089606640 ,  -4.914726597
H8,,         -1.545390643 ,   2.035286250 ,  -3.818951212
} !Geometry is in Bohr

basis=def2-svp

! Step 1: Perform DFT on full system.
! Level shift added for improved convergence for this system.
!  The system is a butene where the double bond is twisted by 90 degrees.
{rks,b-lyp,shifta=-0.3}

! Step 2: Define the procedure that will be used for the embedded calculation.
!         Here, we define a standard MRCI calculation.
proc embedded_mrci
    ! Step 2a: Generate initial orbitals for MCSCF.
    {rhf}

    ! Step 2b: Perform MCSCF calculation.
    ! When performing an embedded MCSCF calculation, always specify the active space.
    {mcscf;
    config,csf;
    OCC,9;
    closed,7;
    frozen,2;}

    ! Step 2c: Add dynamic correlation
    {mrci}
endproc

! Step 3: Call embedding code and specify the atoms corresponding to the twisted double bond. Provide the procedure for the embedded calcuation.
{embed,atoms=[C1,C4,H6,H7,H8],highproc=embedded_mrci}

Density Threshold Truncation

This example performs a CCSD(T)-in-B3LYP calculation on a water dimer, using the density-threshold method with a truncate value of 0.001.

examples/embed_trunc2.inp
!water dimer
memory,150,m

{gthresh,orbital=1e-8,coeff=1e-7}

geometry={
O1,,          2.872302451  , -0.739837226  ,  2.587473481
H11,,         1.114872273  , -1.362093804  ,  2.841067168
H12,,         3.830013761  , -1.316171569  ,  4.018958028
O2,,          1.378439935  ,  3.695314083  , -0.113721828
H21,,         2.070708975  ,  2.651580550  ,  1.227452707
H22,,         2.064068477  ,  2.891071211  , -1.605148488
} !Geometry is in Bohr

Basis,avtz

! Step 1: Perform an LDA calculation on the full system
{df-ks}

! Step 2: Perform orbital localization.
{ibba;save,3100.2}

! Step 3: Define the procedure that will be used for the embedded calculation.
!         Here, we define a standard CCSD(T) calculation.
proc embedded_ccsdt
    ! Step 3a: Perform a HF-in-LDA calculation
    {df-hf;noenest}

    ! Step 3b: Perform a CCSD(T)-in-LDA calculation
    {df-ccsd(t)}
endproc

! Step 4: Run the embedding calculation
! using truncation with the density threshold method.
{df-embed,N_orbitals=5,truncate=0.0001,atoms=[O1,H11,H12],highproc=embedded_ccsdt;orbital,3100.2,local}

Open Shell System

This example performs a RCCSD(T)-in-LDA calculation on the methoxy radical.

examples/embed_open.inp
!Methoxy Radical
memory,100,M
zsymel=nosym  ! Embedding is not implemented with symmetry

{gthresh,grid=1e-12,orbital=1e-8,coeff=1e-7}

Geometry={
O1,,          2.328578323  , -0.516253350  , -0.241152468
H1,,         -3.930899693  ,  0.818199944  ,  0.136592671
C1,,         -2.292181882  , -0.433746123  , -0.018359872
H2,,         -2.391620554  , -1.802771746  ,  1.530291417
C2,,          0.160581025  ,  1.056429748  ,  0.096479104
H3,,          0.284618712  ,  2.116540479  ,  1.878957579
H4,,          0.248705545  ,  2.426749515  , -1.444737877
H5,,         -2.407046505  , -1.464278195  , -1.802909809
}! Geometry is in Bohr

basis,def2-svp

! Step 1: Perform a DFT calculation on the full system.
{ks,lda;wf,25,1,1}

! Step 2: Perform orbital localization.
{ibba}

! Step 3: Define the procedure that will be used for the embedded calculation.
!         Here, we define a standard CCSD(T) calculation.
proc embedded_ccsdt
    ! Step 3a: Call HF. In this procedure,  only subsystem A will be calculated.
    ! This corresponds to a HF-in-LDA calculation.
    {hf}

    ! Step 3b: Call CCSD(T). This will perform a RCCSD(T)-in-LDA calculation
    {rccsd(t)}
endproc

! Step 4: Perform the embedded calculation.
{embed,atoms=[O1],highproc=embedded_ccsdt}

CCSD(T)-in-LDA Gradient Calculation

This example performs a CCSD(T)-in-LDA gradient calculation on ethanol.

examples/embed_grad_ccsd_t-in-lda.inp
!strained_ethanol
memory,200,M
symmetry,nosym ! Embedding is not implemented with symmetry

{gthresh,orbital=1.0d-07}

! Geometry in bohr
geometry={
O1,,     -6.287074 ,  5.128142 ,  0.368844
H2,,     -4.331335 ,  5.076727 ,  0.829762
C3,,     -7.553333 ,  1.646329 ,  0.142339
C4,,     -5.608539 ,  2.964896 ,  0.000000
H5,,     -9.487950 ,  1.458382 , -0.789351
H6,,     -8.648844 ,  0.140577 ,  1.411194
H7,,     -7.564298 , -0.311707 , -0.867451
H8,,     -4.507431 ,  2.577084 , -1.650962
H9,,     -4.726024 ,  2.654073 ,  1.792503
}

basis={default,sto-3g}

{grid,name=NEESE,neese_index=4}

! Step 1: Perform a LDA calculation on the full system.
{rks,LDA;core,0}

! Step 2: Perform orbital localization.
!         Embedding gradients currently only supports Pipek-Mezey localized orbitals
!         such that the core and valence are localized together.
!         Therefore, core must be specified to be 0 so core and valence orbitals are
!         localized together.
{locali,pipek;core,0}

! Step 3: Define the procedure that will be used for the embedded calculation.
!         Here, we define a standard CCSD(T) gradient calculation (the {force} call
!         must be present).
!         HF_COR should be set to 0 when doing gradients so the projector correction
!         is not added.
!         The frozen-core approximation for correlated calculations can be turned on
!         by specifying the correct number of core orbitals in Subsystem A.
!         (1 core orbital for the Oxygen in this case)
proc hi_proc
   {HF,hf_cor=0;core,1}
   {CCSD(T);core,1;cphf,thrmin=1.d-7}
   {force}
endproc

! Step 4: Perform the embedding calculation. This requires:
!         Specification of the option grad to tell the embedding code gradients
!         are requested.
!         Specification of embedded atoms.
!         The procedure to run on the embedded subsystem A.
!         Any embedding options such as cpksthr which controls the level
!         of convergence of the CPKS solutions.
{embed,grad,highproc=hi_proc,atoms=[O1,H2],cpksthr=1.d-7}

! Step 5: The force command must specified after the embedding command so the
!         embedding gradient is calculated.
{force,gridgrad=1}

CCSD(T)-in-LDA Gradient Calculation w/ Density Threshold Truncation

This example performs a CCSD(T)-in-LDA gradient calculation with density threshold truncation on ethanol.

examples/embed_grad_ccsd_t-in-lda_ao_trunc.inp
!strained_ethanol
memory,200,M
symmetry,nosym ! Embedding is not implemented with symmetry

{gthresh,orbital=1.0d-07}

! Geometry in bohr
geometry={
O1,,     -6.287074 ,  5.128142 ,  0.368844
H2,,     -4.331335 ,  5.076727 ,  0.829762
C3,,     -7.553333 ,  1.646329 ,  0.142339
C4,,     -5.608539 ,  2.964896 ,  0.000000
H5,,     -9.487950 ,  1.458382 , -0.789351
H6,,     -8.648844 ,  0.140577 ,  1.411194
H7,,     -7.564298 , -0.311707 , -0.867451
H8,,     -4.507431 ,  2.577084 , -1.650962
H9,,     -4.726024 ,  2.654073 ,  1.792503
}

basis={
default,def2-svp
}

{grid,name=NEESE,neese_index=4}

! Step 1: Perform a LDA calculation on the full system.
{rks,LDA;core,0}

! Step 2: Perform orbital localization.
!         Embedding gradients currently only supports Pipek-Mezey localized orbitals
!         such that the core and valence are localized together.
!         Therefore, core must be specified to be 0 so core and valence orbitals are
!         localized together.
{locali,pipek;core,0}

! Step 3: Define the procedure that will be used for the type-in-type correction.
!         This procedure needs to call the same DFT method used in Step 1.
!         Except now the {force} call must be present so the type-in-type
!         contributions to the embedding gradient can be added.
proc low_proc
    {rks,LDA;core,0}
    {force,gridgrad=1}
endproc

! Step 4: Define the procedure that will be used for the embedded calculation.
!         Here, we define a standard CCSD(T) gradient calculation (the {force} call
!         must be present).
!         HF_COR should be set to 0 when doing gradients so the projector correction
!         is not added.
!         The frozen-core approximation for correlated calculations can be turned on
!         by specifying the correct number of core orbitals in Subsystem A.
!         (1 core orbital for the Oxygen in this case)
proc hi_proc
   {HF,hf_cor=0;core,1}
   {CCSD(T);core,1;cphf,thrmin=1.d-7}
   {force}
endproc

! Step 5: Perform the embedding calculation. This requires:
!         1) Specification of the option grad to tell the embedding code gradients
!         are requested.
!         2) Specification of embedded atoms.
!         3) The procedure to run on the embedded subsystem A.
{embed,grad,highproc=hi_proc,atoms=[O1,H2],lowproc=low_proc,truncate=1e-4,ao_per_atom}

! Step 6: The force command must specified after the embedding command so the
!         embedding gradient is calculated.
{force,gridgrad=1}

MP2-in-LDA Geometry Optimization

This example performs a MP2-in-LDA geometry optimization on ethanol.

examples/embed_grad_opt.inp
!ethanol
memory,200,M

{gthresh,grid=1d-12,orbital=1.0d-07}

! Original geometry
symmetry,nosym
geometry={
O1         -3.3322316672        2.7486714246        0.4634782128
H2         -2.4357007108        3.2239754945        0.4835540980
C3         -4.1946894574        0.5410193913       -0.0217296115
C4         -2.9359507722        1.4310420137        0.0489661919
H5         -4.9218723840        0.9675783876       -0.7301767123
H6         -4.6706714295        0.4744353658        0.9687956993
H7         -3.9246181899       -0.4721658034       -0.3558474595
H8         -2.4446716957        1.4212399976       -0.9554187232
H9         -2.2101426141        0.9539365224        0.7529059549
}

basis={
default,sto-3g
}

proc hi_proc1
   {HF,hf_cor=0;core,0;start,2102.2;save,2102.2}
   {MP2;core,0;cphf,thrmin=1.d-7}
   {force}
endproc

proc runembed
    {rks,LDA;core,0;start,2101.2;save,2101.2}
    {locali,pipek;core,0}
    {embed,grad,highproc=hi_proc1,orbs=[9,10,11,12,13]}
endproc

{optg,procedure=runembed}

MP2-in-LDA Geometry Optimization w/ Density Threshold Truncation

This example performs a MP2-in-LDA geometry optimization with density threshold truncation on ethanol.

examples/embed_grad_opt_ao_trunc.inp
!ethanol
memory,200,M

{gthresh,grid=1d-12,orbital=1.0d-07}

! Original geometry
symmetry,nosym
geometry={
O1         -3.3265927677        2.7508371961        0.4722099522
H2         -2.4299560718        3.2260275232        0.4808555391
C3         -4.1955133281        0.5390054989       -0.0225174814
C4         -2.9373667449        1.4340964566        0.0502923154
H5         -4.9209868031        0.9595827880       -0.7359982137
H6         -4.6761386977        0.4760110751        0.9660347037
H7         -3.9228360337       -0.4750976765       -0.3498843604
H8         -2.4504497951        1.4273894050       -0.9556425922
H9         -2.2107086788        0.9518805278        0.7491777878
}

basis={
default,sto-3g
}

! This works without specifying the start record for each calculation.
! However, I would highly recommend specifying the start and save dump records
! to ensure that it always works.
proc low_proc
    {rks,LDA,hf_cor=0;core,0;start,2103.2;save,2103.2}
    {force}
endproc

proc hi_proc1
    {HF,hf_cor=0;core,0;start,2104.2;save,2104.2}
    {MP2;core,0;cphf,thrmin=1.d-7}
    {force}
endproc

proc runembed
    {rks,LDA;core,0;start,2100.2;save,2100.2}
    {locali,pipek,thrpip=1.0d-14;core,0}
    {embed,grad,highproc=hi_proc1,lowproc=low_proc,truncate=0.01,orbs=[9,10,11,12,13]}
endproc

{optg,procedure=runembed}