Skip to content

Commit 6c88e3e

Browse files
ifilotifilot
committed
Modifying examples and adding interface to PyMoDia
1 parent f453a80 commit 6c88e3e

13 files changed

Lines changed: 145 additions & 102 deletions

‎examples/.gitignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
11
*.txt
22
*.pdf
33
*.png
4+
*.svg

‎examples/benzene.py‎

Lines changed: 4 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,8 @@
11
# -*- coding: utf-8 -*-
2-
import os,sys
3-
4-
# add a reference to load the module
5-
ROOT = os.path.dirname(__file__)
6-
sys.path.insert(1, os.path.join(ROOT, '..'))
7-
82
from pydft import MoleculeBuilder, DFT
93

10-
#
11-
# Example: Calculate total electronic energy for CO using standard
12-
# settings.
13-
#
14-
15-
CO = MoleculeBuilder().from_name("benzene")
16-
dft = DFT(CO, basis='sto3g')
17-
en = dft.scf(1e-4, verbose=True)
18-
print("Total electronic energy: %f Ht" % en)
4+
mol = MoleculeBuilder().from_name("benzene")
5+
dft = DFT(mol, basis='sto3g')
6+
res = dft.scf(1e-4, verbose=True)
7+
print("Total electronic energy: %f Ht" % res['energy'])
198
dft.print_time_statistics()

‎examples/benzene_atomic_cells.py‎

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,6 @@
11
# -*- coding: utf-8 -*-
22
import numpy as np
33
import matplotlib.pyplot as plt
4-
import os, sys
5-
6-
# add a reference to load the module
7-
ROOT = os.path.dirname(__file__)
8-
sys.path.insert(1, os.path.join(ROOT, '..'))
9-
104
from pydft import MoleculeBuilder, MolecularGrid
115

126
def main():

‎examples/bonding_analysis_co.py‎

Lines changed: 2 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,9 @@
11
import matplotlib.pyplot as plt
22
import numpy as np
3-
from mpl_toolkits.axes_grid1 import make_axes_locatable
4-
import os,sys
5-
6-
# add a reference to load the module
7-
ROOT = os.path.dirname(__file__)
8-
sys.path.insert(1, os.path.join(ROOT, '..'))
9-
3+
from mpl_toolkits.axes_grid import make_axes_locatable
104
from pydft import MoleculeBuilder, DFT
115

12-
mol_builder = MoleculeBuilder()
13-
mol = mol_builder.from_name('co')
6+
mol = MoleculeBuilder().from_name('co')
147

158
# construct dft object
169
dft = DFT(mol, basis='sto3g')

‎examples/co.py‎

Lines changed: 2 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,8 @@
11
# -*- coding: utf-8 -*-
2-
import os,sys
3-
4-
# add a reference to load the module
5-
ROOT = os.path.dirname(__file__)
6-
sys.path.insert(1, os.path.join(ROOT, '..'))
7-
82
from pydft import MoleculeBuilder, DFT
93

10-
#
11-
# Example: Calculate total electronic energy for CO using standard
12-
# settings.
13-
#
14-
15-
CO = MoleculeBuilder().from_name("CO")
16-
dft = DFT(CO, basis='sto3g')
4+
mol = MoleculeBuilder().from_name("CO")
5+
dft = DFT(mol, basis='sto3g')
176
res = dft.scf(1e-4, verbose=True)
187
print("Total electronic energy: %f Ht" % res['energy'])
198
dft.print_time_statistics()

‎examples/co_fb.py‎

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
# -*- coding: utf-8 -*-
2+
from pydft import MoleculeBuilder, DFT
3+
from pyqint import FosterBoys, PyQInt
4+
import matplotlib.pyplot as plt
5+
import numpy as np
6+
import os
7+
8+
def main():
9+
mol = MoleculeBuilder().from_name("CO")
10+
dft = DFT(mol, basis='sto3g')
11+
res = dft.scf(1e-4, verbose=False)
12+
13+
cgfs = res['cgfs']
14+
orbe = res['orbe']
15+
coeff = res['orbc']
16+
17+
# visualize orbitals
18+
fig, ax = plt.subplots(2,5, figsize=(12, 5), dpi=300)
19+
sz = 3
20+
for i in range(0,2):
21+
for j in range(0,5):
22+
dens = plot_wavefunction(cgfs, coeff[:,i*5+j], sz=sz)
23+
limit = max(abs(np.min(dens)), abs(np.max(dens)) )
24+
im = ax[i,j].contourf(dens, origin='lower',
25+
extent=[-sz, sz, -sz, sz], cmap='PiYG', vmin=-limit, vmax=limit,
26+
levels=15)
27+
im = ax[i,j].contour(dens, origin='lower', colors='black',
28+
extent=[-sz, sz, -sz, sz], vmin=-limit, vmax=limit,
29+
levels=15)
30+
ax[i,j].set_xlabel('x [Bohr]')
31+
ax[i,j].set_ylabel('z [Bohr]')
32+
ax[i,j].set_aspect('equal', adjustable='box')
33+
ax[i,j].set_xticks(np.linspace(-3,3, 7))
34+
ax[i,j].set_yticks(np.linspace(-3,3, 7))
35+
ax[i,j].grid(linestyle='--', alpha=0.5)
36+
ax[i,j].set_title(r'$\psi_{%i}$ ($\epsilon_{%i}$ = %6.4f Ht)' %
37+
(i*5+j+1, i*5+j+1, orbe[i*5+j]))
38+
plt.tight_layout()
39+
plt.savefig(os.path.join(os.path.dirname(__file__), "co_contour_canonical.png"))
40+
41+
resfb = FosterBoys(res).run()
42+
cgfs = resfb['cgfs']
43+
orbe = resfb['orbe']
44+
coeff = resfb['orbc']
45+
46+
# visualize orbitals
47+
fig, ax = plt.subplots(2,5, figsize=(12, 5), dpi=300)
48+
sz = 3
49+
for i in range(0,2):
50+
for j in range(0,5):
51+
dens = plot_wavefunction(cgfs, coeff[:,i*5+j], sz=sz)
52+
limit = max(abs(np.min(dens)), abs(np.max(dens)) )
53+
im = ax[i,j].contourf(dens, origin='lower',
54+
extent=[-sz, sz, -sz, sz], cmap='PiYG', vmin=-limit, vmax=limit,
55+
levels=15)
56+
im = ax[i,j].contour(dens, origin='lower', colors='black',
57+
extent=[-sz, sz, -sz, sz], vmin=-limit, vmax=limit,
58+
levels=15)
59+
ax[i,j].set_xlabel('x [Bohr]')
60+
ax[i,j].set_ylabel('z [Bohr]')
61+
ax[i,j].set_aspect('equal', adjustable='box')
62+
ax[i,j].set_xticks(np.linspace(-3,3, 7))
63+
ax[i,j].set_yticks(np.linspace(-3,3, 7))
64+
ax[i,j].grid(linestyle='--', alpha=0.5)
65+
ax[i,j].set_title(r'$\psi_{%i}$ ($\epsilon_{%i}$ = %6.4f Ht)' %
66+
(i*5+j+1, i*5+j+1, orbe[i*5+j]))
67+
plt.tight_layout()
68+
plt.savefig(os.path.join(os.path.dirname(__file__), "co_contour_localized.png"))
69+
70+
def plot_wavefunction(cgfs, coeff, sz=3.5):
71+
# build integrator
72+
integrator = PyQInt()
73+
74+
# build grid
75+
x = np.linspace(-sz, sz, 251)
76+
z = np.linspace(-sz, sz, 251)
77+
xx, zz = np.meshgrid(x,z)
78+
yy = np.zeros(len(x) * len(z))
79+
grid = np.vstack([xx.flatten(), yy, zz.flatten()]).reshape(3,-1).T
80+
res = integrator.plot_wavefunction(grid, coeff, cgfs).reshape((len(z), len(x)))
81+
82+
return res
83+
84+
if __name__ == '__main__':
85+
main()

‎examples/co_mo.py‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
import os
2+
from pymodia import MoDia, MoDiaData, autobuild_from_pyqint, MoDiaSettings
3+
from pydft import MoleculeBuilder, DFT
4+
5+
# Perform PyQInt calculations for CO and its localization
6+
mol = MoleculeBuilder().from_name('co')
7+
res = DFT(mol, 'sto3g').scf()
8+
9+
# adjust settings
10+
settings = MoDiaSettings()
11+
settings.orbc_color = '#555555'
12+
settings.arrow_color = '#CC0000'
13+
settings.core_cutoff = -9
14+
15+
# attempt to automatically create mol and fragments from calculation
16+
mol, f1, f2 = autobuild_from_pyqint(res, name='co')
17+
18+
# we make here a small adjustment to the height of the 5σ orbital to avoid
19+
# overlap with the 2x2π MO
20+
moe = res['orbe']
21+
moe[6] += 0.1
22+
23+
# build data object
24+
data = MoDiaData(mol, f1, f2)
25+
data.set_moe(moe)
26+
27+
diagram = MoDia(data, draw_level_labels=True, level_labels_style='mo_ao',
28+
mo_labels=['1σ', '2σ', '3σ', '4σ', '1π', '1π', '5σ', '2π', '2π', '6σ'],
29+
settings=settings)
30+
diagram.export_svg(os.path.join(os.path.dirname(__file__), "mo_co_canonical.svg"))

‎examples/co_screening.py‎

Lines changed: 5 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,31 +1,22 @@
11
# -*- coding: utf-8 -*-
2-
import os,sys
32
import numpy as np
43
import time
5-
6-
# add a reference to load the module
7-
ROOT = os.path.dirname(__file__)
8-
sys.path.insert(1, os.path.join(ROOT, '..'))
9-
104
from pydft import MoleculeBuilder, DFT
115

12-
#
13-
# Example: Calculate total electronic energy for CO using standard
14-
# settings.
15-
#
16-
17-
CO = MoleculeBuilder().from_name("CO")
6+
mol = MoleculeBuilder().from_name("CO")
187

198
rshells = [8,16,32,64,128]
20-
ang = [6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810]
9+
ang = [6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, 302, 350, 434,
10+
590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890,
11+
4334, 4802, 5294, 5810]
2112

2213
energies = np.zeros((len(rshells), len(ang)))
2314
timestats = np.zeros((len(rshells), len(ang)))
2415

2516
for i,nr in enumerate(rshells):
2617
for j,a in enumerate(ang):
2718
start_time = time.time()
28-
dft = DFT(CO, basis='sto3g', nshells=nr, nangpts=a)
19+
dft = DFT(mol, basis='sto3g', nshells=nr, nangpts=a)
2920
en = dft.scf(1e-4)
3021
end_time = time.time()
3122
elapsed_time = end_time - start_time

‎examples/co_xc.py‎

Lines changed: 7 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,12 @@
11
# -*- coding: utf-8 -*-
22
from pydft import MoleculeBuilder,DFT
33

4-
#
5-
# Example: Calculate total electronic energy for CO using different
6-
# XC-functionals.
7-
#
4+
mol = MoleculeBuilder().from_name("CO")
85

9-
CO = MoleculeBuilder().from_name("CO")
10-
dft = DFT(CO, basis='sto3g', functional='svwn5')
11-
en = dft.scf(1e-4)
12-
print("Total electronic energy (SVWN5): %f Ht" % en)
6+
dft = DFT(mol, basis='sto3g', functional='svwn5')
7+
res = dft.scf(1e-4)
8+
print("Total electronic energy (SVWN5): %f Ht" % res['energy'])
139

14-
CO = MoleculeBuilder().from_name("CO")
15-
dft = DFT(CO, basis='sto3g', functional='pbe')
16-
en = dft.scf(1e-4)
17-
print("Total electronic energy (PBE): %f Ht" % en)
10+
dft = DFT(mol, basis='sto3g', functional='pbe')
11+
res = dft.scf(1e-4)
12+
print("Total electronic energy (PBE): %f Ht" % res['energy'])

‎examples/dft_vs_hf.py‎

Lines changed: 0 additions & 16 deletions
This file was deleted.

0 commit comments

Comments
 (0)