User Tools

Site Tools


exercises:2016_uzh_cmest:geometry_optimization

Electronic structure calculation using DFT

In this exercise, you will perform geometry optimization using DFT.

1. Step: Single point energy calculation with separate coordinate file

In the previous exercises we initially specified all parameters – pseudopotential and basis set coefficients as well as atom coordinates – in the input file. Later we used the pseudpotentials and basis sets from a file provided by CP2K.

Now we go further and also factor out the atomic structure to make it easier to automate different calculations for the same structure or the same calculation for different structures. The format used for this is the same you will get for trajectories for example.

First create two files with (different) coordinates for Ethane C2H8 (do not confuse with Ethene C2H6 from before):

ethane1.xyz
       8

C     0.750           0.000           0.000
C    -0.750           0.000           0.000
H    -1.050           0.000          -0.850
H    -1.050           0.736           0.425
H    -1.050          -0.736           0.425
H     1.050           0.000          -0.850
H     1.050           0.736           0.425
H     1.050          -0.736           0.425
ethane2.xyz
       8

C     0.750           0.000           0.000
C    -0.750           0.000           0.000
H    -1.050           0.000           0.850
H    -1.050           0.736          -0.425
H    -1.050          -0.736          -0.425
H     1.050           0.000          -0.850
H     1.050           0.736           0.425
H     1.050          -0.736           0.425

The input file looks almost the same as the one for the previous calculations:

ethane.inp
&GLOBAL
  PROJECT ethane
  RUN_TYPE ENERGY
  PRINT_LEVEL MEDIUM
&END GLOBAL

&FORCE_EVAL
  METHOD Quickstep              ! Electronic structure method (DFT,...)
  &DFT
    BASIS_SET_FILE_NAME  BASIS_MOLOPT
    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                        ! Parameters needed to compute the electronic exchange potential 
      &XC_FUNCTIONAL PBE
      &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
      COORD_FILE_FORMAT xyz
      COORD_FILE_NAME  ./ethane1.xyz
    &END
    &KIND H
      ELEMENT H
      BASIS_SET DZVP-MOLOPT-GTH
      POTENTIAL GTH-PBE-q1
    &END KIND
    &KIND C
      ELEMENT C
      BASIS_SET DZVP-MOLOPT-GTH
      POTENTIAL GTH-PBE-q4
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

and should give the following energy once you run it:

ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):              -14.746153797151329

You can also directly open a XYZ file in VMD:

$ vmd ethane1.xyz

2. Step: Optimizing the geometry

The only thing you have to change to get a geometry optimization instead of a single point energy calculation is the following:

ethane1_opt.inp
&GLOBAL
  PROJECT ethane1_opt
  RUN_TYPE GEO_OPT
  PRINT_LEVEL MEDIUM
&END GLOBAL
[...]

Note the different RUN_TYPE and the changed PROJECT name. The latter is not strictly necessary but recommended, since CP2K automatically creates additional files using this project name as a prefix.

After running this, you will have the following files:

$ ls ethane1_opt*
ethane1_opt-1.restart        ethane1_opt-1.restart.bak-3  ethane1_opt.out          ethane1_opt-RESTART.wfn.bak-1
ethane1_opt-1.restart.bak-1  ethane1_opt-BFGS.Hessian     ethane1_opt-pos-1.xyz    ethane1_opt-RESTART.wfn.bak-2
ethane1_opt-1.restart.bak-2  ethane1_opt.inp              ethane1_opt-RESTART.wfn  ethane1_opt-RESTART.wfn.bak-3

Take a look at the output file, especially the following section (repeated the number of cycles it took to reach convergence):

 --------  Informations at step =     1 ------------
  Optimization Method        =                 BFGS
  Total Energy               =       -14.9417142787
  Real energy change         =        -0.1955604816
  Predicted change in energy =        -0.1885432833
  Scaling factor             =         0.0000000000
  Step size                  =         0.2677976891
  Trust radius               =         0.4724315332
  Decrease in energy         =                  YES
  Used time                  =               19.018

  Convergence check :
  Max. step size             =         0.2677976891
  Conv. limit for step size  =         0.0030000000
  Convergence in step size   =                   NO
  RMS step size              =         0.1458070233
  Conv. limit for RMS step   =         0.0015000000
  Convergence in RMS step    =                   NO
  Max. gradient              =         0.0287243359
  Conv. limit for gradients  =         0.0004500000
  Conv. for gradients        =                   NO
  RMS gradient               =         0.0180771987
  Conv. limit for RMS grad.  =         0.0003000000
  Conv. for gradients        =                   NO
 ---------------------------------------------------

For each convergence criterion you see the value which is used to check whether convergence is reached and convergence is only reached if all of them are satisfied simultaneously.

From the output file, extract the following data and generate 3 plots with the values vs the iteration number:

  • Total FORCE_EVAL ( QS ) energy
  • Max. gradient .. the maximal force (seen over all atoms)
  • Max. step size .. the maximal displacement (seen over all atoms)

3. Step: Optimizing the geometry with an alternative geometry

Now change the used coordinate file to ethane2.xyz and update the PROJECT name to not overwrite the results from the previous simulation and run it again.

  • Compare the final energy reached for both structures and the total number of optimization steps required
  • Open the two output XYZ files (<PROJECT>-pos-1.xyz) in VMD and compare them. Is there a difference? What is it?
  • Which one is likely to be more stable and why?

4. Step: Visualize the geometry optimization

Append the following section to your input file (does not matter for which structure) and run the simulation again.

&MOTION
  &PRINT
    &TRAJECTORY
      LOG_PRINT_KEY T
      FORMAT XYZ
      ADD_LAST NUMERIC
    &END TRAJECTORY
  &END PRINT
&END MOTION

If you check the output XYZ file now (<PROJECT>-pos-1.xyz) and compare it to the input structure, you will notice that the same molecule is now specified multiple times as so-called frames. Checking the CP2K output file you will also notice the following entry:

[...]
 Writing TRAJECTORY 1_1 to ethane2_opt_traj-pos-1.xyz
[...]

Open this new XYZ file again with VMD, choose an appropriate drawing method (Licorice or CPK), hit the play button and enjoy:

exercises/2016_uzh_cmest/geometry_optimization.txt · Last modified: 2016/10/20 17:53 by tmueller