NMR shieldings in Frozen Density Embedding (FDE) scheme with London atomic orbitals (LAOs)

In this tutorial we will show how to calculate the NMR shielding tensor with London atomic orbitals within the FDE scheme.

NMR shielding tensor in FDE scheme - theory overview

In order to calculate the NMR shielding tensor with LAOs, one must take into account modifications of the property gradient due to the use of perturbation-dependent basis sets (see [Ilias2009]). In FDE scheme, further modifications are needed, since the property gradient includes not only terms from isolated subsystems (here marked as “I” and “II”), but also a contribution from the interaction energy between these subsystems (see [Hofener2013]).

For instance, these new contributions to the property gradient of subsystem “I” can be written as (see :cite: olejniczak2017):

(1)\[E_{[1],emb}^{[B]} = - \int v_{emb}^I \Omega_{ia}^{B,I} - \iint w_{emb}^{I,I} \Omega_{ia}^{I} \Omega_{jj}^{B,I} - \iint w_{emb}^{I,II} \Omega_{ia}^{I} \Omega_{jj}^{B,II}\]

where \(\Omega_{pq}^{B}\) is the first-order perturbed overlap distribution, which in a basis of London orbitals consists of two terms

(2)\[\Omega_{pq}^{B} = \frac{i}{2} (R_{MN} \times r) \Omega_{pq} + (T_{pt}^{B*} \Omega_{tq} + \Omega_{pt}T_{tq}^B),\]

“direct” LAO term (the first term) and “reorthonormalization” term (in parenthesis). The \(v_{emb}^I\) is the embedding potential of subsystem I, \(w_{emb}^{I,I}\) is the embedding kernel involving subsystem I and \(w_{emb}^{I,II}\) is the embedding kernel coupling subsystem I and subsystem II.

In this tutorial we will show how to include each of these contributions to the property gradient in the calculations of the NMR shielding tensor. We will also demonstrate how to calculate and plot the first-order perturbed current density induced by magnetic field and how to use it in analysis of the NMR shielding tensor.

A sample calculation of NMR shielding tensor within FDE

We will consider a simple system of SeH2- H2O, in which SeH2constitutes an active subsystem (subsystem I) and H2Ois an environment (subsystem II):

../../_images/h2se_h2o.jpg

The computational protocol is different depending on whether the coupling between two subsystems is included or not (through the term dependent on \(w_{emb}^{I,II}\)). To avoid confusion, we will describe these two protocols separately.

We will need two molecular input files. File h2se.mol, corresponding to an active subsystem I, reads:

DIRAC


C   2    0         A
       34.    1
Se1      -0.003989   -0.034416    0.000005
LARGE BASIS cc-pVDZ
        1.    2
H2        0.089169    1.434145   -0.001364
H3        1.458406   -0.164893    0.001741
LARGE BASIS cc-pVDZ
FINISH

and file h2o.mol, standing for a subsystem II (‘environment’):

DIRAC


C   2    0         A
        1.    2
H4       -0.257088    4.261903    0.765747
H6       -0.259151    4.266296   -0.762573
LARGE BASIS cc-pVDZ
        8.    1
O5        0.036856    3.761402   -0.000251
LARGE BASIS cc-pVDZ
FINISH

In this example we will use the four-component relativistic DC Hamiltonian and DFT/LDA with small grid and small basis sets.

Calculations without the FDE coupling terms

Here we demonstrate how to include the first two terms of Eq. (1) in the property gradient: terms dependent on the embedding potential and the uncoupled embedding kernel.

In this tutorial, we use a coarse grid generated in DIRAC in preceeding calculations on SeH2- H2O ‘supermolecule’, which was then exported on a numerical_grid file. In DIRAC the numerical_grid file, which can be retrieved from DFT calculations, is divided into batches of points. However, FDE module is able to read only the continuous list of grid points, therefore such numerical_grid file has to be reformatted. This can be done by the script which can be found in utils directory of DIRAC, dft_to_fde_grid_convert.py, which should be called from the directory where the numerical_grid file is:

# make sure that the "utils/fde-mag" in DIRAC's directory is in your PATH, then:

./dft_to_fde_grid_convert.py

This script will overwrite the numerical_grid and save its original version as numerical_grid.ORIG.

Step 1: SCF calculation on the subsystem II

The first step is the calculation on subsystem II, the aim of which is to export its density and electrostatic potential on a text file.

The menu file getfrozen.inp reads:

**DIRAC
.WAVE F
**HAMILTONIAN
.DFT
 LDA
.FDE
*FDE
.EXONLY
FILEEX
.GRIDOU
GRIDOUT
.WRFRMT
TXT
**GRID
.IMPORT
numerical_grid
**INTEGRALS
.NUCMOD
1
*READIN
.UNCONTRACT
**WAVE FUNCTIONS
.SCF
*SCF
.EVCCNV
1.0E-08  1.0E-06
**ANALYZE
*END OF

and we make sure to get the text file GRIDOUT (related to the .GRIDOU keyword). We should also remember to supply two files, related to keywords .GRIDOU (GRIDOUT) and .EXONLY (FILEEX), which are copies of a numerical grid file:

cp numerical_grid GRIDOUT
cp numerical_grid FILEEX

pam --inp=getfrozen.inp --mol=h2o.mol --put="numerical_grid FILEEX GRIDOUT" --get="GRIDOUT" --outcmo

mv DFCOEF   DFCOEF.h2o
mv GRIDOUT  GRIDOUT.h2o

After this step, the GRIDOUT file will be updated: additional columns on that file correspond to the electrostatic potential and the density of subsystem II calculated in each point of a grid.

Step 2: SCF and response calculation on the subsystem I

Now we can use the information about the subsystem II written on the GRIDOUT file for constructing the embedding potential, \(v_{emb}^{I}\), and embedding kernel, \(w_{emb}^{I,I}\). In order to do that, we use the GRIDOUT file saved from Step 1 and copy it to a FROZEN file, which is given as an argument to the .FRDENS keyword.

We use the following input file for the SCF step:

**DIRAC
.WAVE F
**HAMILTONIAN
.DFT
 LDA
.FDE
*FDE
.FRDENS
FROZEN
.RDFRMT
TXT
.UPDATE
.GRIDOU
GRIDOUT
.WRFRMT
TXT
**GRID
.IMPORT
numerical_grid
**INTEGRALS
.NUCMOD
1
*READIN
.UNCONTRACT
**WAVE FUNCTIONS
.SCF
*SCF
.EVCCNV
1.0E-08  1.0E-06
*END OF

and the pam command is:

cp GRIDOUT.h2o FROZEN
cp GRIDOUT.h2o GRIDOUT

pam --inp=updatefrozen.inp --mol=h2se.mol --put="numerical_grid FROZEN GRIDOUT" --get="GRIDOUT" --outcmo

mv GRIDOUT GRIDOUT.h2se
mv DFCOEF  DFCOEF.h2se

To calculate the NMR shielding tensor, we use the following input file:

**DIRAC
.PROPERTIES
**GENERAL
.RKBIMP
**GRID
.IMPORT
numerical_grid
**HAMILTONIAN
.URKBAL
.DFT
 LDA
.FDE
*FDE
.FRDENS
FROZEN
.RDFRMT
TXT
.UPDATE
.RSP
.LAO11
**INTEGRALS
.NUCMOD
1
*READIN
.UNCONTRACT
**WAVE FUNCTIONS
.SCF
*SCF
.EVCCNV
1.0E-08  1.0E-05
**PROPERTIES
.PRINT
2
.SHIELDING
*NMR
.LONDON
.SYMCON
*LINEAR RESPONSE
.PRINT
2
.THRESH
1.0E-07
.MAXITR
100
*END OF

in which the keyword .LAO11 denotes that the contributions from the embedding potential, \(v_{emb}^{I}\), and the uncoupled embedding kernel, \(w_{emb}^{I,I}\), will be added to the property gradient. The corresponding pam command is:

cp GRIDOUT.h2o  FROZEN
cp DFCOEF.h2se  DFCOEF

pam --inp=shield_v_w11.inp --mol=h2se.mol --put="numerical_grid FROZEN" --incmo

This is the advised option. However, as the calculations of the uncoupled embedding kernel terms in the property gradient may be expensive, one may want to restrict the FDE-LAO contributions to the part dependent on the embedding potential only. In this case, the input file reads:

**DIRAC
.PROPERTIES
**GENERAL
.RKBIMP
**GRID
.IMPORT
numerical_grid
**HAMILTONIAN
.URKBAL
.DFT
 LDA
.FDE
*FDE
.FRDENS
FROZEN
.RDFRMT
TXT
.UPDATE
.RSP
.LDPT
**INTEGRALS
.NUCMOD
1
*READIN
.UNCONTRACT
**WAVE FUNCTIONS
.SCF
*SCF
.EVCCNV
1.0E-08  1.0E-05
**PROPERTIES
.PRINT
2
.SHIELDING
*NMR
.LONDON
.SYMCON
*LINEAR RESPONSE
.PRINT
2
.THRESH
1.0E-07
.MAXITR
100
*END OF

and the corresponding pam command is:

cp GRIDOUT.h2o  FROZEN
cp DFCOEF.h2se  DFCOEF

pam --inp=shield_v.inp --mol=h2se.mol --put="numerical_grid FROZEN" --incmo

Additionally, we may want to save the PAMXVC and TBMO files, if we are interested in visualization of the NMR shielding density and/or magnetically-induced current density. The calculation on our test molecule give the following results for the isotropic part of the shielding tensor of Se atom:

terms included in property gradient

value

None

2469.8091

v

2370.3941

v+w11

2370.2488

None in the table denotes calculations performed with no “new” FDE-LAO terms in the property gradient (the gauge origin is in the center of mass of SeH2). Also, in this test case the value obtained with the embedding potential and the uncoupled embedding kernel in the property gradient (v+w11) is very close to the result of calculations with the embedding potential only (v), but for heavier nuclei the difference between these two setups may be significant.

Calculations with the FDE coupling terms

Here we will demonstrate how to include the coupling between two subsystems in the property gradient (term dependent on \(w_{emb}^{I,II}\)).

Step 1: SCF and response calculation on the subsystem II

In the first step we perform the calculations on subsystem II, the aim of which is to get its density and electrostatic potential written on a text file. In addition, we will export the perturbed density of subsystem II. We will need the DFCOEF file of subsystem II, which can be taken from getfrozen.inp run with the following input:

**DIRAC
.WAVE F
**HAMILTONIAN
.DFT
 LDA
.FDE
*FDE
.EXONLY
FILEEX
.GRIDOU
GRIDOUT
.WRFRMT
TXT
**GRID
.IMPORT
numerical_grid
**INTEGRALS
.NUCMOD
1
*READIN
.UNCONTRACT
**WAVE FUNCTIONS
.SCF
*SCF
.EVCCNV
1.0E-08  1.0E-06
**ANALYZE
*END OF

with the pam command:

cp numerical_grid GRIDOUT
cp numerical_grid FILEEX

pam --inp=getfrozen.inp --mol=h2o.mol --put="numerical_grid FILEEX GRIDOUT" --get="GRIDOUT" --outcmo

mv DFCOEF   DFCOEF.h2o
mv GRIDOUT  GRIDOUT.h2o

(exactly as done in calculations without the FDE coupling terms presented above). From this step we save the GRIDOUT.h2o file. In order to get the density perturbed by an external magnetic field in LAO basis, the following input can be used:

**DIRAC
.PROPERTIES
**GENERAL
.RKBIMP
**HAMILTONIAN
.URKBAL
.DFT
 LDA
**GRID
.IMPORT
numerical_grid
**INTEGRALS
.NUCMOD
1
*READIN
.UNCONTRACT
**WAVE FUNCTIONS
.SCF
*SCF
.EVCCNV
 1.0D-08 1.0D-7
**PROPERTIES
.PRINT
2
.SHIELDING
*NMR
.LONDON
.SYMCON
.EXPPED
*LINEAR
.PRINT
2
.THRESH
1.0D-08
*END OF

with the keyword .EXPPED asking the program to save the contributions to the perturbed density on two files pertden_direct_lao.FINAL and pertden_reorth_lao.FINAL. The corresponding pam command is:

cp DFCOEF.h2o DFCOEF
pam --inp=getfrozenpert.inp --mol=h2o.mol --incmo --put="numerical_grid" --get="pertden*"

NOTE: This step has to be run in serial mode for the moment.

Step 2: Response calculation on the subsystem I

Now we can finally calculate the FDE contributions to the NMR shielding tensor of an active subsystem (H2Se) which depend on the embedding kernel including the coupling terms between two subsystems. If the SCF step on an active subsystem was not done before, we should run the following input:

**DIRAC
.WAVE F
**HAMILTONIAN
.DFT
 LDA
.FDE
*FDE
.FRDENS
FROZEN
.RDFRMT
TXT
.UPDATE
.GRIDOU
GRIDOUT
.WRFRMT
TXT
**GRID
.IMPORT
numerical_grid
**INTEGRALS
.NUCMOD
1
*READIN
.UNCONTRACT
**WAVE FUNCTIONS
.SCF
*SCF
.EVCCNV
1.0E-08  1.0E-06
*END OF

with the pam command:

cp GRIDOUT.h2o FROZEN
cp GRIDOUT.h2o GRIDOUT

pam --inp=updatefrozen.inp --mol=h2se.mol --put="numerical_grid FROZEN GRIDOUT" --get="GRIDOUT" --outcmo

mv GRIDOUT GRIDOUT.h2se
mv DFCOEF  DFCOEF.h2se

The response input is then the following:

**DIRAC
.PROPERTIES
**GENERAL
.RKBIMP
**GRID
.IMPORT
numerical_grid
**HAMILTONIAN
.URKBAL
.DFT
 LDA
.FDE
*FDE
.FRDENS
FROZEN
.RDFRMT
TXT
.UPDATE
.RSP
.LAO11
.LAOFRZ
.PERTIM
pertden_direct_lao.FINAL
pertden_reorth_lao.FINAL
**INTEGRALS
.NUCMOD
1
*READIN
.UNCONTRACT
**WAVE FUNCTIONS
.SCF
*SCF
.EVCCNV
1.0E-08  1.0E-05
**PROPERTIES
.PRINT
2
.SHIELDING
*NMR
.LONDON
.SYMCON
*LINEAR RESPONSE
.PRINT
2
.THRESH
1.0E-07
.MAXITR
100
*END OF

with two new keywords: .LAOFRZ and .PERTIM. The corresponding pam command is:

cp DFCOEF.h2se  DFCOEF
cp GRIDOUT.h2o  FROZEN

pam --inp=shield_v_w11_w12_nocoulomb.inp  --mol=h2se.mol --put="numerical_grid FROZEN pertden_direct_lao.FINAL pertden_reorth_lao.FINAL" --incmo

This input does not include the computationally expensive Coulomb term in the embedding kernel, only the nonadditive exchange-correlation and kinetic parts of the embedding kernel. If the Coulomb term is to be calculated then an additional keyword is required in the input, .LFCOUL. Instead of using the .LFCOUL and the .LAOFRZ keywords together, it is possible to use a single .LAO12 keyword, which represents both contributions to the coupled embedding kernel:

**DIRAC
.PROPERTIES
**GENERAL
.RKBIMP
**GRID
.IMPORT
numerical_grid
**HAMILTONIAN
.URKBAL
.DFT
 LDA
.FDE
*FDE
.FRDENS
FROZEN
.RDFRMT
TXT
.UPDATE
.RSP
.LAO11
.LAO12
.PERTIM
pertden_direct_lao.FINAL
pertden_reorth_lao.FINAL
**INTEGRALS
.NUCMOD
1
*READIN
.UNCONTRACT
**WAVE FUNCTIONS
.SCF
*SCF
.EVCCNV
1.0E-08  1.0E-05
**PROPERTIES
.PRINT
2
.SHIELDING
*NMR
.LONDON
.SYMCON
*LINEAR RESPONSE
.PRINT
2
.THRESH
1.0E-07
.MAXITR
100
*END OF

with an analogous pam command:

cp DFCOEF.h2se  DFCOEF
cp GRIDOUT.h2o  FROZEN

pam --inp=shield_v_w11_w12.inp  --mol=h2se.mol --put="numerical_grid FROZEN pertden_direct_lao.FINAL pertden_reorth_lao.FINAL" --incmo

This input contains all the necessary contributions to the property gradient. Additionally, we may want to save the PAMXVC and TBMO files, if we are interested in visualization of the NMR shielding density and/or magnetically-induced current density. With this setup the calculations on the test system gave the following results:

terms included in property gradient

value

v+w11+w12, no Coulomb

2370.2536

v+w11+w12, ALL terms

2370.2536

Here the two results are the same within 0.0001 ppm accuracy, as the contributions from the Coulomb embedding contribution to the perturbed Fock matrix are very small (can be printed in the output with the print level set at least to 5 under **PROPERTIES section).