1 /*
2  *                This source code is part of
3  *
4  *                     E  R  K  A  L  E
5  *                             -
6  *                       DFT from Hel
7  *
8  * Written by Susi Lehtola, 2010-2011
9  * Copyright (c) 2010-2011, Susi Lehtola
10  *
11  * This program is free software; you can redistribute it and/or
12  * modify it under the terms of the GNU General Public License
13  * as published by the Free Software Foundation; either version 2
14  * of the License, or (at your option) any later version.
15  */
16 
17 
18 
19 #ifndef ERKALE_LEBEDEV
20 #define ERKALE_LEBEDEV
21 
22 #include <vector>
23 
24 /// Orders of supported rules (integrates spherical harmonics exactly up to given order)
25 const int lebedev_orders[] ={3,  5,  7,  9, 11, 13, 15,  17,  19,  21,  23,  25,  27,  29,  31,  35,  41,  47,  53,   59,   65,   71,   77,   83,   89,   95,  101,  107,  113,  119,  125,  131};
26 /// The number of points used by the supported rules
27 const int lebedev_degrees[]={6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810};
28 
29 /// Coordinates and weight for a point in Lebedev quadrature
30 typedef struct {
31   /// x coordinate of point
32   double x;
33   /// y coordinate of point
34   double y;
35   /// z coordinate of point
36   double z;
37   /// Angular weight of point
38   double w;
39 } lebedev_point_t;
40 
41 /// Get a Lebedev sphere of the wanted order.
42 std::vector<lebedev_point_t> lebedev_sphere(int order);
43 /// Determine the next order of quadrature, if one is supported.
44 int next_lebedev(int order);
45 /// Get closest rule
46 int closest_lebedev(int order);
47 
48 /// Worker routine - get a Lebedev sphere with the wanted amount of points
49 std::vector<lebedev_point_t> getLebedevSphere(int npoints);
50 /// Worker routine - add points to the sphere structure
51 void getLebedevReccurencePoints(int type, int & start, double a, double b, double v, std::vector<lebedev_point_t> & leb);
52 
53 #endif
54