1.1.1 Cartesian Phonon Calculation: Displacing the Atoms
Here we will learn to compute the dynamical matrix and normal-mode / phonon frequencies of a CO2 molecule using pyEPFD and qbox. The first step is to optimize the geometry using qbox. This could be done using a qbox alone. The optimized structure must be saved as a XYZ file. Here we used B3LYP functional to optimize the geometry. We must use the same functional to compute the normal-mode frequencies. Please consult the qbox documentation: http://qboxcode.org/doc/html/usage/intro.html for more details on geometry optimization.
Let us first look at the XYZ file of the optimize geometry.
[1]:
opt_xyz_file = open("co2_b3lyp_opt.xyz","r").read()
print(opt_xyz_file)
3
## CELL(abcABC): 20.000000 20.000000 20.000000 90.0000 90.0000 90.0000 Step: 000 Bead: 0 positions{angstrom} cell{atomic_unit}
O 0.000000 0.000000 1.157132295528
O 0.000000 0.000000 -1.157130707997
C 0.000000 0.000000 -1.587531e-06
The second line of the file defines the cell-parameters: cell lengths and angles. Here cell lengths are in Bohr, however they can be given also in angstrom. In that case one has to use the directive cell{angstrom}. This style is adapted from ipi code: https://ipi-code.org/.
The next step is to perform the finite displacement (fd) moves along cartesian coordinates. For that purpose first we would import coordinate utilities tools from pyepfd.
[2]:
from pyepfd.coord_util import *
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.0
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2023-05-09 19:45:59.045534
Now we will read the given optimized coordinate file using xyz class available in pyepfd.coord_util and store the information in opt_xyz object. Here io = ‘r’ means xyz class is set to reading mode. We will store the cell-parameters and cartesian cordinates in cell & opt_coord objects, respectively.
[3]:
opt_xyz = xyz(file_path = 'co2_b3lyp_opt.xyz', io = 'r')
cell = opt_xyz.cell[0]
opt_coord = opt_xyz.coords[0]
Time spent on xyz class: 0.0001537799835205078 s.
Now we will use the finite displacement moves using the fdmoves tools. The mode = ‘FD’ means we will use finite displacement along cartesian and deltax = 0.005 means we will move each atoms by +0.005 and -0.005 Bohr.
[4]:
fdmoves = ionic_mover( atoms = opt_xyz.atoms, \
opt_coord = opt_coord,\
mode = 'FD', deltax = 0.005)
Time spent on ionic_mover class: 0.0003070831298828125 s.
After this, all the displaced coordinates are saved onto fdmoves object. Now we save those coordinates into an xyz file using the xyz tool again but changing the io = ‘w’ which tell xyz class would be on writing mode.
[5]:
out_qb_inp = xyz(file_path = 'fd_phonon.xyz',\
io ='w', atoms = opt_xyz.atoms)
ndisp = fdmoves.disp_coord.shape[1] #Total number of displacements
### Writing each displaced configuration using a for loop
for i in range(ndisp):
out_qb_inp.write(cell = cell,coord = fdmoves.disp_coord[:,i])
Time spent on xyz class: 0.0001811981201171875 s.
This will create an xyz file with total 19 displaced coordinates (6N+1, where N is the number of atoms) where the first configuration is the geometry optimized configuration.
[6]:
# Deleting the output object to finish the file writing process
# so that we can view the xyz file in the next step
del out_qb_inp
[7]:
%%bash
cat fd_phonon.xyz
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 0 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0 0 -1.1571307
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 1 positions{angstrom} cell{angstrom}
O 0.0026458862 0 1.1571323
O 0 0 -1.1571307
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 2 positions{angstrom} cell{angstrom}
O -0.0026458862 0 1.1571323
O 0 0 -1.1571307
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 3 positions{angstrom} cell{angstrom}
O 0 0.0026458862 1.1571323
O 0 0 -1.1571307
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 4 positions{angstrom} cell{angstrom}
O 0 -0.0026458862 1.1571323
O 0 0 -1.1571307
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 5 positions{angstrom} cell{angstrom}
O 0 0 1.1597782
O 0 0 -1.1571307
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 6 positions{angstrom} cell{angstrom}
O 0 0 1.1544864
O 0 0 -1.1571307
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 7 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0.0026458862 0 -1.1571307
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 8 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O -0.0026458862 0 -1.1571307
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 9 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0 0.0026458862 -1.1571307
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 10 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0 -0.0026458862 -1.1571307
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 11 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0 0 -1.1544848
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 12 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0 0 -1.1597766
C 0 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 13 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0 0 -1.1571307
C 0.0026458862 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 14 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0 0 -1.1571307
C -0.0026458862 0 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 15 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0 0 -1.1571307
C 0 0.0026458862 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 16 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0 0 -1.1571307
C 0 -0.0026458862 -1.587531e-06
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 17 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0 0 -1.1571307
C 0 0 0.0026442987
3
# CELL(abcABC): 10.58354 10.58354 10.58354 90.00000 90.00000 90.00000 PyEPFD-Step: 18 positions{angstrom} cell{angstrom}
O 0 0 1.1571323
O 0 0 -1.1571307
C 0 0 -0.0026474738
Next step is to convert this xyz file into a series of Qbox input files, that would be the topic of tutorial 1.1.2.