by Ronaldo Rodrigues Pela for exciting neon
Purpose: In this tutorial, you will learn how to run a molecular dynamics calculation with real-time time-dependent density-functional theory (RT-TDDFT) to evolve the electrons wavefunctions (also called Ehrenfest Dynamics). We consider cubic BN to illustrate it.
|
Table of Contents
|
0. Define relevant environment variables
Read the following paragraphs before starting with the rest of this tutorial!
Before starting, be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts. Here is a list of the scripts which are relevant for this tutorial with a short description.
- EXECUTE-single.sh: (Bash) shell script for running a single exciting calculation.
- PLOT-multitask.py: Python script for making plots from output files of RT-TDDFT.
From now on the symbol $ will indicate the shell prompt.
Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries.
Important note: All input parameters that will appear will be given in atomic units! For this tutorial, it can be useful to remind the conversion between a.u. and the SI for times: 1 a.u. = 0.02418884254 fs, and for forces 1 Ha/bohr = 51.42208 eV/A.
1. Theoretical background for Ehrenfest Dynamics
Please refer to Ref.[1] for a better comprehension about how Ehrenfest Dynamics has been implemented in exciting.
For more details about RT-TDDFT in exciting, please check Ref.[2] and the tutorial Real-time TDDFT.
2. Ground-state calculation
i) Preparation of the input file
The first step is to create a directory for the system that we want to investigate, BN-Ehrenfest, and to enter it.
$ mkdir BN-Ehrenfest
$ cd BN-EhrenfestAs initial condition for the time evolution done in RT-TDDFT, we need the electron density and potential from a ground-state calculation. Therefore, we create the file input.xml that could look like the following:
<input> <title>BN Ehrenfest</title> <structure speciespath="$EXCITINGROOT/species/"> <crystal scale="6.833"> <basevect>0.0 0.5 0.5 </basevect> <basevect>0.5 0.0 0.5 </basevect> <basevect>0.5 0.5 0.0 </basevect> </crystal> <species speciesfile="B.xml" rmt="1.32"> <atom coord="0.0 0.0 0.0" velocity="0.0 0.0 0.0" /> </species> <species speciesfile="N.xml" rmt="1.28"> <atom coord="0.25 0.25 0.25" velocity="0.0 0.0 0.0"/> </species> </structure> <groundstate ngridk="6 6 6" rgkmax="4.0" gmaxvr="12.0" epsengy="1e-5" do="fromscratch" outputlevel="high" nempty="4" xctype="LDA_PW"> </groundstate> </input>
N. B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT. You can do this by directly editing the input.xml file or by using the following command:
$ SETUP-excitingroot.shii) Running the ground-state calculation
You can start the calculation by invoking the script EXECUTE-single.sh.
$ EXECUTE-single.sh GSAfter the calculation is finished, move to the subdirectory GS.
$ cd GSYou can check the output files, especially the main output file INFO.OUT, for general information. In the case of a successfully finished calculation, the last lines of the INFO.OUT will contain the message
...
================================================================================
| EXCITING NEON stopped =
================================================================================Check if this is indeed the case and if the EFERMI.OUT and STATE.OUT files are present. They contain the Fermi level and the converged electron density and potential, respectively, and are the starting point of the RT-TDDFT calculation.
Please note: To obtain reliable ground-state density and potential, it is necessary to perform convergence tests with respect to the k-point mesh (parameter ngridk) and the size of the basis set (parameter rgkmax). For details see Simple convergence tests.
3. Performing an "Ehrenfest Dynamics" calculation
Move to the previous directory, create a new one named ehrenfest, and enter it:
$ cd ../
$ mkdir Ehrenfest
$ cd EhrenfestYou can prepare the input file (input.xml) of a starting from the existing ground-state input file. Copy the following files to the current directory. You can perform all of this with the following commands:
$ cp ../GS/{input.xml,EFERMI.OUT,STATE.OUT} .Edit the file input.xml and avoid the explicit ground-state calculation with the attribute do = "skip" inside the element groundstate.
... <groundstate do="skip" ... > ... </groundstate> ...
To perform the "Ehrenfest Dynamics" calculation, add the elements xs and MD (as shown below) to the input file inside the input element.
... <xs xstype ="RT-TDDFT" ngridk="4 4 4" vkloff="0.0 0.0 0.0" nosym="true" reducek="false" nempty="2"> <realTimeTDDFT propagator="SE" timeStep="0.2d0" endTime="16.0d0" readPmatBasis="false" printTimingGeneral="true" printTimingDetailed="true" calculateNExcitedElectrons="false" printAfterIterations="1" vectorPotentialSolver="improvedeuler"> <laser fieldType="total"> <sinSq amplitude="10.d0" omega="1.d0" phase="0.d0" t0="0.d0" pulseLength="16.d0" direction="x"/> <sinSq amplitude="10.d0" omega="1.d0" phase="0.d0" t0="0.d0" pulseLength="16.d0" direction="y"/> <sinSq amplitude="10.d0" omega="1.d0" phase="0.d0" t0="0.d0" pulseLength="16.d0" direction="z"/> </laser> </realTimeTDDFT> </xs> <MD type="Ehrenfest" printAllForces="true" timeStep="1.0d0" integrationAlgorithm="HeunSimplified"> </MD> ...
The following figure illustrates $A_x(t)$, the x component of the vector potential:
In this example, we are assuming $A_x(t)= A_y(t) = A_z(t)$, i.e., one laser pulse applied along the [111] direction.
For a detailed explanation of the element realTimeTDDFT, please check the tutorial Real-time TDDFT. We go here through all attributes that appear in MD to clarify their meaning:
- type: Type of MD to be performed. In Ehrenfest MD, the outputs are printed out every n steps, where n is given by printAfterIterations.
- printAllForces: When set to "true", the following components of the total force acting on an ion are printed out: external, Hellmann-Feynman, core and valence corrections..
- timeStep: Time step $\Delta t_N$ of the MD calculation. When doing Ehrenfest MD, if $\Delta t_N$ is not a multiple of the time step $\Delta t_e$ used in the electronic evolution (RT-TDDFT evolution), it will be automatically rounded down to $n \Delta t_e$, where $n$ is the largest integer lower or equal the ratio $\Delta t_N / \Delta t_e$.
- integrationAlgorithm: Algorithm employed for the trajectory of the nuclei
All relevant parameters for an "Ehrenfest Dynamics" calculation are found in Input Reference.
ii) Running exciting
Having modified the input file as described above, you can launch the Ehrenfest Dynamics calculation inside the subdirectory Ehrenfest using
$ export OMP_STACKSIZE=2G
$ time exciting_smpBe aware that, depending on the computer you are using, this calculation may take several minutes.
iii) Output files
Apart from the output files related to RT-TDDFT (explained in the tutorial Real-time TDDFT), an Ehrenfest Dynamics calculation also generates, for each atom, the output files:
- ATOM_xyzw.OUT: contains 10 columns, where
- the first one has the time $t$;
- columns 2-4 have the cartesian coordinates of the atom position $x(t)$, $y(t)$ and $z(t)$ (at the given time $t$);
- columns 5-7 have the cartesian components of the atom velocity $v_x(t)$, $v_y(t)$ and $v_z(t)$;
- columns 8-10 have the cartesian components of the total force acting on that atom $F_x(t)$, $F_y(t)$ and $F_z(t)$.
- FEXT_xywz.OUT, FHF_xywz.OUT, FVAL_xywz.OUT and FCR_xywz.OUT. These files are printed out only if printAllForces is set to "true". In each file, there are four columns, with the first one containing the time, and the other three, the cartesian components ($x$, $y$ and $z$) of the considered force:
- FEXT_xywz.OUT: external force due to the applied electric field;
- FHF_xywz.OUT: Hellman-Feynman force;
- FVAL_xywz.OUT: valence corrections to the Hellman-Feynman force;
- FCR_xywz.OUT: core corrections to the Hellman-Feynman force.
In the previous list of outputs, x, y, z and w are integer numbers 0-9 and the sequence xyzw refers to the number of a specific atom as given in input.xml. For instance 0001 refers to the 1st atom, 0002, to the second, etc.
4. Visualizing the outputs
To visualize the $x$ component of the forces acting on boron (1st atom in input.xml), you can use the script PLOT-multitask.py [3] as follows :
$ PLOT-multitask.py --ylabel '$F_x$ [a.u.]' --xlabel 'Time [a.u.]' -f ATOM_0001.OUT -k 0 7 -f FEXT_0001.OUT -k 0 1 -f FHF_0001.OUT -k 0 1 -f FVAL_0001.OUT -k 0 1 -f FCR_0001.OUT -k 0 1 -c 'Total' 'Ext' 'HF' 'Valence corr.' 'Core corr.' --xlim 0. 16. --legend_position 'upper left'The resulting figure should look like
To check the time evolution of the position of boron, you can run the command:
$ PLOT-multitask.py --ylabel 'Position [bohr]' --xlabel 'Time [a.u.]' -f ATOM_0001.OUT -k 0 1 -f ATOM_0001.OUT -k 0 2 -f ATOM_0001.OUT -k 0 3 -c '$x(t)$' $y(t)$' '$z(t)$' --xlim 0. 16. --legend_position 'upper left'to obtain a figure like
5. Converging the results
From the previous two figures, you can conclude that the test case analyzed here should be considered as a pedagogical example. Following remarks can be done
- The time step $\Delta t_N$ in the MD element needs to be smaller, to obtain a more smooth curve of the forces or the positions over time $t$. There is an optimal ratio between $\Delta t_N/\Delta t_e$ that allows for a good compromise between the computational time and accuracy (see discussion in Ref.[1]).
- endTime is quite small in the example. For reallistic MD calculations, this should be on the order of fs or ps (1 fs = 41.34 a.u. and 1 ps = 41341 a.u.).
- The frequency and the duration of the laser pulse here are not intended to reproduce experimental settings. Typically, the laser pulse could last for some fs. Moreover, a laser frequency closer to the phonon frequencies would lead to more pronounced atom movements.
- For better converged calculations, the following parameters should be carefully checked: