Description of VASP interface to Cc4s

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.

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.

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

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.

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