Charge Displacement
*******************

ADF Bonding Analysis Charge Displacement

The Charge-Displacement (CD) analysis is a simple method able of
analyzing and quantifying the amount of charge transfer that occurs
upon a bond formation between two interacting fragments. This tutorial
shows three different examples how to calculate the Charge-
Displacement function.

The CD function Δq(z) is defined as the partial integration along a
chosen axis, z, of the electron density rearrangement Δρ(x,y,z’)
taking place upon the bond formation:

   \Delta q(z) = \int_{-\infty}^\infty dx \int_{-\infty}^\infty dy
   \int_{-\infty}^z dz' \Delta \rho (x,y,z')

It can be also used in combination with the NOCV theory (where the
electron density rearrangement is defined with respect to the
orthonormalized non-interacting fragments). More information about
Natural Orbitals for Chemical Valence can be found at ETS-NOCV and in
the EDA-NOCV tutorial. In this case the deformation density can be
expressed as a sum of contributions (Δρ_k) defined by the NOCV pairs
of complementary eigenfunctions (Ψ_-k and Ψ_k) corresponding to
eigenvalues -ν_-k and ν_k with the same absolute value but opposite
sign:

   \Delta \rho_k(z) = v_k ( -| \Psi _{-k}|^2 + | \Psi _{+k}|^2 )

and the associated CD function (NOCV-CD) is simply defined as
following:

   \Delta q_k(z) = \int_{-\infty}^\infty dx \int_{-\infty}^\infty dy
   \int_{-\infty}^z dz' \Delta \rho_k (x,y,z')

This tutorial shows how to build the chemical integrands for
evaluating the CD function in three different relevant examples:

1. Quantifying the CT contribution in closed-shell interacting
   fragments: case of noble metal-noble gas interaction in the Xe-AuF
   complex. [1]

2. Quantifying of the Dewar-Chatt-Duncanson bonding components
   (σ-donation and π-backdonation) in transition metal complex
   applying the CD function in combination of the NOCV theory (NOCV-
   CD). [2]

3. Application to open-shell systems to reveal competing pathways of
   proton-coupled electron transfer in oxidation catalysis (hydrogen-
   atom transfer, concerted proton-coupled electron transfer). [3]

In all cases the basic idea is to evaluate the different integrands
(Δρ, Δρ_k or Δρ_α/β) using *densf* or directly from AMS-GUI which
allow for their efficient mapping of a grid (CUBE  format).


Setting Up PYCUBESCD
====================

The CD function is evaluated numerically on the grid. This requires
that you first download the PYCUBESCD suite of programs which is
freely available (LGPL-3.0 License at GitHub).

You can download PYCUBESCD using the git command line interface (if
available), by running the following command in the desired install
location:

   git clone https://github.com/BERTHA-4c-DKS/pycubescd.git

Alternatively, you can manually download the zip file and extract it.

After downloading PYCUBESCD, set the environment variable "$PYCUBESCD"
to point to the directory where the PYCUBESCD program is located. This
will be called "pycubescd" if downloaded via the git CLI, or
"pycubescd-master" if extracted from the zip file.

   export PYCUBESCD=/path/to/pycubescd

You are now ready to use PYCUBESCD throughout this tutorial.


Case 1: CT in closed-shell interacting fragments: Xe—AuF
========================================================

The molecular structures needed for this tutorial are given "here".
The molecular structures are properly aligned with the molecular axis
along the z-axis. In general, one should optimize the structures first
and align with the z-axis coinciding with the direction chosen for
evaluating the CD function.

Note:

  The alignment is important. The z-axis defines the axis on which is
  performed the partial integration for evaluating the CD function.
  The fragments (Xe and AuF) have the geometry they possess within the
  complex.


1. Fragment Calculation
-----------------------

At first, a fragment calculation must be performed. The following
steps from the Fragment Analysis tutorial can be followed. A script
including the fragment calculation (step 1) and step 2 will be
provided in step 2.

   Start AMSinput and insert the "Xe-AuF structure".
   In the Main ADF Panel, select the **Single Point** task
   Select the **GGA:BLYP** XC functional
   Make sure the Relativity is set at **Scalar**
   Select a **TZ2P** Basis set, a **Small** Frozen core
   Select a **Good** Numerical quality
   In the Panel bar, select **Model** → **Regions**
   Select Xe and click the [image: AddButton][image] to add it as region
   Select the Au and F and add them as a second region
   In the Panel bar, select **MultiLevel** → **Fragments**
   Check the **Use fragment** box
   Save and run the calculation

   [image]

The ".rkf" file of the final fragment analysis (the "adf.rkf" in the
"XeAuF.results" folder) is needed in the following step.


2. Generate CD function
-----------------------

-[ CD function via DENSF and PYCUBESCD ]-

The *densf* program must be used, which is an auxiliary program
distributed with AMS. See the documentation of Densf: Volume Maps to
generate a numerical representation of the electron densities
associated with the XeAuF complex and the sum of the isolated Xe and
AuF fragments. There, details for the definition of the grid can be
found. We will map these quantities in the CUBE format. The only
information necessary from the previous fragment calculation (Step 1)
is the rkf file.

The following script can be used to generate a CUBE file using *densf*
for the SCF electron density (**Density SCF**):

   "$AMSBIN"/densf << eor
   InputFile XeAuF.results/adf.rkf
   Density SCF                      ! We require the SCF density to be mapped
   CUBOutput XeAuF                  ! Output is the file: XeAuF%SCF%Density.cub
   Grid
      -10.0 -10.0 -10.0
      100  100  100
      1.0 0.0 0.0 15.03
      0.0 1.0 0.0 15.03
      0.0 0.0 1.0 15.03
   End
   eor

A CUBE file can also be generated using *densf* for the sum of the SCF
electron densities of the isolated Xe and AuF fragments (**Density
SumFrag**):

   "$AMSBIN"/densf << eor
   InputFile XeAuF.results/adf.rkf
   CUBInput
   Density SumFrag                  ! We require the Sum of SCF densities of the isolated Fragments
   CUBOutput XeAuF                  ! Output is the file: XeAuF%SumFrag%Density.cub
   Grid
      -10.0 -10.0 -10.0
      100  100  100
      1.0 0.0 0.0 15.03
      0.0 1.0 0.0 15.03
      0.0 0.0 1.0 15.03
   End
   eor

Note:

  All densities **must** be mapped on the same numerical grid. This
  has been explicitly defined after the grid keyword in the densf
  input. This is important if you use the PYCUBESCD suite of programs
  to calculate the CD function in Step 3.

Now we use the PYCUBESCD suite of programs to make operations between
Cube files and to generate the CD function (see set up instructions).

Run the following command in order to generate the file
"diff_XeAuF.cub", which contains the numerical representation of the
electron density difference between the XeAuF molecule (cube file:
"XeAuF%SCF%Density.cub") and the sum of the electron density of the Xe
and AuF fragments (cube file: "XeAuF%SumFrag%Density.cub") obtained in
the Step 2:

   amspython $PYCUBESCD/PY3/pysub_cube.py -f1 XeAuF%SCF%Density.cub -f2 XeAuF%SumFrag%Density.cub -o  diff_XeAuF.cub

This file ("diff_XeAuF.cub") can be easily visualized within AMS-GUI:

   Open [image: SCMMenu][image] → **View**
   Open the "diff_XeAuF.cub" file (**File** → **Open...**)
   Select **Add** → **Isosurface: With Phase**
   Choose **Select Field...** → **Other...** → **generated**
   Change the isovalue to 0.0005 e/au^3

   [image]

The CD function can be generated by using the "pycd_simple.py" program
available in PYCUBESCD suite (see set up instructions):

   amspython $PYCUBESCD/PY3/pycd_simple.py -f diff_XeAuF.cub

The result is a text format file ("diff_XeAuF.cub_cdz.txt") which
contains the numerical representation of the CD function. It can be
easily visualized with a 2D-plotting software.

   [image]

-[ CD function via AMSview Line graph ]-

   Make sure the XeAuF calculation you just ran is selected in AMSjobs
   Use [image: SCMMenu][image] → **View**
   Use **Fields** → **Grid** → **Fine**
   Use **Fields** → **Calculated**
   Set the left field as **Select Field...** → **Density** → **Density SCF**
   Leave the operation at **-**
   Set the right field as **Select Field...** → **Density** → **Density Sum Frag**
   **Add** → **Isosurface: With Phase**
   Choose the calculated field **C-1** from **Select Field...** → **Other...**
   Change the isovalue to 0.0005 e/au^3
   Make the atom spheres smaller (use **View** → **Molecule** → **Smaller** three times)

   [image]

To get the CD function we need to define a line through the molecule.
For each point along this line all values in a plane through point
perpendicular to the line should be integrated. And finally the values
along the line should be integrated. This functionality is part of the
LineGraph in AMSview:

   Use **Add** → **Graph Along Line**
   At the top of the window next to the (empty) graph window, select **Field** → **Other...** → **C-1**
   Return to the main AMSview window
   Hide the isosurface (using the **Show** check box for the Isosurface: With Phase controls), so we can select atoms
   Select the Xe and F atoms
   Return to the Line Graph window
   In the right-hand panel click the button **Use Two Atoms** to define the line
   Extend the line by 4 Angstrom (click the **Do it** button on the extend line twice)
   Uncheck the show F(r) option for the Line Graph
   Check the show double integral F(r) dSdr  option

In the graph window you will see the result:

   [image]

Note that the grid used is not large enough to provide data for all
points sampled by the planes over which is integrated. You can see
this in the lower left. To fix this, one can extent the grid:

   Use **Fields** → **Grid** → **Details** command, and change the extend value in all directions from 4 to 5

This action has no effect on the graph, as the extra grid points have
no significant contribution to the CD function.

In the Line Graph details window there are more ways to control the
line, which is not very relevant in this case. Another option is to
show the integration planes, which visually gives you an idea of the
size of the planes over which the grid is sampled to calculate the
integrals over the planes.

   Check the **Show (some) Planes** check button in the Details window

   [image]

You can adjust the size of the planes, and the number of points along
the line and the planes using the controls in the Details window.


Case 2: Dewar-Chatt-Duncanson bonding components in a transition metal complex (NOCV-CD)
========================================================================================

This tutorial shows how to generate the NOCV deformation densities
(Δρ_k) and export them on a numerical grid (CUBE file format) directly
using AMS-GUI. The associated CD functions (Δq_k) are evaluated
numerically using the PYCUBESCD suite of programs (see set up
instructions).

The transition state complex needed for this tutorial is given "here".

Note:

  The alignment of the molecule is important. The z-axis defines the
  axis on which is performed the partial integration for evaluating
  the CD function. In this case the integration axis is defined as the
  axis passing for the Ni atom and the midpoint of the CC triple bond
  of the coordinated alkyne.


1. Fragment Calculation
-----------------------

   Start AMSinput and insert the "transition state complex"
   In the Main ADF Panel, select the **Single Point** task
   Select the **GGA:BLYP** XC functional
   Make sure the Relativity is set at **Scalar**
   Select a **TZ2P** Basis set and a **Good** Numerical quality
   In the Panel bar, select **Model** → **Regions**
   Select the C_2 H_2 and click the [image: AddButton][image] to add it as region
   In the menu bar, click **Select** → **Invert Selection** and add them as a second region
   In the Panel bar, select **MultiLevel** → **Fragments**
   Check the **Use fragment** box
   In the Panel bar, select **Properties** → **ETS-NOCV**
   Select the **Closed-Shell** at ETS-NOCV analysis
   Save and run the calculation

   [image]


2. Visualize NOCV deformation densities
---------------------------------------

The NOCV deformation densities can be visualized by performing the
following steps:

   Open [image: SCMMenu][image] → **View** of the Fragment Calculation
   In the menu bar, select **Add** → **Isosurface: With Phase**
   Choose **Select Field...** → **NOCV Def Densities**
   Select a deformation density, for instance Δρ_1
   Change the Iso Value to 0.003 e/au^3

   [image]

In the above figure is shown the deformation density (Δρ_1) associated
with the largest NOCV eigenvalue. It clearly represents the metal-to-
substrate π-backdonation (red/blue colored isodensity surface defines
electron density depletion/accumulation).


3. Generate a CUBE file
-----------------------

One starts with adding all the NOCV deformation densities you want to
map in Cube File format. For instance, in the below example we chose
those associated with the three largest eigenvalues, namely 0.8278,
0.3615 and 0.2194. The isodensity surface in the figure is related to
Δρ_2, which is clearly related to the σ-substrate-to-metal donation
component. You add the NOCV deformation densities with:

   **Add** → **Isosurface: With Phase**
   Choose **Select Field ...** → **NOCV Def Densities**

   [image]

The CUBE file can be generated now:

   Improve the Coarse Grid by selecting **Fields** → **Grid** → **Medium** (or **Fine**)
   In the menu bar, select **File** → **Export Fields as Cube Files**

Note:

  Typically the Medium setting for the Grid  is sufficient to obtain
  accurate CD functions. Clearly, the larger the grid is, the larger
  the size of the generated Cube files will be. For instance using
  Medium Grid the three Cube Files associated at Δρ_1,Δρ_2,Δρ_3 have a
  size of 4.3MB (the latter  increases to 34.3MB whether one sets Fine
  Grid to generate Cube Files).


4. CD functions for the NOCV deformation densities
--------------------------------------------------

The CD functions are generated using the "pycd_simple.py" program
available in  PYCUBESCD suite (see set up instructions):

Give the following command line in the directory where the Cube Files
associated at Δρ_1, Δρ_2, and Δρ_3 have been generated (in the
*.results* folder):

   amspython $PYCUBESCD/PY3/pycd_simple.py   -f  215--Medium-0.cube%NOCV%SumDensities_A_1.cub
   amspython $PYCUBESCD/PY3/pycd_simple.py   -f  218--Medium-0.cube%NOCV%SumDensities_A_2.cub
   amspython $PYCUBESCD/PY3/pycd_simple.py   -f  221--Medium-0.cube%NOCV%SumDensities_A_3.cub

The results are three text format files (*.cub_cdz.txt*) which
contains the numerical representation of the CD functions associated
with Δρ_1, Δρ_2, and Δρ_3. They can be easily visualized with a 2D-
plot software. The picture below can be  obtained combining the 3D
isosurface pictures of the NOCV density deformations.

   [image]


Case 3: Open-Shell CD in the HAT mechanism of the TauD-J intermediate
=====================================================================

This tutorial shows how to generate the *spin-density difference* and
the *total-density difference* (Δρ_α,Δρ_β,Δρ_tot) and export them on a
numerical grid (CUBE file format) using *densf*. The associated CD
functions (Δq_α,Δq_β,Δq_tot) are then evaluated numerically using the
PYCUBESCD suite of programs (see set up instructions).


1. Unrestricted Calculations
----------------------------

First, we need to perform a set of unrestricted ADF calculations on
the transition state structure of the TauD-J__C_2 H_6 model system.
The TauD-J intermediate features a Fe^IV =O unit which is the
H-abstracting species. The two constituting fragments to consider are
the Fe=O moiety of TauD-J and the C_2 H_6 substrate, in the position
they occupy in the transition state geometry. The molecular structures
needed for this tutorial can be downloaded "here".

Note:

  The alignment of the molecule is important. The z-axis defines the
  axis on which is performed the partial integration for evaluating
  the CD function. In this case the integration axis is defined as the
  axis through the Fe=O bond.

Three independent calculations will be performed: (i) TauD-J__C_2 H_6
model system (FeO_CH) (ii) TauD-J fragment (FeO) and (iii) ethane
fragment.

   1. Start AMSinput and insert the "TauD-J__C2H6 system"
   2. In the Main ADF Panel, select the **Single Point** task
   3. Tick the **Unrestricted** calculation box
   4. Fill in for the **Spin polarization** "4.0"
   5. Select the **GGA-D:S12g-D3** XC functional
   6. Make sure the Relativity is set at **Scalar**
   7. Select a **TZ2P** Basis set, a **Small** Frozen core
   8. Select a **Good** Numerical quality
   9. Save and run the calculation

   [image]

Now, the other two single point calculations can be performed:

   Insert the "TauD-J fragment" in a new AMSinput
   Repeat step 2-9 from above

   Insert the "ethane fragment" in a new AMSinput
   Select the **Single Point** task
   Tick the **Unrestricted** calculation box, but set a spin polarization of "0.0"
   Repeat step 5-9 from above

The "adf.rkf" files of these calculations are needed for the following
step.


2. Generate CUBE files
----------------------

We use densf to generate a numerical representation of the electron
spin-densities associated with the unrestricted calculations performed
on the (i) TauD-J__C_2 H_6 model system (FeO_CH) (ii) TauD-J fragment
(FeO) and (iii) ethane fragment (CH). We map these quantities in the
CUBE format. The only information necessary from previous single point
calculations are the "adf.rkf" files.

The following script can be used to generate the CUBE files for the α-
and β-spin electron densities of FeO, CH and FeO_CH.

   for job in "FeO" "CH" "FeO_CH"; do
   $AMSBIN/densf << eor
   InputFile $job.results/adf.rkf
   Density SCF
   CUBOutput $job
   Grid
      -10.0 -10.0 -10.0
      300  300  350
      1.0 0.0 0.0 20.0
      0.0 1.0 0.0 20.0
      0.0 0.0 1.0 20.0
   End
   eor
   done

Note:

  As in any unrestricted calculation, we will get the α- and β-spin
  electron density as a result.


3. Obtain density differences
-----------------------------

We now use the PYCUBESCD suite (see set up instructions) of programs
to make operations between Cube files and to obtain the *spin-density
differences* (Δρ_α, Δρ_β) and the *total density difference* (Δρ).
These density differences will be then integrated along the z-axis to
retrieve the corresponding CD functions.

Use the command lines below in the directory where your just generated
cube files are placed, and when you have named the cuboutput
differently (than FeO, CH and FeO_CH) change the names of the input
cube files as well.

   amspython $PYCUBESCD/PY3/pyadd_cube.py -f1 FeO%SCF%Density_A.cub -f2 CH%SCF%Density_A.cub -o SumFrag_Density_A.cub
   amspython $PYCUBESCD/PY3/pyadd_cube.py -f1 FeO%SCF%Density_B.cub -f2 CH%SCF%Density_B.cub -o SumFrag_Density_B.cub
   amspython $PYCUBESCD/PY3/pysub_cube.py -f1 FeO_CH%SCF%Density_A.cub -f2 SumFrag_Density_A.cub -o diff_FeO_CH_A.cub
   amspython $PYCUBESCD/PY3/pysub_cube.py -f1 FeO_CH%SCF%Density_B.cub -f2 SumFrag_Density_B.cub -o diff_FeO_CH_B.cub
   amspython $PYCUBESCD/PY3/pyadd_cube.py -f1 diff_FeO_CH_A.cub -f2 diff_FeO_CH_B.cub -o diff_FeO_CH.cub

The generated "diff_FeO_CH_A.cub", "diff_FeO_CH_B.cub" and
"diff_FeO_CH.cub" files contain the numerical representation of the
electron (spin) density difference between the FeO_CH system and the
sum of the electron (spin) densities of the FeO and CH fragments.

These files can be easily visualized within AMS-GUI in the same manner
as for Case1 and Case2:

   Open [image: SCMMenu][image] → View
   Open a "diff_FeO_CH.cub" file (**File** → **Open...**)
   **Add** → **Isosurface: With Phase**
   Choose **Select Field...** → **Other...** → **generated**
   Change the Iso Value to 0.01 e/au^3

   [image]


4: Generate the CD functions
----------------------------

The CD functions are finally generated using the pycd_simple.py
program available in PYCUBESCD suite (see set up instructions).

   amspython  $PYCUBESCD/PY3/pycd_simple.py -f diff_FeO_CH_A.cub
   amspython  $PYCUBESCD/PY3/pycd_simple.py -f diff_FeO_CH_B.cub
   amspython  $PYCUBESCD/PY3/pycd_simple.py -f diff_FeO_CH.cub

This results in three text formal files ( *.cub_cdz.txt*) which
contains the numerical representation of the CD functions associated
with Δρ_α, Δρ_β and Δρ_tot. They can be easily visualized with a 2D-
plot software. The picture below can be  obtained combining the
isosurface picture of the total density difference and the three CD
functions.

   [image]


References
==========

[1] L. Belpassi, I. Infante, F. Tarantelli, L. Visscher, *The Chemical
    Bond between Au(I) and the Noble Gases. Comparative Study of NgAuF
    and NgAu+ (Ng = Ar, Kr, Xe) by Density Functional and Coupled
    Cluster Methods*, Journal of the American Chemical Society 130,
    1048 (2008).

[2] G. Bistoni, S. Rampino, F. Tarantelli, L. Belpassi, *Charge-
    displacement analysis via natural orbitals for chemical valence:
    Charge transfer effects in coordination chemistry*, Journal of
    Chemical Physics 142, 084112 (2015).

[3] L. D’Amore, L. Belpassi, J. E. M. N. Klein, M. Swart, *Spin-
    resolved charge displacement analysis as an intuitive tool for the
    evaluation of cPCET and HAT scenarios* Chemical Communications 56,
    12146 (2020).
