Differences
This shows you the differences between two versions of the page.
Next revision | Previous revision | ||
projection-based_wf-in-dft_embedding [2020/06/11 18:17] – external edit 127.0.0.1 | projection-based_wf-in-dft_embedding [2024/01/08 13:24] (current) – external edit 127.0.0.1 | ||
---|---|---|---|
Line 1: | Line 1: | ||
+ | ====== 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, | ||
+ | |||
+ | 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, [[https:// | ||
+ | |||
+ | T. A. Barnes, J. D. Goodpaster, F. R. Manby, and T. F. Miller III, [[https:// | ||
+ | |||
+ | 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, [[https:// | ||
+ | |||
+ | S. J. R. Lee, F. Ding, F. R. Manby, T. F. Miller III, arXiv e-prints, (2019), arXiv: | ||
+ | |||
+ | 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 '' | ||
+ | |||
+ | * '' | ||
+ | |||
+ | The KS calculation is performed on the full system. The '' | ||
+ | |||
+ | Alternatively to setting the '' | ||
+ | |||
+ | * '' | ||
+ | |||
+ | 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 '' | ||
+ | |||
+ | '' | ||
+ | |||
+ | '' | ||
+ | |||
+ | '' | ||
+ | |||
+ | '' | ||
+ | |||
+ | ===== Options ===== | ||
+ | |||
+ | ==== EMBED Options ==== | ||
+ | |||
+ | The following options may be specified on the '' | ||
+ | |||
+ | * '' | ||
+ | * '' | ||
+ | * '' | ||
+ | * '' | ||
+ | * '' | ||
+ | * '' | ||
+ | * '' | ||
+ | * '' | ||
+ | + E_{\text{DFT-in-DFT}}^{\text{whole basis}} | ||
+ | - E_{\text{DFT-in-DFT}}^{\text{trunc. basis}}$$ It is important not to specify the number of electrons or spin state in this procedure, as this is determined automatically from the density partioning specified through other options. | ||
+ | * '' | ||
+ | * '' | ||
+ | * '' | ||
+ | |||
+ | ==== SCF options related to EMBED ==== | ||
+ | |||
+ | * '' | ||
+ | |||
+ | ==== Density fitting ==== | ||
+ | |||
+ | Density fitting may be enabled by specifying '' | ||
+ | |||
+ | '' | ||
+ | |||
+ | Further density fitting options may be supplied to the '' | ||
+ | |||
+ | ==== Saving orbitals and densities ==== | ||
+ | |||
+ | The following options may be used to save densities and orbitals for later analysis. If AO truncation is employed ('' | ||
+ | |||
+ | * '' | ||
+ | * '' | ||
+ | * '' | ||
+ | |||
+ | ===== 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 '' | ||
+ | (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, [[https:// | ||
+ | |||
+ | ===== 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:** | ||
+ | |||
+ | - Only closed-shell gradients are supported right now. | ||
+ | - 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. | ||
+ | - Embedding gradients do not include the projector correction so the option '' | ||
+ | - If density threshold truncation ([[projection-based WF-in-DFT embedding# | ||
+ | - Density fitting embedding gradients are not supported right now. | ||
+ | - Please note the following embedding gradient methods have not yet been implemented: | ||
+ | |||
+ | Example input files for gradients are located under the examples section ([[projection-based WF-in-DFT embedding# | ||
+ | |||
+ | 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: | ||
+ | |||
+ | ==== Embedding Gradient Options ==== | ||
+ | |||
+ | * '' | ||
+ | * '' | ||
+ | |||
+ | ==== Embedding Gradient Best Practices ==== | ||
+ | |||
+ | * We have found the best agreement between numerical gradients and analytical gradients when using the Neese (NEESE) grids. | ||
+ | * Grid weight derivatives should be calculated for all DFT quantities. | ||
+ | |||
+ | ===== 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, | ||
+ | |||
+ | <code - examples/ | ||
+ | !methanol | ||
+ | memory, | ||
+ | symmetry, | ||
+ | |||
+ | r1=2.648205499 ! Bohr | ||
+ | r2=1.780518852 ! Bohr | ||
+ | r3=2.064847045 ! Bohr | ||
+ | a1=110.61344 | ||
+ | d1=0 | ||
+ | d2=119.9573 | ||
+ | |||
+ | {gthresh, | ||
+ | |||
+ | geometry= {C; | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | |||
+ | basis, | ||
+ | |||
+ | ! 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; | ||
+ | |||
+ | ! 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, | ||
+ | orbital, | ||
+ | </ | ||
+ | |||
+ | ==== Embedded Multiconfigurational Calculations ==== | ||
+ | |||
+ | This example performs an MCSCF-in-DFT calculation on butene. | ||
+ | |||
+ | <code - examples/ | ||
+ | !butene | ||
+ | symmetry, | ||
+ | |||
+ | {gthresh, | ||
+ | |||
+ | geometry={ | ||
+ | C1,, 0.835432384 , | ||
+ | C2,, | ||
+ | C3,, | ||
+ | H1,, | ||
+ | H2,, 0.942690000 , -2.802115823 , | ||
+ | H3,, | ||
+ | H4,, | ||
+ | H5,, 1.778149207 , | ||
+ | C4,, | ||
+ | H6,, 2.782544497 , | ||
+ | H7,, | ||
+ | H8,, | ||
+ | } !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, | ||
+ | |||
+ | ! 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, | ||
+ | {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, | ||
+ | </ | ||
+ | |||
+ | ==== 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. | ||
+ | |||
+ | <code - examples/ | ||
+ | !water dimer | ||
+ | memory, | ||
+ | |||
+ | {gthresh, | ||
+ | |||
+ | geometry={ | ||
+ | O1,, 2.872302451 | ||
+ | H11,, | ||
+ | H12,, | ||
+ | O2,, 1.378439935 | ||
+ | H21,, | ||
+ | H22,, | ||
+ | } !Geometry is in Bohr | ||
+ | |||
+ | Basis,avtz | ||
+ | |||
+ | ! Step 1: Perform an LDA calculation on the full system | ||
+ | {df-ks} | ||
+ | |||
+ | ! 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: Perform a HF-in-LDA calculation | ||
+ | {df-hf; | ||
+ | |||
+ | ! 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, | ||
+ | </ | ||
+ | |||
+ | ==== Open Shell System ==== | ||
+ | |||
+ | This example performs a RCCSD(T)-in-LDA calculation on the methoxy radical. | ||
+ | |||
+ | <code - examples/ | ||
+ | !Methoxy Radical | ||
+ | memory, | ||
+ | zsymel=nosym | ||
+ | |||
+ | {gthresh, | ||
+ | |||
+ | Geometry={ | ||
+ | O1,, 2.328578323 | ||
+ | H1,, | ||
+ | C1,, | ||
+ | H2,, | ||
+ | C2,, 0.160581025 | ||
+ | H3,, 0.284618712 | ||
+ | H4,, 0.248705545 | ||
+ | H5,, | ||
+ | }! Geometry is in Bohr | ||
+ | |||
+ | basis, | ||
+ | |||
+ | ! Step 1: Perform a DFT calculation on the full system. | ||
+ | {ks, | ||
+ | |||
+ | ! 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, | ||
+ | ! 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, | ||
+ | </ | ||
+ | |||
+ | ==== CCSD(T)-in-LDA Gradient Calculation ==== | ||
+ | |||
+ | This example performs a CCSD(T)-in-LDA gradient calculation on ethanol. | ||
+ | |||
+ | <code - examples/ | ||
+ | !strained_ethanol | ||
+ | memory, | ||
+ | symmetry, | ||
+ | |||
+ | {gthresh, | ||
+ | |||
+ | ! Geometry in bohr | ||
+ | geometry={ | ||
+ | O1,, | ||
+ | H2,, | ||
+ | C3,, | ||
+ | C4,, | ||
+ | H5,, | ||
+ | H6,, | ||
+ | H7,, | ||
+ | H8,, | ||
+ | H9,, | ||
+ | } | ||
+ | |||
+ | basis={default, | ||
+ | |||
+ | {grid, | ||
+ | |||
+ | ! Step 1: Perform a LDA calculation on the full system. | ||
+ | {rks, | ||
+ | |||
+ | ! Step 2: Perform orbital localization. | ||
+ | ! | ||
+ | ! such that the core and valence are localized together. | ||
+ | ! | ||
+ | ! | ||
+ | {locali, | ||
+ | |||
+ | ! 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). | ||
+ | ! | ||
+ | ! 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 | ||
+ | | ||
+ | | ||
+ | | ||
+ | endproc | ||
+ | |||
+ | ! Step 4: Perform the embedding calculation. This requires: | ||
+ | ! | ||
+ | ! are requested. | ||
+ | ! | ||
+ | ! 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, | ||
+ | |||
+ | ! Step 5: The force command must specified after the embedding command so the | ||
+ | ! | ||
+ | {force, | ||
+ | </ | ||
+ | |||
+ | ==== 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. | ||
+ | |||
+ | <code - examples/ | ||
+ | !strained_ethanol | ||
+ | memory, | ||
+ | symmetry, | ||
+ | |||
+ | {gthresh, | ||
+ | |||
+ | ! Geometry in bohr | ||
+ | geometry={ | ||
+ | O1,, | ||
+ | H2,, | ||
+ | C3,, | ||
+ | C4,, | ||
+ | H5,, | ||
+ | H6,, | ||
+ | H7,, | ||
+ | H8,, | ||
+ | H9,, | ||
+ | } | ||
+ | |||
+ | basis={ | ||
+ | default, | ||
+ | } | ||
+ | |||
+ | {grid, | ||
+ | |||
+ | ! Step 1: Perform a LDA calculation on the full system. | ||
+ | {rks, | ||
+ | |||
+ | ! Step 2: Perform orbital localization. | ||
+ | ! | ||
+ | ! such that the core and valence are localized together. | ||
+ | ! | ||
+ | ! | ||
+ | {locali, | ||
+ | |||
+ | ! 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. | ||
+ | ! | ||
+ | ! | ||
+ | proc low_proc | ||
+ | {rks, | ||
+ | {force, | ||
+ | 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). | ||
+ | ! | ||
+ | ! 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 | ||
+ | | ||
+ | | ||
+ | | ||
+ | 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, | ||
+ | |||
+ | ! Step 6: The force command must specified after the embedding command so the | ||
+ | ! | ||
+ | {force, | ||
+ | </ | ||
+ | |||
+ | ==== MP2-in-LDA Geometry Optimization ==== | ||
+ | |||
+ | This example performs a MP2-in-LDA geometry optimization on ethanol. | ||
+ | |||
+ | <code - examples/ | ||
+ | !ethanol | ||
+ | memory, | ||
+ | |||
+ | {gthresh, | ||
+ | |||
+ | ! Original geometry | ||
+ | symmetry, | ||
+ | geometry={ | ||
+ | O1 | ||
+ | H2 | ||
+ | C3 | ||
+ | C4 | ||
+ | H5 | ||
+ | H6 | ||
+ | H7 | ||
+ | H8 | ||
+ | H9 | ||
+ | } | ||
+ | |||
+ | basis={ | ||
+ | default, | ||
+ | } | ||
+ | |||
+ | proc hi_proc1 | ||
+ | | ||
+ | | ||
+ | | ||
+ | endproc | ||
+ | |||
+ | proc runembed | ||
+ | {rks, | ||
+ | {locali, | ||
+ | {embed, | ||
+ | endproc | ||
+ | |||
+ | {optg, | ||
+ | </ | ||
+ | |||
+ | ==== MP2-in-LDA Geometry Optimization w/ Density Threshold Truncation ==== | ||
+ | |||
+ | This example performs a MP2-in-LDA geometry optimization with density threshold truncation on ethanol. | ||
+ | |||
+ | <code - examples/ | ||
+ | !ethanol | ||
+ | memory, | ||
+ | |||
+ | {gthresh, | ||
+ | |||
+ | ! Original geometry | ||
+ | symmetry, | ||
+ | geometry={ | ||
+ | O1 | ||
+ | H2 | ||
+ | C3 | ||
+ | C4 | ||
+ | H5 | ||
+ | H6 | ||
+ | H7 | ||
+ | H8 | ||
+ | H9 | ||
+ | } | ||
+ | |||
+ | basis={ | ||
+ | default, | ||
+ | } | ||
+ | |||
+ | ! 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, | ||
+ | {force} | ||
+ | endproc | ||
+ | |||
+ | proc hi_proc1 | ||
+ | {HF, | ||
+ | {MP2; | ||
+ | {force} | ||
+ | endproc | ||
+ | |||
+ | proc runembed | ||
+ | {rks, | ||
+ | {locali, | ||
+ | {embed, | ||
+ | endproc | ||
+ | |||
+ | {optg, | ||
+ | </ | ||