Geometry optimisation converges to nonsense

[hfruchtl]

Folks,

I am trying to do my first baby steps with CP2K, after several aborted attempts, so please tell me if I do something very stupid.

One thing that worries me after first tests is that geometry optimisations look converged, when they clearly aren't. In the example below, I start something that should be a water dimer at a very wrong geometry. When I plot the energy changes, they become smaller as expected, and after less than 30 steps the program finishes. But the geometry in the xyz file looks nothing like I expect (H2O)2 to look. When I optimise this structure with the program that must not be named, using the same functional, I get a proper dimer. If I then plug this into CP2K, it remains close to this geometry (good!) with a much lower energy than the first attempt (also good). I understand that BFGS may not be good from a bad guess, but I would expect CG to get there in the end, or at least not look like it converged. I'd be very surprised if this was a local minimum. In short: how do I know if I can trust an optimised geometry?

Input below. Thanks in advance,

  Herbert

&GLOBAL
  PROJECT water-dimer
  RUN_TYPE GEO_OPT
  PRINT_LEVEL high
&END GLOBAL

&MOTION
  &GEO_OPT
    MAX_ITER 2000
    MINIMIZER CG
  &END GEO_OPT
  &PRINT
    &TRAJECTORY
      LOG_PRINT_KEY T
      FORMAT xyz
      UNIT angstrom
      &EACH
        GEO_OPT 1
      &END EACH
      ADD_LAST NUMERIC
    &END TRAJECTORY
  &END PRINT
&END MOTION

&FORCE_EVAL
  METHOD Quickstep
  &DFT
    BASIS_SET_FILE_NAME /usr/local/CP2K/BASIS_MOLOPT
    POTENTIAL_FILE_NAME /usr/local/CP2K/GTH_POTENTIALS
    &MGRID
      CUTOFF 100
    &END MGRID
    &QS
    &END QS
    &SCF
      SCF_GUESS ATOMIC
    &END SCF
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
    &END XC
  &END DFT

  &SUBSYS
    &CELL
      A    15.0 0.0 0.00
      B    0.0 15.0 0.0
      C    0.0 0.0 15.0
      PERIODIC NONE
    &END CELL
    &COORD
      O         0.0006535865        0.0348037292        6.5225572751
      H        -0.0003839912       -0.0241840602        7.5017980482
      H        -0.0019131905        1.0517623171        6.4318839506
      O         0.0244126401        0.0102938023        4.0045222106
      H        -0.0164725514       -0.0290451673        5.4795600635
      H         1.0326311856        0.2923983937        4.3857121035
    &END COORD

    &KIND O
      BASIS_SET DZVP-MOLOPT-GTH-q6
      POTENTIAL GTH-PBE-q6
    &END KIND
    &KIND H
      BASIS_SET DZVP-MOLOPT-GTH-q1
      POTENTIAL GTH-PBE-q1
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

[Matt W]

Hi Herbert,

excitingly there is actually a tutorial available now, that I think addresses you problem

http://www.cp2k.org/howto:converging_cutoff

Please give feedback! 

I will say that the settings in the current version are probably "military grid" standard. If energies are converged to something sensible like microhartrees, I think your issues will be solved.

Matt

[Jörg Saßmannshausen]

Hi Herbert,

I would guess that your convergence criteria are set too loose and hence it
'converges' but the results are nonsense.

You can check that from your output file.
Just search for:
"Informations at step"

and have a look at the printed values here.
You should see something like that:
Optimization Method = BFGS
Total Energy = -160.5982073285
Real energy change = -0.0000000071
Predicted change in energy = -0.0000000070
Scaling factor = 0.8736246364
Step size = 0.0001454725
Trust radius = 0.1000000000
Decrease in energy = YES
Used time = 39.292

Convergence check :
Max. step size = 0.0000532434
Conv. limit for step size = 0.0003000000
Convergence in step size = YES
RMS step size = 0.0000171249
Conv. limit for RMS step = 0.0001500000
Convergence in RMS step = YES
Max. gradient = 0.0000112942
Conv. limit for gradients = 0.0000450000
Conv. in gradients = YES
RMS gradient = 0.0000036405
Conv. limit for RMS grad. = 0.0000300000
Conv. in RMS gradients = YES

Some parameters which might help here are:
In the force_eval section I got these settings:
&DFT
CHARGE ${CHARGE}
&MGRID
CUTOFF ${CUTOFF}
NGRIDS 5
&END MGRID
&QS
EPS_DEFAULT 1.0E-16
&END QS
&SCF
SCF_GUESS random
EPS_SCF ${SCF_CONV}
MAX_SCF 30
&PRINT
[ ... ]

with
SCF_CONV 1.0E-8

I can send you the whole file off-list if you want to.

Let u know how you are getting on here.

All the best from dull London!

Jörg

--
*************************************************************
Jörg Saßmannshausen
University College London
Department of Chemistry
Gordon Street
London
WC1H 0AJ

email: j.sassma...@ucl.ac.uk
web: http://sassy.formativ.net

Please avoid sending me Word or PowerPoint attachments.
See http://www.gnu.org/philosophy/no-word-attachments.html

[Florian Schiffmann]

Hi Herbert,

this is clearly a cutoff problem (100 is way too low for whatever you want to do). This breaks translational invariance with respect to the grid and artificial minima arise simply from the atomic positions relative to the grid which explains the strange "minimal energy" geometry you observe.
Looks to me like you created your input from one of the regtest files. Keep in mind there, the regtest are tuned for speed not for meaningful results. They are simply consistency checks for which convergence doesn't matter. For more meaningful inputs have a look at the benchmarks but even there they are not meant for high quality results.
Flo

[Michael Banck]

Hi,

On Wed, Feb 19, 2014 at 11:29:02PM -0800, Florian Schiffmann wrote:
> this is clearly a cutoff problem (100 is way too low for whatever you want

> to do). This leads to translational invariance with respect to the grid and

> artificial minima arise simply from the atomic positions relative to the
> grid which explains the strange "minimal energy" geometry you observe.
> Looks to me like you created your input from one of the regtest files. Keep
> in mind there, the regtest are tuned for speed not for meaningful results.
> They are simply consistency checks for which convergence doesn't matter.
> For more meaningful inputs have a look at the benchmarks but even there
> they are not meant for high quality results.

While this should be clear to all users, it might make sense to have a
few example files which are production-ready for people to base their
ideas on, clearly seperating them from the regtests, e.g. in a top-level
examples/ directory.


Michael

[hfruchtl]

Just to report back: The cutoff value was indeed the culprit. With 200 it converges to something that at least consists of two H2O molecules, and 300 gives pretty much the result I expected. It takes a very long time though, and I was drawn to CP2K because of its famed speed.

Anyway, I will move on to systems closer to what I really intend to do (materials and surfaces). Using CASTEP or VASP at the moment, but some of the systems are sparse, so using a non-PW basis should improve things. Yes, I know, no K-points and therefore no metals...

Thanks for your help,

  Herbert

On Thursday, February 20, 2014 7:29:02 AM UTC, Florian Schiffmann wrote:

Hi Herbert,

this is clearly a cutoff problem (100 is way too low for whatever you want to do). This breaks translational invariance with respect to the grid and artificial minima arise simply from the atomic positions relative to the grid which explains the strange "minimal energy" geometry you observe.

Looks to me like you created your input from one of the regtest files. Keep in mind there, the regtest are tuned for speed not for meaningful results. They are simply consistency checks for which convergence doesn't matter. For more meaningful inputs have a look at the benchmarks but even there they are not meant for high quality results.

Flo

[Florian Schiffmann]

Hi Herbert,

just a bit about how CP2K works which shines light on your disappointing experience:

CP2K still uses a plane wave basis to compue the electrostatic energy (because it is easy and cheap in PW). This automatically means make your box larger and your cutoff higher will increase the cost of the calculation. Hence, vacuum is not for free in CP2K (not quite sure how any of your other codes react when you put in such a cell). Nevertheless, using Gaussians as a primary basis for the other contributions (localized basis sets are cheap for the rest, see QUICKSTEP paper) at least the size wavefunction coefficients matrix is independent of the system size.
If you have a look at the timing reports you will see all time is spend in the FFT. I can tell you from experience going to real system the FFT becomes very soon neglegible.
The system you use for testing is not really good as a benchmark for speed. As I said above every aproach has its ups and downs. Take that system as an overall performance benchmark and you will find PW codes (VASP, CPMD, CP2K,...) are the worst DFT codes ever while Gaussian based codes are as close to heaven as you can get. And indeed, for (smaller) molecular systems that is the truth but not if you go to larger periodic systems.
If you want to get a feeling for the speed of different codes give them a chance to perform at what they are good in. Try once the 32/64/... water benchmark (better done in parallel but could work in serial as well). That gives a bit a feeling for how the code performs on real dense periodic systems.

Another point you wrote is:

> but some of the systems are sparse, so using a non-PW basis should improve things

As written above that doesn't matter too much for CP2K. True, vacuum only affects the FFT part and if that is neglegible in the overall timings, vacuum is for free. But also dense systems benefit from the mixed basis aproach (try the benchmark).

Another point is that CP2K is meant to be a parallel code. We can't help it that systems of a given size become somewhat expensive. What we can do is try to make everything as parallel as possible. If you go parallel you can 'easiliy' go up to systems with 1000 atoms and get your answer by tomorrow (or after a coffee break, as always it depends on what you are doing and on how much coffee you drink).

That's good about CP2K. Its worst bit you have already said no k-points. That means you won't get away with small systems.
As always which code to use depends on what you are interested in. For crystals and surfaces (and other things for sure but I am not an expert) VASP is neat, it has heaps of features (properties it can compute,...) and in combination with k-points allows for smaller unit cells,... If thats what you want, great choice of code.

If you want to study adsorption or interfaces (which can become very big) you might want to use different codes (OK, I admit I meant CP2K here) which allow you to simulate systems of the necessary size.

Finally the thing I wouldn't do is a mixed use of codes. Yes, it is DFT in both cases but every code has its special bits. Basis set dependece of the results, pseudos,... you really have to make sure that your results are transferable with the settings used in the two codes. And that usually ends up doing the same study with two codes, which means twice the ammount of work. Don't get me wrong, it is definitely a good idea if you get strange or unexpected stuff to check it with a different code/setup. I only say in my oppinion it is better to decide in the beginning for a code and stick to it so your results are consistent (again there can be exceptions).

Flo

PS: Just to make sure, I wrote this to indicate what CP2K is good/meant for as it seems there was a bit of a misunderstanding and not to start a this code is better han that code war.