Electronic structure calculation using DFT

In this exercise, you will perform again an electronic structure calculation (of Ethene), but this time using Density Functional Theory and different functionals.

1. Step: Running a DFT calculation

Create a new directory for this exercise and create an input input file using the following content:

ethene_LDA.inp 
&GLOBAL
  PROJECT ethene
  RUN_TYPE ENERGY
  PRINT_LEVEL MEDIUM
&END GLOBAL

&FORCE_EVAL
  METHOD Quickstep              ! Electronic structure method (DFT,...)
  &DFT
    BASIS_SET_FILE_NAME  BASIS_SET
    POTENTIAL_FILE_NAME  POTENTIAL

    &POISSON                    ! Solver requested for non periodic calculations
      PERIODIC NONE
      PSOLVER  WAVELET          ! Type of solver
    &END POISSON
    &SCF                        ! Parameters controlling the convergence of the scf. This section should not be changed. 
      SCF_GUESS ATOMIC
      EPS_SCF 1.0E-6
      MAX_SCF 300
    &END SCF
    &XC                        ! Parametes needed to compute the electronic exchange potential 
      &XC_FUNCTIONAL PADE
      &END XC_FUNCTIONAL
    &END XC
  &END DFT

  &SUBSYS
    &CELL
      ABC 10 10 10
      PERIODIC NONE              ! Non periodic calculations. That's why the POISSON section is needed 
    &END CELL
    &TOPOLOGY                    ! Section used to center the atomic coordinates in the given box. Useful for big molecules
      &CENTER_COORDINATES
      &END
    &END
    &COORD
    C         -2.15324        3.98235        0.00126
    C         -0.83403        4.16252       -0.00140
    H         -0.25355        3.95641        0.89185
    H         -0.33362        4.51626       -0.89682
    H         -2.65364        3.62861        0.89669
    H         -2.73371        4.18846       -0.89198
    &END COORD
    &KIND H
      ELEMENT H
      BASIS_SET DZVP-GTH-PADE
      POTENTIAL GTH-PADE-q1
    &END KIND
    &KIND C
      ELEMENT C
      BASIS_SET DZVP-GTH-PADE
      POTENTIAL GTH-PADE-q4
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

Comparing this input file to the one from the previous exercise, we notice a couple of things:

As you have seen in the lecture, one has to select a basis set for doing calculations efficiently. Furthermore we approximate the core electrons of an atom by a common pseudopotential instead of calculating them explicitly, reducing the computational complexity even further.

CP2K comes with a number of files, specifying the respective coefficients. Which one of those files is going to be used to lookup the basis sets/pseudopotentials is defined using the BASIS_SET_FILE_NAME and POTENTIAL_FILE_NAME inside the &DFT section. In the &KIND section then one has only to specify an identifier for an entry in the respective file.

The files are located in $CP2K_DATA_DIR. Change to this directory and take a look at the file: 

$ cd $CP2K_DATA_DIR
$ less BASIS_SET
$ less POTENTIAL

In the basis sets you will find entries like: 

[...]
H DZVP-GTH-PADE
  2
  1  0  0  4  2
        8.3744350009  -0.0238943732   0.0000000000
        1.8058681460  -0.1397943259   0.0000000000
        0.4852531032  -0.2530970874   0.0000000000
        0.1658235797  -0.6955307423   1.0000000000
  2  1  1  1  1
        0.7000000000   1.0000000000
[...]

while the pseudopotentials file contains something like: 

H GTH-PADE-q1 GTH-LDA-q1 GTH-PADE GTH-LDA
    1
     0.20000000    2    -4.18023680     0.72507482
    0

Now return to the previous (exercise) directory ($ cd - brings you back to the previous directory) and run the simulation. Compare the energy calculated using DFT-LDA ($ grep 'ENERGY|' ethene_LDA.out) to the one calculated in the last exercise using Hartree-Fock. What do you observe? What does that mean?

2. Step: Selecting different functionals

LDA is fast, but seldom very accurate. This is why there are improved functionals and one class are the Generalized Gradient Approximation (GGA) functionals.

Two prominent implementations of GGA functionals are:

When changing the functional, one usually also has to change the pseudopotential to one optimized for the usage with this functional, while the right basis set is selected based on the selected pseudopotential.

Change the input file given above for LDA to run the same calculation once using PBE and once using BLYP:

  1. Change the parameter for the XC_FUNCTIONAL section to either PBE or BLYP
  2. Update the POTENTIAL to GTH-PBE-q? or GTH-BLYP-q? for each KIND
  3. Set the BASIS_SET to DZVP-MOLOPT-GTH for each KIND (the same for PBE and BLYP)
  4. Since the so-called MOLOPT basis sets are in a separate file, you also have to set the value for the BASIS_SET_FILE_NAME to BASIS_MOLOPT

Using the following command you can measure the total time for a simulation in minutes:seconds: 

$ TIME=%E time cp2k.sopt -i ethene_LDA.inp -o ethene_LDA.out
0:09.34

Measure the time and the energy for LDA, PBE and BLYP. How do they compare?

3. Step: Convergence with regard to the basis set

As has been said in the lectures, the methods we are employing here are in principal exact. But to make calculations feasible we have to truncate the number of basis functions. This is why when you make calculations which should give accurate numbers you will have to prove that you have employed a basis set large enough such that the result does not change anymore when going to a larger basis.

Update the input file from the previous for PBE to run the simulation and collect the timings once for the following basis sets (for both kinds):

For each output, look for the Number of independent orbital functions entry (besides the energy) in the output file.

Create two plots: