nose thermostat

[Diederica]

Hello all,

I did a molecular dynamics simulation in the NVT ensemble using a nose
thermostat, but I always observed drift of the conserved quantity
after few ps, that is just accumulating during longer time
simulations. Is this a normal phenomenon and can it be avoided?

I tested the effect of some parameters: scf convergence, timestep,
noseconstant, temperature and level of theory. We thought that the
temperature and the scf-convergence would be of major importance, but
the effect wasn't that big.
Below a summary is given of these different runs:

drift drift
kJ/mol kJ/mol
level scf conv timestep fs nosect fs Temp (K) steps time (fs) per
step per fs
pm3 1.E-05 0.4 100 2000 13288 5315.2 input1 0.00764 0.01910
pm3 1.E-09 0.8 100 2000 6644 5315.2 input1 0.00951 0.01189
pm3 1.E-05 0.8 1000 2000 6644 5315.2 input2 0.01204 0.01505
pm3 1.E-05 0.8 100 2000 6644 5315.2 input3 0.00537 0.00671
pm3 1.E-05 0.8 100 1500 6644 5315.2 input4 0.01037 0.01296
pm3 1.E-09 0.8 100 1500 6644 5315.2 input4 0.00944 0.01181
pm3 1.E-09 0.8 100 1000 6644 5315.2 input4 0.00752 0.00940
pm3 1.E-05 0.8 100 300 6644 5315.2 input3 0.00214 0.00267
blyp 5.E-07 0.8 100 2000 6644 5315.2 input5 0.00156 0.00195
CHARMM 0.8 100 2000 6644 5315.2 input6
0.00106 0.00132

If you can give any comment on this it's appreciated a lot.

Best greetings,
Diederica

For the pm3 calculations I used inputs analogue to the one given
below:
&FORCE_EVAL
METHOD QS
&DFT
&QS
METHOD PM3
&END QS
&SCF
EPS_SCF 1.0E-5
SCF_GUESS ATOMIC
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 25. 25. 25.
UNIT ANGSTROM
&END CELL
&TOPOLOGY
CONNECTIVITY GENERATE
&GENERATE
CREATE_MOLECULES
&END GENERATE
COORDINATE XYZ
COORD_FILE_NAME input.xyz
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
FFT_LIB FFTSG
PRINT_LEVEL LOW
PROJECT md
RUN_TYPE MD
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS 40000
TIMESTEP 0.8
TEMPERATURE 1500
#&THERMOSTAT
&NOSE
TIMECON 100
&END NOSE
#&END THERMOSTAT
&END MD
&PRINT
&RESTART_HISTORY SILENT
EACH 1000
&END RESTART_HISTORY
&END PRINT
&END MOTION

I used an older version of cp2k (2007-04-27) for the pm3 simulations
as the molecule fell apart using the more recent version of
2008-01-22.
The beginning of the output was:
CP2K| Program compiled at Mon Apr 30 17:13:25
CEST 2007
CP2K| Program compiled
on moldyn48
CP2K| Program compiled for Linux-
x86-64-intel
CP2K| Last CVS entry
CP2K| Input file
name md.inp

GLOBAL| Force Environment
number 1
GLOBAL| Basis set file name
BASIS_SET
GLOBAL| Potential file name
POTENTIAL
GLOBAL| MM Potential file name
MM_POTENTIAL
GLOBAL| Restart file
name RESTART
GLOBAL| Coordinate file name
input.xyz
GLOBAL| Method
name CP2K
GLOBAL| Project
name md
GLOBAL| Default FFT
library FFTSG
GLOBAL| Run
type MD
GLOBAL| All-to-all communication in single
precision F
GLOBAL| Global print
level 1
GLOBAL| Total number of message passing
processes 1
GLOBAL| Number of threads for this
process 1
GLOBAL| This output is from
process 0

GENERATE| Preliminary Number of Bonds
generated: 71
GENERATE| Molecules names have been generated. Now reordering
particle set in order to have
GENERATE| atoms belonging to the same molecule in a sequential order.
GENERATE| Achieved consistency in connectivity generation.
GENERATE| Number of Bonds
generated: 71
REORDER |
REORDER | Reordering Molecules. The Reordering of molecules will
override all
REORDER | information regarding molecule names and residue names.
REORDER | New ones will be provided on the basis of the
connectivity!
REORDER | Number of unique molecules found: 1.
REORDER |
GENERATE| Preliminary Number of Bonds
generated: 71
GENERATE| Achieved consistency in connectivity generation.
GENERATE| Number of Bonds
generated: 71
GENERATE| Preliminary Number of Bends
generated: 130
GENERATE| Number of Bends
generated: 130
GENERATE| Number of UB
generated: 130
GENERATE| Preliminary Number of Torsions
generated: 184
GENERATE| Number of Torsions
generated: 184
GENERATE| Number of Impropers
generated: 13
GENERATE| Number of 1-4 interactions
generated: 184

...
...

WARNING| Integral interaction cutoff has been
redefined
WARNING| Old value
[a.u.] 24.000
WARNING| New value
[a.u.] 23.622



SCF PARAMETERS Density
guess: ATOMIC

--------------------------------------------------------

max_scf: 50

max_scf_history: 0

max_diis: 4

--------------------------------------------------------

eps_scf: 1.00E-05

eps_scf_history: 0.00E+00

eps_diis: 1.00E-01

eps_eigval: 1.00E-05

eps_jacobi: 0.00E+00

--------------------------------------------------------

p_mix: 0.40
G-space mixing
a: 1.00
G-space mixing
b: 0.00

work_syevx: 1.00
level_shift
[a.u.]: 0.00
smear
[a.u.]: 0.00
added
MOs 0 0

--------------------------------------------------------
No outer SCF

MD| Molecular Dynamics Protocol
MD| Ensemble
Type nvt
MD| Number of Time
Steps 40000
MD| Time Step
[fs] 0.80
MD| Temperature
[K] 1500.00
MD| Temperature tolerance
[K] 0.00
MD| Nose-Hoover-Chain parameters
MD| Nose-Hoover-Chain
length 3
MD| Nose-Hoover-Chain time constant
[ fs] 100.00
MD| Order of Yoshida
integrator 3
MD| Number of multiple time
steps 2
MD| Print MD information every
1 step(s)
MD| File type Print frequency[steps]
File names
MD| Coordinates 1 md-
pos-1.xyz
MD| Velocities 1 md-
vel-1.xyz
MD| Energies 1
md-1.ener
MD| Dump 1
md-1.restart


Calculation of degrees of freedom
Number of
atoms: 70
Number of Intramolecular
constraints: 0
Number of Intermolecular
constraints: 0
Invariants(translation +
rotations): 3
Degrees of
freedom: 207


Restraints Information
Number of Intramolecular
restraints: 0
Number of Intermolecular
restraints: 0
********************** begin of velocity initialization
***********************
Initial Temperature
1500.00 K
Centre of mass velocity in direction x:
0.000000000000
Centre of mass velocity in direction y:
0.000000000000
Centre of mass velocity in direction z:
0.000000000000
*********************** end of velocity initialization
************************


Number of electrons: 198
Number of occupied orbitals: 99
Number of orbital functions: 178

Number of independent orbital functions: 178

Extrapolation method: initial_guess


SCF WAVEFUNCTION OPTIMIZATION

Step Update method Time Convergence
Total energy

-----------------------------------------------------------------------------
1 Mixing/Diag. 0.40E+00 0.06 0.8217339710
-193.7823266620
2 Mixing/Diag. 0.40E+00 0.07 0.5394381165
-205.4804849700
3 Mixing/Diag. 0.40E+00 0.06 0.3790033200
-214.8393334760
4 Mixing/Diag. 0.40E+00 0.06 0.2639343339
-221.5026105018
5 Mixing/Diag. 0.40E+00 0.06 0.1968783765
-225.9578659528
6 Mixing/Diag. 0.40E+00 0.06 0.1475535176
-228.8329037865
7 Mixing/Diag. 0.40E+00 0.06 0.1089769279
-230.6497393989
8 Mixing/Diag. 0.40E+00 0.06 0.0799099309
-231.7833560627
9 DIIS/Diag. 0.77E-02 0.06 0.0996616672
-232.4851722626
10 DIIS/Diag. 0.13E-02 0.06 0.0042508339
-233.6056702709
11 DIIS/Diag. 0.19E-02 0.06 0.0050778455
-233.6054467721
12 DIIS/Diag. 0.35E-02 0.06 0.0077720309
-233.6050396940
13 DIIS/Diag. 0.78E-03 0.06 0.0019990584
-233.6057751594
14 DIIS/Diag. 0.38E-03 0.06 0.0015106555
-233.6058004225
15 DIIS/Diag. 0.88E-04 0.06 0.0003069758
-233.6058067479
16 DIIS/Diag. 0.38E-04 0.06 0.0002154762
-233.6058071041
17 DIIS/Diag. 0.17E-04 0.06 0.0000801785
-233.6058072187
18 DIIS/Diag. 0.30E-05 0.06 0.0000143691
-233.6058072374
19 DIIS/Diag. 0.13E-05 0.06 0.0000038178
-233.6058072381

*** SCF run converged in 19 steps ***


Total electronic density (r-space): 0.0000000000
198.0000000000
Total core charge density (r-space): 0.0000000000
-198.0000000000
Total charge density (r-space):
0.0000000000
Total charge density (g-space):
0.0000000000

Core-core repulsion energy [eV]:
79331.42453161157027
Core Hamiltonian energy [eV]:
-9973.29763062268285
Two-electron integral energy [eV]:
-151429.72820039524231
Electronic energy [eV]:
-85688.16173082029854

Total energy [eV]:
-6356.73719920873555

Atomic reference energy [eV]:
6354.13592960308870
Heat of formation [kcal/mol]:
-59.98670797360201

MD_ENERGIES| Initialization proceeding


******************************** GO CP2K GO!
**********************************
INITIAL POTENTIAL ENERGY[hartree] =
-0.233605807238E+03
INITIAL KINETIC ENERGY[hartree] =
0.491648061274E+00
INITIAL TEMPERATURE[K]
= 1500.000
INITIAL VOLUME[bohr^3] =
0.105442728034E+06
INITIAL CELL LNTHS[bohr] = 0.4724315E+02 0.4724315E+02
0.4724315E+02
INITIAL CELL ANGLS[deg] = 0.9000000E+02 0.9000000E+02
0.9000000E+02
******************************** GO CP2K GO!
**********************************

Number of electrons: 198
Number of occupied orbitals: 99
Number of orbital functions: 178

Number of independent orbital functions: 178

Extrapolation method: previous_wf



For blyp we used the following input:
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME /home/didi/GTH_BASIS_SETS
POTENTIAL_FILE_NAME /home/didi/GTH_POTENTIALS
&MGRID
CUTOFF 320
&END MGRID
&SCF
EPS_SCF 5.0E-7
SCF_GUESS RESTART
MIXING 0.4
MAX_SCF 40
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
ENERGY_GAP 0.001
&END OT
&END SCF
&QS
METHOD GPW
EXTRAPOLATION PS
EXTRAPOLATION_ORDER 3
&END QS
&XC
&XC_FUNCTIONAL BLYP
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 20.0 20.0 20.0
UNIT ANGSTROM
&END CELL
&COORD
N -0.0660564363 -1.9984035901 -2.5831179889
...
...
H -4.0615750468 0.6446391816 3.2422118808
&END COORD
&KIND C
BASIS_SET TZVP-GTH
POTENTIAL GTH-BLYP-q4
&END KIND
&KIND O
BASIS_SET TZVP-GTH
POTENTIAL GTH-BLYP-q6
&END KIND
&KIND H
BASIS_SET TZVP-GTH
POTENTIAL GTH-BLYP-q1
&END KIND
&KIND N
BASIS_SET TZVP-GTH
POTENTIAL GTH-BLYP-q5
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PREFERRED_FFT_LIBRARY FFTESSL
EXTENDED_FFT_LENGTHS T
PROJECT md
RUN_TYPE MD
PRINT_LEVEL LOW
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVT
STEPS 5000
TIMESTEP 0.8
TEMPERATURE 2000.0
&THERMOSTAT
&NOSE
LENGTH 3
TIMECON 100.0
&END NOSE
&END THERMOSTAT
&END MD
&END MOTION

[Teodoro Laino]

Dear Diederica,

Let's go step by step:

(1) First and most important.. there should not be differences
between old and new versions of cp2k. If this is the case it may
easily be
that we introduced a bug during the changes. I would highly
appreciate if you could send us the whole input file including geometry.
In this way we can identify the problem and fix it.

(2) The drift that you observe depends very probably on the semi-
empirical (Ben can comment on that.. we spent a couple of
days trying to identify something similar but for NVE runs). Are you
doing a condensed phase? or an isolated molecule?
How large is your molecule compared to the box?
Most important, if you switch to NVE do you also observe the drift?

(3) Nose has been tested extensively and it conserves the energy. So
must necessarily be the quality of your forces.

For the next time, can you please put tables and input files as
attachment (maybe .tgz) so that the formatting of mail readers are
not changed and the input files not digested ? Thanks..

Ciao,
teo

[Diederica Claeys]

Hello!

 

On Mon, Mar 10, 2008 at 8:53 PM, Diederica <diederic...@gmail.com> wrote:

---------- Forwarded message ----------
From: Teodoro Laino <teodoro.la...@gmail.com>
Date: 10 mrt, 19:03
Subject: nose thermostat
To: cp2k


Dear Diederica,

Let's go step by step:

(1) First and most important.. there should not be differences  
between old and new versions of cp2k. If this is the case it may  
easily be
that we introduced a bug during the changes. I would highly  
appreciate if you could send us the whole input file including
geometry.
In this way we can identify the problem and fix it.

 

These are added, as an attachment this time

 

(2) The drift that you observe depends very probably on the semi-
empirical (Ben can comment on that.. we spent a couple of
days trying to identify something similar but for NVE runs). Are you  
doing a condensed phase? or an isolated molecule?
How large is your molecule compared to the box?
Most important, if you switch to NVE do you also observe the drift?

 

I'm modelling an isolated molecule, it should fit in the box. Using NVE without temperature rescaling, with an analogue structure, an energy drift was observed, that got worse after about 20 ps, though no non-convergence was observed. The input file and an energy-plot are added.

 

(3) Nose has been tested extensively and it conserves the energy. So  
must necessarily be the quality of your forces.

 

The potential and kinetic energy are conserved indeed, but the conserved quantity is not

 

For the next time, can you please put tables and input files as  
attachment (maybe .tgz) so that the formatting of mail readers are  
not changed and the input files not digested ? Thanks..

Ciao,
teo

 

Thanks for helping!

 

Greetz,

Diederica

 

[Teodoro Laino]

Ciao Diederica,

>
> These are added, as an attachment this time

Thanks.. will look into this problem ASAP.

> I'm modelling an isolated molecule, it should fit in the box. Using
> NVE without temperature rescaling, with an analogue structure, an
> energy drift was observed, that got worse after about 20 ps, though
> no non-convergence was observed. The input file and an energy-plot
> are added.

If you have a drift with NVE this means that energy and forces are
not consistent.. or (what I assume) that the cutoff for the coulomb
interactions is not large enough..

Please put all of the RC_COULOMB and RC_INTERACTION and RC_RANGE to
50.0.
You've the problem that in your cell the standard range (24 bohr)
does not fully cover the extension of your molecule (more than 26 bohr).
Use a larger box (30.0 Angstrom) and set the above keywords to 50.0
bohr.
You should not see any jump in the energy (or drift)..

> The potential and kinetic energy are conserved indeed, but the
> conserved quantity is not

Sorry to disprove you, but it is not possible that in an NVT you have
kinetic and potential conserved, since the thermostat
is continuously exchanging energy with the system. I't s a matter of
concept.. we may discuss about the magnitude of the
thermostat energy.. but that's another story..

Cheers
teo

>
>
>
> For the next time, can you please put tables and input files as
> attachment (maybe .tgz) so that the formatting of mail readers are
> not changed and the input files not digested ? Thanks..
>
> Ciao,
> teo
>
> Thanks for helping!
>
> Greetz,
> Diederica
>
>
>
>
> >

> <mailcp2k.tgz>

[Nichols A. Romero]

To elaborate on what Teo said, (he can correct me if I am wrong)

In NVT, the total energy of the *real sytem + thermostat* should be conserved.
(I don't remember if CP2K outputs this quantity, but most codes should).

The total energy of just the real system (atoms,molecules,etc.) should have a mean value (after equilibration) with fluctuations proportional to 1/sqrt(N) (or something like that).

If the total energy of the real system increases or decreases without bounds, then there is something wrong.
Or maybe you are near some phase transition.

--
Nichols A. Romero, Ph.D.
DoD User Productivity Enhancement and Technology Transfer (PET) Group
High Performance Technologies, Inc.
Reston, VA
443-567-8328 (C)
410-278-2692 (O)

[Teodoro Laino]

Ciao Nick,

On 10 Mar 2008, at 22:28, Nichols A. Romero wrote:

To elaborate on what Teo said, (he can correct me if I am wrong)

In NVT, the total energy of the *real sytem + thermostat* should be conserved.
(I don't remember if CP2K outputs this quantity, but most codes should).

Exactly.. and this quantity is what we call CONSERVED QNTY or you can find the same number in the 6th column of the .ener file.

The total energy of just the real system (atoms,molecules,etc.) should have a mean value (after equilibration) with fluctuations proportional to 1/sqrt(N) (or something like that).

If the total energy of the real system increases or decreases without bounds, then there is something wrong.
Or maybe you are near some phase transition.

Absolutely correct. In this case I think that the forces computed have something wrong.. 

Let me explain better once for all. The SCF module is still not fully completed. The most important part missing is 

the ewald for the coulomb part (we have a project going on for that, so it should be hopefully fixed one day ;-) ).

This means that at the moment the coulomb is evaluated through a neighbor lists. 

This creates ENORMOUS problem if you condensed phases (same problem of using a cutoff to truncate the electrostatic interactions 

in force field based calculations) but also if you do big molecules (and this is her case) and don't set properly

the parameter I wan mentioning in my last e-mail.

What happens here is that during the MD you've interactions turned on and off at the boundary and the contribution is such that

you observe a drift. 

Of course I could be wrong.. but since the nose-hoover thermostat has been so extensively tested i would rarely look for a bug in it ;-)

Teo

[Diederica Claeys]

hi

 

thanks a lot for all the explanations! The MD is running well now

 

greetz,

Diederica

[Teodoro Laino]

Dear Diederica,

I wanted to have a look at the problem you mentioned (i.e. different energies/forces for two different versions).

I downaloaded a cp2k version dated 2007-04-27 and this does not run your example. It gives a segmentation fault.

Can you please post the line of the output file containing the "Last CVS Entry" .. I want just to be sure that it was really 27 April 2007.

(From what I reproduced it cannot be..)

Thanks

teo

[Diederica Claeys]

Dear Teo,

 

Sorry for giving the wrong date. I didn't install anything myself...

Unfortunately there is nothing on that line, I added the head of the output and the CVS/Entries file I found in the folder with the codes.

 

greetings,

Diederica

[Teodoro Laino]

No problem Diederica,

Did you check if with the changes I told you to do (RC_*) a today cvs snapshot works or not.

Sorry but that CVS/Entries is useless.. I need the cp2k/src/CVS/Entries

The one you sent me looks like the cp2k/CVS/Entries.

Thanks

teo

<pm3_oldversion.tgz>

[Diederica Claeys]

Hello,

 

I just checked it and the molecule is still flying apart

 

Greetz,

[Teodoro Laino]

Hi Diederica,

it took a while but I got it.. there's a bug in the ANALYTICAL_DERIVATIVES..

No idea if the bug has always been there or if it was introduced later..

I'm working now to isolate the bug and fix it.. In the meantime you can use the new version (please do it! it

is very important in terms of bug reporting that people always uses the latest and the best cvs snapshot) setting 

SE%ANALYTICAL_DERIVATIVES  to FALSE

Will send here a message when the bug has been fixed.

Cheers

Teo

<Entries.tgz><testrecent.tgz>

[Teodoro Laino]

Ciao Diederica,

Just to let you know that the bug has been fixed. Thanks a lot for your help!

Please download a new version (the newest the best ;-) ) and no need to use the numerical derivatives.

Analytical should work great now.

Thanks again!

Teo

On 12 Mar 2008, at 10:12, Diederica Claeys wrote:

<Entries.tgz><testrecent.tgz>