# Description of VASP interface to Cc4s

## Table of Contents

This section explains the `VASP`

workflow to calculate occupied and unoccupied one-electron orbitals, as well as all necessary input files for `Cc4s`

.
Basic knowledge with VASP is assumed.
In addition, we refer to the Wiki page of VASP.

An example `bash`

script performing all following steps for a periodic BN crystal is available here here.

**Step 1**: Hartree-Fock: Occupied orbitals

In this and all following steps we assume that the required `VASP`

input files (`INCAR`

, `POSCAR`

, `POTCAR`

and `KPOINTS`

) are present in the same working directory.
We assume that the `KPOINTS`

file samples one \(k\) -point of the first Brillouin zone only.
We initialize the Hartree-Fock (HF) step by a preceding density functional theory (DFT) ground state calculation using the following `INCAR`

file.
The individual flags will be explained below.

ENCUT = $encut NCORE = $ncore ISMEAR = 0 SIGMA = 0.001 ISYM = -1

Although not essential, the DFT preparation usually increases the stability of the subsequent HF iteration procedure in `VASP`

.
The produced `WAVECAR`

file is then read in the HF step for which we use the following `INCAR`

file.

ENCUT = $encut NCORE = $ncore ISMEAR = 0 SIGMA = 0.001 ISYM = -1 ALGO = A EDIFF = 1E-7 LHFCALC = T AEXX = 1.0

The tag `ENCUT`

defines the plane-wave basis set size by a kinetic energy cutoff in eV.
It is advisable to set `ENCUT`

to a value which is at least 20% larger than the largest `ENMAX`

entry in the `POTCAR`

file.
For optimal performance on multi-core machines set `NCORE`

to the number of cores per socket.
The `ISMEAR`

and `SIGMA`

ensure integer occupancies of the one-electron orbitals and can safely be used for all insulating systems with a finite band gap.
We recommend to use the conjugated gradient algorithm for the iterative HF procedure by setting `ALGO=A`

.
For properly converged Hartree–Fock orbital energies, we choose a break condition of \(10^{-7}\,\text{eV}\) by setting `EDIFF=1E-7`

.
Further information on the used tags can be found in the VASP wiki.

**Step 2**: Hartree-Fock: Unoccupied orbitals

Coupled cluster theory approximates the many-electron wave function using excited Slater determinants contructed from occupied and unoccupied orbitals.
We compute all unoccupied HF orbitals in the given plane-wave basis as described below.
The total number of HF orbitals/bands is equal to the number of plane-waves and can be found in the `OUTCAR`

file of the HF ground state calculation in the following way.

nbands=`awk '/number of plane-waves:/ {print $5*2 -1} ' < OUTCAR # OUTCAR file from HF ground state calculation using gamma-only code #nbands=`awk '/number of plane-waves:/ {print $5} ' < OUTCAR # OUTCAR file from HF ground state calculation using complex code

Please note, that `VASP`

automatically adjusts the number of bands to be a multiple of the number of MPI ranks in `VASP`

calculations for reasons of parallelization.
Note that `nbands`

should be determined in a different manner when using the complex version of `VASP`

.

Ensure that the `WAVECAR`

file contains all converged occupied HF orbitals as described above and use the following `INCAR`

file to calculate the unoccupied HF orbitals:

ENCUT = $encut NCORE = $ncore ISMEAR = 0 SIGMA = 0.001 ISYM = -1 ALGO = Exact LHFCALC = .TRUE. AEXX = 1.0 NBANDS = $nbands NELM = 1

Here, an exact diagonalisation (`ALGO=Exact`

) of the Fock matrix with dimension `nbands`

\(\times\) `nbands`

is performed.
This is a single-shot step, i.e. no iterations are necessary: `NELM=1`

.
Save the `WAVECAR`

file containing all unoccupied (virtual) orbitals, e.g. `mv WAVECAR WAVECAR.HFV`

.

**Optional Step 2-(a)**: MP2 for focal-point basis-set correction

For the BasisSetCorrection algorithm the MP2 pair energies matrix
is needed. This input file can be obtained from VASP using automated extrapolation techniques for MP2 pair energies.
Make sure to employ the `WAVECAR`

file from the previous calculation (`WAVECAR.HFV`

) that contains all orbitals.
In dependence of the system size, two VASP algorithms are available.
For small systems (e.g. less than 50 occupied orbitals with `NBANDS`

< 20 * occupied orbitals) the following `INCAR`

can be used

ENCUT = $encut ISMEAR = 0 SIGMA = 0.001 ISYM = -1 ALGO = MP2 LHFCALC = .TRUE. AEXX = 1.0 NBANDS = $nbands LSFACTOR = .TRUE.

For larger systems, the low-scaling algorithm (Schäfer, Ramberger, and Kresse 2017) might be faster and less memory consuming.

ENCUT = $encut ISMEAR = 0 SIGMA = 0.001 ISYM = -1 LMP2LT = .TRUE. ALGO = ACFDTRK NOMEGA = 6 LDUMPMP2IJ = .TRUE. LHFCALC = .TRUE. AEXX = 1.0 NBANDS = $nbands

This algorithm is based on a Laplace transformed (LT) MP2 formulation using `NOMEGA = 6`

sampling points for the Laplace integration. This setting is sufficient for large gap systems. For small gap systems try with `NOMEGA = 8`

or even higher.
For an efficient parallelization, use the `KPAR`

flag and set it to a divisor of the number of mpi-ranks (optimally half of the number of mpi-ranks).
In case of memory issues, reduce `KPAR`

.

NOTE: The basis-set extrapolation procedures between the different MP2 algorithms in VASP differ and for consistency reasons one should not compare extrapolated correlation energy estimates between different algorithms.

NOTE: This low-scaling algorithm heavily relies on FFTs. For efficiency reasons, please use `PRECFOCK = Fast`

consistently in all post-HF or post-DFT VASP runs.
For further details, please check the LTMP2 - Tutorial in the VASP wiki.

At the end of this `VASP`

calculation the following input files needed by `Cc4s`

are written to disk.

- Mp2PairEnergies (\(\epsilon_{ij}\)) :
`Mp2PairEnergies.yaml`

,`Mp2PairEnergies.elements`

**Optional Step 2-(b)**: Approximate natural orbitals

The convergence of the electronic correlation energy is very slow when using canonical Hartree–Fock orbitals.
Approximate natural orbitals allow for achieving a more rapid correlation energy convergence to the complete basis set limit.
Here, we compute these natural orbitals as described below and using Eq.2 from Ref.(Grüneis et al. 2011).
However, alternative approaches such as `ALGO=RPANO`

exist in `VASP`

.
Make sure to employ the `WAVECAR`

file containing all orbitals (`WAVECAR.HFV`

).

ENCUT = $encut ISMEAR = 0 SIGMA = 0.001 ISYM = -1 ALGO = MP2NO LHFCALC = .TRUE. AEXX = 1.0 NBANDS = $nbands LAPPROX = .TRUE.

This `VASP`

calculation will produce the `WAVECAR.FNO`

file containing all unoccupied natural orbitals.
We choose to work with a small subset of these natural orbitals. We recommend to use 10 unoccupied natural orbitals per occupied orbital in combination with the
BasisSetCorrection algorithm of `Cc4s`

. The corresponding number of orbitals can be obtained
using the following command.

nbands=`awk <OUTCAR "/NELEC/ { print $3/2 * 11 }"`

All `Cc4s`

algorithms are currently based on canonical formulations.
To this end we need to re-canonicalize the subset of natural orbitals by performing another `VASP`

calculation.
Make sure to use the `WAVECAR.FNO`

file as input, `cp WAVECAR.FNO WAVECAR`

and employ the following `INCAR`

file.

ENCUT = $encut NCORE = $ncore ISMEAR = 0 SIGMA = 0.001 ISYM = -1 ALGO = sub LHFCALC = .TRUE. AEXX = 1.0 NBANDS = $nbands NBANDSHIGH = $nbands NELM = 1

The `NBANDSHIGH`

tag makes sure that exactly the number of orbitals specified by `NBANDS`

will be used regardless of the number of MPI ranks.
Save the `WAVECAR`

file containing the re-canonicalized unoccupied natural orbitals, e.g. `mv WAVECAR WAVECAR.CNO`

.

**Step 3**: Computing `Cc4s`

input files

In the final step we call `VASP`

using the `WAVECAR`

file with the desired choice of unoccupied orbitals (e.g. `cp WAVECAR.HFV WAVECAR`

or `cp WAVECAR.CNO WAVECAR`

) and
the following `INCAR`

file.

ENCUT = $encut NCORE = $ncore ISMEAR = 0 SIGMA = 0.001 ISYM = -1 ALGO = CC4S EDIFF = 1E-5 NBANDS = $nbands NBANDSHIGH = $nbands ENCUTGW = $encutgw ENCUTGWSOFT = $encutgw ISYM = -1

This step produces the following input files needed by `Cc4s`

- Eigenenergies (\(\epsilon_{p}\)) :
`EigenEnergies.yaml`

,`EigenEnergies.elements`

- CoulombVertex (\(\Gamma^{pG}_{q}\)) :
`CoulombVertex.yaml`

,`CoulombVertex.elements`

- GridVectors (\(\vec G\)) :
`GridVectors.yaml`

,`GridVectors.elements`

- CoulombPotential (\(V(\vec G)\)) :
`CoulombPotential.yaml`

,`CoulombPotential.elements`

- DeltaIntegrals (\(\delta^{ab}_{ij}\)) :
`DeltaPPHH.yaml`

,`DeltaPPHH.elements`

- DeltaIntegrals (\(\delta_{ij}\)) :
`DeltaHH.yaml`

,`DeltaHH.elements`

- CoulombVertexSingularVetors (\(U_{G}^{F}\)) :
`CoulombVertexSingularVectors.yaml`

,`CoulombVertexSingularVectors.elements`

## Literature

*Journal of Chemical Theory and Computation*7 (9): 2780–85. https://doi.org/10.1021/ct200263g.

*The Journal of Chemical Physics*146 (10): 104101. doi:10.1063/1.4976937.