[Pw_forum] How to adjust emass,dt and cutoff in cp.x calculation?

Wed Apr 12 21:30:25 CEST 2006

Hi,

Nicola and Sandro give good advice and I
just want to add to Sandro's list:

5) Check the forces to make sure that they are close
to the ground state forces. If they're not, reduce emass
further until they are. The difference between
the Car-Parrinello forces and the ground state
forces is roughly proportional to emass so you only have
to check the forces once...and then reduce emass to get
your desired level of accuracy.

6) If you reduce emass by a factor of X you reduce
the errors in the forces by a factor of X but you need
to reduce the time step by roughly a factor of sqrt(X)
so that you can still accurately integrate the
equations of motion of the orbitals.

Step 5 is necessary because, in theory, Car-Parrinello
dynamics does not give you the correct forces (even on average)
except in the limit emass --> 0.
Silicon was checked many years ago and it
gave good forces and I think the assumption was made that
you could therefore get good forces for all materials with
emasses of the same order of magnitude.
This is not true, and silicon is actually a very special case
because its valence wavefunctions are unusually smooth (you
can describe them using a very low plane-wave cutoff) and so
they don't change quickly as the ions move.
For some other systems that have been tested (oxides),
emasses that have been chosen by simply using the
prescription of Pastore et al. give you horrendous
forces.

By "horrendous" I mean: forces that have average errors
(with respect to the ground state forces) with a magnitude
that is only a factor of between 3 and 20 smaller than
the (r.m.s) magnitude of the forces themselves.
These errors are 2 or 3 orders of magnitude larger than what
you would get in most Born-Oppenheimer MD simulations.
Very few tests have been published so I'm just going
by those that have.

Most people either don't know about this issue or use the
ostrich approach to computer simulation and skip step 5 but
the tests that have been published suggest that this is
a mistake.

This issue is not discussed in review papers or in Richard
Martin's book which predate the work of Sandro and I:
See P. Tangney, S. Scandolo JCP 116, 14 (2002) and
P. Tangney JCP 124, 044111 (2006). You can download them
here: http://civet.berkeley.edu/tangney/JCPs.html

You could also just use Born-Oppenheimer MD, which is
always faster for a chosen level of accuracy on the forces
or the Kohn-Sham energy.
It doesn't conserve energy as well, but the energy that
is conserved in Car-Parrinello MD is a physically meaningless
quantity anyway, so who cares?
If temperature drifts too much, attach a weak thermostat.

Regards,

Paul

What Nicola says is entirely correct, and you should pay special
attention when studying metallic systems with CP (or seriously
considering using extensions of the CP method like the ones pioneered by
Nicola and by others).

However, should you still be willing to study Na with CP as a test case,
or should you consider studying "easier" (say, nonmetallic) systems in
the future, the recipe for setting emass, emass_cutoff, and dt is
roughly as follows:

1) set emass to a physical sound value. The value of emass determines
how much the fictitiuous dynamics of the electrons couples with the real
dynamics of the ions. A small emass guarantees decoupling. How small is
small is typically determined by the excitation gap E_gap of your
system. The minimum frequency of the fictitious electronic dynamics
scales like sqrt(E_gap/emass) [see, e.g., Pastore et al, PRA 44, 6334
(1991)], and this has to be much higher than the maximum frequency of
the ion dynamics. Typically a factor of three larger is enough. As you
see, metallic systems, where E_gap=0 are a bit of a nightmare for CP,
unless the finite size of your system introduces a finite gap, but then
you have to be extra-careful...

2) start playing with emass_ecutoff, by fixing it to a starting value,
and checking what is the maximum dt that alllows you to integrate the
equations of motion. The error you got is precisely a sign that the dt
you used is too large for the integrator to converge.

3) choose the value of emass_ecutoff that allows you to use the largest dt.

4) Set dt accordingly.

Points (2-4) are described in detail in Tassone, Mauri, Car, PRB 50,
10561 (1994).

Regards,
Sandro

Ruijuan Xiao wrote:

 >Dear all,
 >I met some problems when I do some calculations by cp.x. What I 
calculated is a cubic box in which 54 Na atoms exist. When I use 
dt=5.0d0, emass=400.0d0, and emass_cutoff=3.0d0, it can run without any 
problems. When I put dt=5.0d0,emass=400.0d0, then change the value of 
emass_cutoff, the calculation also runs without problems.
 >
 >But when I put emass=400.0d0 and emass_cutoff=3.0d0,then change dt 
into 10.0d0, the following error message appears and the program stopped.
 >-----------------------------------
 >%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 >     from ortho : error #        21
 >     max number of iterations exceeded
 > 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 >
 >     stopping ...
 >------------------------------------
 >
 >When I put dt=5.0d0 and emass_cutoff=3.0d0, but change emass into 
200.0d0, the program stopped after 12 iterations. The message is :
 >------------------------------------
 >  nfi    ekinc  temph  tempp        etot      enthal       econs 
  econt    vnhh   xnhh0    vnhp   xnhp0
 >   10 ********    0.0    0.0    41.82500    41.82500    41.82500 
146.23636  0.0000  0.0000  0.0000  0.0000
 >   11  NaN        0.0    0.0    58.16717    58.16717    58.16717  NaN 
         0.0000  0.0000  0.0000  0.0000
 >   12  NaN        0.0    0.0  NaN         NaN         NaN         NaN 
         0.0000  0.0000  0.0000  0.0000
 >
 >   MAIN:          EKINC   (thr)          DETOT   (thr)       MAXFORCE 
   (thr)
 >   MAIN:  NaN            0.1D-03 NaN            0.1D-08  0.000000D+00 
  0.1D+11
 >   MAIN: convergence achieved for system relaxation
 >
 >
 >              averaged quantities :
 >         ekinc          ekin          epot          etot     tempp
 > NaN           NaN           NaN           NaN                 0.0
 >-----------------------------------
 >It seems that the calculation is sensitive to these parameters. Since 
on the page 37 of the maunal, it says that "unless you are already 
experienced with the system you are studying or with the code internals, 
usually you need to tune some
 >input parameters, like emass, dt, and cut-offs.", so these parameters 
must be important for calculations, but now I am puzzled about how to 
adjust emass,dt and cutoff in cp calculations.Would you like to give me 
any information or any suggestion about that?
 >
 >Thank you very much.
 >
 >
 >Best regards,
 >
 >Sincerely,
 >Ruijuan Xiao
 >Institute of Physics,
 >Chinese Academy of Sciences
 >
 >
 >
 >_______________________________________________
 >Pw_forum mailing list
 >Pw_forum at pwscf.org
 >http://www.democritos.it/mailman/listinfo/pw_forum
 >
 >
 >

-- 
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
Dr. Paul Tangney
Theory of Nanostructured Materials Facility
The Molecular Foundry
Lawrence Berkeley National Lab.         E-mail: PTTangney at lbl.gov
1 Cyclotron Road, Bldg 66                  Phone: (510) 642-2635
Berkeley, CA 94720                         Fax :  (510) 643-9345
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo