-
Notifications
You must be signed in to change notification settings - Fork 22
Tutorial TDDFT
Authors: Xavier Andrade and Alfredo Correa
In this tutorial you will learn how to do an Ehrenfest-TDDFT molecular-dynamics simulation with Qball.
We will be simulating a particle, a proton in this case, entering a material, Aluminum. As the particle moves through the material, it feels a drag force. At low speed, this force is caused by the interaction with the ions of the material. At high speeds, however, the force is caused by the interaction of the particle with the electrons of the material. This process is known as electronic stopping.
In order to accurately simulate the interaction of the proton with the electrons, we need to use non-adiabatic dynamics where we allow the electrons to get excited. So in this example we will be using TDDFT to model the electrons.
First, we generate a file with the coordinates of an aluminum supercell with 32 atoms (for actual calculations you probably need a larger cell). Call it al32.sys and copy the following content:
set cell 15.306 0.0 0.0 0.0 15.306 0.0 0.0 0.0 15.306 bohr
species aluminum Al.xml
atom Al1 aluminum 0.000000 0.000000 0.000000 bohr
atom Al2 aluminum 3.826500 3.826500 0.000000 bohr
atom Al3 aluminum 3.826500 0.000000 3.826500 bohr
atom Al4 aluminum 0.000000 3.826500 3.826500 bohr
atom Al5 aluminum -7.653000 0.000000 0.000000 bohr
atom Al6 aluminum -3.826500 3.826500 0.000000 bohr
atom Al7 aluminum -3.826500 0.000000 3.826500 bohr
atom Al8 aluminum -7.653000 3.826500 3.826500 bohr
atom Al9 aluminum 0.000000 -7.653000 0.000000 bohr
atom Al10 aluminum 3.826500 -3.826500 0.000000 bohr
atom Al11 aluminum 3.826500 -7.653000 3.826500 bohr
atom Al12 aluminum 0.000000 -3.826500 3.826500 bohr
atom Al13 aluminum 0.000000 0.000000 -7.653000 bohr
atom Al14 aluminum 3.826500 3.826500 -7.653000 bohr
atom Al15 aluminum 3.826500 0.000000 -3.826500 bohr
atom Al16 aluminum 0.000000 3.826500 -3.826500 bohr
atom Al17 aluminum 0.000000 -7.653000 -7.653000 bohr
atom Al18 aluminum 3.826500 -3.826500 -7.653000 bohr
atom Al19 aluminum 3.826500 -7.653000 -3.826500 bohr
atom Al20 aluminum 0.000000 -3.826500 -3.826500 bohr
atom Al21 aluminum -7.653000 0.000000 -7.653000 bohr
atom Al22 aluminum -3.826500 3.826500 -7.653000 bohr
atom Al23 aluminum -3.826500 0.000000 -3.826500 bohr
atom Al24 aluminum -7.653000 3.826500 -3.826500 bohr
atom Al25 aluminum -7.653000 -7.653000 0.000000 bohr
atom Al26 aluminum -3.826500 -3.826500 0.000000 bohr
atom Al27 aluminum -3.826500 -7.653000 3.826500 bohr
atom Al28 aluminum -7.653000 -3.826500 3.826500 bohr
atom Al29 aluminum -7.653000 -7.653000 -7.653000 bohr
atom Al30 aluminum -3.826500 -3.826500 -7.653000 bohr
atom Al31 aluminum -3.826500 -7.653000 -3.826500 bohr
atom Al32 aluminum -7.653000 -3.826500 -3.826500 bohr
Make sure you understand what the commands in the file do. You will need the Al.xml pseudopotential file.
Now, lets create the input file for a ground state calculation. You could call it gs.i. First, you need to force Qball to use complex wave-functions with the force_complex variable.
set force_complex_wf ON
This is necessary since in TDDFT the wave-functions are always complex.
Then we load the atoms of the aluminum supercell, simply by adding the name of the file
al32.sys
Now we need to add the proton to the system. First you will need a pseudopotential for hydrogen, which you can find here. Place the proton at coordinates 0.00000 1.91325 1.91325 using the atom command. To add an atom without electrons we have to append a + sign to the name of the atom, for example call it "+proton".
Now add the standard commands to perform a ground state calculation, using an energy cutoff (ecut) of 30 Rydberg, a preconditioner cutoff (ecutprec) of 4 Rydberg, and a threshold (threshold_scf) of 1e-6 for 5 steps. If you don't know how to set the necessary variables, we recommend you have a look the ground state tutorial.
Since the system is metallic we have to adjust a few parameters. First, we will use the PSDA diagonalizer (wf_dyn variable). We need to add some empty states with
set nempty 12
and now add some smearing with an electronic temperature of 1000 K. This is done with the commands
set smearing fermi
set smearing_width 1000 kelvin
Finally we reduce the amount of charge mixing a bit with
set charge_mix_coeff 0.6
(the default is 0.5).
We now can randomize the wave-functions and issue the run command. This time we tell qball to do 5 optimization steps per SCF iteration by passing a third value to run as this
run 0 200 5
We still have to add a few commands to save the results after the calculation is done. We can then use the results as a starting point for the real-time calculation. First we save the coordinates using the savesys command
savesys md.sys
This will generate an file in the Qball input format that we can load from other input files. Then, we save the wave-functions with the save command
save -states restart/md
This command saves several files, one per state, in the restart/ directory.
You can now go ahead and tell Qball to execute the newly created input file. Carefully check the output for things that might have gone wrong. Qball tends to keep going even when something failed (we are working to change that, meanwhile grep qball -A 1 can help spotting warning messages).
After the run is successful, I recommend you have a look at md.sys to see what is saved there.
Now that we have the set of ground-state orbitals, we can set up the real-time calculation. Create a new input file for this called td.i.
As before, the first thing to do is to set the force_complex_wf variable to force Qball to work with complex states. We also need to set the energy cutoff (ecut) to the value we used before. After that, we load the saved geometry in the input file with
md.sys
and then we load the wave-functions that we obtained from the ground state calculation with
load -states restart/md
If you had a look at the md.sys you probably noticed that after the positions for each atom, there are three zeros. These zeros are the components of the velocities. So to get the proton moving in our real-time calculation, we need to give it some velocity. We do this with the set_velocity command. For example, to set a velocity of 0.1 atomic units in the x direction we add
set_velocity +proton 0.1 0.0 0.0 atomicvelocity
to our input file. Of course, we could have set the velocity when we added the proton while setting up the ground state calculation. However, using set_velocity is more convenient since we can change the velocity without having to redo the ground state calculation (which does not depend on the ionic velocities in any case).
Now we need to tell Qball how to propagate the electrons and the ions. For the electrons, we set the variable wf_dyn. In real-time TDDFT this variable selects the propagator we want to use. In this case we want the enforced time-reversal symmetry propagator, so we set wf_dyn to ETRS (fourth-order Runge-Kutta, FORKTD, is another reasonable alternative). We set the propagation for the ions with the variable atoms_dyn. For example, if we want to do molecular dynamics we set atoms_dyn to MD.
Finally, we need to set the time-step for the time integration. In Qball this is controlled by two variables, dt for the ions, and TD_dt for the electrons. Both must have the same value. We need to do some work to find the optimum value for the time-step. we will do this in the next step, for the moment set both time-step variables to 0.1 atomictime.
We now can add the run command. For TDDFT the first number will determine the number of time-steps we want to do and the second number should always be 1. So, to perform 100 iteration we would use
run 100 1
With this, we should have a complete Qball input file for a TDDFT run. Check that it works and that there are no errors in the execution.
For TDDFT simulations we want to use a time-step that is as large as possible, as this makes our calculation faster. However, for values that are too large, the propagation will become unstable.
So we need to do a few runs of trial and error to determine what is the maximum time-step that we can use. For example, set the time-step to 1.0, and run Qball grepping the output for econst like this
qball td.i | grep econst
(If you are using a non-interactive queue to run you can save the output to a file and use grep to search that file.) The value of econst is the total energy of the electrons plus the ionic kinetic energy. We want a time step for which econst is conserved, but you probably see its value becoming very large just after a few iterations. This is because the time-step is too large. Try different values of the time-step to find the largest one that does not make the calculation explode (do not expect a very good energy conservation in this case).
The maximum time-step depends on the cutoff that you are using in your calculation. If you increase the cutoff you need to decrease your time-step.
Now that we have our optimal time step, we can start the electronic stopping calculation.
We are going to make a small change to our previous dynamics. Instead of using regular molecular dynamics for the ions we are going to use IMPULSIVE. This forces the ions to move without changing their velocities. Since electronic stopping is a very fast process we can neglect the loss of speed of the proton, or the movement of the ions of the material.
Now, set the initial velocity to 0.5 a.u. in the x direction, and propagate for enough time for the proton to move along the whole unit cell. Save the output of Qball to a file. We will assume you call it log_0.5.
The quantity we want to calculate is the stopping "power", which can be understood as the drag force that the proton feels as it moves through the material. We can directly obtain this force from the Qball output file. Since the output is XML we can write a parser routine using an XML library...
Or we can just use grep. The following grep command will extract the force from the output file and add it to a file called "force_0.5"
grep -A 3 proton log_0.5 | grep force > force_0.5
Now go to a plotting tool and have a look at minus the force in the x direction as a function of time. This is the instantaneous stopping power. Why do you think there are oscillations in the force?
Now calculate the average stopping power. You can do this from the command line using the awk utility with the following command
awk '{ sum += $2; n++ } END { if (n > 0) print -sum / n; }' force_0.5
Does the number you get looks reasonable? What is the sign of the force? Compare it with the plot you just made.
Now, we want to see how the electronic stopping changes with respect to the velocity of the particle. So repeat your calculations for larger velocities, I suggest you do calculations for the values 0.5, 1, 1.5, 2, 3, 4, and 5 a.u.. If the calculations run fast in your system, you can add more points. You can also consider that for 0 velocity the stopping power is 0. Note that you have to adjust the simulation time for each velocity so that in each case the particle travels through the unit cell once.
Plot the stopping power as a function of the velocity. What is the location of the maximum? Can you guess why the curve has that shape?
Have a look at the value of econst during your propagation. You will see that the energy is not conserved. Can you guess why? You can also get the electronic stopping from the slope of this curve.
In order to make things run quickly, the parameters we use in this tutorial are not converged. If you want to make predictions you should properly converge the calculation with respect to the energy cutoff and supercell size. You might also need to include semicore states in your aluminum pseudopotential, specially for high electronic velocities.
If you want to know more details about electronic stopping calculations, you can check the following paper.