Problem with 2nd generation CPMD

[Onofrio Tau]

Dear CP2K community,

I'm new to CP2K and I'm trying to set up a second-generation Car-Parrinello (SGCP) MD simulation. I have started with a geometrically optimized system of toluene adsorbed on a 3-layer Cu slab (8x6). I would to know what sections/keywords in the input file must be used to correctly set up a SGCPMD simulation and, possibly, the best settings to improve the speed of such simulations since I have not found any related tutorial (but only few things in the cp2k google group).
If I well understood, it has to be turn on the Langevin dynamics section in combination with the OT method. For example:

&FORCE_EVAL
  ...
  &DFT
    ...
    &QS
      METHOD GPW
      EXTRAPOLATION ASPC
      EXTRAPOLATION_ORDER 3
      EPS_DEFAULT 1.0E-10
    &END QS
    ...
    &SCF
    ...
      &OT
        PRECONDITIONER FULL_SINGLE_INVERSE   ! or FULL_KINETIC
        MINIMIZER DIIS
        N_DIIS 7
      &END OT
      &OUTER_SCF
        MAX_SCF 3
        EPS_SCF 1.0E-05
      &END OUTER_SCF    
    ...
    &END SCF
    ...
  &END DFT
  ...
&END FORCE_EVAL

&MOTION
  &MD
    ENSEMBLE LANGEVIN
    ...
    &LANGEVIN
      GAMMA 0.001
    &END LANGEVIN
    ...
  &END MD
  ...
&END MOTION

Is it correct? To follow what I just mentioned, I prepared an input file (attached below) based on sections/keywords defined just above and I got from the output file (attached below) only
    "Leaving inner SCF loop after reaching   200 steps" 

(200 steps because I set MAX_SCF 200)
resulting in a very bad convergence.
As a result, it seems that I have not been running any SGCPMD simulation but a simple bad BOMD. Am I missing something?
How can I correctly set up my input file for SGCPMD?

Any help would be greatly appreciated.

Thanks.

Onofrio

[Sam Broderick]

Hi Onofrio

Question: have you set this up and run a more standard simulation to work out the convergence first?

If I guess, it's probably Cu as a transition metal with d-orbitals that is difficult (I am suffering with Au).

Kind Regards

Sam 

[Onofrio Tau]

Hi Sam,

yes, of course. I set up more standard BOMD with ENSEMBLE NVT and DIAGONALIZATION with smear (look at my input file) and worked fine: it took more or less 30 steps to reach convergence every AIMD step.

Then, I tried with ENSEMBLE NVT and OT (without smear obviously) and I got (MAX_SCF 200):

 *** WARNING in qs_scf.F:598 :: SCF run NOT converged ***
 *** WARNING in qs_scf.F:598 :: SCF run NOT converged ***
  *** SCF run converged in    79 steps ***
  *** SCF run converged in    72 steps ***
  *** SCF run converged in    84 steps ***
  *** SCF run converged in    53 steps ***
  *** SCF run converged in    49 steps ***
  *** SCF run converged in    48 steps ***
  *** SCF run converged in    38 steps ***
  *** SCF run converged in    40 steps ***
  *** SCF run converged in    38 steps ***
  *** SCF run converged in    30 steps ***
  *** SCF run converged in    23 steps ***
  *** SCF run converged in    26 steps ***
  *** SCF run converged in    27 steps ***
  *** SCF run converged in    26 steps ***
  *** SCF run converged in    26 steps ***
  *** SCF run converged in    23 steps ***
  *** SCF run converged in    16 steps ***
  *** SCF run converged in    16 steps ***
  *** SCF run converged in    16 steps ***
  *** SCF run converged in    18 steps ***

   ...

Though SCF run not converged for the first two AIMD steps, then it worked sometimes fine (SCF run converged in 16 steps, but not always!). Could be a possible solution.

N.B.: setting MAX_SCF below 200 (i.e., 50) makes SCF run not converged more than two initial AIMD steps.

Anyway, I'm looking for correctly set up SGCPMD as that found in the following works:

- Lan, Jinggang, Jürg Hutter, and Marcella Iannuzzi. "First-principles simulations of an aqueous CO/Pt (111) interface." The Journal of Physical Chemistry C 122.42 (2018): 24068-24076.

- Musso, Tiziana, et al. "Second generation Car-Parrinello MD: application to the h-BN/Rh (111) nanomesh." The European Physical Journal B 91.7 (2018): 1-8.

However, I'm struggling with and I kindly need your help.

Kind regards.

Onofrio

[Marcella Iannuzzi]

As you may have read in previous works, to run SGCPMD a preliminary study is needed to set the proper GAMMA and NOISY_GAMMA parameters, which allow for the conservation of energy, for a Gaussian distribution of the error in the forces, and for a correct equipartition of the kinetic energy among the degrees of freedom. 

Moreover the parameters to tune the convergence of the electronic structure at each MD step need also to be properly adjusted to the system. 

This means that a series of preliminary simulations are needed, before the production runs are started. 

This said, my suggestions are the following:

- Always start sgcpmd from a well converged wfn obtained along a standard BOMD run with diagonalisation and smearing

-  Let the first steps run on many scf iterations if needed, i.e., set MAX_SCF large enough

-  But use a relative large eps_scf. The wfn is not going to be optimised as in the BOMD, energy convergence  around 10^-5 is desirable

- Use FULL_ALL preconditioner + ASPC

- Use a small time step

- Tune properly GAMMA and NOISY_GAMMA

- Monitor the conservation of energy and the temperature per kind or per region : TEMP_KIND

Regards

Marcella

[Onofrio Tau]

Dear Marcella,

your suggestions were so useful! I was able to reduce to 3-4 iterations per MD step. Sorry for my late reply but I played a lot with GAMMA and NOISE_GAMMA to equilibrate the MD run by monitoring the temperature per each kind. However, I have gathered some unusual behaviour (following attached files represent Temp per kind vs MD steps):

1. I first assumed GAMMA 0.001 fs^-1 (value typically used in other works) and gave different NOISE_GAMMA values for each &DEFINE_REGION (region1: Cu slab ; region2: one toluene molecule). In all cases, for the first 100 MD steps the temperature rapidly drops to very lower values (see attached gamma0.001.png file as an example), though it partially riequilibrates to the target T (900 K) only for sufficiently larger NOISE_GAMMA. Must be the initial temperature drop avoided or is it sufficient only adjusting the NOISE_GAMMA to larger values for T compensation (thus fixing GAMMA to 0.001)?

2. To avoid initial temperature drops, I have tried to directly increase GAMMA. For 0.1 fs^-1 it seems to work well (attached gamma0.1.png file), but it acts like an 'aggressive' thermostat producing large T fluctuations, especially for the toluene region (similar situations for GAMMA = 0.01 fs^-1, attached gamma0.01.png file). I'm aware of the few degrees of freedom of the single toluene molecule causing large temperature fluctuations; therefore it is better to assume GAMMA 0.001, is it?

3. Most important: the number of iterations per MD step costantly increase from 3-4 to 13-16 iterations during MD runs as it can be seen in the attached output file taken as an example (Cu_toluene.out) with

    grep 'SCF run' Cu_toluene.out
To restore the convergence with 3-4 iterations per MD step I sometimes did external EXIT and then restart the run with .RESTART.wfn. It initially works, but then the number of iterations increases again.
Is it a normal behaviour or is it due to a bad instability of the Langevin dynamics? How to eventually avoid it?

Many thanks!

N.B.: No external thermostat and no velocity rescaling (TEMP_TOL) were used during simulations.

Onofrio

[Marcella Iannuzzi]

Dear Onofrio

The system seems not to be very well equilibrated, indeed

The fact that the number of SCF iterations tends to increase is also not a very good sign

Probably I would try and run longer BOMD before starting sgcp, but it happens that finding the good parameters is anyway cumbersome 

Regards

Marcella