CP2K Open Source Molecular Dynamics

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
======= 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:

<code - 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
</code>

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

  * the ''HF'' section is missing: this is obvious since we are doing a pure DFT calculation now
  * the parameter for the ''XC_FUNCTIONAL'' section is set to ''PADE'' (which is a synonym für ''LDA'') instead of ''NONE'', meaning that we are going to do a DFT calculation using the Local Density Approximation
  * the specifications for the ''BASIS'' and the ''POTENTIAL'' for the different atom ''KIND''s has changed

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:

<code>
$ cd $CP2K_DATA_DIR
$ less BASIS_SET
$ less POTENTIAL
</code>

In the basis sets you will find entries like:

<code>
[...]
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
[...]
</code>

while the pseudopotentials file contains something like:

<code>
H GTH-PADE-q1 GTH-LDA-q1 GTH-PADE GTH-LDA
    1
     0.20000000    2    -4.18023680     0.72507482
    0
</code>

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:

  * Perdew-Burke-Ernzerhof (**PBE**)
  * Becke-Lee-Yang-Parr (**BLYP**)

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:

  - Change the parameter for the ''XC_FUNCTIONAL'' section to either PBE or BLYP
  - Update the ''POTENTIAL'' to ''GTH-PBE-q?'' or ''GTH-BLYP-q?'' for each ''KIND''
  - Set the ''BASIS_SET'' to ''DZVP-MOLOPT-GTH'' for each ''KIND'' (the same for PBE and BLYP)
  - 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//:

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

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):

  * SZV-MOLOPT-GTH
  * DZVP-MOLOPT-GTH
  * TZVP-MOLOPT-GTH
  * TZV2P-MOLOPT-GTH
  * TZV2PX-MOLOPT-GTH

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

Create two plots:

  * total energy vs number of independent orbital functions
  * time in seconds vs number of independent orbital functions