The equation-of-motion coupled-cluster method
We want to use this tutorial to show you how you can perform calculations of excited states using the FHI-aims/CC-aims/Cc4s infrastructure. In 1993, Stanton and Bartlett developed the equation-of-motion coupled-cluster method (EOM-CC), which constitutes an extension of ground-state coupled-cluster (CC) theory to address excited states. Depending on the type of excitation one is interested in, there are different EOM-CC methods, which includes EE-EOM-CC for neutral electronic excitations, IP-EOM-CC for ionization processes and EA-EOM-CC for electron attachment processes. In this tutorial we will focus on IP- and EA-EOM-CC theory.
The EOM-CCSD implementation in Cc4s
The development version of Cc4s does offer the calculation of charged excited states on the EOM-CC level of theory with single and double excitations (IP-EOM-CCSD and EA-EOM-CCSD). These methods are publicly available on the open-source repository of Cc4s,
General EOM-CCSD workflow
As you will see, most of the workflow to perform Cc4s-calculations with FHI-aims does not change for excited states. Indeed, the FHI-aims and CC-aims part of the calculation remains exactly the same, since the EOM-CCSD method makes use of the same quantities as the ground-state CC methods (e.g MP2, CCSD) do.
What does change is the Cc4s part: First of all, as you can read in the original paper the starting-point of a EOM-CCSD calculation is a ground-state CCSD calculation. As a consequence, our new workflow will involve one more step on the Cc4s side. Overall, we will perform the following calculations
- FHI-aims Hartree-Fock calculation involving call to CC-aims library (same as before)
- CCSD calculation using Cc4s (same as in the second tutorial)
- IP- and/or EA-EOM-CCSD calculation
In this tutorial we will show how to perform the third step. Depending on your application, you may want to perform either an IP-EOM-CCSD or an EA-EOM-CCSD calculation or both. In molecules, both the ionization potential (IP) and the electron affinity (EA) are well-defined quantities, and you can use either of the two EOM-CCSD methods to obtain those. For molecular applications, the EOM-CCSD method is considered to be one of the most accurate available methods for excited states and can routinely achieve accuracies of ~100 meV. For periodic solids, the IP and EA is only well-defined relative to a given shape and surface termination of the Born-von-Karman cell, so that a comparison to experiment would require exact knowledge of the experimental setup. However, the difference EA(k)-IP(k) for a given momentum vector k in reciprocal space is a well-defined quantity in periodic systems as well and yields the band gap of the material at that k-point. Note, that due to the sign convention used in Cc4s, the band gap in this tutorial will be calculated by the sum EA+IP.
Currently, there are two limitations on performing EOM-CC band gap calculations using Cc4s with FHI-aims. The first one has been mentioned briefly in the second tutorial and is due to the fact that Cc4s does only perform super cell calculations. Hence, it is not possible for the EOM-CCSD algorithm to pick a specific k-point at which the band gap should be computed. Instead, we will need to provide a super cell with an appropriate shift in reciprocal space. Note, that we can generate Cc4s input via FHI-aims using either a super cell or a k-grid. However, the Cc4s calculation will internally always convert the input to its super cell equivalent. Due to this constraint, we currently can only access the minimal band gap for a given super cell, while entire band structures would require a k-point aware treatment of the EOM-CC equations in Cc4s, which is currently under development. The second limitation stems from FHI-aims: At the time of writing, shifting of a k-grid for a Hartree-Fock calculation is not supported by FHI-aims. However, we have developed a workaround: Instead of directly shifting the k-grid or the super cell in reciprocal space, we can perform a Hartree-Fock calculation on a very dense k-grid and subsequently select a shifted subset of that grid to perform the Cc4s calculations.
In this tutorial, we will perform both calculations which can be performed on a regular desktop computer or laptop (IP of the fluorine molecule, and the 1x1x2 super cell of trans-polyacetylene, tPA) and those for which access to a HPC infrastructure is highly advisable (IP of uracil).
First ionization potential of F2
The first two steps involving the FHI-aims calculation and the CCSD ground-state calculation are identical to the
previously shown examples.
As in the previous cases, start by creating two directories aims and CC4S
mkdir aims
mkdir CC4S
and enter the aims directory
cd aims
where we will create the control.in and geometry.in input files for
the FHI-aims calculation.
In the control.in file we specify
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
To determine the ionization potential of F2, we use the NAO-VCC-3Z basis sets
which you can find in the species_defaults directory of FHI-aims.
Complete the control.in file by concatenating the computational parameters specified above
with the species_defaults basis set files for F.
cat control.in /path/to/FHI-aims/species_defaults/NAO-VCC-nZ/NAO-VCC-3Z/06_C_default > control_F.in
mv control_F.in control.in
As in the previous examples, add auxiliary basis functions
for_aux hydro 4 f 1.0
for_aux hydro 5 g 1.0
for_aux hydro 6 h 1.0
to the basis set specification of every element (here: only F).
Your final control.in file should look like this
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
################################################################################
#
# FHI-aims code project
# Igor Ying Zhang and Xinguo Ren, Fritz Haber Institute Berlin, 2012
#
# Suggested "NAO-VCC-3Z" defaults for F atom (to be pasted into control.in file)
#
# Please cite:
#
# Igor Ying Zhang, Xinguo Ren, Patrick Rinke, Volker Blum, and Matthias Scheffler,
# "Numeric atom-centered-orbital basis sets with valence-correlation consistency from H to Ar"
# New Journal of Physics 15, 123033 (2013).
# http://dx.doi.org/10.1088/1367-2630/15/12/123033
#
################################################################################
species F
# set basis set reference for FHI-aims output
cite_reference NAO-VCC-2013
# global species definitions
nucleus 9
mass 18.9984032
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 0e-4
#
radial_base 37 7.0
radial_multiplier 2
angular_grids specified
division 0.4014 110
division 0.5291 194
division 0.6019 302
division 0.6814 434
division 0.7989 590
# division 0.8965 770
# division 1.3427 974
# outer_grid 974
outer_grid 770
# outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 2 s 2.
valence 2 p 5.
# ion occupancy
ion_occ 2 s 1.
ion_occ 2 p 4.
################################################################################
basis_acc 1e-04
#============================================
# Optimization in F for NAO-VCC-3Z
#============================================
for_aux hydro 4 f 1.0
for_aux hydro 5 g 1.0
for_aux hydro 6 h 1.0
# (sp) correlation set
hydro 1 s 2.79571831
hydro 1 s 3.59883615
hydro 2 p 2.79844642
hydro 2 p 7.10137981
# polarization set
hydro 3 d 9.19074920
hydro 3 d 11.25968007
hydro 4 f 16.05919112
# (sp) enhanced minimal set
hydro 1 s 9.62761974
hydro 1 s 11.10594321
hydro 2 p 7.75321023
For the geometry, we use the bond length reported in this paper, resulting in
the following geometry.in file
atom 0.0000000000000000 0.0000000000000000 0.0000000000000000 F
atom 0.0000000000000000 0.0000000000000000 1.4119300000000000 F
With the input files for the FHI-aims/CC-aims calculation complete, run it for example via
mpirun -np 4 /path/to/FHI-aims-executable > out.aims
which should take a minute at most. Make sure that the calculation finished successfully by
checking that "Have a nice day" is printed at the end of the output file and by
finding some new files most importantly EigenEnergies.elements and CoulombVertex.elements
in your calculation directory.
This concludes the first step of the 3-step workflow outlined above.
It's time for the second step, the ground-state CCSD calculation. Leave the aims directory and
create a new sub-directory ccsd in the CC4S directory
cd ../CC4S
mkdir ccsd
cd ccsd
To perform the CCSD calculation you will need to either copy or link the previously generated input
.yaml and .elements files to this directory.
ln -s ../../aims/*elements .
ln -s ../..//aims/*yaml .
Finally you need to provide a yaml file, which communicates to Cc4s what kind of calculation we
are interested in. In the present case, we want to run a CCSD calculation, which in the end
saves the amplitudes. The following ccsd.yaml file
- name: Read
in:
fileName: "EigenEnergies.yaml"
out:
destination: EigenEnergies
- name: Read
in:
fileName: "CoulombVertex.yaml"
out:
destination: CoulombVertex
- 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: 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
will do exactly that. Run the CCSD calculation via
mpirun -np 4 /path/to/cc4s/executable -i ccsd.yaml > out.ccsd
In contrast to the previous tutorials, we are not really interested in the CCSD correlation energy
but in the wave function, which has been saved in the files Amplitudes.components.ph.elements
and Amplitudes.components.pphh.elements.
To compute the ionization potential of F2, we create one final calculations directory
inside CC4S/and enter it
cd ..
mkdir -p eom/IP
cd eom/IP
Again, we will need the .yaml and .elements files
ln -s ../../ccsd/*yaml .
ln -s ../../ccsd/*elements .
and a Cc4s input file, eom.yaml, to perform the IP-EOM-CCSD calculation:
- 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: Amplitudes
- 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: EquationOfMotion
in:
initialAmplitudes: Amplitudes
slicedEigenEnergies: EigenEnergies
coulombIntegrals: CoulombIntegrals
slicedCoulombVertex: CoulombVertex
numberExcitations: 1
initialGroundstateCalculation: ccsd
evaluateStructureFactor: 0
coulombPotential: CoulombPotential
eigensolver:
type: Davidson
maxIterations: 20
maxBasisSize: 50
evaluator:
type: IP
groundstate: ccsd
out:
rightAmplitudes: RightAmplitudes
As you can see, this input is identical to the CCSD(T) input from the previous tutorial,
except for the fact that we replaced the PerturbativeTriples algorithm at the end
by the EquationOfMotion method. This method has a couple of new keywords, in particular
numberExcitations: Specifies how many excited states you want to compute.initialGroundstateCalculation: Specifies our starting point for the EOM-CC calculation. Currently only the most widely usedccsdstarting point is implemented.evaluateStructureFactor: Specifies if you want to compute the correlation EOM-CC structure factor after the calculation is completed. This is a potentially useful quantity for periodic calculations and gives insights into short- and long-range electronic correlations. This keyword can be either switched on (1) or off (0). For this tutorial, we won't need to compute the EOM-CC structure factor, hence we will always setevaluateStructureFactor: 0.eigensolver: The EOM-CCSD excitation energies are determined by partially diagonalizing the similarity transformed Hamiltonian, which generally is a of very impractical size to store. As a consequence, we need to use an indirect eigensolver, which does not construct the entire Hamiltonian at any point. The most widely used algorithm for this purpose is theDavidsonalgorithm, which is also the one specified in the input file.maxIterations: Number of iterations to be performed in the Davidson algorithmmaxBasisSize: With every iteration the Davidson algorithm computes new basis vectors, which span the subspace in which the eigenvectors of the Hamiltonian are searched. To avoid memory problems, you can specify a maximum number of basis vectors. Once this maximum is reached, the Davidson algorithm only keeps the currently best guess of the solution eigen vectors and deletes the rest. The remaining vectors then serve as the starting point to construct the new subspace.
evaluator: The instance, which defines what type of EOM-CC calculation is to be performed. It consists of two parts:type: What kind of excitations are to be calculated? One can currently choose betweenIPandEAgroundstate: This is not identical to theinitialGroundstateCalculationkeyword. Here, we specify to which order of correlation we want to perform the EOM-CC calculation. The combination ofgroundstate: ccsdwithtype: IPis equivalent to saying that we will perform a IP-EOM-CCSD calculation. Currently, onlygroundstate: ccsdis supported.
Run the EOM-CCSD calculation via
mpirun -np 4 /path/to/cc4s/executable -i eom.yaml > out.ip
This calculation will determine the IPs by performing the Davidson algorithm which iteratively determines the lowest lying
eigenvalue and the corresponding eigenvector of the similiarity-transformed Hamiltonian. After each completed iteration,
the intermediate result along with the convergence of the energy and the excitation wave function will be printed to out.ip.
In the case of our F2 calculation you should get something very similar to this:
grep -A1 "Davidson iteration" out.ip
Davidson iteration 1 completed:
0.75131426 7.5131e-01 5.1900e-01
--
Davidson iteration 2 completed:
0.57670638 1.7461e-01 1.0742e-01
--
Davidson iteration 3 completed:
0.57559085 1.1155e-03 1.5648e-02
--
Davidson iteration 4 completed:
0.57563070 3.9854e-05 4.6677e-03
--
Davidson iteration 5 completed:
0.57556909 6.1612e-05 7.8613e-04
--
Davidson iteration 6 completed:
0.57556327 5.8213e-06 2.8441e-04
--
Davidson iteration 7 completed:
0.57556137 1.9006e-06 7.9887e-05
--
Davidson iteration 8 completed:
0.57556128 9.0691e-08 1.3592e-05
--
Davidson iteration 9 completed:
0.57556173 4.4898e-07 2.0038e-06
--
Davidson iteration 10 completed:
0.57556180 6.9361e-08 3.7864e-07
where the first value is the IP value in Hartree, the second value is the change of the IP compared to the previous iteration
and the third value quantifies the change of the wave function compared to the previous iteration. We find that we obtain
an IP of 0.57556180 Hartree = 15.66 eV, which is in good agreement with the aug-cc-pVTZ result from the previously mentioned
paper
of 15.61 eV.
Electronic band gap of trans-polyacetylene
To perform a periodic EOM-CCSD calculation all you have to do is adapt the FHI-aims input to perform a periodic calculation. Since the Cc4s code is agnostic to the fact if the system is a molecule or under periodic boundary conditions, the CCSD and the EOM-CCSD calculation inputs will remain unchanged.
We again start in a new, empty directory, which we call aims.
The computational settings we set in the control.in file look slightly different this time:
xc hf
frozen_core_postSCF 1
RI_method LVL
k_grid 1 1 1
symmetry_reduced_k_grid .false.
output cc4s
cc4s_screen_thresh 1.0d-3
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
If you have done the previous tutorial on periodic ground-state CC calculations, these should look familiar.
Since the present calculation should be performed under periodic boundary conditions, we need to specify
the k-grid (k_grid 1 1 1) and that we don#t want to exploit the symmetry of the the Brillouin zone
(symmetry_reduced_k_grid .false.). Moreover, if we employ PBC in FHI-aims we need to use the more
memory-efficient RI-LVL scheme (RI_method LVL) instead of the RI-V one. The last change is optional
but notably decreases the computation time of the subsequent CC steps and involves making the
Coulomb vertex more compact via a principal component analysis (cc4s_screen_thresh 1.0d-3). Details on this
keyword can be found in the previous tutorials.
The second half of the control.in file is comprised of the basis set specification, for which we will use the
loc-NAO-VCC-2Z, which can be found in the species_defaults directory of FHI-aims. Don't forget to manually
add the three for_aux auxiliary basis functions for each element as we did before. Note, that
you will first have to remove the instances of the for_aux hydro 4 f 3.0 line first.
The complete control.in file looks like this
xc hf
frozen_core_postSCF 1
RI_method LVL
k_grid 1 1 1
symmetry_reduced_k_grid .false.
output cc4s
cc4s_screen_thresh 1.0d-3
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
################################################################################
#
# FHI-aims code project
#
# Igor Ying Zhang, Fritz Haber Institute Berlin, 2014
#
# Suggested "loc-NAO-VCC-nZ" defaults for C atom (to be pasted into control.in file)
#
################################################################################
species C
# global species definitions
nucleus 6
mass 12.0107
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 0e-4
#
radial_base 34 7.0
radial_multiplier 2
angular_grids specified
division 0.2187 50
division 0.4416 110
division 0.6335 194
division 0.7727 302
division 0.8772 434
division 0.9334 590
# division 0.9924 770
# division 1.0230 974
# division 1.5020 1202
# outer_grid 974
outer_grid 770
# outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 2 s 2.
valence 2 p 2.
# ion occupancy
ion_occ 2 s 1.
ion_occ 2 p 1.
################################################################################
basis_acc 1e-05
#================================================
# Optimization in C for loc-NAO-VCC-2Z
# RI-V prodbas generation parameter in Opt.:
# 1) "prodbas_threshold 1.e-5"
# 2) "default_prodbas_acc 1.e-4"
# RPA total energy : -xxxx.xxxx eV
#================================================
for_aux hydro 4 f 1.0
for_aux hydro 5 g 1.0
for_aux hydro 6 h 1.0
# (sp) correlation set
hydro 1 s 0.63593881
hydro 2 p 3.98166864
# polarization set
hydro 3 d 6.01255805
################################################################################
#
# FHI-aims code project
#
# Igor Ying Zhang, Fritz Haber Institute Berlin, 2014
#
# Suggested "loc-NAO-VCC-nZ" defaults for H atom (to be pasted into control.in file)
#
################################################################################
species H
# global species definitions
nucleus 1
mass 1.00794
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 0e-4
#
radial_base 24 7.0
radial_multiplier 2
angular_grids specified
division 0.1930 50
division 0.3175 110
division 0.4293 194
division 0.5066 302
division 0.5626 434
division 0.5922 590
# division 0.6227 974
# division 0.6868 1202
outer_grid 770
# outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 1 s 1.
# ion occupancy
ion_occ 1 s 0.5
################################################################################
basis_acc 1e-05
#================================================
# Optimization in H for loc-NAO-VCC-2Z
# RI-V prodbas generation parameter in Opt.:
# 1) "prodbas_threshold 1.e-5"
# 2) "default_prodbas_acc 1.e-4"
# RPA total energy : -xxxx.xxxx eV
#================================================
for_aux hydro 4 f 1.0
for_aux hydro 5 g 1.0
for_aux hydro 6 h 1.0
# (sp) correlation set
hydro 1 s 1.23411101
hydro 2 p 2.94900728
# polarization set
To determine the band gap of trans-polyacetylene (tPA) in the bulk-limit, one needs to perform
a series of calculations with increasing chain length. To keep the example simple, we will perform the band gap calculation
for the 1x1x2 super cell, i.e 2 repeating monomer units. We will use the structural parameters of one tPA geometry discussed
in the paper by Windom et al.. In that study the convergence
of the tPA band gap with respect to the chain length was investigated by performing molecular EOM-CCSD calculations.
We will perform periodic calculations using the structural paramaters of the tPA3 molecule, as it is referred to in that
paper.
The geometry.in file for the 1x1x2 super cell of the tPA3 structure contains
lattice_vector 4.898485 0.000000 0.000000
lattice_vector 0.000000 90.000000 0.000000
lattice_vector 0.000000 0.000000 90.000000
atom 1.854683 5.339047 0.000000 C
atom 3.033575 4.660953 0.000000 C
atom 1.862432 6.429020 0.000000 H
atom 3.025826 3.570980 0.000000 H
atom 4.303925 5.339047 0.000000 C
atom 0.584333 4.660953 0.000000 C
atom 4.311674 6.429020 0.000000 H
atom 0.576584 3.570980 0.000000 H
where the lattice vectors in the directions perpendicular to the chain were chosen to be sufficiently long (90A) to avoid any interactions between chains.
Run the FHI-aims/CC-aims calculation again via
mpirun -np 4 /path/to/FHI-aims-executable > out.aims
and check if the calculation completed without problems. The calculation should take less than 10 min.
Using the same ccsd.yaml file as in the F2, you can now leave the aims directory, create a CC4S directory, inside
of which you should make a ccsd sub-directory. Copy your previous ccsd.yaml file to this directory and link or copy
the *.elements and *.yaml files from your FHI-aims calculation there.
Then launch the CCSD calculation.
cd ..
mkdir -p CC4S/ccsd
cd CC4S/ccsd
cp /path/to/F2-calculation/CC4S/ccsd/ccsd.yaml .
ln -s ../../aims/*elements .
ln -s ../../aims/*yaml .
mpirun -np 4 /path/to/cc4s/executable -i ccsd.yaml > out.ccsd
Finally, we want to obtain the IP and EA of tPA, so that this time we need to perform
two calculations using the saved amplitudes of the CCSD run.
For that purpose, leave the ccsd directory and create to more directories
cd ..
mkdir -p eom/IP
mkdir -p eom/EA
While you can use the same eom.yaml file for the IP calculation, for the EA calculation
the only thing you need to do is to replace type: IP by type: EA, so that
the eom.yaml file should look like
- 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: Amplitudes
- 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: EquationOfMotion
in:
initialAmplitudes: Amplitudes
slicedEigenEnergies: EigenEnergies
coulombIntegrals: CoulombIntegrals
slicedCoulombVertex: CoulombVertex
numberExcitations: 1
initialGroundstateCalculation: ccsd
evaluateStructureFactor: 0
coulombPotential: CoulombPotential
eigensolver:
type: Davidson
maxIterations: 20
maxBasisSize: 50
evaluator:
type: EA
groundstate: ccsd
out:
rightAmplitudes: RightAmplitudes
After running the IP
cd eom/IP
cp /path/to/F2-calculation/CC4S/eom/IP/eom.yaml .
ln -s ../../ccsd/*elements .
ln -s ../../ccsd/*yaml .
mpirun -np 4 /path/to/cc4s/executable -i eom.yaml > out.ip
and the EA calculation
cd ../EA
cp /path/to/F2-calculation/CC4S/eom/IP/eom.yaml .
sed -i 's/IP/EA/g' eom.yaml
ln -s ../../ccsd/*elements .
ln -s ../../ccsd/*yaml .
mpirun -np 4 /path/to/cc4s/executable -i eom.yaml > out.ea
you should obtain an IP value of 0.25341246 Ha and an EA value of -0.00063671 Ha. The band gap, IP+EA then equates to 0.25277575 Ha = 6.88 eV. This is over 1 eV more than the band gap value of 5.07 eV in the paper. This is, however, not surprising since we only performed the calculation for the 1x1x2 super cell with a small loc-NAO-VCC-2Z basis set. To obtain a value comparable to the paper, one needs to perform calculations with increasing super cell size and a bigger basis set (loc-NAO-VCC-3Z or more). Note, that the band gap of tPA is at the (0.0, 0.0, 0.5) point in the Brillouin zone. Therefore, one should use super cells of shape 1x1xn where n is even. By performing EOM-CCSD calculations with n up to 16 and extrapolating appropriately to the bulk-limit, one obtains a value of 4.91 eV as shown in Moerman et al.. This paper also contains a thorough discussion how EOM-CCSD band gaps should be extrapolated to the bulk-limit.
First ionization potential of uracil
To compute the ionization potential of uracil, you can mostly reuse the input files from the F2
calculation. You only need to adapt the appended basis set information
in the control.in file to include the elements of the uracil molcule (C,H,O,N) and, naturally, the
geometry.in file should contain the structure of the uracile molecule.
Hence the FHI-aims input in your aims directory should look like this.
geometry.in:
lattice_vector 200.00000000 0.00000000 0.00000000
lattice_vector 0.00000000 200.00000000 0.00000000
lattice_vector 0.00000000 0.00000000 200.00000000
atom 1.17147361362677 0.98429541867005 0.00000252776412 N
atom 0.00418000814178 1.70325956333876 0.00002379122680 C
atom 0.12200942793589 2.77883208151210 0.00003862819434 H
atom -1.20010382983883 1.10592023715235 0.00002550709315 C
atom -2.11732948576409 1.67253822411493 0.00004253934579 H
atom -1.28705241926252 -0.34683196232309 0.00000331557903 C
atom -2.30808895933210 -1.00329830035103 0.00000680904104 O
atom -0.03294018824560 -0.98675259692590 -0.00001521543628 N
atom -0.04578026288494 -1.99822647280132 -0.00002719195527 H
atom 1.21790842532829 -0.40444914969017 -0.00001370630409 C
atom 2.25771571106274 -1.02564176155144 -0.00002744944291 O
atom 2.06814395923261 1.44368171885476 0.00000744489429 H
control.in:
xc hf
frozen_core_postSCF 1
RI_method LVL
k_grid 1 1 1
symmetry_reduced_k_grid .false.
output cc4s
cc4s_screen_thresh 1.0d-3
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
################################################################################
#
# FHI-aims code project
# Igor Ying Zhang and Xinguo Ren, Fritz Haber Institute Berlin, 2012
#
# Suggested "NAO-VCC-3Z" defaults for C atom (to be pasted into control.in file)
#
# Please cite:
#
# Igor Ying Zhang, Xinguo Ren, Patrick Rinke, Volker Blum, and Matthias Scheffler,
# "Numeric atom-centered-orbital basis sets with valence-correlation consistency from H to Ar"
# New Journal of Physics 15, 123033 (2013).
# http://dx.doi.org/10.1088/1367-2630/15/12/123033
#
################################################################################
species C
# set basis set reference for FHI-aims output
cite_reference NAO-VCC-2013
# global species definitions
nucleus 6
mass 12.0107
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 0e-4
#
radial_base 34 7.0
radial_multiplier 2
angular_grids specified
division 0.2187 50
division 0.4416 110
division 0.6335 194
division 0.7727 302
division 0.8772 434
division 0.9334 590
# division 0.9924 770
# division 1.0230 974
# division 1.5020 1202
# outer_grid 974
outer_grid 770
# outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 2 s 2.
valence 2 p 2.
# ion occupancy
ion_occ 2 s 1.
ion_occ 2 p 1.
################################################################################
basis_acc 1e-04
#============================================
# Optimization in C for NAO-VCC-3Z
#============================================
for_aux hydro 4 f 1.0
for_aux hydro 5 g 1.0
for_aux hydro 6 h 1.0
# (sp) correlation set
hydro 1 s 1.80690099
hydro 1 s 3.25966754
hydro 2 p 1.97404286
hydro 2 p 3.16616778
# polarization set
hydro 3 d 5.89063938
hydro 3 d 6.67716137
hydro 4 f 10.14967498
# (sp) enhanced minimal set
hydro 1 s 6.77306992
hydro 1 s 13.85913422
hydro 2 p 3.64746093
################################################################################
#
# FHI-aims code project
# Igor Ying Zhang and Xinguo Ren, Fritz Haber Institute Berlin, 2012
#
# Suggested "NAO-VCC-3Z" defaults for H atom (to be pasted into control.in file)
#
# Please cite:
#
# Igor Ying Zhang, Xinguo Ren, Patrick Rinke, Volker Blum, and Matthias Scheffler,
# "Numeric atom-centered-orbital basis sets with valence-correlation consistency from H to Ar"
# New Journal of Physics 15, 123033 (2013).
# http://dx.doi.org/10.1088/1367-2630/15/12/123033
#
################################################################################
species H
# set basis set reference for FHI-aims output
cite_reference NAO-VCC-2013
# global species definitions
nucleus 1
mass 1.00794
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 0e-4
#
radial_base 24 7.0
radial_multiplier 2
angular_grids specified
division 0.1930 50
division 0.3175 110
division 0.4293 194
division 0.5066 302
division 0.5626 434
division 0.5922 590
# division 0.6227 974
# division 0.6868 1202
outer_grid 770
# outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 1 s 1.
# ion occupancy
ion_occ 1 s 0.5
################################################################################
basis_acc 1e-04
#============================================
# Optimization in H for NAO-VCC-3Z
#============================================
for_aux hydro 4 f 1.0
for_aux hydro 5 g 1.0
for_aux hydro 6 h 1.0
# (sp) correlation set
hydro 1 s 1.21957223
hydro 1 s 1.79015338
hydro 2 p 3.43862914
hydro 2 p 3.48684723
# polarization set
hydro 3 d 6.20828695
# (sp) enhanced minimal set
hydro 1 s 2.62768292
################################################################################
#
# FHI-aims code project
# Igor Ying Zhang and Xinguo Ren, Fritz Haber Institute Berlin, 2012
#
# Suggested "NAO-VCC-3Z" defaults for O atom (to be pasted into control.in file)
#
# Please cite:
#
# Igor Ying Zhang, Xinguo Ren, Patrick Rinke, Volker Blum, and Matthias Scheffler,
# "Numeric atom-centered-orbital basis sets with valence-correlation consistency from H to Ar"
# New Journal of Physics 15, 123033 (2013).
# http://dx.doi.org/10.1088/1367-2630/15/12/123033
#
################################################################################
species O
# set basis set reference for FHI-aims output
cite_reference NAO-VCC-2013
# global species definitions
nucleus 8
mass 15.9994
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 0e-4
#
radial_base 36 7.0
radial_multiplier 2
angular_grids specified
division 0.1817 50
division 0.3417 110
division 0.4949 194
division 0.6251 302
division 0.8014 434
division 0.8507 590
# division 0.8762 770
# division 0.9023 974
# division 1.2339 1202
# outer_grid 974
outer_grid 770
# outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 2 s 2.
valence 2 p 4.
# ion occupancy
ion_occ 2 s 1.
ion_occ 2 p 3.
################################################################################
basis_acc 1e-04
#============================================
# Optimization in O for NAO-VCC-3Z
#============================================
for_aux hydro 4 f 1.0
for_aux hydro 5 g 1.0
for_aux hydro 6 h 1.0
# (sp) correlation set
hydro 1 s 1.84877109
hydro 1 s 3.09577821
hydro 2 p 2.44895734
hydro 2 p 5.85722578
# polarization set
hydro 3 d 8.53701653
hydro 3 d 8.96275978
hydro 4 f 14.09619628
# (sp) enhanced minimal basis
hydro 1 s 3.62165559
hydro 1 s 5.97272514
hydro 2 p 6.29404180
################################################################################
#
# FHI-aims code project
# Igor Ying Zhang and Xinguo Ren, Fritz Haber Institute Berlin, 2012
#
# Suggested "NAO-VCC-3Z" defaults for N atom (to be pasted into control.in file)
#
# Please cite:
#
# Igor Ying Zhang, Xinguo Ren, Patrick Rinke, Volker Blum, and Matthias Scheffler,
# "Numeric atom-centered-orbital basis sets with valence-correlation consistency from H to Ar"
# New Journal of Physics 15, 123033 (2013).
# http://dx.doi.org/10.1088/1367-2630/15/12/123033
#
################################################################################
species N
# set basis set reference for FHI-aims output
cite_reference NAO-VCC-2013
# global species definitions
nucleus 7
mass 14.0067
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 0e-4
#
radial_base 35 7.0
radial_multiplier 2
angular_grids specified
division 0.1841 50
division 0.3514 110
division 0.5126 194
division 0.6292 302
division 0.6939 434
division 0.7396 590
# division 0.7632 770
# division 0.8122 974
# division 1.1604 1202
# outer_grid 974
outer_grid 770
# outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 2 s 2.
valence 2 p 3.
# ion occupancy
ion_occ 2 s 1.
ion_occ 2 p 2.
################################################################################
basis_acc 1e-04
#============================================
# Optimization in N for NAO-VCC-3Z
#============================================
for_aux hydro 4 f 1.0
for_aux hydro 5 g 1.0
for_aux hydro 6 h 1.0
# (sp) correlation set
hydro 1 s 1.20973156
hydro 1 s 2.56620771
hydro 2 p 2.07057135
hydro 2 p 4.62007052
# polarization set
hydro 3 d 7.34425305
hydro 3 d 7.73693275
hydro 4 f 12.37741826
# (sp) enhanced minimal set
hydro 1 s 3.47271250
hydro 1 s 4.95564021
hydro 2 p 4.71939479
The geometry that we use was taken from the IP-EOM-CCSD benchmark paper by Tripathi and Dutta.
To make use of the more efficient RI-LVL scheme, we perform the calculation in a cubic box of 200A edge length.
Alternatively, you can remove the lattice_vector lines and specify the RI-V scheme in the control.in file
(as was the case in the F2). This will,
however most likely, be more computationally expensive without affecting the end result significantly.
Having run the FHI-aims/CC-aims calculation, the remaining steps including the necessary input files are identical to the previous example of the fluorine molecule.
If you run the final IP-EOM-CCSD calculation, you should obtain a first ionization potential of 0.3540 Ha = 9.63 eV, which is in good agreement with the value of 9.54 eV obtained by Tripathi and Dutta. The remaining deviation can be attributed to the basis set incompleteness error, which is the consequence of using only a 3Z basis (cc-pVTZ basis in the paper).