Energy conservation NVT/NVE simulation

alber...@virgilio.it

unread,

Sep 5, 2012, 3:32:26 AM9/5/12

to cp...@googlegroups.com

Dear all,

I'm carrying out some tests to reproduce the VDOS of water and I used as

starting points the input files 'H2O-32.inp' and 'H2O-64.inp' in

tests/QS/benchmark directory.

In both cases I ran a NVT simulation at 300K for 15 ps followed by an NVE

simulation on the equilibrated system.

I'm puzzled about the results in the energies: in both NVT and NVE the plot of

"Cons Qty" show a constant and linear drift while temperature, kinetic and

potential energy are equilibrated and stable.

Since "Cons Qty" should be conserved and this is not the case I think that

there are some problems in the setup of the simulation.

I did not modify the parameters of the benchmark inputs but I just added the

thermostat section in NVT and I modified only the printing option.

Thank you very much for your help and suggestions,

Alberto

hut...@pci.uzh.ch

unread,

Sep 5, 2012, 3:56:13 AM9/5/12

to cp...@googlegroups.com

Hi

the main problem is

&QS
EPS_DEFAULT 1.0E-12
WF_INTERPOLATION PS
EXTRAPOLATION_ORDER 3
&END QS

change to WF_INTERPOLATION ASPC (the default) and maybe
increase the EXTRAPOLATION_ORDER slightly.
See
J. Hutter, Car-Parrinello molecular dynamics
Wiley Interdisciplinary Reviews: Computational Molecular Science 2 (4) , pp. 604-612
for an explanation.

You should also change to a more reasonable DFT functional.
LDA (PADE) will lead to an extremly overstructured system.
See
J. Schmidt et al.,
Isobaric-isothermal molecular dynamics simulations utilizing density functional theory: An assessment of the structure and density of water at near-ambient conditions
Journal of Physical Chemistry B Volume 113, Issue 35, 3 September 2009, Pages 11959-11964

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: "alber...@virgilio.it"
Sent by: cp...@googlegroups.com
Date: 09/05/2012 09:44AM
Subject: [CP2K:4010] Energy conservation NVT/NVE simulation

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To view this discussion on the web visit https://groups.google.com/d/msg/cp2k/-/LgSDjVEpz0YJ.
To post to this group, send email to cp...@googlegroups.com.
To unsubscribe from this group, send email to cp2k+uns...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/cp2k?hl=en.

alber...@virgilio.it

unread,

Sep 11, 2012, 3:58:15 AM9/11/12

to cp...@googlegroups.com

Dear Juerg,

thank you very much.

I changed the functional to BLYP and I used ASPC interpolator.

I carried out different tests where I tried different values of the

EXTRAPOLATION_ORDER (up to 6), I increased the EPS_SCF and used defaults

for most of the other parameters.

Unfortunately, I still observe either drifts or jumps in the total energy.

I think that there is still a problem in the input, reported below.

Many thanks for your help.

Alberto

************

&FORCE_EVAL

METHOD QS

&DFT

BASIS_SET_FILE_NAME GTH_BASIS_SETS

POTENTIAL_FILE_NAME GTH_POTENTIALS

&MGRID

CUTOFF 300

REL_CUTOFF 40

&END MGRID

&QS

EPS_DEFAULT 1.0E-12

WF_INTERPOLATION ASPC

EXTRAPOLATION_ORDER 3

&END QS

&SCF

SCF_GUESS ATOMIC

EPS_SCF 1.0E-6

&OT T

&END OT

&PRINT

&RESTART OFF

&END

&END SCF

&XC

&XC_FUNCTIONAL BLYP

&END XC_FUNCTIONAL

&END XC

&PRINT

&MOMENTS MEDIUM

&END

&END PRINT

&END DFT

&SUBSYS

&CELL

ABC 9.8528 9.8528 9.8528

&END CELL

# 32 H2O (TIP5P,1bar,300K) a = 9.8528

&COORD

O 2.280398 9.146539 5.088696

O 1.251703 2.406261 7.769908

O 1.596302 6.920128 0.656695

O 2.957518 3.771868 1.877387

O 0.228972 5.884026 6.532308

O 9.023431 6.119654 0.092451

O 7.256289 8.493641 5.772041

O 5.090422 9.467016 0.743177

O 6.330888 7.363471 3.747750

O 7.763819 8.349367 9.279457

O 8.280798 3.837153 5.799282

O 8.878250 2.025797 1.664102

O 9.160372 0.285100 6.871004

O 4.962043 4.134437 0.173376

O 2.802896 8.690383 2.435952

O 9.123223 3.549232 8.876721

O 1.453702 1.402538 2.358278

O 6.536550 1.146790 7.609732

O 2.766709 0.881503 9.544263

O 0.856426 2.075964 5.010625

O 6.386036 1.918950 0.242690

O 2.733023 4.452756 5.850203

O 4.600039 9.254314 6.575944

O 3.665373 6.210561 3.158420

O 3.371648 6.925594 7.476036

O 5.287920 3.270653 6.155080

O 5.225237 6.959594 9.582991

O 0.846293 5.595877 3.820630

O 9.785620 8.164617 3.657879

O 8.509982 4.430362 2.679946

O 1.337625 8.580920 8.272484

O 8.054437 9.221335 1.991376

H 1.762019 9.820429 5.528454

H 3.095987 9.107088 5.588186

H 0.554129 2.982634 8.082024

H 1.771257 2.954779 7.182181

H 2.112148 6.126321 0.798136

H 1.776389 7.463264 1.424030

H 3.754249 3.824017 1.349436

H 3.010580 4.524142 2.466878

H 0.939475 5.243834 6.571945

H 0.515723 6.520548 5.877445

H 9.852960 6.490366 0.393593

H 8.556008 6.860063 -0.294256

H 7.886607 7.941321 6.234506

H 7.793855 9.141028 5.315813

H 4.467366 9.971162 0.219851

H 5.758685 10.102795 0.998994

H 6.652693 7.917443 3.036562

H 6.711966 7.743594 4.539279

H 7.751955 8.745180 10.150905

H 7.829208 9.092212 8.679343

H 8.312540 3.218330 6.528858

H 8.508855 4.680699 6.189990

H 9.742249 1.704975 1.922581

H 8.799060 2.876412 2.095861

H 9.505360 1.161677 6.701213

H 9.920117 -0.219794 7.161006

H 4.749903 4.186003 -0.758595

H 5.248010 5.018415 0.403676

H 3.576065 9.078451 2.026264

H 2.720238 9.146974 3.273164

H 9.085561 4.493058 9.031660

H 9.215391 3.166305 9.749133

H 1.999705 2.060411 1.927796

H 1.824184 0.564565 2.081195

H 7.430334 0.849764 7.438978

H 6.576029 1.537017 8.482885

H 2.415851 1.576460 8.987338

H 2.276957 0.099537 9.289499

H 1.160987 1.818023 4.140602

H 0.350256 2.874437 4.860741

H 5.768804 2.638450 0.375264

H 7.221823 2.257514 0.563730

H 3.260797 5.243390 5.962382

H 3.347848 3.732214 5.988196

H 5.328688 9.073059 5.982269

H 5.007063 9.672150 7.334875

H 4.566850 6.413356 3.408312

H 3.273115 7.061666 2.963521

H 3.878372 7.435003 6.843607

H 3.884673 6.966316 8.283117

H 5.918240 3.116802 5.451335

H 5.355924 2.495093 6.711958

H 5.071858 7.687254 10.185667

H 6.106394 7.112302 9.241707

H 1.637363 5.184910 4.169264

H 0.427645 4.908936 3.301903

H 9.971698 7.227076 3.709104

H 10.647901 8.579244 3.629806

H 8.046808 5.126383 2.213838

H 7.995317 4.290074 3.474723

H 1.872601 7.864672 7.930401

H 0.837635 8.186808 8.987268

H 8.314696 10.115534 2.212519

H 8.687134 8.667252 2.448452

&END COORD

&KIND H

BASIS_SET TZV2P-GTH

POTENTIAL GTH-BLYP-q1

&END KIND

&KIND O

BASIS_SET TZV2P-GTH

POTENTIAL GTH-BLYP-q6

&END KIND

&END SUBSYS

&END FORCE_EVAL

&GLOBAL

PROJECT H2O-32

RUN_TYPE MD

PRINT_LEVEL MEDIUM

&TIMINGS

THRESHOLD 0.000001

&END

&END GLOBAL

&MOTION

&MD

ENSEMBLE NVT

STEPS 30000

TIMESTEP 0.5

TEMPERATURE 300.0

&THERMOSTAT

TYPE CSVR

&CSVR

TIMECON 100.0

&END CSVR

&END THERMOSTAT

&END MD

&PRINT

&VELOCITIES LOW

&END

&END PRINT

&END MOTION

**************

Matt W

unread,

Sep 11, 2012, 4:39:07 AM9/11/12

to cp...@googlegroups.com

Hi,

as a test, try running with PBE instead. BLYP can be numerically difficult.

Matt

hut...@pci.uzh.ch

unread,

Sep 11, 2012, 5:03:56 AM9/11/12

to cp...@googlegroups.com

Hi

it is difficult to give good advice without more information.

Possible problems might be:

- Restarts: do you correctly restart all variables
- Thermostat: I don't know if the setting you have is reasonable
and adequate for this timestep.
Do you see a drift also for NVE?
- How big is the drift and how long is your trajectory?

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: "alber...@virgilio.it"
Sent by: cp...@googlegroups.com

Date: 09/11/2012 10:06AM
Subject: Re: [CP2K:4024] Energy conservation NVT/NVE simulation

To view this discussion on the web visit https://groups.google.com/d/msg/cp2k/-/yrA7tbiv8B4J.

alber...@virgilio.it

unread,

Sep 11, 2012, 5:35:42 AM9/11/12

to cp...@googlegroups.com

Hi,

for the restart I used the last restart file of the previous run, modifying the parameters in $QS or $SCF to test the new ones.

I observed a drift also for NVE even if most of the tests have been done for NVT.

I will carry out more tests for NVE to reduce the number of variables and avoid also the problem of thermostat setting.

About the drift, after few ps it is about 10^-3 au, depending on the simulations.

Changing the parameters to increase the accuracy, in some cases I observed no drift but sudden jumps in the total energy.

I will keep you updated about further NVE tests.

Thanks

alberto

garold

unread,

Sep 12, 2012, 6:38:09 AM9/12/12

to cp...@googlegroups.com

Hi Alberto,

May I suggest use of smoothing, i.e., including a section:

&XC_GRID

XC_SMOOTH_RHO NN10 ! (or NN50)

XC_DERIV SPLINE2_SMOOTH

&END XC_GRID

You can also see the previous discussion on this forum:

https://groups.google.com/forum/?hl=it&fromgroups=#!topic/cp2k/dc4kzNqNw-8

where a similar discussion occurred.

Best,

Garold

marco masia

unread,

Sep 12, 2012, 9:38:52 AM9/12/12

to cp...@googlegroups.com

Dear all,

recently I have noticed the same problem with FIST simulations. I took the H2O-32_ewald.inp input, equilibrated it for 1 ns and run a 1ns NVE simulation from the thermalized configuration. I have tried to refine many simulation parameters (e.g. time step, RCUT, Ewald pars) and also (thinking that there was a bad implementation of constraints) I run the simulations for a flexible model. I always get a non-negligible drift (see attachment); the graph shows block averaged energies of 1 ns simulations performed with cp2k and obeli.x (a program of mine). The values are shifted for a better comparison. I have been looking into this, thinking that the problem was in the input file...but since other people are experiencing it with QS, I wonder if there is some bug in the velocity verlet routines. The results are platform independent (I tried both on mac and linux).

best

marco

energy-conservation.pdf

alber...@virgilio.it

unread,

Sep 17, 2012, 5:52:58 AM9/17/12

to cp...@googlegroups.com

Hi,

I have solved my problems either by using PBE functional or using BLYP with the smoothing parameters suggested by Garold.

BLYP seeems indeed numerically quite difficult and delicate while PBE is better from this point of view.

Concerning the other parameters, in particular in &QS section, I used those suggested above by Juerg.

Thank you all for your suggestions,

Alberto

hut...@pci.uzh.ch

unread,

Sep 17, 2012, 7:09:51 AM9/17/12

to cp...@googlegroups.com

Hi

I did some tests with the inputs from
../tests/Fist
and get much better energy conservation that you show.
Probably your input use a not optimal setting for the
cutoffs and Ewald parameters (I was using SPME).

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: marco masia
Sent by: cp...@googlegroups.com
Date: 09/12/2012 03:39PM
Subject: [CP2K:4029] Re: Energy conservation NVT/NVE simulation

--
You received this message because you are subscribed to the Google Groups "cp2k" group.

To view this discussion on the web visit https://groups.google.com/d/msg/cp2k/-/sAiVZ2dob1sJ.

To post to this group, send email to cp...@googlegroups.com.
To unsubscribe from this group, send email to cp2k+uns...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/cp2k?hl=en.

[attachment "energy-conservation.pdf" removed by Jürg Hutter/at/UZH]

paolo nicolini

unread,

Sep 21, 2012, 11:34:56 AM9/21/12

to cp...@googlegroups.com

Hi

I tried to use the H2O-32_SPME.inp in the tests/Fist directory by changing only the length of the simulation (1 ns) and I got also a significant energy drift (see attachment). Then I tried to reduce the integration time step to 1.0 fs. The drift become lower but it do not vanish. I did also 1ns of NVE dynamics (starting from the last restart of the previous simulations) and here the drift become even more clear.
The inputs that I used are below.
Best regards
Paolo

NVT input:
&FORCE_EVAL
METHOD FIST
&MM
    &FORCEFIELD
      &BEND
        ATOMS H O H
        K 0.
        THETA0 1.8
      &END BEND
      &BOND
        ATOMS O H
        K 0.
        R0 1.8
      &END BOND
      &CHARGE
        ATOM O
        CHARGE -0.8476
      &END CHARGE
      &CHARGE
        ATOM H
        CHARGE 0.4238
      &END CHARGE
      &NONBONDED
        &LENNARD-JONES
          atoms O O
          EPSILON 78.198
          SIGMA 3.166
          RCUT 11.4
        &END LENNARD-JONES
        &LENNARD-JONES
          atoms O H
          EPSILON 0.0
          SIGMA 3.6705
          RCUT 11.4
        &END LENNARD-JONES
        &LENNARD-JONES
          atoms H H
          EPSILON 0.0
          SIGMA 3.30523
          RCUT 11.4
        &END LENNARD-JONES
      &END NONBONDED
    &END FORCEFIELD
    &POISSON
      &EWALD
        EWALD_TYPE spme
        ALPHA .5
        GMAX 21
        O_SPLINE 6
      &END EWALD
    &END POISSON
&END MM
&SUBSYS
    &CELL
      ABC 9.865 9.865 9.865
    &END CELL
    &COORD
    O                  -4.583   5.333   1.560   H2O
    H                  -3.777   5.331   0.943   H2O
    H                  -5.081   4.589   1.176   H2O
...
    &END COORD

&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT H2O

RUN_TYPE MD
&END GLOBAL
&MOTION
&CONSTRAINT
    &G3X3
      DISTANCES 1.8897268 1.8897268 3.0859239
      MOLECULE 1
      ATOMS 1 2 3
    &END G3X3
&END CONSTRAINT
&MD
    ENSEMBLE NVT
    STEPS 1000000
    TIMESTEP 1.0
    TEMPERATURE 300.0
    &THERMOSTAT
      REGION MOLECULE
      &NOSE
        LENGTH 3
        YOSHIDA 3
        TIMECON 1000
        MTS 2
      &END NOSE
    &END
&END MD
&END MOTION

NVE input:
&FORCE_EVAL
METHOD FIST
&MM
    &FORCEFIELD
      &BEND
        ATOMS H O H
        K 0.
        THETA0 1.8
      &END BEND
      &BOND
        ATOMS O H
        K 0.
        R0 1.8
      &END BOND
      &CHARGE
        ATOM O
        CHARGE -0.8476
      &END CHARGE
      &CHARGE
        ATOM H
        CHARGE 0.4238
      &END CHARGE
      &NONBONDED
        &LENNARD-JONES
          atoms O O
          EPSILON 78.198
          SIGMA 3.166
          RCUT 11.4
        &END LENNARD-JONES
        &LENNARD-JONES
          atoms O H
          EPSILON 0.0
          SIGMA 3.6705
          RCUT 11.4
        &END LENNARD-JONES
        &LENNARD-JONES
          atoms H H
          EPSILON 0.0
          SIGMA 3.30523
          RCUT 11.4
        &END LENNARD-JONES
      &END NONBONDED
    &END FORCEFIELD
    &POISSON
      &EWALD
        EWALD_TYPE spme
        ALPHA .5
        GMAX 21
        O_SPLINE 6
      &END EWALD
    &END POISSON
&END MM
&END FORCE_EVAL
&GLOBAL
PROJECT H2O
RUN_TYPE MD
&END GLOBAL
&MOTION
&CONSTRAINT
    &G3X3
      DISTANCES 1.8897268 1.8897268 3.0859239
      MOLECULE 1
      ATOMS 1 2 3
    &END G3X3
&END CONSTRAINT
&MD
    ENSEMBLE NVE
    STEPS 1000000
    TIMESTEP 1.0
    TEMPERATURE 300.0
&END MD
&END MOTION
&EXT_RESTART
RESTART_FILE_NAME H2O-1_1000000.restart
&END EXT_RESTART

all.pdf

hut...@pci.uzh.ch

unread,

Sep 24, 2012, 4:36:04 AM9/24/12

to cp...@googlegroups.com

Hi

my guess is that you can get better results by tuning the Ewald
parameters.
I see the following sources of errors that can lead to drifts.
1) the cutoff for the LJ potential (there is no smoothing at the cutoff used)
2) An inadequate combination of length cutoff (real space), alpha
and grid cutoff
3) timestep

In CP2K the two cutoff in 1) and 2) are the same. A rule of thumb
says alpha*rc should be at least 7 (0.5*11.4=5.7 in your case).

To get better conservation of energy I would:

- shorten the timestep (2fs probably)
- increase the LJ cutoff (14)
- increase alpha (1.) and the grid cutoff (more points)
Warning: the grid cutoff is given as number of points, if you
go to bigger systems you have to scale the number of points with
the box size.

I don't do long time simulations with classical force fields,
so the values above are only guesses.

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: paolo nicolini
Sent by: cp...@googlegroups.com
Date: 09/21/2012 05:45PM
Subject: Re: [CP2K:4045] Re: Energy conservation NVT/NVE simulation

To view this discussion on the web visit https://groups.google.com/d/msg/cp2k/-/nmvFu5yyMKQJ.

To post to this group, send email to cp...@googlegroups.com.
To unsubscribe from this group, send email to cp2k+uns...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/cp2k?hl=en.

[attachment "all.pdf" removed by Jürg Hutter/at/UZH]

Matt W

unread,

Sep 24, 2012, 11:00:18 AM9/24/12

to cp...@googlegroups.com

Hi,

Taking the H2O-32_SPME.inp file and reducing the timestep to 2 fs I get no drift with NVT.

In NVE there is still a little drift, but decreasing the SHAKE tolerance to 1.0E-07 this seems to be removed (see attached, changed the SHAKE tolerance at ~2.1 ns.

In all cases the drift seems to be smaller than being reported, I think. Which version did you use? Mine is serial compiled under cygwin with internal FFTSG.

"CP2K| version string: CP2K version 2.2.356 (Development Version)
CP2K| is freely available from http://cp2k.berlios.de/
CP2K| Program compiled at Tue Sep 27 14:10:08 GMTDT 2011
CP2K| Program compiled on palliata
CP2K| Program compiled for Cygwin-i686-gfortran"

Matt

NVE_conserve.png

NVT_conserve.png

paolo nicolini

unread,

Sep 25, 2012, 5:45:03 AM9/25/12

to cp...@googlegroups.com

Hi

thank you for the quick answers! I will try to follow your suggestions...

My version is:

CP2K| version string:                    CP2K version 2.4 (Development Version)
CP2K| source code revision number:                                        12363
CP2K| is freely available from                             http://www.cp2k.org/
CP2K| Program compiled at                         Fri Sep 7 12:38:53 CEST 2012
CP2K| Program compiled on                                              titania1
CP2K| Program compiled for                                     Linux-i686-intel

Best
Paolo

paolo nicolini

unread,

Oct 10, 2012, 5:48:07 AM10/10/12

to cp...@googlegroups.com

Hi,

I tried some test to assess the suggestions from Jorg and Matt. It follows the results.
First the Matt's suggestions (all_matt_50ps.png):
- only decreasing the time step to 2.0 fs don't change the effect (red lines in NVT and NVE runs);
- putting the SHAKE_TOLERANCE key allows to get better results in the NVT runs but in the NVE ones it is appreciable a drift (blue and green lines).
However, the outcome is that the results with time step 2.0 and 2.5 are quite indistinguishable.

Regard to Jorg's suggestions (all_jorg_50ps.png):
- only changing the ewald parameter we get a drift even larger (red lines);
- by combining ewald changes with an enlarged shake tolerance I get the correct behaviour (blue lines).

To confirm it, I performed also a run with flexible molecules (all_flexible_50ps_cut.png) and I got almost flat profiles.

I tried also to change the real space cutoff value (that is not the same of the LJ interactions) but I got identical results as the one that I showed.

Best
Paolo

all_matt_50ps.png

all_jorg_50ps.png

all_flexible_50ps_cut.png

Reply all

Reply to author

Forward