Steps to ensure compilation of atm.f¶
Ensure that you are in QXMD_Course/src/atm folder
Way to ensure this is
pwdin your current folder and look for QXMD_Course/src/atm as relative addressRun
make ifortfrom this directory
In [5]:
! make ifort
Compilation of atm.f continued¶
- Compile the atm.f code by executing
make atmcommand. - This will compile the code with all it's dependencies and move the
./atmexecutable at root level
In [6]:
! make atm
Edit in7.dat if necessary (can be done within the jupyter notebook)¶
- Define the atom by the atomic number
- Change the valence electron distribution
- Change reference energies and cutoff radii for pseudo-wavefuncions
Executing atm executable¶
- This code reads in7.dat and constructs the PAW potential
- Run the exectuable by typing
./atm
In [7]:
! ./atm
Compare all electron and pseudo wavefunctions¶
In [8]:
import numpy as np
import matplotlib.pyplot as plt
# figure size (inch)
figx, figy = 16, 4
# y range
dymin, dymax = -0.0000015, 0.000004
symin, symax = -0.006, 0.006
pymin, pymax = -0.00002, 0.00005
# cutoff
dcut = 2
scut = 2.6
pcut = 2.7
dae0, dae1, dae2 = [], [], []
dps0, dps1, dps2 = [], [], []
sae0, sae1, sae2 = [], [], []
sps0, sps1, sps2 = [], [], []
pae0, pae1, pae2 = [], [], []
pps0, pps1, pps2 = [], [], []
##### 5d orbital
# AE wavefunctions
ffile1 = open('D_Pae.dat', 'r')
# PS wavefunctions
ffile2 = open('D_Pus.dat', 'r')
lines1 = ffile1.readlines()
lines2 = ffile2.readlines()
for line in lines1:
dae0.append(float(line.split()[0]))
dae1.append(float(line.split()[1]))
dae2.append(float(line.split()[2]))
for line in lines2:
dps0.append(float(line.split()[0]))
dps1.append(float(line.split()[1]))
dps2.append(float(line.split()[2]))
ffile1.close()
ffile2.close()
##### 6s orbital
# AE wavefunctions
ffile1 = open('S_Pae.dat', 'r')
# PS wavefunctions
ffile2 = open('S_Pus.dat', 'r')
lines1 = ffile1.readlines()
lines2 = ffile2.readlines()
for line in lines1:
sae0.append(float(line.split()[0]))
sae1.append(float(line.split()[1]))
sae2.append(float(line.split()[2]))
for line in lines2:
sps0.append(float(line.split()[0]))
sps1.append(float(line.split()[1]))
sps2.append(float(line.split()[2]))
ffile1.close()
ffile2.close()
##### 6p orbital
# AE wavefunctions
ffile1 = open('P_Pae.dat', 'r')
# PS wavefunctions
ffile2 = open('P_Pus.dat', 'r')
lines1 = ffile1.readlines()
lines2 = ffile2.readlines()
for line in lines1:
pae0.append(float(line.split()[0]))
pae1.append(float(line.split()[1]))
pae2.append(float(line.split()[2]))
for line in lines2:
pps0.append(float(line.split()[0]))
pps1.append(float(line.split()[1]))
pps2.append(float(line.split()[2]))
ffile1.close()
ffile2.close()
fig = plt.figure(figsize=(figx,figy))
# d-orbital
p1 = fig.add_subplot(131)
p1.plot(dps0, dps1, "blue")
p1.plot(dps0, dps2, "red")
p1.plot(dae0, dae1, "blue", linestyle='dashed')
p1.plot(dae0, dae2, "red", linestyle='dashed')
p1.vlines([dcut], dymin, dymax, linestyle='dotted')
p1.set_ylim(dymin, dymax)
p1.ticklabel_format(style='sci',axis='y',scilimits=(0,0))
p1.set_ylabel('Wavefunction', size=15)
# p1.set_xlabel('Radius (bohr)', size=15)
p1.legend(('PP_E=ae', 'PP_E=ref', 'AE_E=ae', 'AE_E=ref'), loc='upper right')
p1.set_title('5d orbital', size = 15)
# s-orbital
p2 = fig.add_subplot(132)
p2.plot(sps0, sps1, "blue")
p2.plot(sps0, sps2, "red")
p2.plot(sae0, sae1, "blue", linestyle='dashed')
p2.plot(sae0, sae2, "red", linestyle='dashed')
p2.vlines([scut], symin, symax, linestyle='dotted')
p2.set_ylim(symin, symax)
p2.ticklabel_format(style='sci',axis='y',scilimits=(0,0))
# p2.set_ylabel('Wavefunction', size=15)
p2.set_xlabel('Radius (bohr)', size=15)
p2.set_title('6s orbital', size = 15)
# p-obital
p3 = fig.add_subplot(133)
p3.plot(pps0, pps1, "blue")
p3.plot(pps0, pps2, "red")
p3.plot(pae0, pae1, "blue", linestyle='dashed')
p3.plot(pae0, pae2, "red", linestyle='dashed')
p3.vlines([pcut], pymin, pymax, linestyle='dotted')
p3.set_ylim(pymin, pymax)
p3.ticklabel_format(style='sci',axis='y',scilimits=(0,0))
# p3.set_ylabel('Wavefunction', size=15)
# p3.set_xlabel('Radius (bohr)', size=15)
p3.set_title('6p orbital', size = 15)
plt.show()
Compare logarithmic derivatives of all electron and pseudo wavefunctions¶
In [9]:
import numpy as np
import matplotlib.pyplot as plt
# figure size (inch)
figx, figy = 16, 4
# y range
ymin, ymax = -4, 4
dae0, dae1, dae2 = [], [], []
dps0, dps1, dps2 = [], [], []
sae0, sae1, sae2 = [], [], []
sps0, sps1, sps2 = [], [], []
pae0, pae1, pae2 = [], [], []
pps0, pps1, pps2 = [], [], []
##### 5d orbital
# AE wavefunctions
ffile1 = open('D_chiae.dat', 'r')
# PS wavefunctions
ffile2 = open('D_chil.dat', 'r')
lines1 = ffile1.readlines()
lines2 = ffile2.readlines()
for line in lines1:
dae0.append(float(line.split()[0]))
dae1.append(float(line.split()[1]))
for line in lines2:
dps0.append(float(line.split()[0]))
dps1.append(float(line.split()[1]))
ffile1.close()
ffile2.close()
##### 6s orbital
# AE wavefunctions
ffile1 = open('S_chiae.dat', 'r')
# PS wavefunctions
ffile2 = open('S_chil.dat', 'r')
lines1 = ffile1.readlines()
lines2 = ffile2.readlines()
for line in lines1:
sae0.append(float(line.split()[0]))
sae1.append(float(line.split()[1]))
for line in lines2:
sps0.append(float(line.split()[0]))
sps1.append(float(line.split()[1]))
ffile1.close()
ffile2.close()
##### 6p orbital
# AE wavefunctions
ffile1 = open('P_chiae.dat', 'r')
# PS wavefunctions
ffile2 = open('P_chil.dat', 'r')
lines1 = ffile1.readlines()
lines2 = ffile2.readlines()
for line in lines1:
pae0.append(float(line.split()[0]))
pae1.append(float(line.split()[1]))
for line in lines2:
pps0.append(float(line.split()[0]))
pps1.append(float(line.split()[1]))
ffile1.close()
ffile2.close()
fig = plt.figure(figsize=(figx,figy))
# d-orbital
p1 = fig.add_subplot(131)
p1.plot(dps0, dps1, "blue")
p1.plot(dae0, dae1, "red", linestyle='dashed')
p1.set_ylim(ymin, ymax)
p1.set_ylabel('Logarithmic derivative', size=15)
# p1.set_xlabel('Energy (Ry)', size=15)
p1.legend(('PP', 'AE'), loc='upper right')
p1.set_title('5d orbital', size = 15)
# s-orbital
p2 = fig.add_subplot(132)
p2.plot(sps0, sps1, "blue")
p2.plot(sae0, sae1, "red", linestyle='dashed')
p2.set_ylim(ymin, ymax)
# p2.set_ylabel('Logarithmic derivative', size=15)
p2.set_xlabel('Energy (Ry)', size=15)
p2.set_title('6s orbital', size = 15)
# p-obital
p3 = fig.add_subplot(133)
p3.plot(pps0, pps1, "blue")
p3.plot(pae0, pae1, "red", linestyle='dashed')
p3.set_ylim(ymin, ymax)
# p3.set_ylabel('Logarithmic derivative', size=15)
# p3.set_xlabel('Energy (Ry)', size=15)
p3.set_title('6p orbital', size = 15)
plt.show()
Estimation of plane-wave cutoff energies for wavefunctions¶
In [10]:
import numpy as np
import matplotlib.pyplot as plt
xmin, xmax = 0, 50
ymin, ymax = 0, 0.002
d0, d1 = [], []
s0, s1 = [], []
p0, p1 = [], []
ffile1 = open('D_delE.dat', 'r')
ffile2 = open('S_delE.dat', 'r')
ffile3 = open('P_delE.dat', 'r')
lines1 = ffile1.readlines()
lines2 = ffile2.readlines()
lines3 = ffile3.readlines()
for line in lines1:
d0.append(float(line.split()[0]))
d1.append(float(line.split()[1]))
for line in lines2:
s0.append(float(line.split()[0]))
s1.append(float(line.split()[1]))
for line in lines3:
p0.append(float(line.split()[0]))
p1.append(float(line.split()[1]))
fig = plt.gcf()
fig.set_size_inches(6.5, 4.5)
p1 = plt.plot(d0, d1, s0, s1, p0, p1)
p1 = plt.hlines([0.001], xmin, xmax, linestyle='dotted')
plt.xlim(xmin, xmax)
plt.ylim(ymin, ymax)
plt.ticklabel_format(style='sci',axis='y',scilimits=(0,0))
plt.ylabel('Relative energy (Ry)', size=15)
plt.xlabel('Energy (Ry)', size=15)
plt.legend(('d orbital', 's orbital', 'p orbital'), loc='upper right')
plt.show()
Estimation of plane-wave cutoff energies for electron densities¶
In [11]:
import numpy as np
import matplotlib.pyplot as plt
# figure size (inch)
figx, figy = 16, 4
sq0, sq1 = [], []
pq0, pq1, pq2 = [], [], []
dq0, dq1, dq2, dq3 = [], [], [], []
ffile1 = open('D_Qbar.dat', 'r')
ffile2 = open('S_Qbar.dat', 'r')
ffile3 = open('P_Qbar.dat', 'r')
lines1 = ffile1.readlines()
lines2 = ffile2.readlines()
lines3 = ffile3.readlines()
for line in lines1:
dq0.append(float(line.split()[0]))
dq1.append(float(line.split()[1]))
dq2.append(float(line.split()[2]))
dq3.append(float(line.split()[3]))
for line in lines2:
sq0.append(float(line.split()[0]))
sq1.append(float(line.split()[1]))
for line in lines3:
pq0.append(float(line.split()[0]))
pq1.append(float(line.split()[1]))
pq2.append(float(line.split()[2]))
ffile1.close()
ffile2.close()
ffile3.close()
fig = plt.figure(figsize=(figx,figy))
# d-orbital
p1 = fig.add_subplot(131)
p1.plot(dq0, dq1, "blue")
p1.plot(dq0, dq2, "red")
p1.plot(dq0, dq3, "orange")
p1.ticklabel_format(style='sci',axis='y',scilimits=(0,0))
p1.set_ylabel('Fourier compornets', size=15)
# p1.set_xlabel('Energy (Ry)', size=15)
p1.legend(('L = 0', 'L = 2', 'L = 4'), loc='upper right')
p1.set_title('5d orbital', size = 15)
# s-orbital
p2 = fig.add_subplot(132)
p2.plot(sq0, sq1, "blue")
p2.ticklabel_format(style='sci',axis='y',scilimits=(0,0))
# p2.set_ylabel('Fourier compornets', size=15)
p2.set_xlabel('Energy (Ry)', size=15)
p2.set_title('6s orbital', size = 15)
# p-obital
p3 = fig.add_subplot(133)
p3.plot(pq0, pq1, "blue")
p3.plot(pq0, pq2, "red")
p3.ticklabel_format(style='sci',axis='y',scilimits=(0,0))
# p3.set_ylabel('Fourier compornets', size=15)
# p3.set_xlabel('Energy (Ry)', size=15)
p3.set_title('6p orbital', size = 15)
plt.show()
