Differences

This shows you the differences between two versions of the page.

--- exercises:2019_conexs_newcastle:ex1 [2019/09/10 12:17] – [Part 2: Running an single-point calculation] abussy
+++ exercises:2019_conexs_newcastle:ex1 [2020/08/21 10:15] (current) – external edit 127.0.0.1
@@ Line 114: / Line 114: @@
 &END COORD
 </code>
-sets up our simulation cell and atomic coordinates respectively. In the &KIND sections,
+sets up our simulation cell and atomic coordinates respectively. By default, CP2K will run your calculations with periodic boundary conditions but since we are using a large simulation box, the molecule does not interact with its own copies and we are effectively looking at it in gas phase. In the &KIND sections,
 <code>
 &KIND H
@@ Line 137: / Line 137: @@
 or
 <code>
-./cp2k.sopt water.inp | tee water.out
+sbatch cp2k.sh
 </code>
-if you want to follow the output as it runs. When the program is finished, you will have two new files in your working directory:
+if you are logged in to the cluster and have modified the //cp2k.sh// accordingly. When the program is finished, you will have two new files in your working directory:
 <code>
 water.out
@@ Line 216: / Line 216: @@
 </code>
 in the &SCF section.
+=====Part 3: Geometry optimization=====
+In our previous calculation, our guess for the geometry of the water molecule, e.g. with a bond angle of 90°, was not very good. It is very often a good idea to optimize the geometry of your system before extracting any information from your calculations. Make another directory and move a copy of our previous input into it. To perform a geometry optimization, we must add an additional section to our input file.
+<code>
+&MOTION
+  &GEO_OPT
+    OPTIMIZER BFGS
+    MAX_DR    3.0E-03
+    MAX_FORCE 4.5E-04
+    RMS_DR    1.5E-03
+    RMS_FORCE 3.0E-04
+  &END
+  &CONSTRAINT
+    &FIXED_ATOMS
+      COMPONENTS_TO_FIX XYZ
+      LIST 1
+    &END
+  &END CONSTRAINT
+&END MOTION
+</code>
+In the &GEO_OPT section we select which kind of optimizer we want as well as the convergence criterion, such as the maximum force component. Then, since we only want to optimize the hydrogen relative to the oxygen, the &CONSTRAINT section allows us to freeze all components of the first atom in our input coordinates (make sure this is the oxygen by examining the &SUBSYS section again). Finally, we must also change the RUN_TYPE to GEO_OPT in the &GLOBAL section.
+<code>
+&GLOBAL
+  PROJECT h2o_opt
+  RUN_TYPE GEO_OPT
+  IOLEVEL LOW
+&END GLOBAL
+</code>
+We are now ready to start the calculation. As before, type e.g.
+<code>
+./cp2k.sopt -i water_opt.inp -o water_opt.out &
+</code>
+or
+<code>
+./cp2k.sopt water.inp | tee water.out
+</code>
+into your terminal to start the calculation. A geometry optimization contains more SCF cycles and so this calculation will take more time than before. Once the system has converged, we can examine the output file to evaluate our calculation. At each step in the optimization, we get the same information as in the last exercise, along with
+<code>
+ --------  Informations at step =     2 ------------
+  Optimization Method        =                 BFGS
+  Total Energy               =       -17.1947218061
+  Real energy change         =        -0.0002999093
+  Predicted change in energy =        -0.0002336556
+  Scaling factor             =         0.0000000000
+  Step size                  =         0.0178822087
+  Trust radius               =         0.4724315332
+  Decrease in energy         =                  YES
+  Used time                  =               22.246
+  Convergence check :
+  Max. step size             =         0.0178822087
+  Conv. limit for step size  =         0.0030000000
+  Convergence in step size   =                   NO
+  RMS step size              =         0.0091540732
+  Conv. limit for RMS step   =         0.0015000000
+  Convergence in RMS step    =                   NO
+  Max. gradient              =         0.0026713842
+  Conv. limit for gradients  =         0.0004500000
+  Conv. for gradients        =                   NO
+  RMS gradient               =         0.0017198827
+  Conv. limit for RMS grad.  =         0.0003000000
+  Conv. for gradients        =                   NO
+ ---------------------------------------------------
+</code>
+describing each step. It is very important to check that your geometry optimization has converged properly, and this can be done by looking at the end of the output file. The last step should look something like this:
+<code>
+ --------  Informations at step =     4 ------------
+  Optimization Method        =                 BFGS
+  Total Energy               =       -17.1947556096
+  Real energy change         =        -0.0000002510
+  Predicted change in energy =        -0.0000002260
+  Scaling factor             =         0.0000000000
+  Step size                  =         0.0009840890
+  Trust radius               =         0.4724315332
+  Decrease in energy         =                  YES
+  Used time                  =               20.295
+  Convergence check :
+  Max. step size             =         0.0009840890
+  Conv. limit for step size  =         0.0030000000
+  Convergence in step size   =                  YES
+  RMS step size              =         0.0006410487
+  Conv. limit for RMS step   =         0.0015000000
+  Convergence in RMS step    =                  YES
+  Max. gradient              =         0.0000198391
+  Conv. limit for gradients  =         0.0004500000
+  Conv. in gradients         =                  YES
+  RMS gradient               =         0.0000096694
+  Conv. limit for RMS grad.  =         0.0003000000
+  Conv. in RMS gradients     =                  YES
+ ---------------------------------------------------
+ *******************************************************************************
+ ***                    GEOMETRY OPTIMIZATION COMPLETED                      ***
+ *******************************************************************************
+</code>
+Make sure that the calculation has converged by all criteria. You should also have several more files in your working directory, including e.g.
+    *h2o_opt-1.restart
+    *h2o_opt-1.restart.bak-1
+    *h2o_opt-pos-1.xyz
+The h2o_opt-1.restart contains all the input required to restart a calculation from the last step in the calculation. Its h2o_opt-1.restart.bak-1 analogue is the backup containing information from the previous step, and so on. The last file h2o_opt-pos-1.xyz contains the atomic trajectory of the optimization.
+Try to open it in e.g. Molden(?) to visualize the steps.
+Here we can see that the bond angle is now 102.8°, which is better but still quite far from the accepted 104.45°. Also comparing the final total energy to our first single-point calculation, we can see that our new geometry is more energetically stable.
+=====Part 4: Computing PDOS=====
+Now that we have obtained a more reasonable geometry, we will know attempt to compute the eigenvalues of our Kohn-Sham orbitals and compare the energy levels to the attached experimental spectra. We will also use this opportunity to try out the option to restart a calculation.
+-0.316280
+To restart our calculation from the previous geometry optimization, we need to copy two files from the old working directory, e.g.
+    h2o_opt-RESTART.wfn
+    h2o_opt-1.restart
+to our current working directory. To make the code read these files, we have to add a section
+<code>
+&EXT_RESTART
+  RESTART_FILE_NAME h2o_opt-1.restart
+&END EXT_RESTART
+</code>
+as well as adding
+<code>
+SCF_GUESS RESTART
+RESTART_FILE_NAME h2o_opt-RESTART.wfn
+</code>
+to the &SCF and &DFT sections respectively. To print the PDOS, we also need to add the &PDOS section under &DFT
+<code>
+&PRINT
+  &PDOS
+     FILENAME ./pdos
+     NLUMO 2
+  &END
+&END
+</code>
+and specify the spectra file names for example if we want to include unoccupied virtual orbitals. The final input can e.g. look like this:
+<code>
+&DFT
+  BASIS_SET_FILE_NAME GTH_BASIS_SETS
+  POTENTIAL_FILE_NAME GTH_POTENTIALS
+  LSD 0
+  RESTART_FILE_NAME h2o_opt-RESTART.wfn
+  &QS
+    METHOD GPW
+    EPS_DEFAULT 1.0E-10
+  &END QS
+  &SCF
+    MAX_SCF    300
+    EPS_SCF    1.0E-06
+    SCF_GUESS RESTART
+    &MIXING
+      METHOD DIRECT_P_MIXING
+      ALPHA   0.6
+    &END MIXING
+    &DIAGONALIZATION
+      ALGORITHM STANDARD
+    &END DIAGONALIZATION
+  &END SCF
+  &MGRID
+    NGRIDS 4
+    CUTOFF 300
+    REL_CUTOFF 60
+  &END
+  &XC
+    &XC_FUNCTIONAL BLYP
+    &END XC_FUNCTIONAL
+  &END XC
+  &PRINT
+    &PDOS
+      FILENAME ./h2o
+      NLUMO 2
+    &END
+  &END
+&END DFT
+</code>
+When the calculation is done, open the file
+    *h2o-k1-1.pdos
+and compare your calculated energies to the experimental measurements reported in the attached paper, e.g. Fig. 1. As a reminder, one 1 a.u. = 27.212 eV. Look in particular on the energy differences between the three valence states.
+{{:exercises:2019_conexs_newcastle:h2o_pes.pdf | Water PES}}
+----
+=====Part 5: (Optional) Testing your parameters=====
+If you followed this exercise closely, chances are that the computed orbital energies in your PDOS files are quite similar to the experimental ones, given an ad hoc shift of the energies. However, the 1b2 orbital will be slightly off. Can we do better by tweaking our input, or can we get as good results with a faster method? In this exercise, we will investigate and try to get a feel for how different parameters affect the accuracy and time consumption of your calculation. Go back to Part 3 and Part 4 and redo the calculation while changing the following parameters
+<code>
+&XC_FUNCTIONAL  {PADE BLYP B3LYP}
+CUTOFF          {100 300 600}
+BASIS_SET       {SZV-GTH DZVP-GTH TZVP-GTH}
+EPS_SCF         {1.0E-04 1.0E-06 1.0E-0.8}
+</code>
+(choose one of the suggested in the brackets). The suggestions are given from less accurate to more accurate. When you feel comfortable, feel free to experiment with other keywords!   Note that using e.g. a larger basis set will make the calculation take longer, so you might want to do something else in the mean time.