nose thermostat
Diederica
diederic... at gmail.com
Mon Mar 10 17:54:15 UTC 2008
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
More information about the CP2K-user
mailing list