[CPMD-list] problems with MD -- conservation of energy and translational drift
Axel Kohlmeyer
axel.kohlmeyer at theochem.ruhr-uni-bochum.de
Tue Feb 15 11:49:06 CET 2005
On Tue, 15 Feb 2005 jernej at cmm.ki.si wrote:
JS> Dear CPMD community,
dear jernej,
JS>
JS> I wonder if somebody can help me with the following two issues regarding CP
JS> molecular dynamics:
yes, i may have one or two suggestions. ;-)
JS> 1) Conservation of EHAM: when starting CPMD simulation from a previously
JS> minimized structure, EHAM is by no means constant and features a significant
JS> drift. HOWEVER, if I kill the dynamics after a few (say 10) steps and restart
JS> it with the 'RESTART WAVEFUNCTION COORDINATES VELOCITIES NOSE ACCUMULATORS
JS> LATEST' and 'QUENCH BO' options, EHAM does not drift anymore.
when you start a CP-MD without velocities, then it'll will take a little
time until the electronic degrees of freedom have picked up some
kinetic energy and are in equilibrium (20-100 steps). if you have
a very large timestep (like your 5 a.u.) chances are, that the ions
move to fast to stay close enough to the BO-surface for a 'good'
CP-dynamics. with the QUENCH BO you re-optimize the wavefunction.
please note, that to restart the nose-hoover chains for the ions,
you have to use RESTART NOSEP (RESTART NOSE will have no effect,
see the code in control.F).
JS> This weird feature does not make any sense to me. It appears with either of the
JS> three popular pseudopotentials I've used so far (i.e., Goedecker,
JS> Trouiller-Martins, and Vanderbilt). I've been using the BLYP functional
JS> throughout. Can somebody tell me if I did something wrong in the simulation
JS> setup, or otherwise explain this feature? Is it possible that this problem has
JS> anything to do with the way in which the code was compiled (I had quite a few
JS> problems with this)?
no, this should not be a compilation problem. what
platform/compiler/libraries do you use and what were
your problems?
JS> 2) Translational drift: even after "solving" the problem with EHAM, as depicted
JS> above, it happens that after running the MD for quite a while, say 20.000 steps
JS> of 5 a.u., the DIS value that was reasonably small (up to 1.) and stable within
JS> the first 20.000 steps, starts increasing considerably -- it reaches a value of
JS> more than 25. after 60.000 steps! Then the whole MD chrases -- it produces
JS> meaningless characters in the output:
JS>
[...]
JS> 1.91
JS> 60610 0.00675 324.7 -221.25750 -221.23313 -221.22638 26.436
JS> 1.97
JS> 60611+++++++++ 324.2 -221.25749 -221.23314++++++++++++++ 26.438
JS> 1.96
JS> 60612??????????????????????????????????????????????????????????? 26.440
this happens, when your wavefunction 'diverges', e.g. when an ion moves
too fast for the electronic degrees of freedom to be able to follow.
JS> I've checked the center of mass of the system and indeed it features a drift --
JS> therefore there is considerable translation in the simulation, especially after
JS> the DIS value had increased. I'm a bit puzzled about the origin of this drift,
the drift is a consequence of discretization and rounding errors
(you are solving the equations of motions numerically for finite
time steps and you keep your electron structure on a grid),
like it happens for all (long-running) MD simulations.
to counter this effect, you can add to your &CPMD section:
SUBTRACT COMVEL
1000
which will cancel the center of mass velocity every 1000 steps (the energy
is conserved by distributing the removed COM kinetic energy to the other
DOFs). please note, that this code does not (yet) take nose-hoover
thermostats into account, so you'll see a (small) discontinuity in the
total energy, if you use it in your example.
JS> for I've made sure that I started my MD from a previously minimized structure!
JS> Can somebody tell me what is the source of this problem and how to avoid it?
JS> Could the drift have anything to do with the weird conservation of EHAM
JS> depicted in 1) ?
i have a few more comments regarding your input file.
JS> **** Below is the input file that I've been using for my simulations.
JS>
JS> &CPMD
JS> MOLECULAR DYNAMICS CP
JS> RESTART WAVEFUNCTION COORDINATES VELOCITIES NOSE ACCUMULATORS LATEST
as stated above, you need NOSEP instead of NOSE here.
JS> TIMESTEP
JS> 5.
a timestep of 5.0 a.u. is pretty large for a system
containing (light and thus fast moving) hydrogen atoms
JS> EMASS
JS> 600.
my guess is you need to increase EMASS to 600 to get the
CP-MD working at all. right?
JS> NOSE IONS
JS> 300 3200
JS> QUENCH BO
you should not do a QUENCH BO regularly, when restarting a
trajectory. this will disrupt the dynamics. unless you also
do a QUENCH ELECTRONS, your electronic DOFs will retain their
old velocities, so you won't gain a lot from it here.
if you are concerned about 'loosing' your electron structure
you should better use a thermostat for the electrons
as well. ideally, you'll thermostat to the average value
of EKINC from a non-thermostatted run, e.g. something like
0.0065 in your example. of course you run into the risk of
'thermostatting away' the real problem, which to me seems
mainly a too large timestep.
JS> MAXSTEP
JS> 500000
JS> STORE
JS> 50
JS> &END
JS> &DFT
JS> FUNCTIONAL BLYP
with gradient corrected functionals, especially when using
a small cutoff (and even more if you run an MD of an
isolated molecule) you should use
GC-CUTOFF
here (e.g. with a value of 1.0e-6 in your case) for a much
improved energy conservation. the BLYP functional due to its
functional form seems very sensitive to this.
JS> &END
JS> &SYSTEM
JS> ANGSTROM
JS> SYMMETRY
JS> TRICLINIC
JS> CELL ABSOLUTE DEGREE
JS> 6.1026 3.4867 11.954 90.0 105.79 90.0
hmmm. this is quite a small cell. did you do any
checks, whether your k-point sampling is sufficient
with a gamma-point calculation?
this should not affect your energy conservation
issues, though.
JS> CUTOFF
JS> 70.0
a plane wave cutoff of 70.0 is very much pushing
the limit for your system, if you use the standard
MT pseudopotentials. a very nice discussion of
what are reasonable parameters for CP-MD is in
Kuo et al., J. Phys. Chem. B 2004, 108, 12990-12998.
JS> POINT GROUP
JS> AUTO
JS> &END
JS> &ATOMS
JS> *O_MT_BLYP KLEINMAN-BYLANDER
JS> LMAX=D
unless you use custom pseudopotentials, you
should most likely use LMAX=P here.
[...]
JS> 2.895549 2.593059 5.334355
JS> 2.763486 2.207081 2.052808
JS> *C_MT_BLYP KLEINMAN-BYLANDER
JS> LMAX=D
same as above.
[...]
JS> *H_MT_BLYP KLEINMAN-BYLANDER
JS> LMAX=S
here KLEINMANN-BYLANDER is not really
needed, but should do nor harm either.
best regards,
axel.
--
=======================================================================
Dr. Axel Kohlmeyer e-mail: axel.kohlmeyer at theochem.ruhr-uni-bochum.de
Lehrstuhl fuer Theoretische Chemie Phone: ++49 (0)234/32-26673
Ruhr-Universitaet Bochum - NC 03/53 Fax: ++49 (0)234/32-14045
D-44780 Bochum http://www.theochem.ruhr-uni-bochum.de/~axel.kohlmeyer/
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.
More information about the CPMD-list
mailing list