Skip to content

System-size convergence in periodic wave function methods

Having introduced the basic workflow, how correlated calculations can be performed using the CC-aims interface and the Cc4s coupled-cluster code for molecules, it is time to move on to periodic systems. All methods implemented in Cc4s work for both molecules and solids. In contrast to molecules, however, an additional source of error needs to be considered for solid-state calculations: The error resulting from using a finite super cell or k-grid, which we will refer to as the finite-size error in this tutorial. Correlated wave function theories, which includes the MP2 and CC methods, tend to converge only slowly with system size, so that one must be careful in correcting that error. As a demonstration of these challenges and possible solutions, we will compute the cohesive energy of the Neon crystal on the coupled-cluster level of theory including single and double excitations (CCSD). Note, that the biggest calculations performed here required over 3000 cores, so that you will need access to a HPC facility to run these yourself.

The FHI-aims input

As is the case for molecules, we first need to choose a basis set that is sufficiently accurate for our purpose. For that purpose, we perform calculations to test the basis set convergence as we did for the paracetamol molecule as well. To compute the cohesive energy of Neon, we need to perform calculations on the Neon atom and the Neon solid.

The Neon atom

The control.in file

We start by writing a control.in file very similiar to the one we used previously for paracetamol. Note, that molecular calculation in the first tutorial was also performed with periodic boundary conditions to make use of the more cost-efficient RI-LVL scheme. To eliminate the error of artifacts, the box size was chosen appropriately big. For the control.in file we write

xc                       hf
frozen_core_postSCF      1
RI_method                V

output cc4s
KS_method                parallel

relativistic             none
charge                   0

override_illconditioning .true.

sc_accuracy_rho          1E-6
sc_accuracy_eev          1E-4
sc_accuracy_etot         1E-6

Note, that the paracetomal calculations of the previous tutorial wew performed with periodic boundary conditions (k_grid 1 1 1) to make use of the more cost-efficient RI-LVL scheme. To eliminate the error of potential artefacts, the box size was chosen appropriately big. For the Neon atom, we will actually perform a purely molecular calculation without any periodic boundary conditions. Since we are only computing one atom, there is no computational overhead in choosing the RI-V scheme (RI_method V). Finally, to obtain a complete control.in file you will need to attach the basis set information below this set of key-value pairs. For the calculations of Neon, we will again use the NAO-VCC-nZ basis sets which you can find in the species_defaults folder of your FHI-aims code. As was the case before, we will manually add the auxiliary basis functions

 for_aux hydro 4 f 1.0
 for_aux hydro 5 g 1.0
 for_aux hydro 6 h 1.0

for each element (in this case only Neon).

The geometry.in file

The geometry.in file for the Neon atom is particularly simple:

atom 0.000000 0.000000 0.000000 Ne

The Neon solid

The control.in file

In comparison to the atomic calculation, the control.in of the Neon crystal contains additional keywords relevant for the periodic boundary conditions.

xc                       hf 
frozen_core_postSCF      1 
RI_method                LVL 
k_grid                   2 2 2
symmetry_reduced_k_grid  .false.
ri_multipole_threshold    1e-7

output cc4s
cc4s_screen_num 210
gamma_cut_coulomb .true.
KS_method                parallel

relativistic             none
charge                   0

override_illconditioning .true.

sc_accuracy_rho          1E-6
sc_accuracy_eev          1E-4
sc_accuracy_etot         1E-6

This examplary input tell FHI-aims to perform a 2x2x2 k-point calculation. To reduce the memory requirements of the Coulomb vertex, we have found that cc4s_screen_num 210 (see the first tutorial for an explanation of the keyword) allows to reduce the size of the Coulomb vertex without notably affecting the resulting CCSD correlation energy. Finally, we have added ri_multipole_threshold 1e-7 since the calculation will otherwise stop during the construction of the auxiliary basis, where a very tight threshold on the multiple moments of the basis functions is used.

The geometry.in file

For the structure of the Neon crystal, we use the FCC structure

lattice_vector 2.320515888050245 2.320515888050245 0.000000000000000
lattice_vector 0.000000000000000 2.320515888050245 2.320515888050245
lattice_vector 2.320515888050245 0.000000000000000 2.320515888050245

atom_frac 0.000000 0.000000 0.000000 Ne

which has also been used in a CC study by Gruber et al..

Running a CC calculation with Cc4s

After performing the FHI-aims and CC-aims calculations, in exactly the same way as was presented in the first tutorial, we now want to perform a CCSD calculation using Cc4s. As a consequence, in comparison to the paracetamol MP2 calculations, we did before, only the Cc4s input file needs to be adapted. Most of the Cc4s input remains identical, we only need to tell Cc4s to perform a CCSD instead of a MP2 calculation. We do so by replacing

- name: SecondOrderPerturbationTheory
  in:
    slicedEigenEnergies: EigenEnergies
    coulombIntegrals: CoulombIntegrals
  out:
    energy: Mp2Energy

from the previous tutorial

by

- name: CoupledCluster
  in:
    method: Ccsd
    integralsSliceSize: 100
    slicedEigenEnergies: EigenEnergies
    coulombIntegrals: CoulombIntegrals
    slicedCoulombVertex: CoulombVertex
    maxIterations: 100
    energyConvergence: 1.0E-4
    mixer:
      type: DiisMixer
  out:
    energy: CcsdEnergy
    amplitudes: Amplitudes

This examplary input tell Cc4s to perform a CCSD calculation. Since the CCSD equations are solved iteratively. While an energy convergence criterion of 1.0E-4 may sound to big, Cc4s only stops when the CC amplitudes (or their residuum) doesn't change by more than 1.0E-6. By that point the energy is generally also sufficiently converged to 1 micro-Hartree or less. However, be always sure to check how well converged your CC calculation is.

Basis set convergence

Now, that we have presented the workflow in its entirety including the relevant input files, we want to determine which basis set to use. In order to allow a fair comparison to the results by Gruber et al., our calculations will be performed using a complete basis set (CBS) extrapolation as introduced at the beginning of the first tutorial. To determine the remaining basis set incompleteness error, we compute the CBS[23] and CBS[45] values for the cohesive energy of the Neon crystal. Since the NAO-VCC-nZ basis set has only been constructed up to n=5, the CBS[45] value will be our reference. We find, that the CCSD correlation contribution to the cohesive energy using the CBS[23] extrapolation is only 8 meV/atom below the CBS[45] value. Note, that at this point we are not including the Hartree-Fock ground-state energy to compute the cohesive energy. Since the Hartree-Fock method is computationally substantially cheaper, we can determine the converged HF contribution to the cohesive energy separately, which we will do at the end.

Convergence to the bulk-limit

As noted in the beginning, one additional error source one needs to account for in periodic calculation is the slow convergence to the bulk-limit. Due to steep scaling (N6) of the CCSD method with respect to system size (N being for example the number of basis functions), only super cells or k-grids of limited size are computationally feasible. Note that at the time of writing, Cc4s is only able to work with super cells. If a k-point based input is provided to Cc4s, as is done in this tutorial, Cc4s interally converts the relevant quantities to their super cell-equivalent. As a consequence, one can not benefit from a speed-up or savings in memory in Cc4s by employing a k-grid instead of a super cell. A block-sparse, k-aware treatment of the CC tensor contractions is currently under development and will allow to access bigger system sizes, making the convergence to the bulk-limit a more managable task.

One way to overcome the slow convergence to the bulk-limit, is by comparing the convergence of the CCSD method with a different, computationally cheaper method, for which the bulk-limit can be reliably determined and which converges with the same rate as the CCSD method. The simplest method, which fulfills these criteria is the random-phase approximation (RPA), which is implemented in FHI-aims. What we will now do to determine the bulk-limit of the CCSD cohesive energy of Neon is

  1. Determine the RPA cohesive energy for a set of k-point grids (2x2x2, 2x2x3, 2x2x4, 2x3x3, 2x3x4, 2x4x4,3x3x3)
  2. Do the same for the CCSD cohesive energy
  3. Determine the bulk-limit of the RPA cohesive energy.
  4. Plot the system size convergence of the cohesive energy of both methods against each other.

To perform an RPA calculation with FHI-aims, you merely need to replace the cc4s-keywords by total_energy_method rpt2. For the control.in file of the Neon crystal that would be

xc                       hf
total_energy_method      rpt2
frozen_core_postSCF      1
RI_method                LVL
k_grid                   2 2 2
symmetry_reduced_k_grid  .false.
ri_multipole_threshold    1e-7

gamma_cut_coulomb .true.
KS_method                parallel

relativistic             none
charge                   0

override_illconditioning .true.

sc_accuracy_rho          1E-6
sc_accuracy_eev          1E-4
sc_accuracy_etot         1E-6

Make sure, that you use the same Coulomb kernel approximation (gamma_cut_coulomb .true.) and the same basis set for both RPA and CCSD. Due to the relatively small computational overhead of the RPA method (in comparison to the CCSD calculations), we can determine the size-converged RPA cohesive energy by brute-force calculation of sufficiently dense k-grids. Indeed, by performing RPA calculation with up to 12x12x12 k-points, we determine a bulk-limit value of -13.16 meV/atom.

To obtain the bulk-limit value of the CCSD method, we now plot the cohesive energies of boths methods against each other yielding the figure below.

Correlation of the RPA and CCSD cohesive energy

First of all, what we see is that all the data points lie on the same line. This demonstrates that both methods converge to the bulk limit with the same rate. It is therefore only reasonable to follow that line to the point where it crosses

\[ E_{\mathrm{coh}}^{RPA} = -13.16\,\mathrm{meV/atom} \]

, which is indicated by the green line. The abscissa value at that intersection (blue line) indicates the bulk-limit value of the CCSD cohesive energy. We find a converged cohesive energy of -18.6 meV/atom. This is in excellent agreement with the -19 meV/atoms, which were obtained by Gruber et al..

The Hartree-Fock contribution

What we should not forget at this point, is the HF contribution to the cohesive energy. Due to the (comparably) negligible cost of the HF method, we can obtain a fully converged HF cohesive energy of the Neon crystal. Indeed, since the Neon crystal is dominated by van-der-Waals correlation interactions, which are not accounted for by the HF method, the HF cohesive energy converges particularly quickly with respect to the system size, so that already with a 3x3x3 k-grid and the NAO-VCC-5Z basis we obtain a value of

\[ E_{\mathrm{coh}}^{HF} = 0.4\,\mathrm{meV/atom.} \]

Beyond the CCSD method : Inclusion of perturbative triple excitations

One more CC method, we haven looked at yet is the CCSD(T) method, which in addition to the CCSD method includes an approximate, perturbative correction from higher-order excitations. In molecular quantum chemistry, the CCSD(T) method is also known as the gold standard of quantum chemistry, since it routinely achieves accuracies of below 1 kcal/mol (43 meV) for thermochemical applications.

To compute the (T)-correction of a CCSD calculation, you need to adapt the present workflow in two ways. First, since the (T)-contribution is based on the CCSD amplitudes,we need to write those to a file. Hence, in addition to the specification of the CCSD method in the Cc4s input file

- name: CoupledCluster
  in:
    method: Ccsd
    integralsSliceSize: 100
    slicedEigenEnergies: EigenEnergies
    coulombIntegrals: CoulombIntegrals
    slicedCoulombVertex: CoulombVertex
    maxIterations: 100
    energyConvergence: 1.0E-4
    mixer:
      type: DiisMixer
  out:
    energy: CcsdEnergy
    amplitudes: Amplitudes

you will need to add a Write-step

- name: CoupledCluster
  in:
    method: Ccsd
    integralsSliceSize: 100
    slicedEigenEnergies: EigenEnergies
    coulombIntegrals: CoulombIntegrals
    slicedCoulombVertex: CoulombVertex
    maxIterations: 100
    energyConvergence: 1.0E-4
    mixer:
      type: DiisMixer
  out:
    energy: CcsdEnergy
    amplitudes: Amplitudes

- name: Write
  in:
    fileName: "Amplitudes.yaml"
    binary: 1
    source: Amplitudes

To perform the (T) calculation it is highly recommended to switch to a different directory. For example, assuming you have performed the CCSD calculation in the CC4S/ccsd subdirectory, create another directory CC4S/pert_T and enter that directory.

Then copy or creat links to the yaml and elements files of the CCSD calculation

ln -s ../ccsd/*yaml .
ln -s ../ccsd/*elements .

To compute the perturbative triples contribution, you will need a Cc4s-input file like this

- name: Read
  in:
    fileName: "EigenEnergies.yaml"
  out:
    destination: EigenEnergies

- name: Read
  in:
    fileName: "CoulombVertex.yaml"
  out:
    destination: CoulombVertex

- name: Read
  in:
    fileName: "Amplitudes.yaml"
  out:
    destination: InitialAmplitudes


- name: DefineHolesAndParticles
  in:
    eigenEnergies: EigenEnergies
  out:
    slicedEigenEnergies: EigenEnergies

- name: SliceOperator
  in:
    slicedEigenEnergies: EigenEnergies
    operator: CoulombVertex
  out:
    slicedOperator: CoulombVertex

- name: VertexCoulombIntegrals
  in:
    slicedCoulombVertex: CoulombVertex
  out:
    coulombIntegrals: CoulombIntegrals

- name: PerturbativeTriples
  in:
    coulombIntegrals: CoulombIntegrals
    amplitudes: InitialAmplitudes
    slicedEigenEnergies: EigenEnergies
  out:
    {}

As you can easily verify, this is essentially the same as the MP2 input file from the previous tutorial, where a Read section for the CC amplitudes was added and the SecondOrderPerturbationTheory has been replaced by the PerturbativeTriples method. Repeating the Neon atom and Neon crystal calculations with CCSD(T) reveals that the (T)-contribution to the cohesive energy is between -10 and -12 meV/atom and largely independent of the system size. Hence, our final CCSD(T) cohesive energy is found to be around 30 meV/atom, which again agrees very well with the results by Gruber et al..