-
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 doing a calculation of electronic stopping of a proton in Aluminum. This is the force that fast particule feels when entering a material.
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
species aluminum Al.xml
atom Al1 aluminum 0.000000 0.000000 0.000000
atom Al2 aluminum 3.826500 3.826500 0.000000
atom Al3 aluminum 3.826500 0.000000 3.826500
atom Al4 aluminum 0.000000 3.826500 3.826500
atom Al5 aluminum -7.653000 0.000000 0.000000
atom Al6 aluminum -3.826500 3.826500 0.000000
atom Al7 aluminum -3.826500 0.000000 3.826500
atom Al8 aluminum -7.653000 3.826500 3.826500
atom Al9 aluminum 0.000000 -7.653000 0.000000
atom Al10 aluminum 3.826500 -3.826500 0.000000
atom Al11 aluminum 3.826500 -7.653000 3.826500
atom Al12 aluminum 0.000000 -3.826500 3.826500
atom Al13 aluminum 0.000000 0.000000 -7.653000
atom Al14 aluminum 3.826500 3.826500 -7.653000
atom Al15 aluminum 3.826500 0.000000 -3.826500
atom Al16 aluminum 0.000000 3.826500 -3.826500
atom Al17 aluminum 0.000000 -7.653000 -7.653000
atom Al18 aluminum 3.826500 -3.826500 -7.653000
atom Al19 aluminum 3.826500 -7.653000 -3.826500
atom Al20 aluminum 0.000000 -3.826500 -3.826500
atom Al21 aluminum -7.653000 0.000000 -7.653000
atom Al22 aluminum -3.826500 3.826500 -7.653000
atom Al23 aluminum -3.826500 0.000000 -3.826500
atom Al24 aluminum -7.653000 3.826500 -3.826500
atom Al25 aluminum -7.653000 -7.653000 0.000000
atom Al26 aluminum -3.826500 -3.826500 0.000000
atom Al27 aluminum -3.826500 -7.653000 3.826500
atom Al28 aluminum -7.653000 -3.826500 3.826500
atom Al29 aluminum -7.653000 -7.653000 -7.653000
atom Al30 aluminum -3.826500 -3.826500 -7.653000
atom Al31 aluminum -3.826500 -7.653000 -3.826500
atom Al32 aluminum -7.653000 -3.826500 -3.826500
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 putting 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, you can find it 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 set a ground state calculation, using an energy cutoff of 30 Rydberg, a preconditioner cutoff of 4 Rydberg, and a threshold of 1e-6 for 5 steps. If you don't know how to set this, we recommend you to 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 (0.006333621 Rydberg). This is done with the commands
set smearing fermi
set smearing_width 0.006333621
Finally we reduce a bit the amount of charge mixing 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 command so after the calculation is done, the results are saved to use them 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).
After the run is successful, I recommend you to have a look to the md.sys to see what it is saved there.
Now that we have the set of ground-state orbitals, we can setup the real-time calculation. Create a new input file for this, call it 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 to the value we used before. After that, we load the saved geometry in the input file with
md.sys
and the 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 are 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
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 (that 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 it to ETRS (fourth-order Runge-Kutta, FORK, is another reasonable alternative). For the ions the variable us called atoms_dyn. If we want to do molecular dynamics we set it to MD.
Finally, we need to set 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 it to 0.1.
We now can add the run command, for TDDFT the first number will determine the number of time-step we want to do and the second number should always be 1. So, to 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 over that file.) The value of econst is the total energy electrons plus the ionic kinetic energy. In this case it should be conserved, but you probably see its value becoming very large just after a few iteration, 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 stopping calculation.
We are going to make a small change to our previous dynamic, instead of using regular molecular dynamics for the ions we are going to use IMPULSIVE, this moves the ions 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 direction x, 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", it 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 put into 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? 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 of up to 5 a.u., storing for each one the average stopping power. Note that faster particles require shorted simulation times.
Plot the stopping power as a function of the velocity, what do you think is happening here?
== Extras
If you want to know more details about electronic stopping calculations, you can check the following paper.