scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

CP user-guide update

This commit is contained in:
Riccardo Bertossa 2023-03-14 11:42:50 +00:00 committed by giannozz
parent ae67887bb2
commit 67e6882a86
5 changed files with 262 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
CPV/Doc/Makefile
View File

@ -9,7 +9,7 @@ AUXS = $(PDFS:.pdf=.aux)
LOGS = $(PDFS:.pdf=.log)
OUTS = $(PDFS:.pdf=.out)
TOCS = $(PDFS:.pdf=.toc)
USER_GUIDES = user_guide.pdf user_guide.html
USER_GUIDES = user_guide.pdf user_guide.html autopilot_guide.pdf autopilot_guide.html
PANDOC_PDF_ENGINE_CMD = --latex-engine
PANDOC_PDF_ENGINE_CMD = --pdf-engine
PANDOC_ADDITIONAL_DEF = --metadata csquotes=yes
@ -20,10 +20,16 @@ pdf: $(PDFS)
user_guide: $(USER_GUIDES)
user_guide.pdf: user_guide.md pandoc_template.tex tp.tex
$(shell sed -e '/[[_TOC_]]/d' -e 's/\$$`/$$/g' -e 's/`\$$/$$/g' $< | $(PANDOC) -N --toc --metadata title='CP user guide' --metadata customtitlepage='tp' --metadata graphics --metadata geometry='top=3cm' --metadata geometry='left=3cm' --metadata geometry='right=3cm' --metadata geometry='bottom=3cm' --metadata hyperrefoptions='linktoc=all' ${PANDOC_PDF_ENGINE_CMD}=xelatex ${PANDOC_ADDITIONAL_DEF} -o $@ --template $(word 2,$^) )
$(shell sed -e '/[[_TOC_]]/d' -e 's/\$$`/$$/g' -e 's/`\$$/$$/g' $< | $(PANDOC) -N --toc --metadata title="CP's user guide" --metadata customtitlepage='tp' --metadata graphics --metadata geometry='top=3cm' --metadata geometry='left=3cm' --metadata geometry='right=3cm' --metadata geometry='bottom=3cm' --metadata hyperrefoptions='linktoc=all' ${PANDOC_PDF_ENGINE_CMD}=xelatex ${PANDOC_ADDITIONAL_DEF} -o $@ --template $(word 2,$^) )
autopilot_guide.pdf: autopilot_guide.md pandoc_template.tex tp.tex
$(shell sed -e '/[[_TOC_]]/d' -e 's/\$$`/$$/g' -e 's/`\$$/$$/g' $< | $(PANDOC) -N --toc --metadata title="CP's AUTOPILOT user guide" --metadata customtitlepage='tp' --metadata graphics --metadata geometry='top=3cm' --metadata geometry='left=3cm' --metadata geometry='right=3cm' --metadata geometry='bottom=3cm' --metadata hyperrefoptions='linktoc=all' ${PANDOC_PDF_ENGINE_CMD}=xelatex ${PANDOC_ADDITIONAL_DEF} -o $@ --template $(word 2,$^) )
user_guide.html: user_guide.md
$(shell sed -e '/[[_TOC_]]/d' -e 's/\$$`/$$/g' -e 's/`\$$/$$/g' $< | $(PANDOC) -s --toc --metadata title='CP user guide' --mathml -o $@ )
$(shell sed -e '/[[_TOC_]]/d' -e 's/\$$`/$$/g' -e 's/`\$$/$$/g' $< | $(PANDOC) -s --toc --metadata title="CP's user guide" --mathml -o $@ )
autopilot_guide.html: autopilot_guide.md
$(shell sed -e '/[[_TOC_]]/d' -e 's/\$$`/$$/g' -e 's/`\$$/$$/g' $< | $(PANDOC) -s --toc --metadata title="CP's AUTOPILOT user guide" --mathml -o $@ )
clean:
- rm -f $(PDFS) $(AUXS) $(LOGS) $(OUTS) $(TOCS) $(USER_GUIDES) *~

114
CPV/Doc/autopilot_guide.md
View File

@ -17,11 +17,21 @@ This documentation, like the software it accompanies, is distributed in the hope
that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY FOR A PARTICULAR PURPOSE.
\
\
\
\
\
\
\
\
[[_TOC_]]
--------------------------------------------------------------------------------
AUTOPILOT DOCUMENTATION
--------------------------------------------------------------------------------
# Autopilot Documentation
The Autopilot Feature Suite is a user level enhancement for directing
Car-Parrinello simulations based on CP.X packaged in ESPRESSO.
@ -33,9 +43,7 @@ The following features are incorporated:
- Autopilot Course Correction (Steering)
--------------------------------------------------------------------------------
Auto Restart Mode
--------------------------------------------------------------------------------
# Auto Restart Mode
Auto Restart Mode is an extension of restart_mode declared in the CONTROL section
of the input file. When restart mode is set to "auto", control determines if the
@ -71,9 +79,7 @@ placing more control with the user.
--------------------------------------------------------------------------------
Autopilot Course Configuration (Dynamic Rules)
--------------------------------------------------------------------------------
# Autopilot Course Configuration (Dynamic Rules)
Autopilot Course Configuration (Dynamic Rules) is a method that allows select
input parameters (Autopilot variables) to change during the course of a
@ -169,9 +175,7 @@ runtime error may occur.
--------------------------------------------------------------------------------
Autopilot Course Correction (Steering)
--------------------------------------------------------------------------------
# Autopilot Course Correction (Steering)
Autopilot Course Correction (Steering) provides a run-time method of changing
Autopilot variables on the fly, after the simulation is underway. Autopilot
@ -282,49 +286,63 @@ will pause.
# List of autopilot variables
--------------------------------------------------------------------------------
TESTING
--------------------------------------------------------------------------------
The Autopilot module supports the following input variables:
The entire Autopilot Feature Suite issues directives to slave nodes under
MPI, with the fewest broadcasted parameters. All features have been tested
under Intel 8.1 with MKL 7.0.1 libraries on a Linux-32 single processor and under
PGI 5.2 with MPI on Linux-64 with 1, 2 and 4 processors.
- isave
- iprint
- tprint
- dt (automatically set also tolp to old dt and electron_velocities/ions_velocities to 'change_step')
- electron_dynamics (if 'cg', this has little side-effects)
- electron_damping
- ion_dynamics
- orthogonalization
- ion_damping
- ion_temperature
- tempw
- nhpcl
- fnosep
--------------------------------------------------------------------------------
ADDING AN AUTOPILOT VARIABLE
--------------------------------------------------------------------------------
# How to add support for a new autopilot variable
See `autopilot.f90` for examples.
- Select the input parameter from the list in file INPUT_CP
- Identify parameter dependencies, initializations, assignments, etc
- Edit autopilot.f90 to add the following,
where VARNAME is the name of the new Autopilot variable:
VARTYPE :: rule_VARNAME(max_event_step) at module scope
LOGICAL :: event_VARNAME(max_event_step) at module scope
* Remember to add to the PUBLIC block as well
event_VARNAME(:) = .false. to init_autopilot subroutine
rule_VARNAME(:) = VARDEFAULT to init_autopilot subroutine
* Import VARNAME with USE to employ_rules subroutine
* In employ_rules, add conditional clause on event_VARNAME to assign VARNAME:
! VARNAME
if (event_VARNAME(event_index)) then
VARNAME = rule_VARNAME(event_index)
CALL init_other_VARNAME_dependent_variables( VARNAME)
write(*,*) 'RULE EVENT: VARNAME', VARNAME
endif
* Import VARNAME with USE to assign_rule subroutine
* In assign_rule, add condition clause matching the VARNAME create rule as so:
ELSEIF ( matches( "VARNAME", var ) ) THEN
read(value, *) VARTYPE_value
rule_VARNAME(event) = VARTYPE_value
event_VARNAME(event) = .true.
* TEST
- Select the input parameter from the list in file INPUT_CP
- Identify parameter dependencies, initializations, assignments, etc
- Edit autopilot.f90 to add the following,
where VARNAME is the name of the new Autopilot variable:
VARTYPE :: rule_VARNAME(max_event_step) !at module scope
LOGICAL :: event_VARNAME(max_event_step) !at module scope
* Remember to add to the PUBLIC block as well
event_VARNAME(:) = .false. !to init_autopilot subroutine
rule_VARNAME(:) = VARDEFAULT !to init_autopilot subroutine
* Import VARNAME with USE to employ_rules subroutine
* In employ_rules, add conditional clause on event_VARNAME to assign VARNAME:
! VARNAME
if (event_VARNAME(event_index)) then
VARNAME = rule_VARNAME(event_index)
CALL init_other_VARNAME_dependent_variables( VARNAME)
write(*,*) 'RULE EVENT: VARNAME', VARNAME
endif
* Import VARNAME with USE to assign_rule subroutine
* In assign_rule, add condition clause matching the VARNAME create rule as so:
ELSEIF ( matches( "VARNAME", var ) ) THEN
read(value, *) VARTYPE_value
rule_VARNAME(event) = VARTYPE_value
event_VARNAME(event) = .true.
* TEST
WARNING: Some Autopilot variables may create "side-effects". For example, the
inclusion of a rule for TEMPW rules invokes a side-effect call to ions_nose_init.

2
CPV/Doc/pandoc_template.tex
View File

@ -450,6 +450,7 @@ $if(beamer)$
\frame{\titlepage}
$else$
$if(customtitlepage)$
\def\intitle{$title$}
\begin{titlepage}
\input{$customtitlepage$}
\end{titlepage}
@ -486,6 +487,7 @@ $if(colorlinks)$
$endif$
\setcounter{tocdepth}{$toc-depth$}
\tableofcontents
\newpage
}
$endif$
$endif$

2
CPV/Doc/tp.tex
View File

@ -7,5 +7,5 @@
% title
\vspace{5.5cm}
\Huge CP User's Guide (v.\version)
\Huge \intitle\ (v.\version)
\end{center}

228
CPV/Doc/user_guide.md
View File

@ -9,7 +9,7 @@ beyond what is provided in this guide, can be found in the directory
`CPV/Doc/`, containing a copy of this guide.
This guide assumes that you know the physics that `CP` describes and the
methods it implements. It also assumes that you have already installed,
methods it implements. A good reference for the topic is the book [D, Marx and J. Hutter, "Ab Initio Molecular Dynamics"](https://doi.org/10.1017/CBO9780511609633) It also assumes that you have already installed,
or know how to install, Quantum ESPRESSO. If not, please read the
general User's Guide for Quantum ESPRESSO, found in directory `Doc/` two
levels above the one containing this guide; or consult the web site:
@ -18,7 +18,7 @@ levels above the one containing this guide; or consult the web site:
People who want to modify or contribute to `CP` should read the
Developer Manual: [gitlab.com/QEF/q-e/-/wikis/home](https://gitlab.com/QEF/q-e/-/wikis/home).
`CP` can perform Car-Parrinello molecular dynamics, including
`CP` can perform [Car-Parrinello molecular dynamics](https://doi.org/10.1103/PhysRevLett.55.2471), including
variable-cell dynamics. The `CP` package is based on the original code
written by Roberto Car
and Michele Parrinello. `CP` was developed by Alfredo Pasquarello (EPF
@ -31,12 +31,12 @@ and others. We quote in particular:
- Sergio Orlandini (CINECA) for completing the CUDA Fortran acceleration
started by Carlo Cavazzoni
- Fabio Affinito and Maruella Ippolito (CINECA) for testing and benchmarking
- Fabio Affinito and Mariella Ippolito (CINECA) for testing and benchmarking
- Ivan Carnimeo and Pietro Delugas (SISSA) for further openACC acceleration
- Riccardo Bertossa (SISSA) for extensive refactoring of ensemble dynamics /
conjugate gradient part
conjugate gradient part, contributions to the documentation
- Federico Grasselli and Riccardo Bertossa (SISSA) for bug fixes,
extensions to Autopilot;
@ -181,6 +181,7 @@ keywords with self-explanatory names:
> ATOMIC\_POSITIONS\
> CELL\_PARAMETERS (optional)\
> OCCUPATIONS (optional)
> AUTOPILOT (optional)
The keywords may be followed on the same line by an option. Unknown
fields are ignored. See the files mentioned above for details on the
@ -246,14 +247,16 @@ The `prefix.for` file, formatted like the previous two, contains the computed
forces, in Hartree atomic units as well. It is written only if a molecular
dynamics calculation is performed, or if `tprnfor = .true.` is set in input.
In the `prefix.str` file you can find the stress tensor (with the ionic kinetic part included), in GPa.
The simulation cell is written in a file named `prefix.cel` with the same header as the previous
described files, and the cell matrix is then listed. NB: **THE CELL MATRIX IN THE
OUTPUT IS TRANSPOSED** that means that if you want to reuse it again for a new input file,
you have to pick the one that you find in `prefix.cel` and write in the input file
after inverting rows and columns.
after inverting rows and columns. In the `prefix.cel` file you will find the cell vectors in the **columns** of the matrix.
The file `prefix.evp` has one line per printed step and contains some
thermodynamical data.
thermodynamics data.
The first line of the file names the columns:
```
# nfi time(ps) ekinc Tcell(K) Tion(K) etot enthal econs econt Volume Pressure(GPa)
@ -264,7 +267,7 @@ where:
- `etot` is the DFT (potential) energy of the system, $`E_{DFT}`$
- `econs` is a physically meaningful constant of motion, $`E_{DFT} + K_{NUCLEI}`$,
in the limit of zero electronic fictitious mass
- `econt` is the constant of motion of the lagrangian$`E_{DFT} + K_{IONS} + K_{ELECTRONS}`$ t.
- `econt` is the constant of motion of the Lagrangian$`E_{DFT} + K_{IONS} + K_{ELECTRONS}`$ t.
If the time step `dt` is small enough this will be up to a very good precision a constant.
It is not a physical quantity, since $`K_{ELECTRONS}`$ has _nothing_ to do with the quantum
kinetic energy of the electrons.
@ -282,6 +285,7 @@ To prepare and run a CP simulation you should first of all define the
system:
> atomic positions\
> atomic velocities (can be zero, read from the input, sampled from Maxwell-Boltzmann)\
> system cell\
> pseudopotentials\
> cut-offs\
@ -354,6 +358,8 @@ An example of input file (Benzene Molecule):
You can find the description of the input variables in file `Doc/INPUT_CP.*`.
The best way to initialize and run the simulation is first doing a few steps of Born-Oppenheimer (BO) molecular dynamics (MD) using the `electron_dynamic = 'cg'` input for minimizing the electronic degrees of freedom each MD step, after sampling the ionic velocities from a Maxwell-Boltzmann distribution. Then usually you run a Nosé thermostat and/or a barostat to reach some temperature and pressure that you choose. Then you can run in the microcanonical ensemble if needed. During the CP MD you must check that the fake electronic kinetic energy does not increase too much, otherwise at some point your simulation may become completely wrong.
Reaching the electronic ground state
------------------------------------
@ -361,9 +367,9 @@ The first run, when starting from scratch, is always an electronic
minimization, with fixed ions and cell, to bring the electronic system
on the ground state (GS) relative to the starting atomic configuration.
This step is conceptually very similar to self-consistency in a
`pw.x` run.
`pw.x` run. The suggested method is to use the `electron_dynamics = ’cg’` to initialize the simulation.
Sometimes a single run is not enough to reach the GS. In this case, you
Sometimes when you do not use the CG routine, a single run is not enough to reach the GS. In this case, you
need to re-run the electronic minimization stage. Use the input of the
first run, changing `restart_mode = ’from_scratch’` to
`restart_mode = ’restart’`.
@ -385,37 +391,66 @@ Usually we consider the system on the GS when `ekin_conv_thr`
$`< 10^{-5}`$. You could check the value of the fictitious kinetic energy
on the standard output (column EKINC).
Different strategies are available to minimize electrons, but the most
frequently used is _damped dynamics_: `electron_dynamics = ’damp’` and
`electron_damping` = a number typically ranging from 0.1 and 0.5.
See the input description to compute the optimal damping factor.
Steepest descent: `electron_dynamics = ’sd’`, is also available but it
is typicallyslower than damped dynamics and should be used only to
Different strategies are available to minimize electrons.
You can choose among the following:
- _conjugate gradient_ (CG): direct minimization of the energy functional `electron_dynamics = ’cg’`
- damped dynamics: `electron_dynamics = ’damp’` and `electron_damping` = a number typically ranging from 0.1 and 0.5.
- Steepest descent: `electron_dynamics = ’sd’`
The CG routine computes also the wavefunction time derivative in the parallel transport gauge.
The time derivative is performed by computing twice the ground state using the atomic positions at the current timestep,
and the atomic positions at the next timestep, that are computed by moving the atoms with their atomic velocities.
This procedure allows to start the dynamics with a wavefunction velocity consistent with the atomic velocities.
This usually results in a smoother start of the simulation.
An other advantage of this technique is that you can apply it also in the middle of the dynamics introducing
a much lower discontinuity in the trajectory than apply a minimization only and settings again to zero the wavefunction
velocity. All the other minimization routines in CP when called set the wavefunction velocity to zero.
At the moment the CG routine works only with the plane wave parallelization scheme, both on CPU and GPU machines. It can run "on the fly" by using the autopilot module.
See the input description to compute the optimal damping factor of the damped dynamics routine.
Steepest descent is also available but it is typically slower than damped dynamics and should be used only to
start the minimization.
Note that the CG routine is used also for doing a Born-Oppenheimer dynamic, or for using the ensemble DFT molecular dynamics. For this reason the meaning of the input variable
`nstep` is different. In the CG minimization `nstep` is the number of ions dynamic step: the conjugate gradient algorithm stops
when the required accuracy is reached, and the maximum number of steps that the CG algorithm performs is the input variable
`maxiter`. In the damped dynamics and the steepest descent the `nstep` input variable is the number of sd/damped dynamics steps performed
to reach the ground state, and ions usually are kept fixed.
Relax the system
----------------
If the input atomic positions are far from the equilibrium position you need to relax the system.
Note that you can consider to do that with the pw.x code that may be faster for this task.
Keep in mind that a recipe that works for every case does not exist, and it is impossible to predict
what will happen in an arbitrary system without running the simulation. The best approach is to try
and see what happens. Here we give some ideas on a possible approach to the simulation.
Once your system is in the GS, depending on how you have prepared the
starting atomic configuration:
1. if you have set the atomic positions \"by hand\" and/or from a
classical code, check the forces on atoms, and if they are large
($`\sim 0.1 \div 1.0`$ atomic units), you should perform an ionic
minimization, otherwise the system could break up during the
dynamics.
minimization or a few step of CG dynamic with a very small timestep
to equilibrate a little the simulation,
otherwise the system could break up during the CP dynamics.
2. if you have taken the positions from a previous run or a previous
ab-initio simulation, check the forces, and if they are too small
($`\sim 10^{-4}`$ atomic units), this means that atoms are already in
equilibrium positions and, even if left free, they will not move.
Then you need to randomize positions a little bit (see below).
You will need to set the atomic velocities from input or let the code
chose them by sampling a Maxwell-Boltzmann distribution at a given temperature.
An alternative approach is also to randomize a little bit the atomic positions.
Let us consider case 1). There are different strategies to relax the
system, but the most used are again steepest-descent or damped-dynamics
Let us consider case 1). Suppose that you decided to minimize the system with the cp.x code.
There are different strategies to relax the
system in cp.x, but the most used are again steepest-descent or damped-dynamics
for ions and electrons. You could also mix electronic and ionic
minimization scheme freely, i.e. ions in steepest-descent and electron
in with damped-dynamics or vice versa.
in with conjugate-gradient.
- suppose we want to perform steepest-descent for ions. Then we should
specify the following section for ions:
@ -459,9 +494,10 @@ in with damped-dynamics or vice versa.
electron and ions together, the ionic forces are not properly
correct, then it is often better to perform a ionic step every N
electronic steps, or to move ions only when electron are in their GS
(within the chosen threshold).
(within the chosen threshold). You can also use the CG dynamics for the electrons,
that take care of fully minimizing the electronic energy.
This can be specified by adding, in the ionic section, the
You can perform a ionic step every N electronic steps done with 'sd' or 'damp', by adding, in the ionic section, the
`ion_nstepe` parameter, then the &IONS namelist become as follows:
&ions
@ -483,10 +519,12 @@ in with damped-dynamics or vice versa.
process thus continues until the forces become smaller than
`forc_conv_thr`.
*Note* that to fully relax the system you need many runs, and
different strategies, that you should mix and change in order to
speed-up the convergence. The process is not automatic, but is
strongly based on experience, and trial and error.
*Note* that to fully relax the system
you need many runs, and different strategies, that you should mix and change in order to
speed-up the convergence. This process is not automatic, but is
strongly based on experience, and trial and error. For this reason
we suggest to use the CG minimization algorithm, that simplifies the
process and usually works well, at least for the electronic part of the problem.
Remember also that the convergence to the equilibrium positions
depends on the energy threshold for the electronic GS, in fact
@ -495,8 +533,65 @@ in with damped-dynamics or vice versa.
on forces could not be satisfied, if you do not require an even
smaller threshold on total energy.
Let us now move to case 2: randomization of positions.
A different approach is to use a small timestep
and run the ion dynamic using ground states computed with the CG routine.
The system, given a small enough timestep and a big enough number of steps,
will thermalize at a defined temperature. Note that if the potential energy of the
initial configuration is too big the final temperature may be too large.
To run some steps of Born-Oppenheimer (BO) molecular dynamics you can set the following input variables:
&control
! ...
nstep = 100, ! number of MD steps. It has to be big enought
dt = 1.0d0, ! put a small number here...
! ...
/
!...
&electrons
electron_dynamics = 'cg'
/
!...
&ions
ion_dynamics = 'verlet' ! you also mix different approaches by changing this to sd
/
!...
An additional feature of the code that you can use to perform the simulation initialization doing less restarts is the autopilot module ( see [the autopilot guide](autopilot_guide.md) ). For example you can play with the timestep _while the simulation is running_ by writing a file called `pilot.mb` in the folder where the simulation is running with the following content:
PILOT
NOW : DT = 0.5
NOW + 10 : DT = 1.0
NOW + 20 : DT = 5.0
ENDRULES
Let us now move to case 2.
If you have forces that are too small in your initial state, you can chose to initialize the simulation with
random initial atomic velocities. The input variables are `ion_velocities = 'random'` and `tempw = 300.0` to sample from a 300K distribution.
In this case the best approach is to perform the ground state calculation and the small initial thermalization in a single run by performing few steps of BO dynamics using CG
&control
! ...
nstep = 50, ! number of MD steps. It has to be big enought
dt = 20.0d0, ! if the forces are small, you don't need a small number...
! ...
/
!...
&electrons
electron_dynamics = 'cg'
/
!...
&ions
ion_dynamics = 'verlet' ! you also mix different approaches by changing this to sd
ion_velocities = 'random'
tempw = 300.d0
/
!...
A different approach is to randomize the atomic positions, and set the ionic velocities to zero.
If you have relaxed the system or if the starting system is already in
the equilibrium positions, then you need to displace ions from the
equilibrium positions, otherwise they will not move in a dynamics
@ -531,41 +626,86 @@ energy plateau.
CP dynamics
-----------
At this point after having minimized the electrons, and with ions
displaced from their equilibrium positions, we are ready to start a CP
dynamics. We need to specify `’verlet’` both in ionic and electronic
At this point after having:
1. minimized the electrons
2. a good initial configuration of the ionic degrees of freedom (good positions and velocities)
3. optionally computed the initial wavefunction velocity with the CG routine
we are ready to start a CP dynamics.
The parameter specific to the CP method that you must choose very carefully is the electron fake mass `emass`. The fake electron mass is the parameter that change how the wavefunctions follow the minimum of the DFT energy during the dynamics. The smaller the better. In the limit where both the fake mass and the timestep are zero we recover a perfect minimization of the wavefunction at each timestep, but this would require an infinite amount of steps to simulate a trajectory of a given length in time. So we need a compromise between efficiency and precision. You should always check that the forces on the ions are well converged with respect to the emass parameter.
What happens is that, the bigger the fake electron mass, the bigger the systematic error that you have on the forces. The error, intuitively, is related to the fact that the electrons have a finite classical mass, that is interacting with the ions through the DFT potential, so they have an inertia that adds up to the mass of the ion, slowing it down. This effect causes that, with a good approximation, the forces calculated with the CP method, on average, are the _true_ ground state forces time a factor between 0.0 and 1.0 that depends on `emass`. You can correct the leading order of this error by changing the ionic mass by the same factor to recover the same inertia of the atom+electron classical object. Remember also that static properties of an equilibrium molecular dynamics simulation do not depend on the mass of the ions. So if only static properties like radial distribution function or lattice parameters are needed this leading order correction is not necessary, but it becomes important if you are computing, for example, a diffusion coefficient. For a detailed analysis of this aspect of the theory see for example [P. Tangney, "On the theory underlying the Car-Parrinello method and the role of the fictitious mass parameter"](https://doi.org/10.1063/1.2162893)([arXiv](https://arxiv.org/abs/cond-mat/0601130)) or the Marx and Hutter book.
In practice, you can safely use `emass` of about 50 atomic units without affecting the dynamical properties too much.
Keep in mind that the lower the `emass` the lower the integration timestep `dt`, because the electronic degrees of freedom will move faster. You have to check
that with your parameters the CP constant of motion in the output is conserved with a good approximation.
To run the CP simulation you need to specify `’verlet’` both in ionic and electronic
dynamics. The threshold in control input section will be ignored, like
any parameter related to minimization strategy. The first time we
perform a CP run after a minimization, it is always better to put
velocities equal to zero, unless we have velocities, from a previous
simulation, to specify in the input file. Restore the proper masses for
the ions. In this way we will sample the microcanonical ensemble. The
input section changes as follow:
any parameter related to minimization strategy.
To start the simulation you must use the restart from the previous step.
Restore the proper masses for the ions. In this way we will sample the microcanonical ensemble.
If you did not use the CG routine (that computes electron velocities),
the input section changes as follow:
&electrons
emass = 400.d0,
emass = 50.d0,
emass_cutoff = 2.5d0,
electron_dynamics = 'verlet',
electron_velocities = 'zero'
electron_velocities = 'zero' ! only if you did NOT use CG
/
&ions
ion_dynamics = 'verlet',
ion_velocities = 'zero'
ion_velocities = 'zero' ! if you did NOT use CG!
/
ATOMIC_SPECIES
C 12.0d0 c_blyp_gia.pp
H 1.00d0 h.ps
If you used a method that computes the ionic velocities and the electronic velocities in a consistent way you can use the following:
If you want to specify the initial velocities for ions, you have to set
&electrons
emass = 50.d0,
emass_cutoff = 2.5d0,
electron_dynamics = 'verlet',
electron_velocities = 'default' ! if you have wfc velocities
/
&ions
ion_dynamics = 'verlet',
ion_velocities = 'default' ! if you have velocities
/
ATOMIC_SPECIES
C 12.0d0 c_blyp_gia.pp
H 1.00d0 h.ps
If you want to change the timestep, for example because you used a big timestep in a small
CG thermalization performed in the relaxation step, you must specify the following additional parameters:
&control
dt = 5.0d0, ! new integration timestep
/
&electrons
electron_velocities = 'change_step' ! if you have wfc velocities
/
&ions
ion_velocities = 'change_step' ! if you have velocities
tolp = 20.d0 ! old integration timestep
/
If you want to specify a new set of initial velocities for ions, you have to set
`ion_velocities =’from_input’`, and add the ATOMIC\_VELOCITIES card,
after the ATOMIC\_POSITION card, with the list of velocities in atomic
units.
after the ATOMIC\_POSITION card, with the list of velocities with time in atomic
units, and length in the same units specified for the atomic positions.
NOTA BENE: in restarting the dynamics after the first CP run, remember
to remove or comment the velocities parameters:
&electrons
emass = 400.d0,
emass = 50.d0,
emass_cutoff = 2.5d0,
electron_dynamics = 'verlet'
! electron_velocities = 'zero'