pyfock.Basis
1# Basis.py 2# Author: Manas Sharma ([email protected]) 3# This is a part of CrysX (https://bragitoff.com/crysx) 4# 5# 6# .d8888b. Y88b d88P 8888888b. 8888888888 888 7# d88P Y88b Y88b d88P 888 Y88b 888 888 8# 888 888 Y88o88P 888 888 888 888 9# 888 888d888 888 888 .d8888b Y888P 888 d88P 888 888 8888888 .d88b. .d8888b 888 888 10# 888 888P" 888 888 88K d888b 8888888P" 888 888 888 d88""88b d88P" 888 .88P 11# 888 888 888 888 888 "Y8888b. d88888b 888888 888 888 888 888 888 888 888 888888K 12# Y88b d88P 888 Y88b 888 X88 d88P Y88b 888 Y88b 888 888 Y88..88P Y88b. 888 "88b 13# "Y8888P" 888 "Y88888 88888P' d88P Y88b 888 "Y88888 888 "Y88P" "Y8888P 888 888 14# 888 888 15# Y8b d88P Y8b d88P 16# "Y88P" "Y88P" 17import numpy as np 18import scipy 19from scipy.special import factorial2, binom 20import re 21from . import Data 22from collections import Counter 23 24import sys 25import os 26 27# print(__file__) 28# print(sys.argv[0]) 29# print(sys.path[0]) 30 31#Class to store basis function properties 32class Basis: 33 """ 34 Class to store and manage basis function properties for a given molecular system. 35 36 This class processes atomic basis set information for each atom in the molecule, 37 including primitive Gaussian functions, shells, and contracted basis functions (AOs). 38 It also supports TURBOMOLE-style shell ordering. 39 """ 40 def __init__(self, mol, basis, tmoleFormat=False): 41 """ 42 Initialize a Basis object for storing basis function properties. 43 44 Args: 45 mol: Mol object containing molecular information 46 basis: Dictionary specifying the basis set to be used for each atom, or None to use mol.basis 47 tmoleFormat (bool): Whether to use TURBOMOLE format ordering for basis functions 48 49 Returns: 50 None 51 52 Raises: 53 None - prints error message if mol is None 54 """ 55 if mol is None: 56 print('Error: Please provide a Mol object.') 57 return None 58 #Basis set name 59 #If no basis set name is provided specifically, 60 #then use the basis set name of the mol. 61 #'basis' is a dictionary specifying the basis set to be used for a particular atom 62 if basis is None: 63 self.basis = mol.basis 64 else: 65 self.basis = basis 66 67 68 self.basisSet = '' 69 """Basis set name or label. 70 Type: str 71 """ 72 73 self.nao = 0 74 """Number of basis functions (atomic orbitals). 75 Type: int 76 """ 77 78 self.nshells = 0 79 """Number of shells. 80 Type: int 81 """ 82 83 self.totalnprims = 0 84 """Total number of primitive Gaussian functions. 85 Type: int 86 """ 87 88 self.nprims = [] 89 """Number of primitives in each shell. 90 Type: List[int] 91 """ 92 93 self.prim_coeffs = [] 94 """Primitive contraction coefficients. 95 Type: List[float] 96 """ 97 #Exponents of primitives 98 self.prim_expnts = [] 99 """Exponents of primitives. 100 Type: List[float] 101 """ 102 #Coords of primitives 103 self.prim_coords = [] 104 """Coordinates of each primitive function. 105 Type: List[List[float]] 106 """ 107 #Normalization factors of primitives 108 self.prim_norms = [] 109 """Normalization constants of primitives. 110 Type: List[float] 111 """ 112 # This list contains the indices of atoms (first = #0) that correspond to a particular primitive function 113 self.prim_atoms = [] #(This list should be of the size of the number of self.totalnprims) 114 """Atom index corresponding to each primitive. 115 Type: List[int] 116 """ 117 # A list of tuples (2D list) that gives the number of primitives belonging to each atomic index 118 self.nprims_atoms = [] # Size is the same as that of no. of atoms. Looks like this: [(0, 12), (1, 5), (2, 5)] 119 """Number of primitives per atom as (atom_index, nprims). 120 Type: List[Tuple[int, int]] 121 """ 122 # This list contains the indices of shells (first = #0) that correspond to a particular primitive function 123 self.prim_shells = [] #(This list should be of the size of the number of self.totalnprims) 124 """Shell index corresponding to each primitive. 125 Type: List[int] 126 """ 127 # A list of tuples (2D list) that gives the number of primitives belonging to each shell index 128 self.nprims_shells = [] # Size is the same as that of no. of shells. Looks like this: [(0, 12), (1, 5), (2, 5)] 129 """Number of primitives per shell as (shell_index, nprims). 130 Type: List[Tuple[int, int]] 131 """ 132 # A special 2d list that contains the angular momentum of shell in the first index and the no. of primitives corresponding to the shell as second index 133 self.nprims_shell_l_list = [] # This will be of the same size as the no. of shells 134 """Angular momentum and primitive count per shell as (l, nprims). 135 Type: List[Tuple[int, int]] 136 """ 137 # A list of size natoms, that contains the values of the largest primitive exponents for all the atoms 138 self.alpha_max = [] 139 """Maximum exponent among primitives for each atom. 140 Type: List[float] 141 """ 142 # A list of size natoms, that contains the list of the smallest primitive exponents for each shell for all the atoms 143 self.alpha_min = [] 144 """Minimum exponent of primitives per shell for each atom. 145 Type: List[List[float]] 146 """ 147 #Number of contractions 148 self.ncontrs = 0 149 """Total number of contracted basis functions. 150 Type: int 151 """ 152 #Normalization factors for contractions 153 self.contrs_norm = [] 154 """Normalization factors for contracted basis functions. 155 Type: List[float] 156 """ 157 #Degeneracy of shells 158 self.shell_degen = [] 159 """Degeneracy of each shell. 160 Type: List[int] 161 """ 162 #Shells 163 self.shellsLabel = [] 164 """Label of each shell (e.g., 'S', 'P', etc.). 165 Type: List[str] 166 """ 167 self.shells = [] 168 """Shell index (0: s, 1: p, 2:d). 169 Type: List[int] 170 """ 171 self.shell_coords = [] 172 """Coordinates of each shell center. 173 Type: List[List[float]] 174 """ 175 #--------------------------------- 176 #Information for basis function 177 #--------------------------------- 178 #Coordinates of basis functions (This list should be of the size of the number of AOs/BFs) 179 self.bfs_coords = [] 180 """Coordinates of basis functions (same as their parent shell). 181 Type: List[List[float]] 182 """ 183 # Number of primitives corresponding to each basis function (This list should be of the size of the number of AOs/BFs) 184 self.bfs_nprim = [] 185 """Number of primitives corresponding to each basis function. 186 Type: List[int] 187 """ 188 # Radius cutoff for each basis function (This list should be of the size of the number of AOs/BFs) 189 # This is used to accelerate the evaluation of AOs on grid points, as those gridpoints farther than the radius cutoff 190 # of the basis function can be discarded for calcualtion. 191 self.bfs_radius_cutoff = [] 192 """Radius cutoff for each basis function to speed up AO evaluation. 193 Type: List[float] 194 """ 195 #A list of the size of number of atomic orbitals (BFs), 196 #that contains lists of primitive contraction coefficients/exponents 197 #corresponding to every basis function 198 self.bfs_expnts = [] 199 """Primitive exponents used in each basis function. 200 Type: List[List[float]] 201 """ 202 203 self.bfs_coeffs = [] 204 """Primitive coefficients for each basis function. 205 Type: List[List[float]] 206 """ 207 # A list of the size of number of AOs (BFs) 208 # that contains the lists of normalization factors for 209 # every primitive corresponding to the basis functions 210 self.bfs_prim_norms = [] 211 """Primitive normalization constants per basis function. 212 Type: List[List[float]] 213 """ 214 # Normalization factor for the contraction of primitives making up a BF (This list should be of the size of the number of AOs/BFs) 215 self.bfs_contr_prim_norms = [] 216 """Normalization for the contraction of primitives forming a BF. 217 Type: List[float] 218 """ 219 # l,m,n indices of BFs (This list should be of the size of the number of AOs/BFs) 220 self.bfs_lmn = [] 221 """Cartesian angular momentum indices (l, m, n) per basis function. 222 Type: List[Tuple[int, int, int]] 223 """ 224 # Angular momentum of a basis function (This list should be of the size of the number of AOs/BFs) 225 self.bfs_lm = [] 226 """Angular momentum `l` for each basis function. 227 Type: List[int] 228 """ 229 # Label of a basis function (This list should be of the size of the number of AOs/BFs) 230 self.bfs_label = [] 231 """Label for each basis function. 232 Type: List[str] 233 """ 234 # Number of BFs in a shell 235 self.bfs_nbfshell = [] 236 """Number of BFs in the shell of each basis function. 237 Type: List[int] 238 """ 239 # Shell index of each bf 240 self.bfs_shell_index =[] 241 """Shell index for each basis function. 242 Type: List[int] 243 """ 244 # BF offset index of each shell (This should be of the size of the number of shells) 245 self.shell_bfs_offset =[] 246 """Offset index of the first basis function in each shell. 247 Type: List[int] 248 """ 249 # Total number of basis functions (BFs) also called Atomic orbitals (AOs) 250 self.bfs_nao = [] 251 """Index of each basis function (redundant with `nao`). 252 Type: List[int] 253 """ 254 # This list contains the indices of atoms (first = #0) that correspond to a particular bf 255 self.bfs_atoms = [] #(This list should be of the size of the number of AOs/BFs) 256 """Atom index corresponding to each basis function. 257 Type: List[int] 258 """ 259 #--------------------------------- 260 #If all data was parsed successfully. 261 self.success = False 262 """Indicates if basis parsing was successful. 263 Type: bool 264 """ 265 #Convert the keys of basis dictionary to lower case to avoid ambiguity 266 self.basisSet = self.createCompleteBasisSet(mol, basis) 267 self.readBasisFunctionInfo(mol, self.basisSet) 268 self.totalnprims = sum(self.nprims) 269 self.nshells = len(self.nprims) 270 # Now lets create information about the basis functions 271 # A basis function is a combination of gaussian primitives, 272 # so we need 273 # number of primitives for a given bf, 274 # the exponents and contraction coefficients of the primitives, 275 # the l,m,n triplet information, 276 # the angular momentum l+m+n, 277 # coords, 278 # normalization factors of primitives, 279 # normalization factor for the contraction. 280 #--------------------------------------- 281 offset = 0 282 ishell = 0 283 # Loop through all the shells 284 for i in range(len(self.nprims)): 285 # Number of primitives in a shell 286 numprim = self.nprims[i] 287 # Angular momentum of a shell 288 lm = self.shells[i] - 1 289 # no of bfs in this shell 290 nbf_shell = Data.shell_degen[lm] 291 self.bfs_nbfshell.append(nbf_shell) 292 # Assign the lmn triplets and other information corresponding to basis functions 293 # by looping through the degeneracy of the shells. 294 # Example: If the degeneracy is 1 then we loop through [0,1) that is we look 295 # for the 0th element in the Data.shell_lmn dictionary as that corresponds to the 's' lmn triplet. 296 # Similarly if the denegeracy was 3 as for 'p' shell, then we loop from [1,3] and get the 297 # lmn triplets corresponding to 'p' shell and so on. 298 offset_shell_labels = Data.shell_lmn_offset[lm] 299 for j in range(offset_shell_labels, nbf_shell + offset_shell_labels): 300 #Assign the label to the basis function 301 if tmoleFormat: 302 label = list(Data.shell_lmn_tmole)[j] 303 else: 304 label = list(Data.shell_lmn)[j] 305 self.bfs_label.append(label) 306 # Assign l,m,n triplet for a basis function 307 if tmoleFormat: 308 lmn = Data.shell_lmn_tmole[label] 309 else: 310 lmn = Data.shell_lmn[label] 311 self.bfs_lmn.append(lmn) 312 # Assign the number of primitives in a basis function 313 self.bfs_nprim.append(numprim) 314 # Assign the the total angular momentum of a basis function 315 self.bfs_lm.append(sum(lmn)) 316 # Assign the exponents and contraction coefficients 317 expnts = self.prim_expnts[offset : offset+numprim ] 318 self.bfs_expnts.append(expnts) 319 coeffs = self.prim_coeffs[offset : offset+numprim ] 320 self.bfs_coeffs.append(coeffs) 321 # Assign Radius Cutoffs 322 radius_cutoff = 0.0 323 # Loop over primitives 324 eeta_log = np.log(1.0E-9) 325 for iprim_ in range(len(expnts)): 326 radius_cutoff = np.maximum(radius_cutoff, np.sqrt(1.0/expnts[iprim_]*(np.log(expnts[iprim_])/2.0 - eeta_log))) 327 self.bfs_radius_cutoff.append(radius_cutoff) 328 # Assign the normalization factors for all the primitives corresponding to a BF 329 norms = [] 330 for k in range(offset,offset+numprim): 331 norms.append(self.normalizationFactor(self.prim_expnts[k],lmn[0],lmn[1],lmn[2])) 332 self.bfs_prim_norms.append(norms) 333 self.bfs_contr_prim_norms.append(self.normalizationFactorContraction(expnts, coeffs, norms, lmn[0],lmn[1],lmn[2], lm)) 334 # Assign coordinates to the basis functions 335 self.bfs_coords.append(self.shell_coords[i]) 336 # Assign shell index to the basis functions 337 self.bfs_shell_index.append(ishell) 338 339 ishell += 1 340 offset = offset + numprim 341 self.bfs_nao = len(self.bfs_label) 342 343 ##### The following was mainly developed to calculated alpha_min and alpha_max for numgrid 344 345 # Calculate the number of primitives corresponding to each atom 346 cnt = Counter(self.prim_atoms) # Gives a Counter object similar to a dictionary 347 # Convert the Counter object to a list 348 self.nprims_atoms = list(cnt.items()) #a list of tuples, each containing a value and its count in the list 349 ### Start processing to calculate alpha_min and alpha_max for grid generation 350 ## Calculate alpha_max (https://github.com/dftlibs/numgrid) 351 self.alpha_max = [] 352 offset = 0 353 for iat in range(mol.natoms): 354 alpha_max = max(self.prim_expnts[offset : offset + self.nprims_atoms[iat][1]]) 355 self.alpha_max.append(alpha_max) 356 offset = offset + self.nprims_atoms[iat][1] 357 358 # Calculate the number of primitives corresponding to each shell 359 cnt = Counter(self.prim_shells) # Gives a Counter object similar to a dictionary 360 # Convert the Counter object to a list 361 self.nprims_shells = list(cnt.items()) #a list of tuples, each containing a value and its count in the list 362 # Make a list that is pretty much the same as nprims_shells, however, it contains the shell angular momentum in the first index 363 for ishell in range(self.nshells): 364 self.nprims_shell_l_list.append((self.shells[ishell], self.nprims_shells[ishell][1])) 365 # Create self.nprims_angmom_list 366 # self.calc_nprim_for_each_angular_momentum_l() 367 ## Calculate alpha_min 368 ishell_offset = 0 369 iprim_offset = 0 370 for iat in range(mol.natoms): 371 alpha_min_dict_iatom = {} 372 nshell_atom = 0 373 sum_prim = 0 374 for ishell in range(ishell_offset, self.nshells): 375 nshell_atom = nshell_atom + 1 376 sum_prim = sum_prim + self.nprims_shells[ishell][1] 377 if sum_prim==self.nprims_atoms[iat][1]: 378 break 379 subset_nprims_shell_l_list = self.nprims_shell_l_list[ishell_offset : ishell_offset + nshell_atom] 380 ishell_offset = ishell_offset + nshell_atom 381 # Create self.nprims_angmom_list 382 subset_nprims_angmom_list = self.calc_nprim_for_each_angular_momentum_l(subset_nprims_shell_l_list) 383 # alpha_min 384 for i in range(len(subset_nprims_angmom_list)): 385 alpha_min_dict_iatom[subset_nprims_angmom_list[i][0] - 1] = min(self.prim_expnts[iprim_offset : iprim_offset + subset_nprims_angmom_list[i][1]]) 386 iprim_offset = iprim_offset + subset_nprims_angmom_list[i][1] 387 388 self.alpha_min.append(alpha_min_dict_iatom) 389 390 # Calculate the BF offset index for a given shell 391 self.shell_bfs_offset = np.cumsum(self.bfs_nbfshell) - self.bfs_nbfshell 392 393 394 395 # NORMALIZATION STUFF 396 #Normalization factor for a primitive gaussian 397 def normalizationFactor(self, alpha, l, m, n): 398 """ 399 Calculate the normalization factor for a primitive Gaussian function. 400 401 Args: 402 alpha (float): Exponent of the Gaussian primitive 403 l (int): Angular momentum quantum number in x-direction 404 m (int): Angular momentum quantum number in y-direction 405 n (int): Angular momentum quantum number in z-direction 406 407 Returns: 408 float: Normalization factor for the primitive Gaussian 409 """ 410 return np.sqrt((2*alpha/np.pi)**(3/2)*(4*alpha)**(l+m+n)/(factorial2(2*l-1)*factorial2(2*m-1)*factorial2(2*n-1))) 411 412 #Normalization factor for a contraction of gaussians 413 def normalizationFactorContraction(self, alphas, coeffs, norms, l, m, n, lm): 414 """ 415 Calculate the normalization factor for a contraction of Gaussian primitives. 416 417 Args: 418 alphas (list): List of exponents for the primitive Gaussians 419 coeffs (list): List of contraction coefficients 420 norms (list): List of normalization factors for individual primitives 421 l (int): Angular momentum quantum number in x-direction 422 m (int): Angular momentum quantum number in y-direction 423 n (int): Angular momentum quantum number in z-direction 424 lm (int): Total angular momentum (l+m+n) 425 426 Returns: 427 float: Normalization factor for the contracted Gaussian 428 """ 429 430 temp = np.pi**(3/2)*factorial2(2*l-1)*factorial2(2*m-1)*factorial2(2*n-1)/(2**lm) 431 sum = 0.0 432 for i in range(len(alphas)): 433 for j in range(len(alphas)): 434 sum = sum + (coeffs[i]*coeffs[j]*norms[i]*norms[j]/((alphas[i]+alphas[j])**(lm+3/2))) 435 factor = np.power(temp*sum,-1/2) 436 return factor 437 438 439 def readBasisSetFromFile(key, filename): 440 """ 441 Read the basis set block corresponding to a particular atom from a TURBOMOLE format file. 442 443 Args: 444 key (str): Atomic species symbol to search for 445 filename (str): Path to the basis set file 446 447 Returns: 448 str or bool: Basis set string for the atom, or False if not found 449 """ 450 #This functions returns the basis set block corresponding to a particular atom (key) 451 #The file to be read from is assumed to follow TURBOMOLE's format 452 basisString = '' 453 lookfor = '' 454 file = open(filename, 'r') 455 fileContentsInString=file.read() 456 file.close() 457 lines = fileContentsInString.splitlines() 458 pattern = re.compile('\*\n(.*)\n\*') 459 result = pattern.findall(fileContentsInString) 460 if len(result)==0: 461 print('The basis set corresponding to atom',key, ' was not found in ',filename) 462 return False 463 for res in result: 464 if res.split()[0].lower()==key.lower(): 465 lookfor = res 466 basisString = '\n*\n'+lookfor+'\n*' 467 468 currentLineNo = 0 469 startReading = False 470 for line in lines: 471 line = line.strip() 472 if line==lookfor: 473 startReading = True 474 currentLineNo = 1 475 if startReading and currentLineNo>=3: 476 if line.strip()=='*': 477 break 478 else: 479 basisString = basisString + '\n'+line 480 currentLineNo = currentLineNo + 1 481 basisString = basisString + '\n*' 482 return basisString 483 484 485 486 487 488 #Loads the complete basis set as a string for a given mol/atom and a standard basis set available in the 489 #CrysX library or same directory as the python script. 490 def load( atom=None, mol=None, basis_name=None): 491 """ 492 Load the complete basis set as a string for a given atom/molecule from the CrysX library. 493 494 Args: 495 atom (str, optional): Atomic species symbol 496 mol (Mol, optional): Mol object containing molecular information 497 basis_name (str, optional): Name of the basis set to load 498 499 Returns: 500 str: Complete basis set string in TURBOMOLE format 501 """ 502 # The CrysX library contains basis sets in the TURBOMOLE format. 503 # The basis sets are downloaded from https://www.basissetexchange.org 504 # BSSE Github: https://github.com/MolSSI-BSE/basis_set_exchange 505 # License of BSSE: BSD 3 506 507 # We (CrysX) have saved the basis sets in the 'BasisSets' directory. 508 # The basis set name should be in lower-case and followed by '.version' and '.tm' extension. 509 # Ex: def2-tzvp.1.tm 510 # We can't expect the user to enter the basis set name with such an extension 511 # or in lower or upper case. So we should process this input name. 512 basis_name = basis_name.strip() 513 #The basis set for atom / mol 514 basisSet = basis_name+'\n' 515 basis_name = basis_name.replace(" ", "") 516 basis_name = basis_name.lower() 517 basis_name = basis_name+'.1.tm' 518 # Here we need to create the path name to the BasisSets directory supplied with crysx 519 # Since, the path should be with respect to hte module's location and not the working directory we can follow 520 # this link for various methods: https://stackoverflow.com/questions/4060221/how-to-reliably-open-a-file-in-the-same-directory-as-a-python-script 521 # Method 1: __file__ 522 # Method 2: sys.argv[0] 523 # Method 3: sys.path[0] (probably not relevant to our purpose) 524 # Currently we use method 1 525 temp = __file__ 526 temp = temp.replace('Basis.py', '') 527 basis_name = temp + 'BasisSets/'+basis_name 528 #basis_name = 'BasisSets/'+basis_name 529 530 if basis_name is None: 531 print('Error: A basis set name is needed!') 532 return basisSet 533 #If mol is provided then load the same basis set to all atoms 534 if not mol is None: 535 for i in range(mol.natoms): 536 basisSet = basisSet + Basis.readBasisSetFromFile( mol.atomicSpecies[i], basis_name) 537 elif atom is not None: 538 basisSet = Basis.readBasisSetFromFile(atom, basis_name) 539 basisSet = basisSet.replace('D+', 'E+') 540 basisSet = basisSet.replace('D-', 'E-') 541 return basisSet 542 543 #Loads the complete basis set as a string for a given mol/atom from a file 544 #specified by the user. 545 def loadfromfile( atom=None, mol=None, basis_name=None): 546 """ 547 Load the complete basis set as a string for a given atom/molecule from a user-specified file. 548 549 Args: 550 atom (str, optional): Atomic species symbol 551 mol (Mol, optional): Mol object containing molecular information 552 basis_name (str, optional): Path to the basis set file 553 554 Returns: 555 str: Complete basis set string in TURBOMOLE format 556 """ 557 #The basis set for atom / mol 558 basisSet = '' 559 if basis_name is None: 560 print('Error: A basis set name is needed!') 561 return basisSet 562 #If mol is provided then load the same basis set to all atoms 563 if not mol is None: 564 for i in range(mol.natoms): 565 basisSet = basisSet + Basis.readBasisSetFromFile( mol.atomicSpecies[i], basis_name) 566 elif atom is not None: 567 basisSet = Basis.readBasisSetFromFile(atom, basis_name) 568 basisSet = basisSet.replace('D+', 'E+') 569 basisSet = basisSet.replace('D-', 'E-') 570 return basisSet 571 572 #Creates the complete basis set for a given molecule along with basis dictionary 573 def createCompleteBasisSet(self, mol, basis): 574 """ 575 Create the complete basis set string for a given molecule using the basis dictionary. 576 577 Args: 578 mol: Mol object containing molecular information 579 basis (dict): Dictionary mapping atom types to basis set strings 580 581 Returns: 582 str: Complete basis set string for all atoms in the molecule 583 """ 584 totBasisSet = '' 585 if basis is None: 586 print('Error: Please provide a basis dictionary.') 587 return totBasisSet 588 else: 589 if 'all' in basis: 590 return basis['all'] 591 else: 592 #TODO Change this loop to only read the basis set for unique atoms 593 #TODO otherwise the basis set contains repeat basis sets for repeated atoms 594 #TODO The current form may(it doesn't yet) cause problem when creating basis function information 595 #TODO For even further generality let users label atoms, so that different basis may be used for same atoms with different labels 596 for i in range(mol.natoms): 597 if mol.atomicSpecies[i].lower() in basis: 598 totBasisSet = totBasisSet+ basis[mol.atomicSpecies[i].lower()]+'\n' 599 totBasisSet = totBasisSet.replace('D+', 'E+') 600 totBasisSet = totBasisSet.replace('D-', 'E-') 601 return totBasisSet 602 totBasisSet = totBasisSet.replace('D+', 'E+') 603 totBasisSet = totBasisSet.replace('D-', 'E-') 604 return totBasisSet 605 606 607 def cart2sph(l, ordering='pyscf'): 608 """ 609 Get the transformation matrix from Cartesian to real spherical harmonics for a given angular momentum. 610 611 Args: 612 l (int): Angular momentum quantum number 613 ordering (str): Ordering convention ('pyscf' or other) 614 615 Returns: 616 numpy.ndarray: Transformation matrix from Cartesian to spherical basis 617 """ 618 # This routine is used to get the transformation matrix from cartesian to real spherical basis for a given l 619 # This assumes that the basis functions are normalized. 620 # TODO: Make it work for both normalized and unnormalized basis functions. 621 622 # REFERENCE: 623 # https://theochem.github.io/horton/2.0.1/tech_ref_gaussian_basis.html#collected-notes-on-gaussian-basis-sets 624 # https://onlinelibrary.wiley.com/doi/epdf/10.1002/qua.560540202 625 # PySCF ordering: https://github.com/pyscf/pyscf/issues/1023 626 # https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics 627 l0 = np.array([[1.0]]) 628 l1 = np.array([[0, 0, 1.0], 629 [1.0, 0, 0], 630 [0, 1.0, 0]]) 631 l1_pyscf = np.array([[1.0, 0, 0], 632 [0, 1.0, 0], 633 [0, 0, 1.0]]) 634 l2 = np.array([[-0.5, 0, 0, -0.5, 0, 1.0], 635 [0, 0, 1.0, 0, 0, 0], 636 [0, 0, 0, 0, 1.0, 0], 637 [0.86602540378443864676, 0, 0, -0.86602540378443864676, 0, 0], 638 [0, 1.0, 0, 0, 0, 0], 639 ]) 640 # I don't remember what this was for (perhaps TURBOMOLE??) 641 l2_alternative =np.array([[-0.5, 0, 0, -0.5, 0, 1], 642 [0, 0, 0.7071067811865475, 0, 0, 0], 643 [0, 0, 0, 0, 0.7071067811865475, 0], 644 [0.6123724356957945, 0, 0, -0.6123724356957945, 0, 0], 645 [0, 0.7071067811865475, 0, 0, 0, 0], 646 ]) 647 l2_pyscf = np.array([[0, 1.0, 0, 0, 0, 0], 648 [0, 0, 0, 0, 1.0, 0], 649 [-0.5, 0, 0, -0.5, 0, 1.0], 650 [0, 0, 1.0, 0, 0, 0], 651 [0.86602540378443864676, 0, 0, -0.86602540378443864676, 0, 0], 652 ]) 653 l3 = np.array([[0, 0, -0.67082039324993690892, 0, 0, 0, 0, -0.67082039324993690892, 0, 1.0], 654 [-0.61237243569579452455, 0, 0, -0.27386127875258305673, 0, 1.0954451150103322269, 0, 0, 0, 0], 655 [0, -0.27386127875258305673, 0, 0, 0, 0, -0.61237243569579452455, 0, 1.0954451150103322269, 0], 656 [0, 0, 0.86602540378443864676, 0, 0, 0, 0, -0.86602540378443864676, 0, 0], 657 [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0], 658 [0.790569415042094833, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0], 659 [0, 1.0606601717798212866, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 660 ]) 661 l3_pyscf = np.array([[0, 1.0606601717798212866, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 662 [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0], 663 [0, -0.27386127875258305673, 0, 0, 0, 0, -0.61237243569579452455, 0, 1.0954451150103322269, 0], 664 [0, 0, -0.67082039324993690892, 0, 0, 0, 0, -0.67082039324993690892, 0, 1.0], 665 [-0.61237243569579452455, 0, 0, -0.27386127875258305673, 0, 1.0954451150103322269, 0, 0, 0, 0], 666 [0, 0, 0.86602540378443864676, 0, 0, 0, 0, -0.86602540378443864676, 0, 0], 667 [0.790569415042094833, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0], 668 ]) 669 l4 = np.array([[0.375, 0, 0, 0.21957751641341996535, 0, -0.87831006565367986142, 0, 0, 0, 0, 0.375, 0, -0.87831006565367986142, 0, 1.0], 670 [0, 0, -0.89642145700079522998, 0, 0, 0, 0, -0.40089186286863657703, 0, 1.19522860933439364, 0, 0, 0, 0, 0], 671 [0, 0, 0, 0, -0.40089186286863657703, 0, 0, 0, 0, 0, 0, -0.89642145700079522998, 0, 1.19522860933439364, 0], 672 [-0.5590169943749474241, 0, 0, 0, 0, 0.9819805060619657157, 0, 0, 0, 0, 0.5590169943749474241, 0, -0.9819805060619657157, 0, 0], 673 [0, -0.42257712736425828875, 0, 0, 0, 0, -0.42257712736425828875, 0, 1.1338934190276816816, 0, 0, 0, 0, 0, 0], 674 [0, 0, 0.790569415042094833, 0, 0, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0, 0], 675 [0, 0, 0, 0, 1.0606601717798212866, 0, 0, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 676 [0.73950997288745200532, 0, 0, -1.2990381056766579701, 0, 0, 0, 0, 0, 0, 0.73950997288745200532, 0, 0, 0, 0], 677 [0, 1.1180339887498948482, 0, 0, 0, 0, -1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0, 0], 678 ]) 679 l4_pyscf = np.array([[0, 1.1180339887498948482, 0, 0, 0, 0, -1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0, 0], 680 [0, 0, 0, 0, 1.0606601717798212866, 0, 0, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 681 [0, -0.42257712736425828875, 0, 0, 0, 0, -0.42257712736425828875, 0, 1.1338934190276816816, 0, 0, 0, 0, 0, 0], 682 [0, 0, 0, 0, -0.40089186286863657703, 0, 0, 0, 0, 0, 0, -0.89642145700079522998, 0, 1.19522860933439364, 0], 683 [0.375, 0, 0, 0.21957751641341996535, 0, -0.87831006565367986142, 0, 0, 0, 0, 0.375, 0, -0.87831006565367986142, 0, 1.0], 684 [0, 0, -0.89642145700079522998, 0, 0, 0, 0, -0.40089186286863657703, 0, 1.19522860933439364, 0, 0, 0, 0, 0], 685 [-0.5590169943749474241, 0, 0, 0, 0, 0.9819805060619657157, 0, 0, 0, 0, 0.5590169943749474241, 0, -0.9819805060619657157, 0, 0], 686 [0, 0, 0.790569415042094833, 0, 0, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0, 0], 687 [0.73950997288745200532, 0, 0, -1.2990381056766579701, 0, 0, 0, 0, 0, 0, 0.73950997288745200532, 0, 0, 0, 0], 688 ]) 689 l5 = np.array([[0, 0, 0.625, 0, 0, 0, 0, 0.36596252735569994226, 0, -1.0910894511799619063, 0, 0, 0, 0, 0, 0, 0.625, 0, -1.0910894511799619063, 0, 1.0], 690 [0.48412291827592711065, 0, 0, 0.21128856368212914438, 0, -1.2677313820927748663, 0, 0, 0, 0, 0.16137430609197570355, 0, -0.56694670951384084082, 0, 1.2909944487358056284, 0, 0, 0, 0, 0, 0], 691 [0, 0.16137430609197570355, 0, 0, 0, 0, 0.21128856368212914438, 0, -0.56694670951384084082, 0, 0, 0, 0, 0, 0, 0.48412291827592711065, 0, -1.2677313820927748663, 0, 1.2909944487358056284, 0], 692 [0, 0, -0.85391256382996653193, 0, 0, 0, 0, 0, 0, 1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0.85391256382996653193, 0, -1.1180339887498948482, 0, 0], 693 [0, 0, 0, 0, -0.6454972243679028142, 0, 0, 0, 0, 0, 0, -0.6454972243679028142, 0, 1.2909944487358056284, 0, 0, 0, 0, 0, 0, 0], 694 [-0.52291251658379721749, 0, 0, 0.22821773229381921394, 0, 0.91287092917527685576, 0, 0, 0, 0, 0.52291251658379721749, 0, -1.2247448713915890491, 0, 0, 0, 0, 0, 0, 0, 0], 695 [0, -0.52291251658379721749, 0, 0, 0, 0, -0.22821773229381921394, 0, 1.2247448713915890491, 0, 0, 0, 0, 0, 0, 0.52291251658379721749, 0, -0.91287092917527685576, 0, 0, 0], 696 [0, 0, 0.73950997288745200532, 0, 0, 0, 0, -1.2990381056766579701, 0, 0, 0, 0, 0, 0, 0, 0, 0.73950997288745200532, 0, 0, 0, 0], 697 [0, 0, 0, 0, 1.1180339887498948482, 0, 0, 0, 0, 0, 0, -1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0, 0, 0], 698 [0.7015607600201140098, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 1.169267933366856683, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 699 [0, 1.169267933366856683, 0, 0, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 0, 0, 0.7015607600201140098, 0, 0, 0, 0, 0], 700 ]) 701 l6 = np.array([[-0.3125, 0, 0, -0.16319780245846672329, 0, 0.97918681475080033975, 0, 0, 0, 0, -0.16319780245846672329, 0, 0.57335309036732873772, 0, -1.3055824196677337863, 0, 0, 0, 0, 0, 0, -0.3125, 0, 0.97918681475080033975, 0, -1.3055824196677337863, 0, 1.0], 702 [0, 0, 0.86356159963469679725, 0, 0, 0, 0, 0.37688918072220452831, 0, -1.6854996561581052156, 0, 0, 0, 0, 0, 0, 0.28785386654489893242, 0, -0.75377836144440905662, 0, 1.3816985594155148756, 0, 0, 0, 0, 0, 0, 0], 703 [0, 0, 0, 0, 0.28785386654489893242, 0, 0, 0, 0, 0, 0, 0.37688918072220452831, 0, -0.75377836144440905662, 0, 0, 0, 0, 0, 0, 0, 0, 0.86356159963469679725, 0, -1.6854996561581052156, 0, 1.3816985594155148756, 0], 704 [0.45285552331841995543, 0, 0, 0.078832027985861408788, 0, -1.2613124477737825406, 0, 0, 0, 0, -0.078832027985861408788, 0, 0, 0, 1.2613124477737825406, 0, 0, 0, 0, 0, 0, -0.45285552331841995543, 0, 1.2613124477737825406, 0, -1.2613124477737825406, 0, 0], 705 [0, 0.27308215547040717681, 0, 0, 0, 0, 0.26650089544451304287, 0, -0.95346258924559231545, 0, 0, 0, 0, 0, 0, 0.27308215547040717681, 0, -0.95346258924559231545, 0, 1.4564381625088382763, 0, 0, 0, 0, 0, 0, 0, 0], 706 [0, 0, -0.81924646641122153043, 0, 0, 0, 0, 0.35754847096709711829, 0, 1.0660035817780521715, 0, 0, 0, 0, 0, 0, 0.81924646641122153043, 0, -1.4301938838683884732, 0, 0, 0, 0, 0, 0, 0, 0, 0], 707 [0, 0, 0, 0, -0.81924646641122153043, 0, 0, 0, 0, 0, 0, -0.35754847096709711829, 0, 1.4301938838683884732, 0, 0, 0, 0, 0, 0, 0, 0, 0.81924646641122153043, 0, -1.0660035817780521715, 0, 0, 0], 708 [-0.49607837082461073572, 0, 0, 0.43178079981734839863, 0, 0.86356159963469679725, 0, 0, 0, 0, 0.43178079981734839863, 0, -1.5169496905422946941, 0, 0, 0, 0, 0, 0, 0, 0, -0.49607837082461073572, 0, 0.86356159963469679725, 0, 0, 0, 0], 709 [0, -0.59829302641309923139, 0, 0, 0, 0, 0, 0, 1.3055824196677337863, 0, 0, 0, 0, 0, 0, 0.59829302641309923139, 0, -1.3055824196677337863, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 710 [0, 0, 0.7015607600201140098, 0, 0, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 0, 0, 1.169267933366856683, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 711 [0, 0, 0, 0, 1.169267933366856683, 0, 0, 0, 0, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7015607600201140098, 0, 0, 0, 0, 0], 712 [0.67169328938139615748, 0, 0, -1.7539019000502850245, 0, 0, 0, 0, 0, 0, 1.7539019000502850245, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.67169328938139615748, 0, 0, 0, 0, 0, 0], 713 [0, 1.2151388809514737933, 0, 0, 0, 0, -1.9764235376052370825, 0, 0, 0, 0, 0, 0, 0, 0, 1.2151388809514737933, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 714 ]) 715 if ordering=='pyscf': 716 c2s_l = [l0,l1_pyscf,l2_pyscf,l3_pyscf,l4_pyscf,l5,l6] # PYSCF ordering of SAOs 717 else: 718 c2s_l = [l0,l1,l2,l3,l4,l5,l6] # HORTON ordering of SAOs 719 return c2s_l[l] 720 721 def cart2sph_basis(self): 722 """ 723 Return the complete Cartesian to spherical harmonic basis transformation matrix for all shells. 724 725 Returns: 726 matrix: Block diagonal transformation matrix for the entire basis set 727 """ 728 # Returns the complete Cartesian to Spherical (real) basis transformation matrix 729 cart2sph = [] 730 for i in range(self.nshells): 731 l = self.shells[i]-1 732 cart2sph.append(Basis.cart2sph(l)) 733 734 return scipy.linalg.block_diag(*cart2sph) 735 736 def sph2cart_basis(self): 737 """ 738 Return the complete spherical to Cartesian basis transformation matrix for all shells. 739 740 Returns: 741 matrix: Block diagonal transformation matrix (pseudoinverse of cart2sph) 742 """ 743 # Returns the complete Cartesian to Spherical (real) basis transformation matrix 744 sph2cart = [] 745 for i in range(self.nshells): 746 l = self.shells[i]-1 747 sph2cart_pseudo = np.linalg.pinv(Basis.cart2sph(l)) 748 sph2cart.append(sph2cart_pseudo) 749 750 return scipy.linalg.block_diag(*sph2cart) 751 752 753 # #Probably won't be used 754 # def readBasisFunctionsfromFile(self, mol, basis_name): 755 # #This function will read and parse the basis sets for the basis functions and other properties. 756 # #Returns true if successful or false if failed. 757 # if mol is None: 758 # print('Error: Please provide a Mol object.') 759 # return False 760 # if basis_name is None: 761 # print('Error: Please provide a basis set name dictionary for Mol atoms.') 762 # return False 763 # #If all atoms are assigned the same basis set 764 # if 'all' in basis_name: 765 # for i in range(mol.natoms): 766 # self.basis_name_local[mol.atomicSpecies[i]] = basis_name['all'] 767 # else: 768 # for i in range(mol.natoms): 769 # #If a particular basis set is assigned to an atom 770 # if mol.atomicSpecies[i] in basis_name: 771 # self.basis_name_local[mol.atomicSpecies[i]] = basis_name[mol.atomicSpecies[i]] 772 # else: 773 # #Or assign the default basis set 774 # self.basis_name_local[mol.atomicSpecies[i]] = basis_name_local['default'] 775 776 def readBasisFunctionInfo(self, mol, basisSet): 777 """ 778 Read and parse the basis set information to populate basis function properties. 779 780 Args: 781 mol: Mol object containing molecular information 782 basisSet (str): Complete basis set string in TURBOMOLE format 783 784 Returns: 785 None - populates internal data structures with basis function information 786 """ 787 #TODO Add error checks, 788 #TODO Like check if the basis set contains the needed atoms or not, and other checks that you can think of 789 #TODO Also, add error messages explaining what is happening 790 791 #Read the basis set and set the basis set information 792 indx_shell = 0 793 for i in range(mol.natoms): 794 lookfor = '' 795 lines = basisSet.splitlines() 796 pattern = re.compile('\*\n(.*)\n\*') 797 result = pattern.findall(basisSet) 798 if len(result)==0: 799 print('The basis set corresponding to atom ',mol.atomicSpecies[i], ' was not found in the provided basis set.') 800 for res in result: 801 if res.split()[0].lower()==mol.atomicSpecies[i].lower(): 802 lookfor = res 803 startReading = False 804 #print(lines) 805 currentLineNo = 0 806 while currentLineNo<len(lines): 807 line = lines[currentLineNo].strip() 808 #print(line) 809 if line==lookfor: 810 startReading = True 811 currentLineNo = currentLineNo + 2 812 if startReading: 813 line = lines[currentLineNo].strip() 814 splitLine = line.split() 815 if '*' in line: 816 startReading = False 817 currentLineNo = currentLineNo + 1 818 break 819 #print(line) 820 #Check if the line looks like ' 3 s' 821 if isinstance(int(splitLine[0]), int) and isinstance(str(splitLine[1]), str): 822 #print('here') 823 #Read the shell information 824 #Read the number of primitives in the given shell Ex: '3 s' 825 self.nprims.append(int(splitLine[0])) #Read '3' 826 #Read the shell type Ex: '3 s' 827 self.shellsLabel.append(str(splitLine[1])) #Read 's' 828 #Get the angular momentum 'l' corresponding to the shell type. 829 #Ex: 's'->1, 'p'->2, 'd'->3, etc. 830 l = Data.shell_dict[splitLine[1]] 831 self.shells.append(l) 832 self.shell_degen.append(Data.shell_degen[l-1]) 833 self.shell_coords.append(mol.coordsBohrs[i]) 834 #create the information for basis functions 835 #for ibf in range(int(l*(l+1)/2.0)): 836 # self.bfcoords.append(mol.coords[i]) 837 #Run a loop over the number of primitives in the current shell 838 #and read the exponents and contraction coefficients 839 for k in range(int(splitLine[0])): 840 currentLineNo = currentLineNo + 1 841 line = lines[currentLineNo].strip() 842 splitLine = line.split() 843 if '*' in line: 844 startReading = False 845 currentLineNo = currentLineNo + 1 846 break 847 splitLine[0] = splitLine[0].replace('D+', 'E+') 848 splitLine[1] = splitLine[1].replace('D-', 'E-') 849 self.prim_expnts.append(float(splitLine[0])) 850 self.prim_coeffs.append(float(splitLine[1])) 851 self.prim_coords.append(mol.coordsBohrs[i]) 852 self.prim_atoms.append(i) 853 self.prim_shells.append(indx_shell) 854 #self.prim_norms.append(float(normalizationFactor())) 855 #print(self.prim_coeffs, self.prim_expnts) 856 #print(currentLineNo) 857 858 indx_shell +=1 859 for ibf in range(int(l*(l+1)/2.0)): 860 self.bfs_atoms.append(i) 861 862 currentLineNo = currentLineNo + 1 863 864 865 866 #Now that the reading of shells and their primitives is done 867 868 #Run a loop over shells, and create information for 869 # basis functions 870 #for i in range(shells 871 872 def calc_nprim_for_each_angular_momentum_l(self, tuple_list): 873 """ 874 Calculate the number of primitives for each angular momentum quantum number. 875 876 Args: 877 tuple_list (list): List of tuples containing (angular_momentum, num_primitives) 878 879 Returns: 880 list: List of tuples with summed primitives for each unique angular momentum 881 """ 882 # Initialize the list of tuples 883 # tuple_list = [(0,12),(0,5),(1,2)] 884 # tuple_list = self.nprims_shells 885 886 # Create a dictionary to store the sums of each common angular momentum 887 sums_dict = {} 888 889 # Loop through each tuple in the list 890 for t in tuple_list: 891 # Get the first element of the tuple 892 key = t[0] 893 894 # Check if the element is already in the dictionary 895 if key in sums_dict: 896 # If it is, add the second element of the tuple to the sum 897 sums_dict[key] += t[1] 898 else: 899 # If it isn't, add the element to the dictionary with the value of the second element of the tuple 900 sums_dict[key] = t[1] 901 902 # Initialize the list of tuples to return 903 result = [] 904 905 # Loop through each key-value pair in the dictionary 906 for key, value in sums_dict.items(): 907 # Append a tuple with the key and value to the result list 908 result.append((key, value)) 909 910 # self.nprims_angmom_list = result 911 return result 912 913
class
Basis:
33class Basis: 34 """ 35 Class to store and manage basis function properties for a given molecular system. 36 37 This class processes atomic basis set information for each atom in the molecule, 38 including primitive Gaussian functions, shells, and contracted basis functions (AOs). 39 It also supports TURBOMOLE-style shell ordering. 40 """ 41 def __init__(self, mol, basis, tmoleFormat=False): 42 """ 43 Initialize a Basis object for storing basis function properties. 44 45 Args: 46 mol: Mol object containing molecular information 47 basis: Dictionary specifying the basis set to be used for each atom, or None to use mol.basis 48 tmoleFormat (bool): Whether to use TURBOMOLE format ordering for basis functions 49 50 Returns: 51 None 52 53 Raises: 54 None - prints error message if mol is None 55 """ 56 if mol is None: 57 print('Error: Please provide a Mol object.') 58 return None 59 #Basis set name 60 #If no basis set name is provided specifically, 61 #then use the basis set name of the mol. 62 #'basis' is a dictionary specifying the basis set to be used for a particular atom 63 if basis is None: 64 self.basis = mol.basis 65 else: 66 self.basis = basis 67 68 69 self.basisSet = '' 70 """Basis set name or label. 71 Type: str 72 """ 73 74 self.nao = 0 75 """Number of basis functions (atomic orbitals). 76 Type: int 77 """ 78 79 self.nshells = 0 80 """Number of shells. 81 Type: int 82 """ 83 84 self.totalnprims = 0 85 """Total number of primitive Gaussian functions. 86 Type: int 87 """ 88 89 self.nprims = [] 90 """Number of primitives in each shell. 91 Type: List[int] 92 """ 93 94 self.prim_coeffs = [] 95 """Primitive contraction coefficients. 96 Type: List[float] 97 """ 98 #Exponents of primitives 99 self.prim_expnts = [] 100 """Exponents of primitives. 101 Type: List[float] 102 """ 103 #Coords of primitives 104 self.prim_coords = [] 105 """Coordinates of each primitive function. 106 Type: List[List[float]] 107 """ 108 #Normalization factors of primitives 109 self.prim_norms = [] 110 """Normalization constants of primitives. 111 Type: List[float] 112 """ 113 # This list contains the indices of atoms (first = #0) that correspond to a particular primitive function 114 self.prim_atoms = [] #(This list should be of the size of the number of self.totalnprims) 115 """Atom index corresponding to each primitive. 116 Type: List[int] 117 """ 118 # A list of tuples (2D list) that gives the number of primitives belonging to each atomic index 119 self.nprims_atoms = [] # Size is the same as that of no. of atoms. Looks like this: [(0, 12), (1, 5), (2, 5)] 120 """Number of primitives per atom as (atom_index, nprims). 121 Type: List[Tuple[int, int]] 122 """ 123 # This list contains the indices of shells (first = #0) that correspond to a particular primitive function 124 self.prim_shells = [] #(This list should be of the size of the number of self.totalnprims) 125 """Shell index corresponding to each primitive. 126 Type: List[int] 127 """ 128 # A list of tuples (2D list) that gives the number of primitives belonging to each shell index 129 self.nprims_shells = [] # Size is the same as that of no. of shells. Looks like this: [(0, 12), (1, 5), (2, 5)] 130 """Number of primitives per shell as (shell_index, nprims). 131 Type: List[Tuple[int, int]] 132 """ 133 # A special 2d list that contains the angular momentum of shell in the first index and the no. of primitives corresponding to the shell as second index 134 self.nprims_shell_l_list = [] # This will be of the same size as the no. of shells 135 """Angular momentum and primitive count per shell as (l, nprims). 136 Type: List[Tuple[int, int]] 137 """ 138 # A list of size natoms, that contains the values of the largest primitive exponents for all the atoms 139 self.alpha_max = [] 140 """Maximum exponent among primitives for each atom. 141 Type: List[float] 142 """ 143 # A list of size natoms, that contains the list of the smallest primitive exponents for each shell for all the atoms 144 self.alpha_min = [] 145 """Minimum exponent of primitives per shell for each atom. 146 Type: List[List[float]] 147 """ 148 #Number of contractions 149 self.ncontrs = 0 150 """Total number of contracted basis functions. 151 Type: int 152 """ 153 #Normalization factors for contractions 154 self.contrs_norm = [] 155 """Normalization factors for contracted basis functions. 156 Type: List[float] 157 """ 158 #Degeneracy of shells 159 self.shell_degen = [] 160 """Degeneracy of each shell. 161 Type: List[int] 162 """ 163 #Shells 164 self.shellsLabel = [] 165 """Label of each shell (e.g., 'S', 'P', etc.). 166 Type: List[str] 167 """ 168 self.shells = [] 169 """Shell index (0: s, 1: p, 2:d). 170 Type: List[int] 171 """ 172 self.shell_coords = [] 173 """Coordinates of each shell center. 174 Type: List[List[float]] 175 """ 176 #--------------------------------- 177 #Information for basis function 178 #--------------------------------- 179 #Coordinates of basis functions (This list should be of the size of the number of AOs/BFs) 180 self.bfs_coords = [] 181 """Coordinates of basis functions (same as their parent shell). 182 Type: List[List[float]] 183 """ 184 # Number of primitives corresponding to each basis function (This list should be of the size of the number of AOs/BFs) 185 self.bfs_nprim = [] 186 """Number of primitives corresponding to each basis function. 187 Type: List[int] 188 """ 189 # Radius cutoff for each basis function (This list should be of the size of the number of AOs/BFs) 190 # This is used to accelerate the evaluation of AOs on grid points, as those gridpoints farther than the radius cutoff 191 # of the basis function can be discarded for calcualtion. 192 self.bfs_radius_cutoff = [] 193 """Radius cutoff for each basis function to speed up AO evaluation. 194 Type: List[float] 195 """ 196 #A list of the size of number of atomic orbitals (BFs), 197 #that contains lists of primitive contraction coefficients/exponents 198 #corresponding to every basis function 199 self.bfs_expnts = [] 200 """Primitive exponents used in each basis function. 201 Type: List[List[float]] 202 """ 203 204 self.bfs_coeffs = [] 205 """Primitive coefficients for each basis function. 206 Type: List[List[float]] 207 """ 208 # A list of the size of number of AOs (BFs) 209 # that contains the lists of normalization factors for 210 # every primitive corresponding to the basis functions 211 self.bfs_prim_norms = [] 212 """Primitive normalization constants per basis function. 213 Type: List[List[float]] 214 """ 215 # Normalization factor for the contraction of primitives making up a BF (This list should be of the size of the number of AOs/BFs) 216 self.bfs_contr_prim_norms = [] 217 """Normalization for the contraction of primitives forming a BF. 218 Type: List[float] 219 """ 220 # l,m,n indices of BFs (This list should be of the size of the number of AOs/BFs) 221 self.bfs_lmn = [] 222 """Cartesian angular momentum indices (l, m, n) per basis function. 223 Type: List[Tuple[int, int, int]] 224 """ 225 # Angular momentum of a basis function (This list should be of the size of the number of AOs/BFs) 226 self.bfs_lm = [] 227 """Angular momentum `l` for each basis function. 228 Type: List[int] 229 """ 230 # Label of a basis function (This list should be of the size of the number of AOs/BFs) 231 self.bfs_label = [] 232 """Label for each basis function. 233 Type: List[str] 234 """ 235 # Number of BFs in a shell 236 self.bfs_nbfshell = [] 237 """Number of BFs in the shell of each basis function. 238 Type: List[int] 239 """ 240 # Shell index of each bf 241 self.bfs_shell_index =[] 242 """Shell index for each basis function. 243 Type: List[int] 244 """ 245 # BF offset index of each shell (This should be of the size of the number of shells) 246 self.shell_bfs_offset =[] 247 """Offset index of the first basis function in each shell. 248 Type: List[int] 249 """ 250 # Total number of basis functions (BFs) also called Atomic orbitals (AOs) 251 self.bfs_nao = [] 252 """Index of each basis function (redundant with `nao`). 253 Type: List[int] 254 """ 255 # This list contains the indices of atoms (first = #0) that correspond to a particular bf 256 self.bfs_atoms = [] #(This list should be of the size of the number of AOs/BFs) 257 """Atom index corresponding to each basis function. 258 Type: List[int] 259 """ 260 #--------------------------------- 261 #If all data was parsed successfully. 262 self.success = False 263 """Indicates if basis parsing was successful. 264 Type: bool 265 """ 266 #Convert the keys of basis dictionary to lower case to avoid ambiguity 267 self.basisSet = self.createCompleteBasisSet(mol, basis) 268 self.readBasisFunctionInfo(mol, self.basisSet) 269 self.totalnprims = sum(self.nprims) 270 self.nshells = len(self.nprims) 271 # Now lets create information about the basis functions 272 # A basis function is a combination of gaussian primitives, 273 # so we need 274 # number of primitives for a given bf, 275 # the exponents and contraction coefficients of the primitives, 276 # the l,m,n triplet information, 277 # the angular momentum l+m+n, 278 # coords, 279 # normalization factors of primitives, 280 # normalization factor for the contraction. 281 #--------------------------------------- 282 offset = 0 283 ishell = 0 284 # Loop through all the shells 285 for i in range(len(self.nprims)): 286 # Number of primitives in a shell 287 numprim = self.nprims[i] 288 # Angular momentum of a shell 289 lm = self.shells[i] - 1 290 # no of bfs in this shell 291 nbf_shell = Data.shell_degen[lm] 292 self.bfs_nbfshell.append(nbf_shell) 293 # Assign the lmn triplets and other information corresponding to basis functions 294 # by looping through the degeneracy of the shells. 295 # Example: If the degeneracy is 1 then we loop through [0,1) that is we look 296 # for the 0th element in the Data.shell_lmn dictionary as that corresponds to the 's' lmn triplet. 297 # Similarly if the denegeracy was 3 as for 'p' shell, then we loop from [1,3] and get the 298 # lmn triplets corresponding to 'p' shell and so on. 299 offset_shell_labels = Data.shell_lmn_offset[lm] 300 for j in range(offset_shell_labels, nbf_shell + offset_shell_labels): 301 #Assign the label to the basis function 302 if tmoleFormat: 303 label = list(Data.shell_lmn_tmole)[j] 304 else: 305 label = list(Data.shell_lmn)[j] 306 self.bfs_label.append(label) 307 # Assign l,m,n triplet for a basis function 308 if tmoleFormat: 309 lmn = Data.shell_lmn_tmole[label] 310 else: 311 lmn = Data.shell_lmn[label] 312 self.bfs_lmn.append(lmn) 313 # Assign the number of primitives in a basis function 314 self.bfs_nprim.append(numprim) 315 # Assign the the total angular momentum of a basis function 316 self.bfs_lm.append(sum(lmn)) 317 # Assign the exponents and contraction coefficients 318 expnts = self.prim_expnts[offset : offset+numprim ] 319 self.bfs_expnts.append(expnts) 320 coeffs = self.prim_coeffs[offset : offset+numprim ] 321 self.bfs_coeffs.append(coeffs) 322 # Assign Radius Cutoffs 323 radius_cutoff = 0.0 324 # Loop over primitives 325 eeta_log = np.log(1.0E-9) 326 for iprim_ in range(len(expnts)): 327 radius_cutoff = np.maximum(radius_cutoff, np.sqrt(1.0/expnts[iprim_]*(np.log(expnts[iprim_])/2.0 - eeta_log))) 328 self.bfs_radius_cutoff.append(radius_cutoff) 329 # Assign the normalization factors for all the primitives corresponding to a BF 330 norms = [] 331 for k in range(offset,offset+numprim): 332 norms.append(self.normalizationFactor(self.prim_expnts[k],lmn[0],lmn[1],lmn[2])) 333 self.bfs_prim_norms.append(norms) 334 self.bfs_contr_prim_norms.append(self.normalizationFactorContraction(expnts, coeffs, norms, lmn[0],lmn[1],lmn[2], lm)) 335 # Assign coordinates to the basis functions 336 self.bfs_coords.append(self.shell_coords[i]) 337 # Assign shell index to the basis functions 338 self.bfs_shell_index.append(ishell) 339 340 ishell += 1 341 offset = offset + numprim 342 self.bfs_nao = len(self.bfs_label) 343 344 ##### The following was mainly developed to calculated alpha_min and alpha_max for numgrid 345 346 # Calculate the number of primitives corresponding to each atom 347 cnt = Counter(self.prim_atoms) # Gives a Counter object similar to a dictionary 348 # Convert the Counter object to a list 349 self.nprims_atoms = list(cnt.items()) #a list of tuples, each containing a value and its count in the list 350 ### Start processing to calculate alpha_min and alpha_max for grid generation 351 ## Calculate alpha_max (https://github.com/dftlibs/numgrid) 352 self.alpha_max = [] 353 offset = 0 354 for iat in range(mol.natoms): 355 alpha_max = max(self.prim_expnts[offset : offset + self.nprims_atoms[iat][1]]) 356 self.alpha_max.append(alpha_max) 357 offset = offset + self.nprims_atoms[iat][1] 358 359 # Calculate the number of primitives corresponding to each shell 360 cnt = Counter(self.prim_shells) # Gives a Counter object similar to a dictionary 361 # Convert the Counter object to a list 362 self.nprims_shells = list(cnt.items()) #a list of tuples, each containing a value and its count in the list 363 # Make a list that is pretty much the same as nprims_shells, however, it contains the shell angular momentum in the first index 364 for ishell in range(self.nshells): 365 self.nprims_shell_l_list.append((self.shells[ishell], self.nprims_shells[ishell][1])) 366 # Create self.nprims_angmom_list 367 # self.calc_nprim_for_each_angular_momentum_l() 368 ## Calculate alpha_min 369 ishell_offset = 0 370 iprim_offset = 0 371 for iat in range(mol.natoms): 372 alpha_min_dict_iatom = {} 373 nshell_atom = 0 374 sum_prim = 0 375 for ishell in range(ishell_offset, self.nshells): 376 nshell_atom = nshell_atom + 1 377 sum_prim = sum_prim + self.nprims_shells[ishell][1] 378 if sum_prim==self.nprims_atoms[iat][1]: 379 break 380 subset_nprims_shell_l_list = self.nprims_shell_l_list[ishell_offset : ishell_offset + nshell_atom] 381 ishell_offset = ishell_offset + nshell_atom 382 # Create self.nprims_angmom_list 383 subset_nprims_angmom_list = self.calc_nprim_for_each_angular_momentum_l(subset_nprims_shell_l_list) 384 # alpha_min 385 for i in range(len(subset_nprims_angmom_list)): 386 alpha_min_dict_iatom[subset_nprims_angmom_list[i][0] - 1] = min(self.prim_expnts[iprim_offset : iprim_offset + subset_nprims_angmom_list[i][1]]) 387 iprim_offset = iprim_offset + subset_nprims_angmom_list[i][1] 388 389 self.alpha_min.append(alpha_min_dict_iatom) 390 391 # Calculate the BF offset index for a given shell 392 self.shell_bfs_offset = np.cumsum(self.bfs_nbfshell) - self.bfs_nbfshell 393 394 395 396 # NORMALIZATION STUFF 397 #Normalization factor for a primitive gaussian 398 def normalizationFactor(self, alpha, l, m, n): 399 """ 400 Calculate the normalization factor for a primitive Gaussian function. 401 402 Args: 403 alpha (float): Exponent of the Gaussian primitive 404 l (int): Angular momentum quantum number in x-direction 405 m (int): Angular momentum quantum number in y-direction 406 n (int): Angular momentum quantum number in z-direction 407 408 Returns: 409 float: Normalization factor for the primitive Gaussian 410 """ 411 return np.sqrt((2*alpha/np.pi)**(3/2)*(4*alpha)**(l+m+n)/(factorial2(2*l-1)*factorial2(2*m-1)*factorial2(2*n-1))) 412 413 #Normalization factor for a contraction of gaussians 414 def normalizationFactorContraction(self, alphas, coeffs, norms, l, m, n, lm): 415 """ 416 Calculate the normalization factor for a contraction of Gaussian primitives. 417 418 Args: 419 alphas (list): List of exponents for the primitive Gaussians 420 coeffs (list): List of contraction coefficients 421 norms (list): List of normalization factors for individual primitives 422 l (int): Angular momentum quantum number in x-direction 423 m (int): Angular momentum quantum number in y-direction 424 n (int): Angular momentum quantum number in z-direction 425 lm (int): Total angular momentum (l+m+n) 426 427 Returns: 428 float: Normalization factor for the contracted Gaussian 429 """ 430 431 temp = np.pi**(3/2)*factorial2(2*l-1)*factorial2(2*m-1)*factorial2(2*n-1)/(2**lm) 432 sum = 0.0 433 for i in range(len(alphas)): 434 for j in range(len(alphas)): 435 sum = sum + (coeffs[i]*coeffs[j]*norms[i]*norms[j]/((alphas[i]+alphas[j])**(lm+3/2))) 436 factor = np.power(temp*sum,-1/2) 437 return factor 438 439 440 def readBasisSetFromFile(key, filename): 441 """ 442 Read the basis set block corresponding to a particular atom from a TURBOMOLE format file. 443 444 Args: 445 key (str): Atomic species symbol to search for 446 filename (str): Path to the basis set file 447 448 Returns: 449 str or bool: Basis set string for the atom, or False if not found 450 """ 451 #This functions returns the basis set block corresponding to a particular atom (key) 452 #The file to be read from is assumed to follow TURBOMOLE's format 453 basisString = '' 454 lookfor = '' 455 file = open(filename, 'r') 456 fileContentsInString=file.read() 457 file.close() 458 lines = fileContentsInString.splitlines() 459 pattern = re.compile('\*\n(.*)\n\*') 460 result = pattern.findall(fileContentsInString) 461 if len(result)==0: 462 print('The basis set corresponding to atom',key, ' was not found in ',filename) 463 return False 464 for res in result: 465 if res.split()[0].lower()==key.lower(): 466 lookfor = res 467 basisString = '\n*\n'+lookfor+'\n*' 468 469 currentLineNo = 0 470 startReading = False 471 for line in lines: 472 line = line.strip() 473 if line==lookfor: 474 startReading = True 475 currentLineNo = 1 476 if startReading and currentLineNo>=3: 477 if line.strip()=='*': 478 break 479 else: 480 basisString = basisString + '\n'+line 481 currentLineNo = currentLineNo + 1 482 basisString = basisString + '\n*' 483 return basisString 484 485 486 487 488 489 #Loads the complete basis set as a string for a given mol/atom and a standard basis set available in the 490 #CrysX library or same directory as the python script. 491 def load( atom=None, mol=None, basis_name=None): 492 """ 493 Load the complete basis set as a string for a given atom/molecule from the CrysX library. 494 495 Args: 496 atom (str, optional): Atomic species symbol 497 mol (Mol, optional): Mol object containing molecular information 498 basis_name (str, optional): Name of the basis set to load 499 500 Returns: 501 str: Complete basis set string in TURBOMOLE format 502 """ 503 # The CrysX library contains basis sets in the TURBOMOLE format. 504 # The basis sets are downloaded from https://www.basissetexchange.org 505 # BSSE Github: https://github.com/MolSSI-BSE/basis_set_exchange 506 # License of BSSE: BSD 3 507 508 # We (CrysX) have saved the basis sets in the 'BasisSets' directory. 509 # The basis set name should be in lower-case and followed by '.version' and '.tm' extension. 510 # Ex: def2-tzvp.1.tm 511 # We can't expect the user to enter the basis set name with such an extension 512 # or in lower or upper case. So we should process this input name. 513 basis_name = basis_name.strip() 514 #The basis set for atom / mol 515 basisSet = basis_name+'\n' 516 basis_name = basis_name.replace(" ", "") 517 basis_name = basis_name.lower() 518 basis_name = basis_name+'.1.tm' 519 # Here we need to create the path name to the BasisSets directory supplied with crysx 520 # Since, the path should be with respect to hte module's location and not the working directory we can follow 521 # this link for various methods: https://stackoverflow.com/questions/4060221/how-to-reliably-open-a-file-in-the-same-directory-as-a-python-script 522 # Method 1: __file__ 523 # Method 2: sys.argv[0] 524 # Method 3: sys.path[0] (probably not relevant to our purpose) 525 # Currently we use method 1 526 temp = __file__ 527 temp = temp.replace('Basis.py', '') 528 basis_name = temp + 'BasisSets/'+basis_name 529 #basis_name = 'BasisSets/'+basis_name 530 531 if basis_name is None: 532 print('Error: A basis set name is needed!') 533 return basisSet 534 #If mol is provided then load the same basis set to all atoms 535 if not mol is None: 536 for i in range(mol.natoms): 537 basisSet = basisSet + Basis.readBasisSetFromFile( mol.atomicSpecies[i], basis_name) 538 elif atom is not None: 539 basisSet = Basis.readBasisSetFromFile(atom, basis_name) 540 basisSet = basisSet.replace('D+', 'E+') 541 basisSet = basisSet.replace('D-', 'E-') 542 return basisSet 543 544 #Loads the complete basis set as a string for a given mol/atom from a file 545 #specified by the user. 546 def loadfromfile( atom=None, mol=None, basis_name=None): 547 """ 548 Load the complete basis set as a string for a given atom/molecule from a user-specified file. 549 550 Args: 551 atom (str, optional): Atomic species symbol 552 mol (Mol, optional): Mol object containing molecular information 553 basis_name (str, optional): Path to the basis set file 554 555 Returns: 556 str: Complete basis set string in TURBOMOLE format 557 """ 558 #The basis set for atom / mol 559 basisSet = '' 560 if basis_name is None: 561 print('Error: A basis set name is needed!') 562 return basisSet 563 #If mol is provided then load the same basis set to all atoms 564 if not mol is None: 565 for i in range(mol.natoms): 566 basisSet = basisSet + Basis.readBasisSetFromFile( mol.atomicSpecies[i], basis_name) 567 elif atom is not None: 568 basisSet = Basis.readBasisSetFromFile(atom, basis_name) 569 basisSet = basisSet.replace('D+', 'E+') 570 basisSet = basisSet.replace('D-', 'E-') 571 return basisSet 572 573 #Creates the complete basis set for a given molecule along with basis dictionary 574 def createCompleteBasisSet(self, mol, basis): 575 """ 576 Create the complete basis set string for a given molecule using the basis dictionary. 577 578 Args: 579 mol: Mol object containing molecular information 580 basis (dict): Dictionary mapping atom types to basis set strings 581 582 Returns: 583 str: Complete basis set string for all atoms in the molecule 584 """ 585 totBasisSet = '' 586 if basis is None: 587 print('Error: Please provide a basis dictionary.') 588 return totBasisSet 589 else: 590 if 'all' in basis: 591 return basis['all'] 592 else: 593 #TODO Change this loop to only read the basis set for unique atoms 594 #TODO otherwise the basis set contains repeat basis sets for repeated atoms 595 #TODO The current form may(it doesn't yet) cause problem when creating basis function information 596 #TODO For even further generality let users label atoms, so that different basis may be used for same atoms with different labels 597 for i in range(mol.natoms): 598 if mol.atomicSpecies[i].lower() in basis: 599 totBasisSet = totBasisSet+ basis[mol.atomicSpecies[i].lower()]+'\n' 600 totBasisSet = totBasisSet.replace('D+', 'E+') 601 totBasisSet = totBasisSet.replace('D-', 'E-') 602 return totBasisSet 603 totBasisSet = totBasisSet.replace('D+', 'E+') 604 totBasisSet = totBasisSet.replace('D-', 'E-') 605 return totBasisSet 606 607 608 def cart2sph(l, ordering='pyscf'): 609 """ 610 Get the transformation matrix from Cartesian to real spherical harmonics for a given angular momentum. 611 612 Args: 613 l (int): Angular momentum quantum number 614 ordering (str): Ordering convention ('pyscf' or other) 615 616 Returns: 617 numpy.ndarray: Transformation matrix from Cartesian to spherical basis 618 """ 619 # This routine is used to get the transformation matrix from cartesian to real spherical basis for a given l 620 # This assumes that the basis functions are normalized. 621 # TODO: Make it work for both normalized and unnormalized basis functions. 622 623 # REFERENCE: 624 # https://theochem.github.io/horton/2.0.1/tech_ref_gaussian_basis.html#collected-notes-on-gaussian-basis-sets 625 # https://onlinelibrary.wiley.com/doi/epdf/10.1002/qua.560540202 626 # PySCF ordering: https://github.com/pyscf/pyscf/issues/1023 627 # https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics 628 l0 = np.array([[1.0]]) 629 l1 = np.array([[0, 0, 1.0], 630 [1.0, 0, 0], 631 [0, 1.0, 0]]) 632 l1_pyscf = np.array([[1.0, 0, 0], 633 [0, 1.0, 0], 634 [0, 0, 1.0]]) 635 l2 = np.array([[-0.5, 0, 0, -0.5, 0, 1.0], 636 [0, 0, 1.0, 0, 0, 0], 637 [0, 0, 0, 0, 1.0, 0], 638 [0.86602540378443864676, 0, 0, -0.86602540378443864676, 0, 0], 639 [0, 1.0, 0, 0, 0, 0], 640 ]) 641 # I don't remember what this was for (perhaps TURBOMOLE??) 642 l2_alternative =np.array([[-0.5, 0, 0, -0.5, 0, 1], 643 [0, 0, 0.7071067811865475, 0, 0, 0], 644 [0, 0, 0, 0, 0.7071067811865475, 0], 645 [0.6123724356957945, 0, 0, -0.6123724356957945, 0, 0], 646 [0, 0.7071067811865475, 0, 0, 0, 0], 647 ]) 648 l2_pyscf = np.array([[0, 1.0, 0, 0, 0, 0], 649 [0, 0, 0, 0, 1.0, 0], 650 [-0.5, 0, 0, -0.5, 0, 1.0], 651 [0, 0, 1.0, 0, 0, 0], 652 [0.86602540378443864676, 0, 0, -0.86602540378443864676, 0, 0], 653 ]) 654 l3 = np.array([[0, 0, -0.67082039324993690892, 0, 0, 0, 0, -0.67082039324993690892, 0, 1.0], 655 [-0.61237243569579452455, 0, 0, -0.27386127875258305673, 0, 1.0954451150103322269, 0, 0, 0, 0], 656 [0, -0.27386127875258305673, 0, 0, 0, 0, -0.61237243569579452455, 0, 1.0954451150103322269, 0], 657 [0, 0, 0.86602540378443864676, 0, 0, 0, 0, -0.86602540378443864676, 0, 0], 658 [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0], 659 [0.790569415042094833, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0], 660 [0, 1.0606601717798212866, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 661 ]) 662 l3_pyscf = np.array([[0, 1.0606601717798212866, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 663 [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0], 664 [0, -0.27386127875258305673, 0, 0, 0, 0, -0.61237243569579452455, 0, 1.0954451150103322269, 0], 665 [0, 0, -0.67082039324993690892, 0, 0, 0, 0, -0.67082039324993690892, 0, 1.0], 666 [-0.61237243569579452455, 0, 0, -0.27386127875258305673, 0, 1.0954451150103322269, 0, 0, 0, 0], 667 [0, 0, 0.86602540378443864676, 0, 0, 0, 0, -0.86602540378443864676, 0, 0], 668 [0.790569415042094833, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0], 669 ]) 670 l4 = np.array([[0.375, 0, 0, 0.21957751641341996535, 0, -0.87831006565367986142, 0, 0, 0, 0, 0.375, 0, -0.87831006565367986142, 0, 1.0], 671 [0, 0, -0.89642145700079522998, 0, 0, 0, 0, -0.40089186286863657703, 0, 1.19522860933439364, 0, 0, 0, 0, 0], 672 [0, 0, 0, 0, -0.40089186286863657703, 0, 0, 0, 0, 0, 0, -0.89642145700079522998, 0, 1.19522860933439364, 0], 673 [-0.5590169943749474241, 0, 0, 0, 0, 0.9819805060619657157, 0, 0, 0, 0, 0.5590169943749474241, 0, -0.9819805060619657157, 0, 0], 674 [0, -0.42257712736425828875, 0, 0, 0, 0, -0.42257712736425828875, 0, 1.1338934190276816816, 0, 0, 0, 0, 0, 0], 675 [0, 0, 0.790569415042094833, 0, 0, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0, 0], 676 [0, 0, 0, 0, 1.0606601717798212866, 0, 0, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 677 [0.73950997288745200532, 0, 0, -1.2990381056766579701, 0, 0, 0, 0, 0, 0, 0.73950997288745200532, 0, 0, 0, 0], 678 [0, 1.1180339887498948482, 0, 0, 0, 0, -1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0, 0], 679 ]) 680 l4_pyscf = np.array([[0, 1.1180339887498948482, 0, 0, 0, 0, -1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0, 0], 681 [0, 0, 0, 0, 1.0606601717798212866, 0, 0, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 682 [0, -0.42257712736425828875, 0, 0, 0, 0, -0.42257712736425828875, 0, 1.1338934190276816816, 0, 0, 0, 0, 0, 0], 683 [0, 0, 0, 0, -0.40089186286863657703, 0, 0, 0, 0, 0, 0, -0.89642145700079522998, 0, 1.19522860933439364, 0], 684 [0.375, 0, 0, 0.21957751641341996535, 0, -0.87831006565367986142, 0, 0, 0, 0, 0.375, 0, -0.87831006565367986142, 0, 1.0], 685 [0, 0, -0.89642145700079522998, 0, 0, 0, 0, -0.40089186286863657703, 0, 1.19522860933439364, 0, 0, 0, 0, 0], 686 [-0.5590169943749474241, 0, 0, 0, 0, 0.9819805060619657157, 0, 0, 0, 0, 0.5590169943749474241, 0, -0.9819805060619657157, 0, 0], 687 [0, 0, 0.790569415042094833, 0, 0, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0, 0], 688 [0.73950997288745200532, 0, 0, -1.2990381056766579701, 0, 0, 0, 0, 0, 0, 0.73950997288745200532, 0, 0, 0, 0], 689 ]) 690 l5 = np.array([[0, 0, 0.625, 0, 0, 0, 0, 0.36596252735569994226, 0, -1.0910894511799619063, 0, 0, 0, 0, 0, 0, 0.625, 0, -1.0910894511799619063, 0, 1.0], 691 [0.48412291827592711065, 0, 0, 0.21128856368212914438, 0, -1.2677313820927748663, 0, 0, 0, 0, 0.16137430609197570355, 0, -0.56694670951384084082, 0, 1.2909944487358056284, 0, 0, 0, 0, 0, 0], 692 [0, 0.16137430609197570355, 0, 0, 0, 0, 0.21128856368212914438, 0, -0.56694670951384084082, 0, 0, 0, 0, 0, 0, 0.48412291827592711065, 0, -1.2677313820927748663, 0, 1.2909944487358056284, 0], 693 [0, 0, -0.85391256382996653193, 0, 0, 0, 0, 0, 0, 1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0.85391256382996653193, 0, -1.1180339887498948482, 0, 0], 694 [0, 0, 0, 0, -0.6454972243679028142, 0, 0, 0, 0, 0, 0, -0.6454972243679028142, 0, 1.2909944487358056284, 0, 0, 0, 0, 0, 0, 0], 695 [-0.52291251658379721749, 0, 0, 0.22821773229381921394, 0, 0.91287092917527685576, 0, 0, 0, 0, 0.52291251658379721749, 0, -1.2247448713915890491, 0, 0, 0, 0, 0, 0, 0, 0], 696 [0, -0.52291251658379721749, 0, 0, 0, 0, -0.22821773229381921394, 0, 1.2247448713915890491, 0, 0, 0, 0, 0, 0, 0.52291251658379721749, 0, -0.91287092917527685576, 0, 0, 0], 697 [0, 0, 0.73950997288745200532, 0, 0, 0, 0, -1.2990381056766579701, 0, 0, 0, 0, 0, 0, 0, 0, 0.73950997288745200532, 0, 0, 0, 0], 698 [0, 0, 0, 0, 1.1180339887498948482, 0, 0, 0, 0, 0, 0, -1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0, 0, 0], 699 [0.7015607600201140098, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 1.169267933366856683, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 700 [0, 1.169267933366856683, 0, 0, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 0, 0, 0.7015607600201140098, 0, 0, 0, 0, 0], 701 ]) 702 l6 = np.array([[-0.3125, 0, 0, -0.16319780245846672329, 0, 0.97918681475080033975, 0, 0, 0, 0, -0.16319780245846672329, 0, 0.57335309036732873772, 0, -1.3055824196677337863, 0, 0, 0, 0, 0, 0, -0.3125, 0, 0.97918681475080033975, 0, -1.3055824196677337863, 0, 1.0], 703 [0, 0, 0.86356159963469679725, 0, 0, 0, 0, 0.37688918072220452831, 0, -1.6854996561581052156, 0, 0, 0, 0, 0, 0, 0.28785386654489893242, 0, -0.75377836144440905662, 0, 1.3816985594155148756, 0, 0, 0, 0, 0, 0, 0], 704 [0, 0, 0, 0, 0.28785386654489893242, 0, 0, 0, 0, 0, 0, 0.37688918072220452831, 0, -0.75377836144440905662, 0, 0, 0, 0, 0, 0, 0, 0, 0.86356159963469679725, 0, -1.6854996561581052156, 0, 1.3816985594155148756, 0], 705 [0.45285552331841995543, 0, 0, 0.078832027985861408788, 0, -1.2613124477737825406, 0, 0, 0, 0, -0.078832027985861408788, 0, 0, 0, 1.2613124477737825406, 0, 0, 0, 0, 0, 0, -0.45285552331841995543, 0, 1.2613124477737825406, 0, -1.2613124477737825406, 0, 0], 706 [0, 0.27308215547040717681, 0, 0, 0, 0, 0.26650089544451304287, 0, -0.95346258924559231545, 0, 0, 0, 0, 0, 0, 0.27308215547040717681, 0, -0.95346258924559231545, 0, 1.4564381625088382763, 0, 0, 0, 0, 0, 0, 0, 0], 707 [0, 0, -0.81924646641122153043, 0, 0, 0, 0, 0.35754847096709711829, 0, 1.0660035817780521715, 0, 0, 0, 0, 0, 0, 0.81924646641122153043, 0, -1.4301938838683884732, 0, 0, 0, 0, 0, 0, 0, 0, 0], 708 [0, 0, 0, 0, -0.81924646641122153043, 0, 0, 0, 0, 0, 0, -0.35754847096709711829, 0, 1.4301938838683884732, 0, 0, 0, 0, 0, 0, 0, 0, 0.81924646641122153043, 0, -1.0660035817780521715, 0, 0, 0], 709 [-0.49607837082461073572, 0, 0, 0.43178079981734839863, 0, 0.86356159963469679725, 0, 0, 0, 0, 0.43178079981734839863, 0, -1.5169496905422946941, 0, 0, 0, 0, 0, 0, 0, 0, -0.49607837082461073572, 0, 0.86356159963469679725, 0, 0, 0, 0], 710 [0, -0.59829302641309923139, 0, 0, 0, 0, 0, 0, 1.3055824196677337863, 0, 0, 0, 0, 0, 0, 0.59829302641309923139, 0, -1.3055824196677337863, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 711 [0, 0, 0.7015607600201140098, 0, 0, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 0, 0, 1.169267933366856683, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 712 [0, 0, 0, 0, 1.169267933366856683, 0, 0, 0, 0, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7015607600201140098, 0, 0, 0, 0, 0], 713 [0.67169328938139615748, 0, 0, -1.7539019000502850245, 0, 0, 0, 0, 0, 0, 1.7539019000502850245, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.67169328938139615748, 0, 0, 0, 0, 0, 0], 714 [0, 1.2151388809514737933, 0, 0, 0, 0, -1.9764235376052370825, 0, 0, 0, 0, 0, 0, 0, 0, 1.2151388809514737933, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 715 ]) 716 if ordering=='pyscf': 717 c2s_l = [l0,l1_pyscf,l2_pyscf,l3_pyscf,l4_pyscf,l5,l6] # PYSCF ordering of SAOs 718 else: 719 c2s_l = [l0,l1,l2,l3,l4,l5,l6] # HORTON ordering of SAOs 720 return c2s_l[l] 721 722 def cart2sph_basis(self): 723 """ 724 Return the complete Cartesian to spherical harmonic basis transformation matrix for all shells. 725 726 Returns: 727 matrix: Block diagonal transformation matrix for the entire basis set 728 """ 729 # Returns the complete Cartesian to Spherical (real) basis transformation matrix 730 cart2sph = [] 731 for i in range(self.nshells): 732 l = self.shells[i]-1 733 cart2sph.append(Basis.cart2sph(l)) 734 735 return scipy.linalg.block_diag(*cart2sph) 736 737 def sph2cart_basis(self): 738 """ 739 Return the complete spherical to Cartesian basis transformation matrix for all shells. 740 741 Returns: 742 matrix: Block diagonal transformation matrix (pseudoinverse of cart2sph) 743 """ 744 # Returns the complete Cartesian to Spherical (real) basis transformation matrix 745 sph2cart = [] 746 for i in range(self.nshells): 747 l = self.shells[i]-1 748 sph2cart_pseudo = np.linalg.pinv(Basis.cart2sph(l)) 749 sph2cart.append(sph2cart_pseudo) 750 751 return scipy.linalg.block_diag(*sph2cart) 752 753 754 # #Probably won't be used 755 # def readBasisFunctionsfromFile(self, mol, basis_name): 756 # #This function will read and parse the basis sets for the basis functions and other properties. 757 # #Returns true if successful or false if failed. 758 # if mol is None: 759 # print('Error: Please provide a Mol object.') 760 # return False 761 # if basis_name is None: 762 # print('Error: Please provide a basis set name dictionary for Mol atoms.') 763 # return False 764 # #If all atoms are assigned the same basis set 765 # if 'all' in basis_name: 766 # for i in range(mol.natoms): 767 # self.basis_name_local[mol.atomicSpecies[i]] = basis_name['all'] 768 # else: 769 # for i in range(mol.natoms): 770 # #If a particular basis set is assigned to an atom 771 # if mol.atomicSpecies[i] in basis_name: 772 # self.basis_name_local[mol.atomicSpecies[i]] = basis_name[mol.atomicSpecies[i]] 773 # else: 774 # #Or assign the default basis set 775 # self.basis_name_local[mol.atomicSpecies[i]] = basis_name_local['default'] 776 777 def readBasisFunctionInfo(self, mol, basisSet): 778 """ 779 Read and parse the basis set information to populate basis function properties. 780 781 Args: 782 mol: Mol object containing molecular information 783 basisSet (str): Complete basis set string in TURBOMOLE format 784 785 Returns: 786 None - populates internal data structures with basis function information 787 """ 788 #TODO Add error checks, 789 #TODO Like check if the basis set contains the needed atoms or not, and other checks that you can think of 790 #TODO Also, add error messages explaining what is happening 791 792 #Read the basis set and set the basis set information 793 indx_shell = 0 794 for i in range(mol.natoms): 795 lookfor = '' 796 lines = basisSet.splitlines() 797 pattern = re.compile('\*\n(.*)\n\*') 798 result = pattern.findall(basisSet) 799 if len(result)==0: 800 print('The basis set corresponding to atom ',mol.atomicSpecies[i], ' was not found in the provided basis set.') 801 for res in result: 802 if res.split()[0].lower()==mol.atomicSpecies[i].lower(): 803 lookfor = res 804 startReading = False 805 #print(lines) 806 currentLineNo = 0 807 while currentLineNo<len(lines): 808 line = lines[currentLineNo].strip() 809 #print(line) 810 if line==lookfor: 811 startReading = True 812 currentLineNo = currentLineNo + 2 813 if startReading: 814 line = lines[currentLineNo].strip() 815 splitLine = line.split() 816 if '*' in line: 817 startReading = False 818 currentLineNo = currentLineNo + 1 819 break 820 #print(line) 821 #Check if the line looks like ' 3 s' 822 if isinstance(int(splitLine[0]), int) and isinstance(str(splitLine[1]), str): 823 #print('here') 824 #Read the shell information 825 #Read the number of primitives in the given shell Ex: '3 s' 826 self.nprims.append(int(splitLine[0])) #Read '3' 827 #Read the shell type Ex: '3 s' 828 self.shellsLabel.append(str(splitLine[1])) #Read 's' 829 #Get the angular momentum 'l' corresponding to the shell type. 830 #Ex: 's'->1, 'p'->2, 'd'->3, etc. 831 l = Data.shell_dict[splitLine[1]] 832 self.shells.append(l) 833 self.shell_degen.append(Data.shell_degen[l-1]) 834 self.shell_coords.append(mol.coordsBohrs[i]) 835 #create the information for basis functions 836 #for ibf in range(int(l*(l+1)/2.0)): 837 # self.bfcoords.append(mol.coords[i]) 838 #Run a loop over the number of primitives in the current shell 839 #and read the exponents and contraction coefficients 840 for k in range(int(splitLine[0])): 841 currentLineNo = currentLineNo + 1 842 line = lines[currentLineNo].strip() 843 splitLine = line.split() 844 if '*' in line: 845 startReading = False 846 currentLineNo = currentLineNo + 1 847 break 848 splitLine[0] = splitLine[0].replace('D+', 'E+') 849 splitLine[1] = splitLine[1].replace('D-', 'E-') 850 self.prim_expnts.append(float(splitLine[0])) 851 self.prim_coeffs.append(float(splitLine[1])) 852 self.prim_coords.append(mol.coordsBohrs[i]) 853 self.prim_atoms.append(i) 854 self.prim_shells.append(indx_shell) 855 #self.prim_norms.append(float(normalizationFactor())) 856 #print(self.prim_coeffs, self.prim_expnts) 857 #print(currentLineNo) 858 859 indx_shell +=1 860 for ibf in range(int(l*(l+1)/2.0)): 861 self.bfs_atoms.append(i) 862 863 currentLineNo = currentLineNo + 1 864 865 866 867 #Now that the reading of shells and their primitives is done 868 869 #Run a loop over shells, and create information for 870 # basis functions 871 #for i in range(shells 872 873 def calc_nprim_for_each_angular_momentum_l(self, tuple_list): 874 """ 875 Calculate the number of primitives for each angular momentum quantum number. 876 877 Args: 878 tuple_list (list): List of tuples containing (angular_momentum, num_primitives) 879 880 Returns: 881 list: List of tuples with summed primitives for each unique angular momentum 882 """ 883 # Initialize the list of tuples 884 # tuple_list = [(0,12),(0,5),(1,2)] 885 # tuple_list = self.nprims_shells 886 887 # Create a dictionary to store the sums of each common angular momentum 888 sums_dict = {} 889 890 # Loop through each tuple in the list 891 for t in tuple_list: 892 # Get the first element of the tuple 893 key = t[0] 894 895 # Check if the element is already in the dictionary 896 if key in sums_dict: 897 # If it is, add the second element of the tuple to the sum 898 sums_dict[key] += t[1] 899 else: 900 # If it isn't, add the element to the dictionary with the value of the second element of the tuple 901 sums_dict[key] = t[1] 902 903 # Initialize the list of tuples to return 904 result = [] 905 906 # Loop through each key-value pair in the dictionary 907 for key, value in sums_dict.items(): 908 # Append a tuple with the key and value to the result list 909 result.append((key, value)) 910 911 # self.nprims_angmom_list = result 912 return result
Class to store and manage basis function properties for a given molecular system.
This class processes atomic basis set information for each atom in the molecule, including primitive Gaussian functions, shells, and contracted basis functions (AOs). It also supports TURBOMOLE-style shell ordering.
Basis(mol, basis, tmoleFormat=False)
41 def __init__(self, mol, basis, tmoleFormat=False): 42 """ 43 Initialize a Basis object for storing basis function properties. 44 45 Args: 46 mol: Mol object containing molecular information 47 basis: Dictionary specifying the basis set to be used for each atom, or None to use mol.basis 48 tmoleFormat (bool): Whether to use TURBOMOLE format ordering for basis functions 49 50 Returns: 51 None 52 53 Raises: 54 None - prints error message if mol is None 55 """ 56 if mol is None: 57 print('Error: Please provide a Mol object.') 58 return None 59 #Basis set name 60 #If no basis set name is provided specifically, 61 #then use the basis set name of the mol. 62 #'basis' is a dictionary specifying the basis set to be used for a particular atom 63 if basis is None: 64 self.basis = mol.basis 65 else: 66 self.basis = basis 67 68 69 self.basisSet = '' 70 """Basis set name or label. 71 Type: str 72 """ 73 74 self.nao = 0 75 """Number of basis functions (atomic orbitals). 76 Type: int 77 """ 78 79 self.nshells = 0 80 """Number of shells. 81 Type: int 82 """ 83 84 self.totalnprims = 0 85 """Total number of primitive Gaussian functions. 86 Type: int 87 """ 88 89 self.nprims = [] 90 """Number of primitives in each shell. 91 Type: List[int] 92 """ 93 94 self.prim_coeffs = [] 95 """Primitive contraction coefficients. 96 Type: List[float] 97 """ 98 #Exponents of primitives 99 self.prim_expnts = [] 100 """Exponents of primitives. 101 Type: List[float] 102 """ 103 #Coords of primitives 104 self.prim_coords = [] 105 """Coordinates of each primitive function. 106 Type: List[List[float]] 107 """ 108 #Normalization factors of primitives 109 self.prim_norms = [] 110 """Normalization constants of primitives. 111 Type: List[float] 112 """ 113 # This list contains the indices of atoms (first = #0) that correspond to a particular primitive function 114 self.prim_atoms = [] #(This list should be of the size of the number of self.totalnprims) 115 """Atom index corresponding to each primitive. 116 Type: List[int] 117 """ 118 # A list of tuples (2D list) that gives the number of primitives belonging to each atomic index 119 self.nprims_atoms = [] # Size is the same as that of no. of atoms. Looks like this: [(0, 12), (1, 5), (2, 5)] 120 """Number of primitives per atom as (atom_index, nprims). 121 Type: List[Tuple[int, int]] 122 """ 123 # This list contains the indices of shells (first = #0) that correspond to a particular primitive function 124 self.prim_shells = [] #(This list should be of the size of the number of self.totalnprims) 125 """Shell index corresponding to each primitive. 126 Type: List[int] 127 """ 128 # A list of tuples (2D list) that gives the number of primitives belonging to each shell index 129 self.nprims_shells = [] # Size is the same as that of no. of shells. Looks like this: [(0, 12), (1, 5), (2, 5)] 130 """Number of primitives per shell as (shell_index, nprims). 131 Type: List[Tuple[int, int]] 132 """ 133 # A special 2d list that contains the angular momentum of shell in the first index and the no. of primitives corresponding to the shell as second index 134 self.nprims_shell_l_list = [] # This will be of the same size as the no. of shells 135 """Angular momentum and primitive count per shell as (l, nprims). 136 Type: List[Tuple[int, int]] 137 """ 138 # A list of size natoms, that contains the values of the largest primitive exponents for all the atoms 139 self.alpha_max = [] 140 """Maximum exponent among primitives for each atom. 141 Type: List[float] 142 """ 143 # A list of size natoms, that contains the list of the smallest primitive exponents for each shell for all the atoms 144 self.alpha_min = [] 145 """Minimum exponent of primitives per shell for each atom. 146 Type: List[List[float]] 147 """ 148 #Number of contractions 149 self.ncontrs = 0 150 """Total number of contracted basis functions. 151 Type: int 152 """ 153 #Normalization factors for contractions 154 self.contrs_norm = [] 155 """Normalization factors for contracted basis functions. 156 Type: List[float] 157 """ 158 #Degeneracy of shells 159 self.shell_degen = [] 160 """Degeneracy of each shell. 161 Type: List[int] 162 """ 163 #Shells 164 self.shellsLabel = [] 165 """Label of each shell (e.g., 'S', 'P', etc.). 166 Type: List[str] 167 """ 168 self.shells = [] 169 """Shell index (0: s, 1: p, 2:d). 170 Type: List[int] 171 """ 172 self.shell_coords = [] 173 """Coordinates of each shell center. 174 Type: List[List[float]] 175 """ 176 #--------------------------------- 177 #Information for basis function 178 #--------------------------------- 179 #Coordinates of basis functions (This list should be of the size of the number of AOs/BFs) 180 self.bfs_coords = [] 181 """Coordinates of basis functions (same as their parent shell). 182 Type: List[List[float]] 183 """ 184 # Number of primitives corresponding to each basis function (This list should be of the size of the number of AOs/BFs) 185 self.bfs_nprim = [] 186 """Number of primitives corresponding to each basis function. 187 Type: List[int] 188 """ 189 # Radius cutoff for each basis function (This list should be of the size of the number of AOs/BFs) 190 # This is used to accelerate the evaluation of AOs on grid points, as those gridpoints farther than the radius cutoff 191 # of the basis function can be discarded for calcualtion. 192 self.bfs_radius_cutoff = [] 193 """Radius cutoff for each basis function to speed up AO evaluation. 194 Type: List[float] 195 """ 196 #A list of the size of number of atomic orbitals (BFs), 197 #that contains lists of primitive contraction coefficients/exponents 198 #corresponding to every basis function 199 self.bfs_expnts = [] 200 """Primitive exponents used in each basis function. 201 Type: List[List[float]] 202 """ 203 204 self.bfs_coeffs = [] 205 """Primitive coefficients for each basis function. 206 Type: List[List[float]] 207 """ 208 # A list of the size of number of AOs (BFs) 209 # that contains the lists of normalization factors for 210 # every primitive corresponding to the basis functions 211 self.bfs_prim_norms = [] 212 """Primitive normalization constants per basis function. 213 Type: List[List[float]] 214 """ 215 # Normalization factor for the contraction of primitives making up a BF (This list should be of the size of the number of AOs/BFs) 216 self.bfs_contr_prim_norms = [] 217 """Normalization for the contraction of primitives forming a BF. 218 Type: List[float] 219 """ 220 # l,m,n indices of BFs (This list should be of the size of the number of AOs/BFs) 221 self.bfs_lmn = [] 222 """Cartesian angular momentum indices (l, m, n) per basis function. 223 Type: List[Tuple[int, int, int]] 224 """ 225 # Angular momentum of a basis function (This list should be of the size of the number of AOs/BFs) 226 self.bfs_lm = [] 227 """Angular momentum `l` for each basis function. 228 Type: List[int] 229 """ 230 # Label of a basis function (This list should be of the size of the number of AOs/BFs) 231 self.bfs_label = [] 232 """Label for each basis function. 233 Type: List[str] 234 """ 235 # Number of BFs in a shell 236 self.bfs_nbfshell = [] 237 """Number of BFs in the shell of each basis function. 238 Type: List[int] 239 """ 240 # Shell index of each bf 241 self.bfs_shell_index =[] 242 """Shell index for each basis function. 243 Type: List[int] 244 """ 245 # BF offset index of each shell (This should be of the size of the number of shells) 246 self.shell_bfs_offset =[] 247 """Offset index of the first basis function in each shell. 248 Type: List[int] 249 """ 250 # Total number of basis functions (BFs) also called Atomic orbitals (AOs) 251 self.bfs_nao = [] 252 """Index of each basis function (redundant with `nao`). 253 Type: List[int] 254 """ 255 # This list contains the indices of atoms (first = #0) that correspond to a particular bf 256 self.bfs_atoms = [] #(This list should be of the size of the number of AOs/BFs) 257 """Atom index corresponding to each basis function. 258 Type: List[int] 259 """ 260 #--------------------------------- 261 #If all data was parsed successfully. 262 self.success = False 263 """Indicates if basis parsing was successful. 264 Type: bool 265 """ 266 #Convert the keys of basis dictionary to lower case to avoid ambiguity 267 self.basisSet = self.createCompleteBasisSet(mol, basis) 268 self.readBasisFunctionInfo(mol, self.basisSet) 269 self.totalnprims = sum(self.nprims) 270 self.nshells = len(self.nprims) 271 # Now lets create information about the basis functions 272 # A basis function is a combination of gaussian primitives, 273 # so we need 274 # number of primitives for a given bf, 275 # the exponents and contraction coefficients of the primitives, 276 # the l,m,n triplet information, 277 # the angular momentum l+m+n, 278 # coords, 279 # normalization factors of primitives, 280 # normalization factor for the contraction. 281 #--------------------------------------- 282 offset = 0 283 ishell = 0 284 # Loop through all the shells 285 for i in range(len(self.nprims)): 286 # Number of primitives in a shell 287 numprim = self.nprims[i] 288 # Angular momentum of a shell 289 lm = self.shells[i] - 1 290 # no of bfs in this shell 291 nbf_shell = Data.shell_degen[lm] 292 self.bfs_nbfshell.append(nbf_shell) 293 # Assign the lmn triplets and other information corresponding to basis functions 294 # by looping through the degeneracy of the shells. 295 # Example: If the degeneracy is 1 then we loop through [0,1) that is we look 296 # for the 0th element in the Data.shell_lmn dictionary as that corresponds to the 's' lmn triplet. 297 # Similarly if the denegeracy was 3 as for 'p' shell, then we loop from [1,3] and get the 298 # lmn triplets corresponding to 'p' shell and so on. 299 offset_shell_labels = Data.shell_lmn_offset[lm] 300 for j in range(offset_shell_labels, nbf_shell + offset_shell_labels): 301 #Assign the label to the basis function 302 if tmoleFormat: 303 label = list(Data.shell_lmn_tmole)[j] 304 else: 305 label = list(Data.shell_lmn)[j] 306 self.bfs_label.append(label) 307 # Assign l,m,n triplet for a basis function 308 if tmoleFormat: 309 lmn = Data.shell_lmn_tmole[label] 310 else: 311 lmn = Data.shell_lmn[label] 312 self.bfs_lmn.append(lmn) 313 # Assign the number of primitives in a basis function 314 self.bfs_nprim.append(numprim) 315 # Assign the the total angular momentum of a basis function 316 self.bfs_lm.append(sum(lmn)) 317 # Assign the exponents and contraction coefficients 318 expnts = self.prim_expnts[offset : offset+numprim ] 319 self.bfs_expnts.append(expnts) 320 coeffs = self.prim_coeffs[offset : offset+numprim ] 321 self.bfs_coeffs.append(coeffs) 322 # Assign Radius Cutoffs 323 radius_cutoff = 0.0 324 # Loop over primitives 325 eeta_log = np.log(1.0E-9) 326 for iprim_ in range(len(expnts)): 327 radius_cutoff = np.maximum(radius_cutoff, np.sqrt(1.0/expnts[iprim_]*(np.log(expnts[iprim_])/2.0 - eeta_log))) 328 self.bfs_radius_cutoff.append(radius_cutoff) 329 # Assign the normalization factors for all the primitives corresponding to a BF 330 norms = [] 331 for k in range(offset,offset+numprim): 332 norms.append(self.normalizationFactor(self.prim_expnts[k],lmn[0],lmn[1],lmn[2])) 333 self.bfs_prim_norms.append(norms) 334 self.bfs_contr_prim_norms.append(self.normalizationFactorContraction(expnts, coeffs, norms, lmn[0],lmn[1],lmn[2], lm)) 335 # Assign coordinates to the basis functions 336 self.bfs_coords.append(self.shell_coords[i]) 337 # Assign shell index to the basis functions 338 self.bfs_shell_index.append(ishell) 339 340 ishell += 1 341 offset = offset + numprim 342 self.bfs_nao = len(self.bfs_label) 343 344 ##### The following was mainly developed to calculated alpha_min and alpha_max for numgrid 345 346 # Calculate the number of primitives corresponding to each atom 347 cnt = Counter(self.prim_atoms) # Gives a Counter object similar to a dictionary 348 # Convert the Counter object to a list 349 self.nprims_atoms = list(cnt.items()) #a list of tuples, each containing a value and its count in the list 350 ### Start processing to calculate alpha_min and alpha_max for grid generation 351 ## Calculate alpha_max (https://github.com/dftlibs/numgrid) 352 self.alpha_max = [] 353 offset = 0 354 for iat in range(mol.natoms): 355 alpha_max = max(self.prim_expnts[offset : offset + self.nprims_atoms[iat][1]]) 356 self.alpha_max.append(alpha_max) 357 offset = offset + self.nprims_atoms[iat][1] 358 359 # Calculate the number of primitives corresponding to each shell 360 cnt = Counter(self.prim_shells) # Gives a Counter object similar to a dictionary 361 # Convert the Counter object to a list 362 self.nprims_shells = list(cnt.items()) #a list of tuples, each containing a value and its count in the list 363 # Make a list that is pretty much the same as nprims_shells, however, it contains the shell angular momentum in the first index 364 for ishell in range(self.nshells): 365 self.nprims_shell_l_list.append((self.shells[ishell], self.nprims_shells[ishell][1])) 366 # Create self.nprims_angmom_list 367 # self.calc_nprim_for_each_angular_momentum_l() 368 ## Calculate alpha_min 369 ishell_offset = 0 370 iprim_offset = 0 371 for iat in range(mol.natoms): 372 alpha_min_dict_iatom = {} 373 nshell_atom = 0 374 sum_prim = 0 375 for ishell in range(ishell_offset, self.nshells): 376 nshell_atom = nshell_atom + 1 377 sum_prim = sum_prim + self.nprims_shells[ishell][1] 378 if sum_prim==self.nprims_atoms[iat][1]: 379 break 380 subset_nprims_shell_l_list = self.nprims_shell_l_list[ishell_offset : ishell_offset + nshell_atom] 381 ishell_offset = ishell_offset + nshell_atom 382 # Create self.nprims_angmom_list 383 subset_nprims_angmom_list = self.calc_nprim_for_each_angular_momentum_l(subset_nprims_shell_l_list) 384 # alpha_min 385 for i in range(len(subset_nprims_angmom_list)): 386 alpha_min_dict_iatom[subset_nprims_angmom_list[i][0] - 1] = min(self.prim_expnts[iprim_offset : iprim_offset + subset_nprims_angmom_list[i][1]]) 387 iprim_offset = iprim_offset + subset_nprims_angmom_list[i][1] 388 389 self.alpha_min.append(alpha_min_dict_iatom) 390 391 # Calculate the BF offset index for a given shell 392 self.shell_bfs_offset = np.cumsum(self.bfs_nbfshell) - self.bfs_nbfshell
Initialize a Basis object for storing basis function properties.
Args:
mol: Mol object containing molecular information
basis: Dictionary specifying the basis set to be used for each atom, or None to use mol.basis
tmoleFormat (bool): Whether to use TURBOMOLE format ordering for basis functions
Returns: None
Raises: None - prints error message if mol is None
nprims_shell_l_list
Angular momentum and primitive count per shell as (l, nprims). Type: List[Tuple[int, int]]
bfs_radius_cutoff
Radius cutoff for each basis function to speed up AO evaluation. Type: List[float]
bfs_contr_prim_norms
Normalization for the contraction of primitives forming a BF. Type: List[float]
bfs_lmn
Cartesian angular momentum indices (l, m, n) per basis function. Type: List[Tuple[int, int, int]]
def
normalizationFactor(self, alpha, l, m, n):
398 def normalizationFactor(self, alpha, l, m, n): 399 """ 400 Calculate the normalization factor for a primitive Gaussian function. 401 402 Args: 403 alpha (float): Exponent of the Gaussian primitive 404 l (int): Angular momentum quantum number in x-direction 405 m (int): Angular momentum quantum number in y-direction 406 n (int): Angular momentum quantum number in z-direction 407 408 Returns: 409 float: Normalization factor for the primitive Gaussian 410 """ 411 return np.sqrt((2*alpha/np.pi)**(3/2)*(4*alpha)**(l+m+n)/(factorial2(2*l-1)*factorial2(2*m-1)*factorial2(2*n-1)))
Calculate the normalization factor for a primitive Gaussian function.
Args:
alpha (float): Exponent of the Gaussian primitive
l (int): Angular momentum quantum number in x-direction
m (int): Angular momentum quantum number in y-direction
n (int): Angular momentum quantum number in z-direction
Returns: float: Normalization factor for the primitive Gaussian
def
normalizationFactorContraction(self, alphas, coeffs, norms, l, m, n, lm):
414 def normalizationFactorContraction(self, alphas, coeffs, norms, l, m, n, lm): 415 """ 416 Calculate the normalization factor for a contraction of Gaussian primitives. 417 418 Args: 419 alphas (list): List of exponents for the primitive Gaussians 420 coeffs (list): List of contraction coefficients 421 norms (list): List of normalization factors for individual primitives 422 l (int): Angular momentum quantum number in x-direction 423 m (int): Angular momentum quantum number in y-direction 424 n (int): Angular momentum quantum number in z-direction 425 lm (int): Total angular momentum (l+m+n) 426 427 Returns: 428 float: Normalization factor for the contracted Gaussian 429 """ 430 431 temp = np.pi**(3/2)*factorial2(2*l-1)*factorial2(2*m-1)*factorial2(2*n-1)/(2**lm) 432 sum = 0.0 433 for i in range(len(alphas)): 434 for j in range(len(alphas)): 435 sum = sum + (coeffs[i]*coeffs[j]*norms[i]*norms[j]/((alphas[i]+alphas[j])**(lm+3/2))) 436 factor = np.power(temp*sum,-1/2) 437 return factor
Calculate the normalization factor for a contraction of Gaussian primitives.
Args: alphas (list): List of exponents for the primitive Gaussians coeffs (list): List of contraction coefficients norms (list): List of normalization factors for individual primitives l (int): Angular momentum quantum number in x-direction m (int): Angular momentum quantum number in y-direction n (int): Angular momentum quantum number in z-direction lm (int): Total angular momentum (l+m+n)
Returns: float: Normalization factor for the contracted Gaussian
def
readBasisSetFromFile(key, filename):
440 def readBasisSetFromFile(key, filename): 441 """ 442 Read the basis set block corresponding to a particular atom from a TURBOMOLE format file. 443 444 Args: 445 key (str): Atomic species symbol to search for 446 filename (str): Path to the basis set file 447 448 Returns: 449 str or bool: Basis set string for the atom, or False if not found 450 """ 451 #This functions returns the basis set block corresponding to a particular atom (key) 452 #The file to be read from is assumed to follow TURBOMOLE's format 453 basisString = '' 454 lookfor = '' 455 file = open(filename, 'r') 456 fileContentsInString=file.read() 457 file.close() 458 lines = fileContentsInString.splitlines() 459 pattern = re.compile('\*\n(.*)\n\*') 460 result = pattern.findall(fileContentsInString) 461 if len(result)==0: 462 print('The basis set corresponding to atom',key, ' was not found in ',filename) 463 return False 464 for res in result: 465 if res.split()[0].lower()==key.lower(): 466 lookfor = res 467 basisString = '\n*\n'+lookfor+'\n*' 468 469 currentLineNo = 0 470 startReading = False 471 for line in lines: 472 line = line.strip() 473 if line==lookfor: 474 startReading = True 475 currentLineNo = 1 476 if startReading and currentLineNo>=3: 477 if line.strip()=='*': 478 break 479 else: 480 basisString = basisString + '\n'+line 481 currentLineNo = currentLineNo + 1 482 basisString = basisString + '\n*' 483 return basisString
Read the basis set block corresponding to a particular atom from a TURBOMOLE format file.
Args: key (str): Atomic species symbol to search for filename (str): Path to the basis set file
Returns: str or bool: Basis set string for the atom, or False if not found
def
load(atom=None, mol=None, basis_name=None):
491 def load( atom=None, mol=None, basis_name=None): 492 """ 493 Load the complete basis set as a string for a given atom/molecule from the CrysX library. 494 495 Args: 496 atom (str, optional): Atomic species symbol 497 mol (Mol, optional): Mol object containing molecular information 498 basis_name (str, optional): Name of the basis set to load 499 500 Returns: 501 str: Complete basis set string in TURBOMOLE format 502 """ 503 # The CrysX library contains basis sets in the TURBOMOLE format. 504 # The basis sets are downloaded from https://www.basissetexchange.org 505 # BSSE Github: https://github.com/MolSSI-BSE/basis_set_exchange 506 # License of BSSE: BSD 3 507 508 # We (CrysX) have saved the basis sets in the 'BasisSets' directory. 509 # The basis set name should be in lower-case and followed by '.version' and '.tm' extension. 510 # Ex: def2-tzvp.1.tm 511 # We can't expect the user to enter the basis set name with such an extension 512 # or in lower or upper case. So we should process this input name. 513 basis_name = basis_name.strip() 514 #The basis set for atom / mol 515 basisSet = basis_name+'\n' 516 basis_name = basis_name.replace(" ", "") 517 basis_name = basis_name.lower() 518 basis_name = basis_name+'.1.tm' 519 # Here we need to create the path name to the BasisSets directory supplied with crysx 520 # Since, the path should be with respect to hte module's location and not the working directory we can follow 521 # this link for various methods: https://stackoverflow.com/questions/4060221/how-to-reliably-open-a-file-in-the-same-directory-as-a-python-script 522 # Method 1: __file__ 523 # Method 2: sys.argv[0] 524 # Method 3: sys.path[0] (probably not relevant to our purpose) 525 # Currently we use method 1 526 temp = __file__ 527 temp = temp.replace('Basis.py', '') 528 basis_name = temp + 'BasisSets/'+basis_name 529 #basis_name = 'BasisSets/'+basis_name 530 531 if basis_name is None: 532 print('Error: A basis set name is needed!') 533 return basisSet 534 #If mol is provided then load the same basis set to all atoms 535 if not mol is None: 536 for i in range(mol.natoms): 537 basisSet = basisSet + Basis.readBasisSetFromFile( mol.atomicSpecies[i], basis_name) 538 elif atom is not None: 539 basisSet = Basis.readBasisSetFromFile(atom, basis_name) 540 basisSet = basisSet.replace('D+', 'E+') 541 basisSet = basisSet.replace('D-', 'E-') 542 return basisSet
Load the complete basis set as a string for a given atom/molecule from the CrysX library.
Args: atom (str, optional): Atomic species symbol mol (Mol, optional): Mol object containing molecular information basis_name (str, optional): Name of the basis set to load
Returns: str: Complete basis set string in TURBOMOLE format
def
loadfromfile(atom=None, mol=None, basis_name=None):
546 def loadfromfile( atom=None, mol=None, basis_name=None): 547 """ 548 Load the complete basis set as a string for a given atom/molecule from a user-specified file. 549 550 Args: 551 atom (str, optional): Atomic species symbol 552 mol (Mol, optional): Mol object containing molecular information 553 basis_name (str, optional): Path to the basis set file 554 555 Returns: 556 str: Complete basis set string in TURBOMOLE format 557 """ 558 #The basis set for atom / mol 559 basisSet = '' 560 if basis_name is None: 561 print('Error: A basis set name is needed!') 562 return basisSet 563 #If mol is provided then load the same basis set to all atoms 564 if not mol is None: 565 for i in range(mol.natoms): 566 basisSet = basisSet + Basis.readBasisSetFromFile( mol.atomicSpecies[i], basis_name) 567 elif atom is not None: 568 basisSet = Basis.readBasisSetFromFile(atom, basis_name) 569 basisSet = basisSet.replace('D+', 'E+') 570 basisSet = basisSet.replace('D-', 'E-') 571 return basisSet
Load the complete basis set as a string for a given atom/molecule from a user-specified file.
Args:
atom (str, optional): Atomic species symbol
mol (Mol, optional): Mol object containing molecular information
basis_name (str, optional): Path to the basis set file
Returns: str: Complete basis set string in TURBOMOLE format
def
createCompleteBasisSet(self, mol, basis):
574 def createCompleteBasisSet(self, mol, basis): 575 """ 576 Create the complete basis set string for a given molecule using the basis dictionary. 577 578 Args: 579 mol: Mol object containing molecular information 580 basis (dict): Dictionary mapping atom types to basis set strings 581 582 Returns: 583 str: Complete basis set string for all atoms in the molecule 584 """ 585 totBasisSet = '' 586 if basis is None: 587 print('Error: Please provide a basis dictionary.') 588 return totBasisSet 589 else: 590 if 'all' in basis: 591 return basis['all'] 592 else: 593 #TODO Change this loop to only read the basis set for unique atoms 594 #TODO otherwise the basis set contains repeat basis sets for repeated atoms 595 #TODO The current form may(it doesn't yet) cause problem when creating basis function information 596 #TODO For even further generality let users label atoms, so that different basis may be used for same atoms with different labels 597 for i in range(mol.natoms): 598 if mol.atomicSpecies[i].lower() in basis: 599 totBasisSet = totBasisSet+ basis[mol.atomicSpecies[i].lower()]+'\n' 600 totBasisSet = totBasisSet.replace('D+', 'E+') 601 totBasisSet = totBasisSet.replace('D-', 'E-') 602 return totBasisSet 603 totBasisSet = totBasisSet.replace('D+', 'E+') 604 totBasisSet = totBasisSet.replace('D-', 'E-') 605 return totBasisSet
Create the complete basis set string for a given molecule using the basis dictionary.
Args: mol: Mol object containing molecular information basis (dict): Dictionary mapping atom types to basis set strings
Returns: str: Complete basis set string for all atoms in the molecule
def
cart2sph(l, ordering='pyscf'):
608 def cart2sph(l, ordering='pyscf'): 609 """ 610 Get the transformation matrix from Cartesian to real spherical harmonics for a given angular momentum. 611 612 Args: 613 l (int): Angular momentum quantum number 614 ordering (str): Ordering convention ('pyscf' or other) 615 616 Returns: 617 numpy.ndarray: Transformation matrix from Cartesian to spherical basis 618 """ 619 # This routine is used to get the transformation matrix from cartesian to real spherical basis for a given l 620 # This assumes that the basis functions are normalized. 621 # TODO: Make it work for both normalized and unnormalized basis functions. 622 623 # REFERENCE: 624 # https://theochem.github.io/horton/2.0.1/tech_ref_gaussian_basis.html#collected-notes-on-gaussian-basis-sets 625 # https://onlinelibrary.wiley.com/doi/epdf/10.1002/qua.560540202 626 # PySCF ordering: https://github.com/pyscf/pyscf/issues/1023 627 # https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics 628 l0 = np.array([[1.0]]) 629 l1 = np.array([[0, 0, 1.0], 630 [1.0, 0, 0], 631 [0, 1.0, 0]]) 632 l1_pyscf = np.array([[1.0, 0, 0], 633 [0, 1.0, 0], 634 [0, 0, 1.0]]) 635 l2 = np.array([[-0.5, 0, 0, -0.5, 0, 1.0], 636 [0, 0, 1.0, 0, 0, 0], 637 [0, 0, 0, 0, 1.0, 0], 638 [0.86602540378443864676, 0, 0, -0.86602540378443864676, 0, 0], 639 [0, 1.0, 0, 0, 0, 0], 640 ]) 641 # I don't remember what this was for (perhaps TURBOMOLE??) 642 l2_alternative =np.array([[-0.5, 0, 0, -0.5, 0, 1], 643 [0, 0, 0.7071067811865475, 0, 0, 0], 644 [0, 0, 0, 0, 0.7071067811865475, 0], 645 [0.6123724356957945, 0, 0, -0.6123724356957945, 0, 0], 646 [0, 0.7071067811865475, 0, 0, 0, 0], 647 ]) 648 l2_pyscf = np.array([[0, 1.0, 0, 0, 0, 0], 649 [0, 0, 0, 0, 1.0, 0], 650 [-0.5, 0, 0, -0.5, 0, 1.0], 651 [0, 0, 1.0, 0, 0, 0], 652 [0.86602540378443864676, 0, 0, -0.86602540378443864676, 0, 0], 653 ]) 654 l3 = np.array([[0, 0, -0.67082039324993690892, 0, 0, 0, 0, -0.67082039324993690892, 0, 1.0], 655 [-0.61237243569579452455, 0, 0, -0.27386127875258305673, 0, 1.0954451150103322269, 0, 0, 0, 0], 656 [0, -0.27386127875258305673, 0, 0, 0, 0, -0.61237243569579452455, 0, 1.0954451150103322269, 0], 657 [0, 0, 0.86602540378443864676, 0, 0, 0, 0, -0.86602540378443864676, 0, 0], 658 [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0], 659 [0.790569415042094833, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0], 660 [0, 1.0606601717798212866, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 661 ]) 662 l3_pyscf = np.array([[0, 1.0606601717798212866, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 663 [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0], 664 [0, -0.27386127875258305673, 0, 0, 0, 0, -0.61237243569579452455, 0, 1.0954451150103322269, 0], 665 [0, 0, -0.67082039324993690892, 0, 0, 0, 0, -0.67082039324993690892, 0, 1.0], 666 [-0.61237243569579452455, 0, 0, -0.27386127875258305673, 0, 1.0954451150103322269, 0, 0, 0, 0], 667 [0, 0, 0.86602540378443864676, 0, 0, 0, 0, -0.86602540378443864676, 0, 0], 668 [0.790569415042094833, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0], 669 ]) 670 l4 = np.array([[0.375, 0, 0, 0.21957751641341996535, 0, -0.87831006565367986142, 0, 0, 0, 0, 0.375, 0, -0.87831006565367986142, 0, 1.0], 671 [0, 0, -0.89642145700079522998, 0, 0, 0, 0, -0.40089186286863657703, 0, 1.19522860933439364, 0, 0, 0, 0, 0], 672 [0, 0, 0, 0, -0.40089186286863657703, 0, 0, 0, 0, 0, 0, -0.89642145700079522998, 0, 1.19522860933439364, 0], 673 [-0.5590169943749474241, 0, 0, 0, 0, 0.9819805060619657157, 0, 0, 0, 0, 0.5590169943749474241, 0, -0.9819805060619657157, 0, 0], 674 [0, -0.42257712736425828875, 0, 0, 0, 0, -0.42257712736425828875, 0, 1.1338934190276816816, 0, 0, 0, 0, 0, 0], 675 [0, 0, 0.790569415042094833, 0, 0, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0, 0], 676 [0, 0, 0, 0, 1.0606601717798212866, 0, 0, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 677 [0.73950997288745200532, 0, 0, -1.2990381056766579701, 0, 0, 0, 0, 0, 0, 0.73950997288745200532, 0, 0, 0, 0], 678 [0, 1.1180339887498948482, 0, 0, 0, 0, -1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0, 0], 679 ]) 680 l4_pyscf = np.array([[0, 1.1180339887498948482, 0, 0, 0, 0, -1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0, 0], 681 [0, 0, 0, 0, 1.0606601717798212866, 0, 0, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0], 682 [0, -0.42257712736425828875, 0, 0, 0, 0, -0.42257712736425828875, 0, 1.1338934190276816816, 0, 0, 0, 0, 0, 0], 683 [0, 0, 0, 0, -0.40089186286863657703, 0, 0, 0, 0, 0, 0, -0.89642145700079522998, 0, 1.19522860933439364, 0], 684 [0.375, 0, 0, 0.21957751641341996535, 0, -0.87831006565367986142, 0, 0, 0, 0, 0.375, 0, -0.87831006565367986142, 0, 1.0], 685 [0, 0, -0.89642145700079522998, 0, 0, 0, 0, -0.40089186286863657703, 0, 1.19522860933439364, 0, 0, 0, 0, 0], 686 [-0.5590169943749474241, 0, 0, 0, 0, 0.9819805060619657157, 0, 0, 0, 0, 0.5590169943749474241, 0, -0.9819805060619657157, 0, 0], 687 [0, 0, 0.790569415042094833, 0, 0, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0, 0], 688 [0.73950997288745200532, 0, 0, -1.2990381056766579701, 0, 0, 0, 0, 0, 0, 0.73950997288745200532, 0, 0, 0, 0], 689 ]) 690 l5 = np.array([[0, 0, 0.625, 0, 0, 0, 0, 0.36596252735569994226, 0, -1.0910894511799619063, 0, 0, 0, 0, 0, 0, 0.625, 0, -1.0910894511799619063, 0, 1.0], 691 [0.48412291827592711065, 0, 0, 0.21128856368212914438, 0, -1.2677313820927748663, 0, 0, 0, 0, 0.16137430609197570355, 0, -0.56694670951384084082, 0, 1.2909944487358056284, 0, 0, 0, 0, 0, 0], 692 [0, 0.16137430609197570355, 0, 0, 0, 0, 0.21128856368212914438, 0, -0.56694670951384084082, 0, 0, 0, 0, 0, 0, 0.48412291827592711065, 0, -1.2677313820927748663, 0, 1.2909944487358056284, 0], 693 [0, 0, -0.85391256382996653193, 0, 0, 0, 0, 0, 0, 1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0.85391256382996653193, 0, -1.1180339887498948482, 0, 0], 694 [0, 0, 0, 0, -0.6454972243679028142, 0, 0, 0, 0, 0, 0, -0.6454972243679028142, 0, 1.2909944487358056284, 0, 0, 0, 0, 0, 0, 0], 695 [-0.52291251658379721749, 0, 0, 0.22821773229381921394, 0, 0.91287092917527685576, 0, 0, 0, 0, 0.52291251658379721749, 0, -1.2247448713915890491, 0, 0, 0, 0, 0, 0, 0, 0], 696 [0, -0.52291251658379721749, 0, 0, 0, 0, -0.22821773229381921394, 0, 1.2247448713915890491, 0, 0, 0, 0, 0, 0, 0.52291251658379721749, 0, -0.91287092917527685576, 0, 0, 0], 697 [0, 0, 0.73950997288745200532, 0, 0, 0, 0, -1.2990381056766579701, 0, 0, 0, 0, 0, 0, 0, 0, 0.73950997288745200532, 0, 0, 0, 0], 698 [0, 0, 0, 0, 1.1180339887498948482, 0, 0, 0, 0, 0, 0, -1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0, 0, 0], 699 [0.7015607600201140098, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 1.169267933366856683, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 700 [0, 1.169267933366856683, 0, 0, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 0, 0, 0.7015607600201140098, 0, 0, 0, 0, 0], 701 ]) 702 l6 = np.array([[-0.3125, 0, 0, -0.16319780245846672329, 0, 0.97918681475080033975, 0, 0, 0, 0, -0.16319780245846672329, 0, 0.57335309036732873772, 0, -1.3055824196677337863, 0, 0, 0, 0, 0, 0, -0.3125, 0, 0.97918681475080033975, 0, -1.3055824196677337863, 0, 1.0], 703 [0, 0, 0.86356159963469679725, 0, 0, 0, 0, 0.37688918072220452831, 0, -1.6854996561581052156, 0, 0, 0, 0, 0, 0, 0.28785386654489893242, 0, -0.75377836144440905662, 0, 1.3816985594155148756, 0, 0, 0, 0, 0, 0, 0], 704 [0, 0, 0, 0, 0.28785386654489893242, 0, 0, 0, 0, 0, 0, 0.37688918072220452831, 0, -0.75377836144440905662, 0, 0, 0, 0, 0, 0, 0, 0, 0.86356159963469679725, 0, -1.6854996561581052156, 0, 1.3816985594155148756, 0], 705 [0.45285552331841995543, 0, 0, 0.078832027985861408788, 0, -1.2613124477737825406, 0, 0, 0, 0, -0.078832027985861408788, 0, 0, 0, 1.2613124477737825406, 0, 0, 0, 0, 0, 0, -0.45285552331841995543, 0, 1.2613124477737825406, 0, -1.2613124477737825406, 0, 0], 706 [0, 0.27308215547040717681, 0, 0, 0, 0, 0.26650089544451304287, 0, -0.95346258924559231545, 0, 0, 0, 0, 0, 0, 0.27308215547040717681, 0, -0.95346258924559231545, 0, 1.4564381625088382763, 0, 0, 0, 0, 0, 0, 0, 0], 707 [0, 0, -0.81924646641122153043, 0, 0, 0, 0, 0.35754847096709711829, 0, 1.0660035817780521715, 0, 0, 0, 0, 0, 0, 0.81924646641122153043, 0, -1.4301938838683884732, 0, 0, 0, 0, 0, 0, 0, 0, 0], 708 [0, 0, 0, 0, -0.81924646641122153043, 0, 0, 0, 0, 0, 0, -0.35754847096709711829, 0, 1.4301938838683884732, 0, 0, 0, 0, 0, 0, 0, 0, 0.81924646641122153043, 0, -1.0660035817780521715, 0, 0, 0], 709 [-0.49607837082461073572, 0, 0, 0.43178079981734839863, 0, 0.86356159963469679725, 0, 0, 0, 0, 0.43178079981734839863, 0, -1.5169496905422946941, 0, 0, 0, 0, 0, 0, 0, 0, -0.49607837082461073572, 0, 0.86356159963469679725, 0, 0, 0, 0], 710 [0, -0.59829302641309923139, 0, 0, 0, 0, 0, 0, 1.3055824196677337863, 0, 0, 0, 0, 0, 0, 0.59829302641309923139, 0, -1.3055824196677337863, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 711 [0, 0, 0.7015607600201140098, 0, 0, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 0, 0, 1.169267933366856683, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 712 [0, 0, 0, 0, 1.169267933366856683, 0, 0, 0, 0, 0, 0, -1.5309310892394863114, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7015607600201140098, 0, 0, 0, 0, 0], 713 [0.67169328938139615748, 0, 0, -1.7539019000502850245, 0, 0, 0, 0, 0, 0, 1.7539019000502850245, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.67169328938139615748, 0, 0, 0, 0, 0, 0], 714 [0, 1.2151388809514737933, 0, 0, 0, 0, -1.9764235376052370825, 0, 0, 0, 0, 0, 0, 0, 0, 1.2151388809514737933, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 715 ]) 716 if ordering=='pyscf': 717 c2s_l = [l0,l1_pyscf,l2_pyscf,l3_pyscf,l4_pyscf,l5,l6] # PYSCF ordering of SAOs 718 else: 719 c2s_l = [l0,l1,l2,l3,l4,l5,l6] # HORTON ordering of SAOs 720 return c2s_l[l]
Get the transformation matrix from Cartesian to real spherical harmonics for a given angular momentum.
Args: l (int): Angular momentum quantum number ordering (str): Ordering convention ('pyscf' or other)
Returns: numpy.ndarray: Transformation matrix from Cartesian to spherical basis
def
cart2sph_basis(self):
722 def cart2sph_basis(self): 723 """ 724 Return the complete Cartesian to spherical harmonic basis transformation matrix for all shells. 725 726 Returns: 727 matrix: Block diagonal transformation matrix for the entire basis set 728 """ 729 # Returns the complete Cartesian to Spherical (real) basis transformation matrix 730 cart2sph = [] 731 for i in range(self.nshells): 732 l = self.shells[i]-1 733 cart2sph.append(Basis.cart2sph(l)) 734 735 return scipy.linalg.block_diag(*cart2sph)
Return the complete Cartesian to spherical harmonic basis transformation matrix for all shells.
Returns: matrix: Block diagonal transformation matrix for the entire basis set
def
sph2cart_basis(self):
737 def sph2cart_basis(self): 738 """ 739 Return the complete spherical to Cartesian basis transformation matrix for all shells. 740 741 Returns: 742 matrix: Block diagonal transformation matrix (pseudoinverse of cart2sph) 743 """ 744 # Returns the complete Cartesian to Spherical (real) basis transformation matrix 745 sph2cart = [] 746 for i in range(self.nshells): 747 l = self.shells[i]-1 748 sph2cart_pseudo = np.linalg.pinv(Basis.cart2sph(l)) 749 sph2cart.append(sph2cart_pseudo) 750 751 return scipy.linalg.block_diag(*sph2cart)
Return the complete spherical to Cartesian basis transformation matrix for all shells.
Returns: matrix: Block diagonal transformation matrix (pseudoinverse of cart2sph)
def
readBasisFunctionInfo(self, mol, basisSet):
777 def readBasisFunctionInfo(self, mol, basisSet): 778 """ 779 Read and parse the basis set information to populate basis function properties. 780 781 Args: 782 mol: Mol object containing molecular information 783 basisSet (str): Complete basis set string in TURBOMOLE format 784 785 Returns: 786 None - populates internal data structures with basis function information 787 """ 788 #TODO Add error checks, 789 #TODO Like check if the basis set contains the needed atoms or not, and other checks that you can think of 790 #TODO Also, add error messages explaining what is happening 791 792 #Read the basis set and set the basis set information 793 indx_shell = 0 794 for i in range(mol.natoms): 795 lookfor = '' 796 lines = basisSet.splitlines() 797 pattern = re.compile('\*\n(.*)\n\*') 798 result = pattern.findall(basisSet) 799 if len(result)==0: 800 print('The basis set corresponding to atom ',mol.atomicSpecies[i], ' was not found in the provided basis set.') 801 for res in result: 802 if res.split()[0].lower()==mol.atomicSpecies[i].lower(): 803 lookfor = res 804 startReading = False 805 #print(lines) 806 currentLineNo = 0 807 while currentLineNo<len(lines): 808 line = lines[currentLineNo].strip() 809 #print(line) 810 if line==lookfor: 811 startReading = True 812 currentLineNo = currentLineNo + 2 813 if startReading: 814 line = lines[currentLineNo].strip() 815 splitLine = line.split() 816 if '*' in line: 817 startReading = False 818 currentLineNo = currentLineNo + 1 819 break 820 #print(line) 821 #Check if the line looks like ' 3 s' 822 if isinstance(int(splitLine[0]), int) and isinstance(str(splitLine[1]), str): 823 #print('here') 824 #Read the shell information 825 #Read the number of primitives in the given shell Ex: '3 s' 826 self.nprims.append(int(splitLine[0])) #Read '3' 827 #Read the shell type Ex: '3 s' 828 self.shellsLabel.append(str(splitLine[1])) #Read 's' 829 #Get the angular momentum 'l' corresponding to the shell type. 830 #Ex: 's'->1, 'p'->2, 'd'->3, etc. 831 l = Data.shell_dict[splitLine[1]] 832 self.shells.append(l) 833 self.shell_degen.append(Data.shell_degen[l-1]) 834 self.shell_coords.append(mol.coordsBohrs[i]) 835 #create the information for basis functions 836 #for ibf in range(int(l*(l+1)/2.0)): 837 # self.bfcoords.append(mol.coords[i]) 838 #Run a loop over the number of primitives in the current shell 839 #and read the exponents and contraction coefficients 840 for k in range(int(splitLine[0])): 841 currentLineNo = currentLineNo + 1 842 line = lines[currentLineNo].strip() 843 splitLine = line.split() 844 if '*' in line: 845 startReading = False 846 currentLineNo = currentLineNo + 1 847 break 848 splitLine[0] = splitLine[0].replace('D+', 'E+') 849 splitLine[1] = splitLine[1].replace('D-', 'E-') 850 self.prim_expnts.append(float(splitLine[0])) 851 self.prim_coeffs.append(float(splitLine[1])) 852 self.prim_coords.append(mol.coordsBohrs[i]) 853 self.prim_atoms.append(i) 854 self.prim_shells.append(indx_shell) 855 #self.prim_norms.append(float(normalizationFactor())) 856 #print(self.prim_coeffs, self.prim_expnts) 857 #print(currentLineNo) 858 859 indx_shell +=1 860 for ibf in range(int(l*(l+1)/2.0)): 861 self.bfs_atoms.append(i) 862 863 currentLineNo = currentLineNo + 1 864 865 866 867 #Now that the reading of shells and their primitives is done 868 869 #Run a loop over shells, and create information for 870 # basis functions 871 #for i in range(shells
Read and parse the basis set information to populate basis function properties.
Args: mol: Mol object containing molecular information basisSet (str): Complete basis set string in TURBOMOLE format
Returns: None - populates internal data structures with basis function information
def
calc_nprim_for_each_angular_momentum_l(self, tuple_list):
873 def calc_nprim_for_each_angular_momentum_l(self, tuple_list): 874 """ 875 Calculate the number of primitives for each angular momentum quantum number. 876 877 Args: 878 tuple_list (list): List of tuples containing (angular_momentum, num_primitives) 879 880 Returns: 881 list: List of tuples with summed primitives for each unique angular momentum 882 """ 883 # Initialize the list of tuples 884 # tuple_list = [(0,12),(0,5),(1,2)] 885 # tuple_list = self.nprims_shells 886 887 # Create a dictionary to store the sums of each common angular momentum 888 sums_dict = {} 889 890 # Loop through each tuple in the list 891 for t in tuple_list: 892 # Get the first element of the tuple 893 key = t[0] 894 895 # Check if the element is already in the dictionary 896 if key in sums_dict: 897 # If it is, add the second element of the tuple to the sum 898 sums_dict[key] += t[1] 899 else: 900 # If it isn't, add the element to the dictionary with the value of the second element of the tuple 901 sums_dict[key] = t[1] 902 903 # Initialize the list of tuples to return 904 result = [] 905 906 # Loop through each key-value pair in the dictionary 907 for key, value in sums_dict.items(): 908 # Append a tuple with the key and value to the result list 909 result.append((key, value)) 910 911 # self.nprims_angmom_list = result 912 return result
Calculate the number of primitives for each angular momentum quantum number.
Args: tuple_list (list): List of tuples containing (angular_momentum, num_primitives)
Returns: list: List of tuples with summed primitives for each unique angular momentum