[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