energy drift in NPT MD simulation

satya

unread,

Sep 19, 2013, 12:13:47 AM9/19/13

to cp...@googlegroups.com

Hi,

I am performing ab initio NPT simulations for periodic system. After few MD steps, the energy became positive suddenly.

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

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

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

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

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

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

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

I enclosed my input below. Could you please suggest possible ways to rectify it.

&FORCE_EVAL

METHOD Quickstep

&DFT

BASIS_SET_FILE_NAME /home1/bala/cp2k/tests/QS/BASIS_MOLOPT

POTENTIAL_FILE_NAME /home1/bala/cp2k/tests/QS/GTH_POTENTIALS

CHARGE = 0

&MGRID

CUTOFF 350

NGRIDS 5

REL_CUTOFF 50

&END MGRID

&QS

METHOD GPW

EPS_DEFAULT 1.0E-11

&END QS

&SCF

SCF_GUESS RESTART

EPS_SCF 1.0E-7

MAX_SCF 50

CHOLESKY OFF

&OT

MINIMIZER DIIS

PRECONDITIONER FULL_ALL

ENERGY_GAP 0.001

# STEPSIZE 0.05

&END OT

&OUTER_SCF

EPS_SCF 1E-7

MAX_SCF 150

&END OUTER_SCF

&END SCF

&XC

&XC_FUNCTIONAL PBE

&END XC_FUNCTIONAL

&vDW_POTENTIAL

POTENTIAL_TYPE PAIR_POTENTIAL

&PAIR_POTENTIAL

TYPE DFTD3

PARAMETER_FILE_NAME /home1/bala/cp2k/tests/QS/dftd3.dat

REFERENCE_FUNCTIONAL PBE

REFERENCE_C9_TERM .TRUE.

CALCULATE_C9_TERM .TRUE.

R_CUTOFF 25.0

&END PAIR_POTENTIAL

&END vDW_POTENTIAL

&XC_GRID

XC_DERIV SPLINE2

XC_SMOOTH_RHO NN50

&END XC_GRID

&END XC

&END DFT

&SUBSYS

&CELL

ABC 21.8620 16.2099 43.6020

ALPHA_BETA_GAMMA 90.0 92.7540 90.0

&CELL_REF

ABC 22.8620 17.2099 44.6020

ALPHA_BETA_GAMMA 90.0 126.0 90.0

&END CELL_REF

&END CELL

&COORD

@INCLUDE coor.xyz

&END COORD

&KIND H

BASIS_SET DZVP-MOLOPT-SR-GTH

POTENTIAL GTH-PBE-q1

&END KIND

&KIND C

BASIS_SET DZVP-MOLOPT-SR-GTH

POTENTIAL GTH-PBE-q4

&END KIND

&KIND N

BASIS_SET DZVP-MOLOPT-SR-GTH

POTENTIAL GTH-PBE-q5

&END KIND

&KIND O

BASIS_SET DZVP-MOLOPT-SR-GTH

POTENTIAL GTH-PBE-q6

&END KIND

&KIND Fe

BASIS_SET DZVP-MOLOPT-SR-GTH

POTENTIAL GTH-PBE-q16

&END KIND

&KIND Zn

BASIS_SET DZVP-MOLOPT-SR-GTH

POTENTIAL GTH-PBE-q12

&END KIND

&KIND Cd

BASIS_SET DZVP-MOLOPT-SR-GTH

POTENTIAL GTH-PBE-q12

&END KIND

&END SUBSYS

STRESS_TENSOR ANALYTICAL

&END FORCE_EVAL

&GLOBAL

PROJECT OUT-GEO-OPT

RUN_TYPE MD

PRINT_LEVEL MEDIUM

&END GLOBAL

&MOTION

&MD

ENSEMBLE NPT_F

STEPS 10000

TIMESTEP 0.5

TEMPERATURE 195.0

&BAROSTAT

PRESSURE 1.0

TIMECON [wavenumber_t] 4000.

&THERMOSTAT

TYPE CSVR

&END THERMOSTAT

&END BAROSTAT

&THERMOSTAT

TYPE NOSE

REGION GLOBAL

&NOSE

TIMECON [wavenumber_t] 4000.0

LENGTH 3

YOSHIDA 3

MTS 2

&END NOSE

&END THERMOSTAT

&END MD

&PRINT

&RESTART

LOG_PRINT_KEY T

&EACH

MD 50

&END EACH

ADD_LAST NUMERIC

&END RESTART

&TRAJECTORY

LOG_PRINT_KEY T

&EACH

MD 1

&END EACH

ADD_LAST NUMERIC

&END TRAJECTORY

&VELOCITIES

LOG_PRINT_KEY T

&EACH

MD 1

&END EACH

ADD_LAST NUMERIC

&END VELOCITIES

&END PRINT

&END MOTION

#&EXT_RESTART

# RESTART_FILE_NAME OUT-GEO-OPT-1.restart

#&END EXT_RESTART

Thank you..

Satyanarayana B

Matt W

unread,

Sep 19, 2013, 4:49:49 AM9/19/13

to cp...@googlegroups.com

Hi,

well your system is clearly exploding. Your cell_ref parameters look inconsistent with your actual cell - maybe the stress tensor is very bad!

Bonakala Satyanarayana

unread,

Sep 19, 2013, 5:01:02 AM9/19/13

to cp...@googlegroups.com

Dear Matt W,

The cell_ref parameters are expected cell parameters, which are observed through experiment. When I had gone through the previous mailing list for reference cell parameters, they suggested that these parameters should be close to the expected cell parameters. So that I have given those parameters.

Thank you

B.Satyanarayana, Ph. D
Molecular modeling lab

JNCASR
Bangalore

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To post to this group, send email to cp...@googlegroups.com.
Visit this group at http://groups.google.com/group/cp2k.
For more options, visit https://groups.google.com/groups/opt_out.

Ari Paavo Seitsonen

unread,

Sep 19, 2013, 5:04:30 AM9/19/13

to cp...@googlegroups.com

Dear Satyanarayana B,

Indeed the energy does not "drift", it explodes. What did happen in an NVT simulation with otherwise the same input, did the simulation work fine?

To me the cut-off energy of 350 Ry looks low, even with the PBE functional (smoother than BLYP); did you test that it is sufficient?

Greetings from Zurich,

apsi

PS Does some one know, are the methods in &XC_GRID ... &END XC_GRID consistent with an NPT simulation?

2013/9/19 Matt W <MattWa...@gmail.com>

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To post to this group, send email to cp...@googlegroups.com.
Visit this group at http://groups.google.com/group/cp2k.
For more options, visit https://groups.google.com/groups/opt_out.

--
-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-
Ari P Seitsonen / Ari.P.S...@iki.fi / http://www.iki.fi/~apsi/
Physikalisch-Chemisches Institut der Universität Zürich
Tel: +41 44 63 55 44 97 / Mobile: +41 79 71 90 935

hut...@pci.uzh.ch

unread,

Sep 19, 2013, 5:09:14 AM9/19/13

to cp...@googlegroups.com

Hi

- Is your system stable in NVT dynamics?
- Have also a look at the forces and stress tensor. Anything special there?
- What about the SCF convergence in the first steps? Is it smooth or
do you get strange behaviour?

regards

Juerg

--------------------------------------------------------------
Juerg Hutter Phone : ++41 44 635 4491
Physical Chemistry Institute FAX : ++41 44 635 6838
University of Zurich E-mail: hut...@pci.uzh.ch
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland
---------------------------------------------------------------

-----cp...@googlegroups.com wrote: -----
To: cp...@googlegroups.com
From: Bonakala Satyanarayana
Sent by: cp...@googlegroups.com
Date: 09/19/2013 11:01AM
Subject: Re: [CP2K:4619] Re: energy drift in NPT MD simulation

Bonakala Satyanarayana

unread,

Sep 20, 2013, 12:41:05 AM9/20/13

to cp...@googlegroups.com

Dear Seitsonen,

Thank you for your reply.

1) I performed NVT simulation with same input, what I have enclosed. It ran well. These are the converged energy values:

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

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

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

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

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

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

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

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

2) I have taken cut-off as per literature (Elucidating the Breathing of the Metal–Organic Framework MIL-53(Sc) with ab Initio Molecular Dynamics Simulations and in Situ X-ray Powder Diffraction Experiments, Duren et al., JACS). Here, they have tested for MIL MOFs and reported optimum cut-off is 350 Ry.

Dear Hutter,

1) Yes, my system is stable in NVT simulation

2) Det( stress tensor) values are

Det(stress tensor) : -7.84528984E-04

Det(stress tensor) : -1.66559235E-02

Det(stress tensor) : -1.44160312E+00

Det(stress tensor) : 2.03287049E+02

Det(stress tensor) : -4.82023324E+03

Det(stress tensor) : 4.36251297E+07

BAROSTAT temperature is increasing:

BAROSTAT TEMP[K] = 0.250067070500E+03

BAROSTAT TEMP[K] = 0.407906043312E+05

BAROSTAT TEMP[K] = 0.373591928441E+06

BAROSTAT TEMP[K] = 0.458891514234E+08

3) SCF convergence is smooth. I enclosed the convergence vs steps

Thank you

B.Satyanarayana, Ph. D
Molecular modeling lab

JNCASR
Bangalore

conver.pdf

Ari Paavo Seitsonen

unread,

Sep 20, 2013, 1:17:05 AM9/20/13

to cp...@googlegroups.com

Dear B.Satyanarayana,

My personal opinions:

1) What is determined as "converged" depends on the quantity observed and accuracy wanted (thus on the person making the choice). For example the pressure tensor, being a second derivative of the energy, converges much slower than the forces (first-order derivative) that in turn converges much slower than the energy (zero-order quantity)

2) In my opinion there are several simulations in the literature that are not converged (and this does not apply to CP2K only). I don't make any statement about the convergence in the article that you cited, just that it is better to test oneself: After that one also understands the convergence, the system, the method etc better

3) The energy on the line 'FORCE_EVAL' is the potential energy, not the quantity that should be conserved. This shows that the energy does not explode, not more (and the sample is very short, plus probably taken at the beginning of the simulation when the things still vary more, due to the rough equilibration needed

4) You see that the value of the determinant of the stress tensor diverges _very_ fast, so there is certainly at least one problem. I was hinted that the time constant that you have taken for the barostat is excessively small: This would basically try to make the cell respond to the fastest oscillations in your system (I hope that there _is_ something oscillating at 4000 cm^-1, it would be necessary for the thermostat chain you have). I don't know if this is the only reason, but I would switch the time constant to something larger, for example of the order picoseconds

5) The convergence in the electronic structure indeed looks smooth on the linear scale, but here one usually uses logarithmic scale in order to see how the energy converges closer to the self-consistency. Another useful check is to calculate the eigenvalues and see that the HOMO-LUMO gap (in molecules) or between the valence and conducation band is more or less that you expected

This after a quick glimpse on your report, maybe I missed again something relevant though.

Good Luck, :)

apsi

2013/9/20 Bonakala Satyanarayana <satyanaray...@gmail.com>

Reply all

Reply to author

Forward