[Pw_forum] cp.x temperature decreases

Mon Jul 21 21:57:53 CEST 2008

On Mon, 21 Jul 2008, Tim Teatro wrote:

TT> Hello All,

hi tim,

there may be one issue with your inputs. it looks as if 
you are using a norm-conserving pseudopotential, but your 
cutoffs (wavefunction/density) are the typical values 
for ultra-soft pseudopotentials. for norm-conserving psp
the wavefunction cutoff has to be larger (my guess would 
be 50-60 rydberg for carbon to be on the safe side), but 
there is no gain from using a larger density cutoff than
the default 4x wavefunction cutoff.

have you checked the convergence wrt. cutoffs?

TT> Thank you in advance for your time and any help you can offer.
TT> 
TT> We are running on HP/XC clusters, Opteron CPUs and pathscale
TT> compilers, MPI, etc. The clusters have version 3.2.3 of the
TT> Quantum-ESPRESSO code.

of course at this time, using version 4.0.1 is highly recommended
(many performance and scaling improvements and bugfixes).

TT> We are experiencing an issue with two MD simulations using cp.x in
TT> certain atomic systems. The two systems of note are a C(111) 2x1 slab
TT> type structure (12 layers with 96 atoms), and the other is a simple
TT> cubic 64 atom Si system. The problem is that over a 40 000 time step
TT> run, there is a very noticeable decrease in average temperature. Of
TT> course, the temperature oscillates over the run, but the C111 system
TT> appears to have decreasing temperature over time (from about 800K to
TT> about 600K for example). The temperature is not controlled and thermal
TT> motion is achieved by initial random displacement of atoms within the

if you don't control the temperature with a thermostat and
you only inject initial potential energy that means your system 
would be equilibrating over time. it would be better to use 
temperature rescaling until your system is equilibrated and then
switch to a thermostat to complete equilibration. this process
is frequently accelerated by using a "massive" nose-hoover scheme
(one thermostat per atom) which at least guarantees equipartitioning
of the kinetic energy throughout the system.

if you have only carbon atoms, you can also tweak the choice of
fictious mass and time step. a value of 400a.u. would be pretty
much at the limit of accuracy of the verlet integrator with a 5a.u.
time step. with only carbon atoms you may be able to crank this up
to 800-1000.0 a.u and perhaps ramp up the timestep a little as well.

OTOH, you have to monitor the size of your band gap and with higher
fictitious mass you run a higher risk of closing it => breakdown of
the cp method. i'm currently running a near metallic system where 
i'm limited to 300a.u. and 1.5a.u. and i actually have to start
with some tricks to slow down the ficitious dynamics after a quench
to the BO surface.

TT> run (using the tranp switch). At the end of this e-mail, I will attach
TT> the set of input files used in the run. It appears that the rate of
TT> energy loss is related to the energy of the system. A system which at
TT> a higher temperature appears to lose energy more quickly. For C111,
TT> over 40 000 time steps, this visually resembles a linear decrease, but
TT> as I mention in a moment, the Si system suggests that the relationship
TT> is exponential.
TT> 
TT> The Si system seems to lose temperature in a way that resembles
TT> (visually) an exponential decrease. The temperature eventually levels
TT> off after falling from some initial value. This only appears to happen
TT> at higher temperatures (I've only seen it in one of my colleagues
TT> runs, and the temperature was above 1500K)

again, that would be quite typical for an equilibration process,
since you are running without a thermostat. the question is how
equilibrated or how close to the XXXX kelvin structure are your
initial coordinates.

TT> Sixty-four atom GaAs, which has the same structure as Si, shows a very
TT> stable temperature.

hmmm... 
hard to tell. could be your choice of pseudopotential
or initial structure has some impact here.

in all cases, you should check whether your total
energy (EHAM) is conserved.

to do an efficient and correct CP dynamics always requires 
some tweaking. if your system is (near) metallic, it is an 
order of magnitude more complicated and requires you to be 
extremely careful and paranoid.

BO dynamics is usually more straightforward to use, but even
though you can use a "classical" MD time step, typically more 
cpu time consuming, since you need to converge wavefunctions 
very tightly to avoid losing (total) energy.

the latter can be countered to a large degree by using a
smart wavefunction extrapolation scheme (i've toyed with 
this recently in a different CP code and was for one specific
system able to do BO-dynamics twice as fast as CP. typically
BO-dynamics is a bit to about 3x slower depending on how
effective the extrapolation schemes work).

hope that helps,
    axel.

TT> Best regards,
TT> 
TT> 
TT> Timothy A.V Teatro
TT> 
TT> --
TT> Timothy A.V. Teatro <timothy.teatro at uoit.ca>
TT> University of Ontario Institute of Technology
TT> Faculty of Science
TT> 2000 Simcoe Street North
TT> Oshawa, ON, CANADA, L1H 7K4
TT> 
TT> 
TT> Input files for C(111) 2x1 (96 atom, 12 layer slab)_________________
TT> **************************************************************
TT> 
TT> &control
TT> title = 'c96 GS',
TT> calculation='cp',
TT> restart_mode='from_scratch',
TT> prefix='c96',
TT> outdir = '/scratch/timtro/c96/'
TT> ndr = 51,
TT> ndw = 51,
TT> dt = 5.0,
TT> nstep = 100,
TT> isave=100,
TT> tprnfor=.TRUE.,
TT> forc_conv_thr=1.0d-3,
TT> etot_conv_thr = 1.d-6,
TT> ekin_conv_thr = 1.d-5,
TT> wf_collect=.TRUE.,
TT> /
TT> &system
TT> ibrav = 8,
TT> celldm(1)=16.510372468,
TT> celldm(2)=0.577350269,
TT> celldm(3)=2.558782523,
TT> nat=96,
TT> ntyp=1,
TT> xc_type = 'BLYP',
TT> ecutwfc = 30.0,
TT> ecutrho = 300.0,
TT> /
TT> &electrons
TT> emass=400.d0,
TT> emass_cutoff=2.5d0,
TT> electron_dynamics='cg',
TT> /
TT> &ions
TT> ion_positions='default',
TT> ion_velocities='zero',
TT> ion_dynamics='none'
TT> /
TT> ATOMIC_SPECIES
TT>  C   12.0107  C.pz-vbc.UPF
TT> ATOMIC_POSITIONS (bohr)
TT> ...
TT> 
TT> 
TT> 
TT> &control
TT> title = 'c96 SD',
TT> calculation='cp',
TT> restart_mode='restart',
TT> prefix='c96',
TT> outdir = '/scratch/timtro/c96/'
TT> ndr = 51,
TT> ndw = 51,
TT> dt = 5.0,
TT> nstep = 500,
TT> isave=100,
TT> tprnfor=.TRUE.,
TT> forc_conv_thr=1.0d-6,
TT> etot_conv_thr = 1.d-6,
TT> ekin_conv_thr = 1.d-5,
TT> wf_collect=.TRUE.,
TT> /
TT> &system
TT> ibrav = 8,
TT> celldm(1)=16.510372468,
TT> celldm(2)=0.577350269,
TT> celldm(3)=2.558782523,
TT> nat=96,
TT> ntyp=1,
TT> xc_type = 'BLYP',
TT> ecutwfc = 30.0,
TT> ecutrho = 300.0,
TT> /
TT> &electrons
TT> emass=400.d0,
TT> emass_cutoff=2.5d0,
TT> electron_dynamics='sd',
TT> electron_velocities='zero',
TT> /
TT> &ions
TT> ion_positions='default',
TT> ion_velocities='zero',
TT> ion_dynamics='none'
TT> /
TT> ATOMIC_SPECIES
TT>  C   12.0107  C.pz-vbc.UPF
TT> ATOMIC_POSITIONS (bohr)
TT> ...
TT> 
TT> 
TT> &control
TT> title = 'c96 RAND',
TT> calculation='cp',
TT> restart_mode='restart',
TT> prefix='c96',
TT> outdir = '/scratch/timtro/c96/'
TT> ndr = 51,
TT> ndw = 51,
TT> dt = 5.0,
TT> nstep = 100,
TT> isave=100,
TT> tprnfor=.TRUE.,
TT> forc_conv_thr=1.0d-6,
TT> etot_conv_thr = 1.d-6,
TT> ekin_conv_thr = 1.d-5,
TT> wf_collect=.TRUE.,
TT> /
TT> &system
TT> ibrav = 8,
TT> celldm(1)=16.510372468,
TT> celldm(2)=0.577350269,
TT> celldm(3)=2.558782523,
TT> nat=96,
TT> ntyp=1,
TT> xc_type = 'BLYP',
TT> ecutwfc = 30.0,
TT> ecutrho = 300.0,
TT> /
TT> &electrons
TT> emass=400.d0,
TT> emass_cutoff=2.5d0,
TT> electron_dynamics='sd',
TT> electron_velocities='zero',
TT> /
TT> &ions
TT> ion_positions='default',
TT> ion_dynamics='none',
TT> tranp(1)=.TRUE.,
TT> amprp(1)=0.42898,
TT> /
TT> ATOMIC_SPECIES
TT>  C   12.0107  C.pz-vbc.UPF
TT> ATOMIC_POSITIONS (bohr)
TT> ...
TT> 
TT> 
TT> &control
TT> title = 'c96 MD',
TT> calculation='cp',
TT> restart_mode='reset_counters',
TT> prefix='c96',
TT> outdir = '/scratch/timtro/c96/'
TT> ndr = 51,
TT> ndw = 51,
TT> dt = 5.0,
TT> nstep = 40000,
TT> isave=200,
TT> tprnfor=.TRUE.,
TT> forc_conv_thr=1.0d-6,
TT> etot_conv_thr = 1.d-6,
TT> ekin_conv_thr = 1.d-5,
TT> wf_collect=.TRUE.,
TT> /
TT> &system
TT> ibrav = 8,
TT> celldm(1)=16.510372468,
TT> celldm(2)=0.577350269,
TT> celldm(3)=2.558782523,
TT> nat=96,
TT> ntyp=1,
TT> xc_type = 'BLYP',
TT> ecutwfc = 30.0,
TT> ecutrho = 300.0,
TT> /
TT> &electrons
TT> emass=400.d0,
TT> emass_cutoff=2.5d0,
TT> electron_dynamics='verlet',
TT> electron_velocities='zero',
TT> /
TT> &ions
TT> ion_positions='default',
TT> ion_dynamics='verlet',
TT> ion_velocities='zero',
TT> /
TT> ATOMIC_SPECIES
TT>  C   12.0107  C.pz-vbc.UPF
TT> ATOMIC_POSITIONS (bohr)
TT> _______________________________________________
TT> Pw_forum mailing list
TT> Pw_forum at pwscf.org
TT> http://www.democritos.it/mailman/listinfo/pw_forum
TT> 

-- 
=======================================================================
Axel Kohlmeyer   akohlmey at cmm.chem.upenn.edu   http://www.cmm.upenn.edu
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.