1+ # -*- coding: utf-8 -*-
2+
3+ from pyqint import MoleculeBuilder , HF , PyQInt
4+ import matplotlib .pyplot as plt
5+ import numpy as np
6+ from copy import deepcopy
7+
8+ def main ():
9+ mol = MoleculeBuilder ().from_name ('ethylene' )
10+ cgfs , nuclei = mol .build_basis ('sto3g' )
11+ from pyqint import CGF , GTO
12+ res = HF ().rhf (mol , 'sto3g' )
13+
14+ # build transformation matrix that casts the original basis onto a
15+ # symmetry adapted basis; no normalization is applied
16+ B = np .zeros ((14 ,14 ))
17+ for i in range (0 ,5 ):
18+ B [i * 2 ,i ] = 1
19+ B [i * 2 ,i + 7 ] = 1
20+ B [i * 2 + 1 ,i ] = 1
21+ B [i * 2 + 1 ,i + 7 ] = - 1
22+
23+ B [10 ,5 ] = 1
24+ B [10 ,6 ] = 1
25+ B [10 ,- 2 ] = 1
26+ B [10 ,- 1 ] = 1
27+
28+ B [11 ,5 ] = 1
29+ B [11 ,6 ] = - 1
30+ B [11 ,- 2 ] = 1
31+ B [11 ,- 1 ] = - 1
32+
33+ B [12 ,5 ] = 1
34+ B [12 ,6 ] = 1
35+ B [12 ,- 2 ] = - 1
36+ B [12 ,- 1 ] = - 1
37+
38+ B [13 ,5 ] = 1
39+ B [13 ,6 ] = - 1
40+ B [13 ,- 2 ] = - 1
41+ B [13 ,- 1 ] = 1
42+
43+ # build yet another transformation basis that re-orders the functions such
44+ # that all irreps belonging to the same symmetry group are consecutive; this
45+ # yields block-diagonal matrices
46+ order = []
47+ overlap = B @ res ['overlap' ] @ B .T
48+ for i in range (len (overlap )):
49+ for j in range (len (overlap )):
50+ if abs (overlap [i ,j ]) > 0.1 :
51+ if j not in order :
52+ order .append (j )
53+ n = len (order )
54+ P = np .zeros ((n , n ), dtype = int )
55+ P [np .arange (n ), order ] = 1
56+
57+ symfuncs = {
58+ 'A$_{g}$' : 4 ,
59+ 'B$_{3u}$' : 4 ,
60+ 'B$_{2g}$' : 1 ,
61+ 'B$_{1u}$' : 1 ,
62+ 'B$_{1g}$' : 2 ,
63+ 'B$_{2u}$' : 2
64+ }
65+ symlabels = []
66+ for k ,v in symfuncs .items ():
67+ for i in range (v ):
68+ symlabels .append (k + '(%i)' % (i + 1 ))
69+ overlap = P @ overlap @ P .T
70+
71+ basislabels = []
72+ counts = {}
73+ for n in nuclei :
74+ if n [1 ] not in counts .keys ():
75+ counts [n [1 ]] = 1
76+ else :
77+ counts [n [1 ]] += 1
78+
79+ if n [1 ] == 6 :
80+ for s in ['1s' , '2s' , '2p_{x}' , '2p_{y}' , '2p_{z}' ]:
81+ basislabels .append ('C$_{%s}^{(%i)}$' % (s , counts [n [1 ]]))
82+ else :
83+ basislabels .append ('H$_{1s}^{(%i)}$' % (counts [n [1 ]]))
84+
85+ overlap = P @ overlap @ P .T
86+
87+ fig , ax = plt .subplots (1 ,1 , dpi = 144 , figsize = (7 ,7 ))
88+ plot_matrix (ax , P @ B , basislabels , symlabels )
89+ plt .show ()
90+
91+ # construct new basis
92+ B = P @ B
93+ cgfs_symad = [CGF () for i in range (len (B ))]
94+ for i in range (len (B )): # loop over new basis functions
95+ for j in range (len (cgfs )): # loop over old basis functions
96+ if abs (B [i ,j ]) > 0.01 : # verify non-negligble contribution
97+ for g in cgfs [j ].gtos :
98+ cgfs_symad [i ].gtos .append (GTO (g .c * B [i ,j ], g .p , g .alpha , g .l , g .m , g .n ))
99+
100+ # re-perform Hartree-Fock calculation using the symmetry adapted basis
101+ res = HF ().rhf (mol , cgfs_symad , verbose = True )
102+ fig , ax = plt .subplots (1 ,1 , dpi = 144 , figsize = (7 ,7 ))
103+ plot_matrix (ax , res ['overlap' ], symlabels , symlabels )
104+ plt .show ()
105+
106+ fig , ax = plt .subplots (1 ,1 , dpi = 144 , figsize = (7 ,7 ))
107+ plot_matrix (ax , res ['fock' ], symlabels , symlabels )
108+ plt .show ()
109+
110+ def plot_matrix (ax , mat , xlabels , ylabels , xlabelrot = 0 ):
111+ """
112+ Produce plot of matrix
113+ """
114+ mv = np .max (np .abs (mat ))
115+ ax .imshow (mat , vmin = - mv , vmax = mv , cmap = 'PiYG' )
116+ for i in range (mat .shape [0 ]):
117+ for j in range (mat .shape [1 ]):
118+ ax .text (i , j , '%.2f' % mat [j ,i ], ha = 'center' , va = 'center' ,
119+ fontsize = 7 , color = 'black' if abs (mat [j ,i ]) < (0.7 * mv ) else 'white' )
120+ ax .set_xticks ([])
121+ ax .set_yticks ([])
122+ ax .hlines (np .arange (1 , mat .shape [0 ])- 0.5 , - 0.5 , mat .shape [0 ] - 0.5 ,
123+ color = 'black' , linestyle = '--' , linewidth = 1 )
124+ ax .vlines (np .arange (1 , mat .shape [0 ])- 0.5 , - 0.5 , mat .shape [0 ] - 0.5 ,
125+ color = 'black' , linestyle = '--' , linewidth = 1 )
126+
127+ # add basis functions as axes labels
128+ ax .set_xticks (np .arange (0 , mat .shape [0 ]))
129+ ax .set_xticklabels (xlabels , rotation = xlabelrot )
130+ ax .set_yticks (np .arange (0 , mat .shape [0 ]))
131+ ax .set_yticklabels (ylabels , rotation = 0 )
132+ ax .tick_params (axis = 'both' , which = 'major' , labelsize = 7 )
133+
134+ if __name__ == '__main__' :
135+ main ()
0 commit comments