5. Tutorial (Shell and Python)
INQ’s command line interface takes a unique approach to run an electronic structure simulation. The basic premise of INQ command interface is that a certain directory/folder corresponds to a specific simulation. All commands that you run, either through a Unix shell or Python functions, are relative to the directory and are persistent (they are stored in files and outlive the run time of the commands).
This approach allows to do away with input files and make it simpler and much less error-prone to run a simulation; each command will give you feedback and at any time you can query the simulation parameters before running. This strategy avoid too common situations like having a simulation job waiting for hours in a supercomputer queue, just to discover that there was a small syntax error in your input file.
INQ commands can be called directly from the Unix shell command line through the inq executable, or from Python using the pinq module. Note that both ways of using the interface modify the underlying meta-data in the current folder, so it is possible to mix shell and python commands for the same simulation. For example you can set up your simulation using a python script and then run it using the shell inq run command.
5.1. Nitrogen molecule
5.1.1. Setting up the run directory
We will start our tutorial by a simple molecular system calculation. We are assuming you are right now running in a Unix terminal and you have INQ installed. First, let’s create a folder for our simulation and get into it:
mkdir nitrogen
cd nitrogen
Note
The basic premise of INQ command interface is that a certain directory/folder corresponds to a specific simulation. All commands that you run, either through a Unix shell or Python functions, are relative to the directory and are persistent (they are stored in files and outlive the run time of the commands).
This tutorial is designed to be done using either the shell command line interface or the Python interface. Instructions for both cases will be given. So just pick which one you will like to use (or you can mix and match, since both modify the persistent data in the current simulation).
If you are going to use Python, start the python interpreter and import the INQ module pinq:
import pinq
You don’t need to do anything at this point if you are going to use the shell command. Just verify that the inq command is available. The first command we will run is clear. This command removes all INQ meta-data in the current folder. Since we just created the directory it should be empty, but let’s run it anyways.
inq clear
pinq.clear()
Note
It is a good idea to run clear each time you start a new simulation (especially if you are making a script) to get rid of all previous information that might exist in the current directory.
5.1.2. Defining the simulation parameters
Now we are ready to start setting up the simulation. First we define a simulation cell of size 10 x 10 x 12 bohr. This is how to do it:
inq cell orthorhombic 10 10 12 bohr
pinq.cell.orthorhombic(10, 10, 12, "bohr")
Note
We need to explicitly tell INQ what units we will be using. This is by design, we want to avoid ambiguities in units that can cause errors (and crash spacecraft ).
Once we defined the cell, we can query what is the cell in our simulation at any time with the following command:
inq cell
pinq.cell.status()
You should get an output like this:
Cell:
Lattice vectors [b] = 10 0 0
0 10 0
0 0 12
Volume [b^3] = 1200
Periodicity = 3d (fully periodic)
Query commands like this exist for most settings and can be used at any time. It is a good practice to use them check to make sure the code is actually going to do what we want.
Now that the cell is defined, let’s add the atoms. This is done with the ions insert command. The following commands will insert two nitrogen atoms at a distance of 1.3 Angstrom.
inq ions insert N 0 0 1.3/2 Angstrom
inq ions insert N 0 0 -1.3/2 Angstrom
pinq.ions.insert("N", [0.0, 0.0, 1.3/2], "Angstrom")
pinq.ions.insert("N", [0.0, 0.0, -1.3/2], "Angstrom")
As you can see the inq shell command has some basic math interpretation capabilities, so we can write 1.3/2 and let INQ do the division.
Now that we defined the ions in the simulation we have to take care of the electrons. The only parameter we really need is the plane-wave cutoff, as for other values we can rely on defaults.
The cutoff is set by the electrons cutoff command. In this case we want to set a value of 40 Hartree that we do by:
inq electrons cutoff 40.0 Hartree
pinq.electrons.cutoff(40.0, "Hartree")
You can check the options that are set for the electrons running the following command:
inq electrons
pinq.electrons.status()
This will give you an idea of the parameters available and their default values.
5.1.3. Running the simulation
We have now set all the parameters needed to do a ground-state calculation for nitrogen. However, the simulation has not been done yet, we need to tell INQ to run it. This is done with the run ground state command:
inq run ground state
pinq.run.ground_state()
This should take a few seconds. At the end of the output INQ should print something like:
SCF ended after 54 iterations with resulting eigenvalues and energies:
st = 1 occ = 2.000 evalue = -0.916321357828 Ha ( -24.934371418951 eV) res = 7e-09
st = 2 occ = 2.000 evalue = -0.523560815064 Ha ( -14.246813862499 eV) res = 3e-08
st = 3 occ = 2.000 evalue = -0.356182318108 Ha ( -9.692213475863 eV) res = 1e-08
st = 4 occ = 2.000 evalue = -0.356180895755 Ha ( -9.692174771677 eV) res = 2e-08
st = 5 occ = 2.000 evalue = -0.351473943512 Ha ( -9.564092091423 eV) res = 5e-08
Energy:
total = -20.705372548423 Ha
kinetic = 12.720356398467 Ha
eigenvalues = -5.007438660535 Ha
hartree = 13.547930959702 Ha
external = -37.150127253812 Ha
non-local = -1.418675597565 Ha
xc = -5.731429377763 Ha
nvxc = -6.254854127029 Ha
exact-exchange = 0.000000000000 Ha
ion = -2.673427677451 Ha
5.1.4. Querying the ground-state results
Now that we finished the calculation we can access the individual simulation results using commands. The command results ground state will give us an overview of the simulation data:
inq results ground state
pinq.results.ground_state.status()
It should print something like this:
Ground-state results:
iterations = 54
dipole [a.u.] = 0 0 0
magnetization = 0 0 0
stress [a.u.] = 0.00151421 6.18129e-09 -1.35831e-07
= 6.18129e-09 0.00151421 -6.06216e-08
= -1.35831e-07 -6.06216e-08 -0.000340571
pressure = -0.000895947 Ha/b^3 | -26.3597 GPa
total energy = -20.70537255 Ha | -563.42188978 eV
We can also use this command to query specific values, this is particularly useful if we want to use scripts to post process the results. For example, to get the total energy we can use the following command
inq results ground state energy total
pinq.results.ground_state.energy.total()
These function print/return a single value that can be captured to use it in scripting. In bash to capture the value you can use backticks to assign it to a variable. In python this is straightforward as the function returns a floating point value that you can directly assign to a variable.
energy=`inq results ground state energy total`
energy = pinq.results.ground_state.energy.total()
You can also query the forces with the command results ground state forces.
inq results ground state forces
pinq.results.ground_state.forces()
This variable will print/return the forces for all atoms and all directions. In the shell if you want to query a specific value, you can add more arguments to the command to specify which one you want. The first option specifies the index of the atom you want the force for (starting for 0) and the second one the component of the forces you want (x, y, or z). In python the forces are returned in an array, so you can just access the element you want. For example, if you want the z component of the force for atom 0 you can use:
inq results ground state forces 0 z
pinq.results.ground_state.forces()[0][2]
This should give you a value around 0.2972.
Based on the results of the simulation would you consider the interatomic distance we used, 1.3 angstrom, corresponds to stable configuration?
5.1.5. Changing the ion positions
As you probably figure out the atoms in our simulation are not at the minimum energy configuration since the force is rather large. The experimental value of the interatomic distance of the N2 molecule is 110 pm. So let’s use that value in our simulation.
INQ has several commands to handle atoms that include remove (to remove a single ions) and ions clear that removes all ions in the system. You can check how to use them using the documentation for ions that can be accessed from the same command interface.
inq help ions
help(pinq.ions)
In this case we are going to use the ions clear command to remove all ions and then insert them again in the new position (make sure you replace NNN by the experimental value in pm):
inq ions clear
inq ions insert N 0 0 NNN/2 pm
inq ions insert N 0 0 -NNN/2 pm
pinq.ions.clear()
pinq.ions.insert("N", [0.0, 0.0, NNN/2], "pm")
pinq.ions.insert("N", [0.0, 0.0, -NNN/2], "pm")
After the command verify that you have two atoms in the correct position. Now run inq again using the run ground state command.
What can you say about the results? Did the energy become smaller? Did the force become close to zero?
5.1.6. Scripting: the potential energy surface of N2
We are now going to explore the power of INQ interface by writing a script that iterates over the interatomic distance calculating the energy for each value.
This is the template of the script we are going to use:
set -e #make the script fail if a command fails
echo "#distance energy" > pes.dat
for dist in `seq 0.8 0.025 1.301`; do
#run the ground state for distance `$dist` here
energy=0.0
echo $dist $energy >> pes.dat
done
import pinq
import numpy
pes = []
for dist in numpy.arange(0.8, 1.301, 0.025):
#run the ground state for distance `dist` here
energy = 0.0
pes.append([dist, energy])
Note
For bash save the script on a file called pes.sh, to run it just type bash pes.sh. For Python name the file pes.py and run it with python3 pes.py.
This scripts iterates over the distance from 0.8 to 1.3 with steps of 0.025 size. For each distance it saves the value of the distance and the energy (that we here just set to zero) For bash the distance and energy values are saved in a file called pes.dat. You can easily plot this file using a program like gnuplot. For Python we store the results in an array. You can use matplotlib to plot them (you might need to transpose the array using numpy.transpose(pes)).
Now it is your job to modify this script to do the actual calculations. Inside the loop iteration you need to set the interatomic distance and run INQ to calculate the new energy. You also need to obtain the total energy value for INQ and store it in the energy variable, we show you how to do this previously in this tutorial. Note how much simpler it is to do get the energy compared to other codes where you would need to parse some text file.
There are two approaches you can take to write your script:
The first one is to rely on the calculation you already setup in the current folder, so you only need to change the coordinates. This makes a simpler script but it depends on the simulation already setup. It won’t run on a different directory.
The second approach is to run a standalone script by adding all the calculation setup (the cell, cutoff). This script is a bit more complex but can run independently of the previous simulation you set. If you decide to do it this way, don’t forget to use the clear command before hand to avoid issues.
Once you have your script, execute it. Now you can plot the potential energy surface with your favorite plotting tool. It should look something like this (hopefully prettier):
-20.2 +-------------------------------------------------------------------+
| + + + + |
| |
-20.3 |*+ +-|
| * |
| * |
| * |
-20.4 |-+ * +-|
| ** |
| * |
-20.5 |-+ ** +-|
| * |
| ** |
-20.6 |-+ ** +-|
| ** |
| ** |
| *** |
-20.7 |-+ **** *****|
| ******* ********** |
| + + ******************** + |
-20.8 +-------------------------------------------------------------------+
0.8 0.9 1 1.1 1.2 1.3
Now modify your script to also get the interatomic force and plot that value too. The result should look something like this:
3 +--------------------------------------------------------------------+
| + + + + |
| |
2.5 |*+ +-|
| * |
| * |
2 |-+ ** +-|
| ** |
| ** |
1.5 |-+ ** +-|
| *** |
1 |-+ ** +-|
| ** |
| **** |
0.5 |-+ *** +-|
| **** |
| ******* |
0 |-+ ****** +-|
| ****************** |
| + + + + ********|
-0.5 +--------------------------------------------------------------------+
0.8 0.9 1 1.1 1.2 1.3
5.2. Benzene optical absorption spectrum
Now we are going to move to the domain of time-dependent density functional theory (TDDFT) to calculate excited states properties.
Note
This tutorial builds upon what you learned before, so make sure to check the previous tutorial if you haven’t done it.
5.2.1. Simulation parameters
We will start by a simple case, the optical absorption of a molecule, in this case benzene. The first thing to do is to create a new folder for our simulation (as we did before). Remember that in INQ you need a separate folder for each simulation.
This time, instead of giving the ion positions one by one we will load it from an XYZ file. Create that file, naming it benzene.xyz with the following content:
12
Geometry of benzene (in Angstrom)
C 0.000 1.396 0.000
C 1.209 0.698 0.000
C 1.209 -0.698 0.000
C 0.000 -1.396 0.000
C -1.209 -0.698 0.000
C -1.209 0.698 0.000
H 0.000 2.479 0.000
H 2.147 1.240 0.000
H 2.147 -1.240 0.000
H 0.000 -2.479 0.000
H -2.147 -1.240 0.000
H -2.147 1.240 0.000
To load the XYZ into INQ we will use the following command.
inq ions file benzene.xyz radius 10 bohr
pinq.ions.file("benzene.xyz", 10.0, "bohr")
This simply tells inq to load all the atoms from an XYZ file. Since a standard XYZ does not contain information about the cell, we need to give this information as an extra argument. We could give INQ a full cell definition as before, but for convenience INQ can also take a distance argument when importing an XYZ. This automatically creates a finite orthorhombic cell where no atom is closer than this distance to the cell boundaries. This is practical for molecular simulations where you need to converge the size of the cell.
Note
INQ can also read the cell and atomic coordinate information from CIF and VASP POSCAR files.
Now we need to setup the parameters for the electrons. This time, however, instead of specifying the energy cutoff we will use the real-space grid spacing as
inq electrons spacing 0.43 bohr
pinq.electrons.spacing(0.43, "bohr")
Note
Using the spacing or the cutoff are equivalent and produce the same result. It’s just a matter of convenience.
5.2.2. Ground state calculation
The first step in a real-time calculation is to calculate the ground state, that will set initial condition for our real-time simulation. So calculate the ground state with the run ground state command.
5.2.3. First time propagation, finding the optimal time-step
Now that we have a ground state we can do a real-time calculation. For this we need to set a few options. The first one is the time-step, this is a crucial quantity in a real-time propagation. If the time step is too large the simulation is unstable, if it is too small the calculation would be slower than necessarily. The time-step depends on the energy cutoff/spacing used for the simulation and possibly other parameters, so we should optimize it for each case. In INQ you set the time-step with the following command:
inq real-time time-step 0.08 atu
pinq.real_time.time_step(0.08, "atu")
Note
As with other quantities, INQ requires units for time magnitudes. In this case atu stands for atomictimeunits. You can also use attoseconds (as), femtoseconds (fs), and picoseconds (ps).
Our first task will be to find the optimal time-step for this case. After you add the propagation line above to your file, run INQ with:
inq run real time
pinq.run.real_time()
You should get an output like this:
[trace] starting real-time propagation
[info] step 0 : t = 0.000 e = -39.547750090759
[info] step 1 : t = 0.080 e = -39.547744740076 wtime = 0.444
[info] step 2 : t = 0.160 e = -39.547742901329 wtime = 0.416
[info] step 3 : t = 0.240 e = -39.547736053030 wtime = 0.427
[info] step 4 : t = 0.320 e = -39.547722024915 wtime = 0.414
[info] step 5 : t = 0.400 e = -39.547611402019 wtime = 0.417
[info] step 6 : t = 0.480 e = -39.546318301288 wtime = 0.438
[info] step 7 : t = 0.560 e = -39.529646335901 wtime = 0.442
[info] step 8 : t = 0.640 e = -39.302794826777 wtime = 0.418
[info] step 9 : t = 0.720 e = -36.145930399491 wtime = 0.419
[info] step 10 : t = 0.800 e = 5.739670746138 wtime = 0.412
...
We are interested in the energy values that start with e = …. Since our system does not have any (time-dependent) perturbation yet the energy should be conserved. But in this case we have deliberately chosen a too-large time step that will make the energy increase rapidly because the simulation is unstable. Let it run for some steps to convince yourself that is the case.
Now modify the time-step in your input and see what happens with the energy values. Repeat this until a few times you find the largest value of the time-step that give you a stable energy. Note that we are looking for stability rather than exact energy conservation, with a proper time-step the energy will oscillate around a value instead of exploding.
5.2.4. Adding a perturbation
Now that you have found the optimal step, let’s do some more physically interesting. Up to now our propagation doesn’t have any perturbation, so it stays in the ground-state just changing the phase of the orbitals. Let’s change that.
INQ defines the concept of a perturbation: an external change in the systems that affects the electrons. Note that a perturbation does not necessarily have to be small, real-time TDDFT can study the response of strongly-driven systems as well.
In this case we will use a kick, this is an electric field perturbation of the form \(\vec{E}(r, t)=\vec{\kappa}\delta(t)\). The advantage of a kick is that it is flat in Fourier space, so it excites all frequencies equally, making it useful for linear response calculations. This is how we define a kick:
inq perturbation kick 0.01 0 0
pinq.perturbation.kick([0.01, 0, 0])
The argument is the direction and intensity of the kick \(\kappa\). Since we are interested in linear response we use a small amplitude to avoid non-linear effects to appear on our simulation (these effects depend at least on \(\kappa^2\)).
5.2.5. Defining the observables
To make our simulation useful we also need to obtain some data out of it. There are many possible observables in a simulation that we want to look at, however calculating them and storing them can be expensive so we only want to calculate what we really need. For this reason we need to tell INQ beforehand what observables we want.
For the calculation of the absorption spectrum of a molecule the quantity we are interested is the dipole. So we tell INQ:
inq real-time observables dipole
pinq.real_time.observables.dipole()
5.2.6. Run the propagation
We are almost ready to run the simulation, the only part we are missing is top specify the length of the simulation. The longer the simulation is, the more well defined the absorption peaks will be, but at the same time the simulation will be more expensive. In practice, the time we need around \(10\) fs (\(~400\) atu) for a converged spectrum. In this tutorial we will use a shorter time of \(5\) fs to avoid waiting for too long. This is set with the command:
inq real-time propagation-time 5 fs
pinq.real_time.propagation_time(5, "fs")
Note
The propagation time command sets the number of steps based on the currently set time-step. If you change the time-step, the total propagation time will change.
Now we can run the time-propagation with INQ. This is going to be longer than the runs we have done so far. If your machine has multiple CPUs or GPUs it is recommended that you run this in parallel. Depending on your system it could take around 10 minutes to an hour. Can estimate how long it will take to run in your machine? (INQ prints the time it take for each iteration as the wtime value.)
Note
Even if you are working in Python you can call the shell version of inq run as long as you are standing in the same folder.
5.2.7. Analyzing the results and obtaining the spectrum
Once the simulation is done we can analyze the results. The command
inq result real-time dipole
pinq.results.real_time.dipole()
will give us the back the time-dependent dipole function.
The shell version will give us back the data to the terminal 4 columns: the time and the x, y, and z components of the dipole. You can directly plot it using gnuplot. First run gnuplot and then type following command:
plot "<inq results real time dipole" w l
In Python things are a bit different, the dipole function will return a Python array of 3D vectors with the dipole for each step. You can get an array with the time for each step using the following command:
times = pinq.results.real_time.time()
You should be able to pass this information to a plotting library.
Since we perturbed our system with an electric field, the normalized Fourier transform of the time-dependent dipole function is the polarizability. (To normalize it you have to divide by \(\kappa\), the intensity of the kick.) In this particular case we will calculate \(\alpha_{11}\), the \(xx\) component of the polarizabilty tensor, as we perturbed the system in the \(x\) direction and we are looking at the \(x\) component of the dipole. To calculate the full tensor we would need 3 different time-propagations with perturbations in each direction.
You could use a third party program to calculate the Fourier transform, however note you need to add a window function to dump the dipole for long times (this is equivalent to adding a broadening in Fourier space). Conveniently, INQ has functionality to do this calculation using the command inq spectrum. This is how we can use it from the shell:
inq results real time dipole > dipole.dat
inq spectrum file dipole.dat > spectrum.dat
Unfortunately the python interface of this function is still under development so we recommend you to use the shell version for now.
The file spectrum.dat contains 7 columns, the first one is the frequency (in Hartree) and the next 5 are the real and imaginary parts of the Fourier transforms of the 3 dipole components. The experimentally relevant quantity is
so it would be interesting to plot that. Remember that the response should be properly normalized by the intensity of the perturbation, the spectrum utility doesn’t do it. Also remember that we only have part of the full polarizability tensor.
The Thomas-Reiche-Kuhn \(f\)-sum rule for the number of electrons, states that
where \(N\) is the number of (valence) electrons. This is a good check to validate your results. Can you write a small code to calculate it?