# 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):

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

“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 SeH_{2}- H_{2}O, in which SeH_{2}constitutes an active subsystem (subsystem I) and H_{2}Ois an environment (subsystem II):

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. :eq:`fde_pg` 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 SeH_{2}- H_{2}O ‘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).