[CP2K:4045] Re: Energy conservation NVT/NVE simulation

hut... at pci.uzh.ch hut... at pci.uzh.ch
Mon Sep 24 08:36:04 UTC 2012


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... at pci.uzh.ch
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland
---------------------------------------------------------------

-----cp... at googlegroups.com wrote: -----
To: cp... at googlegroups.com
From: paolo nicolini 
Sent by: cp... at googlegroups.com
Date: 09/21/2012 05:45PM
Subject: Re: [CP2K:4045] Re: Energy conservation NVT/NVE simulation

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



On Monday, September 17, 2012 1:09:54 PM UTC+2, jgh wrote: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... at pci.uzh.ch 
Winterthurerstrasse 190 
CH-8057 Zurich, Switzerland 
--------------------------------------------------------------- 
 
-----cp... at googlegroups.com wrote: ----- 
To: cp... at googlegroups.com 
From: marco masia  
Sent by: cp... at googlegroups.com 
Date: 09/12/2012 03:39PM 
Subject: [CP2K:4029] Re: Energy conservation NVT/NVE simulation 
 
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 
 
On Wednesday, September 5, 2012 3:32:26 AM UTC-4, albe... at virgilio.it wrote: 
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 
 
   
  --  
 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... at googlegroups.com. 
 To unsubscribe from this group, send email to cp2k+... at 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]  
  -- 
 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/-/nmvFu5yyMKQJ.
  To post to this group, send email to cp... at googlegroups.com.
 To unsubscribe from this group, send email to cp2k+uns... at 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]


More information about the CP2K-user mailing list