pyfock.Grids

  1__all__ = ["Grids"]
  2# Grids.py
  3# Author: Manas Sharma ([email protected])
  4# This is a part of CrysX (https://bragitoff.com/crysx)
  5#
  6#
  7#  .d8888b.                            Y88b   d88P       8888888b.           8888888888                888      
  8# d88P  Y88b                            Y88b d88P        888   Y88b          888                       888      
  9# 888    888                             Y88o88P         888    888          888                       888      
 10# 888        888d888 888  888 .d8888b     Y888P          888   d88P 888  888 8888888  .d88b.   .d8888b 888  888 
 11# 888        888P"   888  888 88K         d888b          8888888P"  888  888 888     d88""88b d88P"    888 .88P 
 12# 888    888 888     888  888 "Y8888b.   d88888b  888888 888        888  888 888     888  888 888      888888K  
 13# Y88b  d88P 888     Y88b 888      X88  d88P Y88b        888        Y88b 888 888     Y88..88P Y88b.    888 "88b 
 14#  "Y8888P"  888      "Y88888  88888P' d88P   Y88b       888         "Y88888 888      "Y88P"   "Y8888P 888  888 
 15#                         888                                            888                                    
 16#                    Y8b d88P                                       Y8b d88P                                    
 17#                     "Y88P"                                         "Y88P"                                       
 18import numpy as np
 19from . import Data
 20from . import Mol
 21from . import Basis
 22import numgrid
 23from joblib import Parallel, delayed
 24#import multiprocessing
 25import os
 26from timeit import default_timer as timer
 27
 28
 29#TODO Change it to return actual cores rather than threads
 30#num_cores = multiprocessing.cpu_count()
 31num_cores = os.cpu_count()
 32"""Number of CPU cores to use for parallel grid generation via Joblib's threading backend."""
 33
 34lebedevOrdering = {
 35    0  : 1   ,
 36    3  : 6   ,
 37    5  : 14  ,
 38    7  : 26  ,
 39    9  : 38  ,
 40    11 : 50  ,
 41    13 : 74  ,
 42    15 : 86  ,
 43    17 : 110 ,
 44    19 : 146 ,
 45    21 : 170 ,
 46    23 : 194 ,
 47    25 : 230 ,
 48    27 : 266 ,
 49    29 : 302 ,
 50    31 : 350 ,
 51    35 : 434 ,
 52    41 : 590 ,
 53    47 : 770 ,
 54    53 : 974 ,
 55    59 : 1202,
 56    65 : 1454,
 57    71 : 1730,
 58    77 : 2030,
 59    83 : 2354,
 60    89 : 2702,
 61    95 : 3074,
 62    101: 3470,
 63    107: 3890,
 64    113: 4334,
 65    119: 4802,
 66    125: 5294,
 67    131: 5810
 68}
 69
 70#         Period    1   2   3   4   5   6   7         # level
 71Mapping = np.array([[11, 15, 17, 17, 17, 17, 17],     # 0        
 72                   [17, 23, 23, 23, 23, 23, 23],      # 1
 73                   [23, 29, 29, 29, 29, 29, 29],      # 2
 74                   [29, 29, 35, 35, 35, 35, 35],      # 3
 75                   [35, 41, 41, 41, 41, 41, 41],      # 4 This is the minimum that user can specify.
 76                   [41, 47, 47, 47, 47, 47, 47],      # 5
 77                   [47, 53, 53, 53, 53, 53, 53],      # 6
 78                   [53, 59, 59, 59, 59, 59, 59],      # 7
 79                   [59, 59, 59, 59, 59, 59, 59],      # 8  
 80                   [65, 65, 65, 65, 65, 65, 65]])     # 9  This is the maximum that user can specify.
 81
 82def max_ang(charge, level):
 83    #Mapping the LEBEDEV order with level of grids
 84    period = int(Data.elementPeriod[charge])
 85    index = Mapping[level+1,period]
 86    return lebedevOrdering[index]
 87
 88
 89def min_ang(charge, level):
 90    #Mapping the LEBEDEV order with level of grids
 91    period = int(Data.elementPeriod[charge])
 92    index = Mapping[level-3,period]
 93    return lebedevOrdering[index]
 94
 95
 96
 97def genGridiNewAPI(radial_precision,proton_charges,center_coordinates_bohr,basis_set_name,center_index,level,alpha_min, alpha_max):
 98    #This function generates the grid for a given atom (center_index)
 99    #This function is used to enable the generation of grids in parallel using joblib library
100    #The idea being, that grids are generated atom-by-atom,
101    #so it would be better to generate them parallely.
102    
103    # min_num_angular_points = 110#110#86#min_ang(proton_charges[center_index], level - 4)
104    # max_num_angular_points = 434#302#max_ang(proton_charges[center_index], level)
105    min_num_angular_points = min_ang(proton_charges[center_index], level)
106    max_num_angular_points = max_ang(proton_charges[center_index], level)
107    hardness = 3
108    # alpha_max = [
109    #                 11720.0,  # O
110    #                 13.01,  # H
111    #                 13.01,  # H
112    #             ]
113    # alpha_min = [
114    #                 {0: 0.3023, 1: 0.2753, 2: 1.185},  # O
115    #                 {0: 0.122, 1: 0.727},  # H
116    #                 {0: 0.122, 1: 0.727},  # H
117    #             ]
118
119    # atom grid using explicit basis set parameters
120    coordinates, weights = numgrid.atom_grid(
121                                    alpha_min[center_index],
122                                    alpha_max[center_index],
123                                    radial_precision,
124                                    min_num_angular_points,
125                                    max_num_angular_points,
126                                    proton_charges,
127                                    center_index,
128                                    center_coordinates_bohr,
129                                    hardness
130                                )
131    
132    # Doesn't seem to be working
133    # The problem is that the REST API request to basis set exchange results in a timeout
134    # coordinates, weights = numgrid.atom_grid_bse(basis_set_name,radial_precision,
135    #                                 min_num_angular_points,
136    #                                 max_num_angular_points,
137    #                                 proton_charges,center_index,
138    #                                 center_coordinates_bohr,
139    #                                 hardness
140    #                                 )
141
142
143    coordinates = np.array(coordinates)
144    x = coordinates[:,0]
145    y = coordinates[:,1]
146    z = coordinates[:,2]
147
148
149    return x,y,z,weights
150
151class Grids:
152    """
153    Class for generating molecular integration grids for DFT and other quantum chemistry calculations.
154
155    This class uses the `numgrid` library to generate atom-centered grids composed of:
156    1. A `level` indicating the fineness of the grid (from 3 to 8, with 0 reserved internally).
157    2. `coords`: a (N, 3) NumPy array of Cartesian coordinates (in Bohrs) for the N grid points.
158    3. `weights`: a NumPy array of length N containing the integration weights for each grid point.
159
160    Grid density depends not only on the level but also on the basis set and the radial precision,
161    making it more flexible and customizable than typical grid generation schemes.
162    """
163
164    #Class for the generation of molecular grids for numerical integration.
165    #The grid has three main components.
166    # 1. The grid 'level' refers to the coarsness or fineness of the grid. level = 4 (coarse) to 9 (fine)
167    # Although, internally CrysX would treat 0 as the minimum limit.
168    # This is to allow us to set a smaller min_num_angular_points and a larger max_num_angular_points.
169    # 2. The grid 'coords' (Bohrs) refer to a (N,3) numpy array with 3D cartesian coordinates 
170    # of N grid points.
171    # 3. The grid 'weights' that is an N-element numpy array with the weights corresponding to the grid points.
172
173    #In order to initialize a Grids object we need to provide the Mol object, Basis object and level (integer).
174    #One very important thing to note here is that, we use the 'numgrid' library to generate grids.
175    # https://github.com/dftlibs/numgrid
176    #Because of the way numgrid is implemented, there are mainly three ways to control the density of grids
177    # 1. Basis set: for radial grid
178    # 2. LEBEDEV ORDER or min/max angular points. We take care of this using the 'level' keyword.
179    # 3. Final the radial_precision.
180    # So, unlike other packages where the grid size would just depend on the level specified,
181    # here even the size of basis set also affect the grid size.
182
183
184    def __init__(self, mol=None, basis=None, level=3, radial_precision=1.0e-13, ncores = os.cpu_count()):
185        """
186        Initialize a Grids object for molecular numerical integration using atom-centered grids.
187
188        Parameters
189        ----------
190        mol : object
191            A Mol object containing atomic coordinates (`mol.coordsBohrs`) and nuclear charges (`mol.Zcharges`).
192            Required for placing atom-centered grid points.
193
194        basis : object or None
195            A Basis object specifying the basis set to use. If None, defaults to 'def2-QZVPP'.
196            Affects the radial distribution of the grid.
197
198        level : int, default=3
199            Grid resolution level. Valid values are 3 (coarse) to 8 (fine).
200            Internally mapped to a min/max angular point scheme used by the `numgrid` library.
201
202        radial_precision : float, default=1.0e-13
203            Precision used in radial grid generation. Lower values increase the number of radial points.
204
205        ncores : int, optional
206            Number of CPU cores to use for parallel grid generation. Defaults to the number of available cores.
207
208        Raises
209        ------
210        ValueError
211            If `mol` is None or `level` is outside the supported range [3, 8].
212
213        Notes
214        -----
215        Uses the `numgrid` library for generating atomic-centered numerical integration grids.
216        The final grid coordinates and weights are stored in the `coords` and `weights` attributes, respectively.
217        """
218        # For level user can enter an integer from 3 to 8.
219
220        # To get the same energies as PySCF (level=5) upto 1e-7 au, use the following settings
221        # radial_precision=1.0e-13
222        # level=3
223        # pruning by density with threshold = 1e-011
224        # alpha_min and alpha_max corresponding to QZVP
225
226
227        # As a default we can use the def2-TZVP basis set, in case no basis set is provided
228        # or a custom basis set is provided, in which case we won't have a specific name to search for in BSE database.
229
230        if basis is None:
231            basis_set_name = 'def2-QZVPP'
232        else:
233            basis_set_name = basis.basisSet.splitlines()[0]  #TODO FIXME But this would only work if the basis set was 
234                                                             #specified from CrysX lib using Basis.load()
235                                                             #If a user specified basis set was given via a string
236                                                             #or a file, then this would fail.
237        
238        if mol is None:
239            print("ERROR: Can't generate grids without molecular information!")
240            return
241
242        if level<3 or level>8:
243            print("ERROR: Enter a valid value for level between 3 to 8!")
244            return
245
246        self.level = level
247        """The chosen grid level (integer between 3 and 8). Controls the angular resolution and density of the integration grid."""
248
249
250        #Create the atomic coordinate arrays needed by numgrid
251        x_coordinates_bohr = []
252        y_coordinates_bohr = []
253        z_coordinates_bohr = []
254        center_coordinates_bohrs = []
255        for i in range (mol.natoms):
256            x_coordinates_bohr.append(mol.coordsBohrs[i][0])
257            y_coordinates_bohr.append(mol.coordsBohrs[i][1])
258            z_coordinates_bohr.append(mol.coordsBohrs[i][2])
259            center_coordinates_bohrs.append((mol.coordsBohrs[i][0],mol.coordsBohrs[i][1],mol.coordsBohrs[i][2]))
260
261        
262        #Get the required values
263        num_centers = mol.natoms
264        proton_charges = mol.Zcharges
265
266        #Now start the grid stuff
267
268        
269        #Create arrays to store the returned grid points and their weights
270        #These will later be turned into numpy arrays
271        x_grid = []
272        y_grid = []
273        z_grid = []
274        w_grid = []
275
276
277        print('Grid generation using '+str(ncores)+' threads for Joblib.')
278        output = Parallel(n_jobs=ncores, backend='threading', require='sharedmem')(delayed(genGridiNewAPI)(radial_precision,proton_charges,center_coordinates_bohrs,basis_set_name,center_index,level,basis.alpha_min, basis.alpha_max) for center_index in range(num_centers))
279        num_total_grid_points = 0
280        for center_index in range(num_centers):
281            num_total_grid_points = num_total_grid_points + len(output[center_index][0])
282            x_grid.extend(output[center_index][0])
283            y_grid.extend(output[center_index][1])
284            z_grid.extend(output[center_index][2])
285            w_grid.extend(output[center_index][3])
286
287
288        #print(num_total_grid_points)
289        #Now create the numpy arrays that store the grid points and weights
290        #self.coords = np.empty((len(num_total_grid_points),3))
291        #self.weights = np.empty((len(num_total_grid_points)))
292        self.coords = np.empty((len(x_grid),3))
293        """A NumPy array of shape (N, 3), storing the Cartesian coordinates of N grid points (in Bohrs)."""
294        self.weights = np.empty((len(x_grid)))
295        """A NumPy array of length N, storing the quadrature weights corresponding to each grid point."""
296        #Convert the Python arrays x_grid, y_grid, z_grid to the Grids.coords numpy arrays
297        self.coords[:,0] = np.array(x_grid)
298        self.coords[:,1] = np.array(y_grid)
299        self.coords[:,2] = np.array(z_grid)
300        # center_coordinates_bohrs = np.array(center_coordinates_bohrs)
301        # self.coords = self.coords[(np.abs(np.subtract.outer(self.coords,center_coordinates_bohrs)) > 1e-5).all(1)]
302        self.weights[:] = np.array(w_grid)
303        return
304
305        
306
307
308
309