[CP2K:760] nose thermostat
Teodoro Laino
teodor... at gmail.com
Mon Mar 10 18:03:14 UTC 2008
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
On 10 Mar 2008, at 18:54, Diederica wrote:
>
> 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