1 /* 2 * This source code is part of 3 * 4 * E R K A L E 5 * - 6 * HF/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 #ifndef ERKALE_LMGRID 19 #define ERKALE_LMGRID 20 21 #include <complex> 22 #include <vector> 23 #include "basis.h" 24 25 /// Index of (l,m) in tables: l^2 + l + m 26 #define lmind(l,m) ( ((size_t) (l))*(size_t (l)) + (size_t) (l) + (size_t) (m)) 27 28 /// Structure for radial integration 29 typedef struct { 30 /// Radius 31 double r; 32 /// Radial weight 33 double w; 34 } radial_grid_t; 35 36 /// Form radial grid 37 std::vector<radial_grid_t> form_radial_grid(int nrad); 38 39 /// Structure for angular integration 40 typedef struct { 41 /// Coordinate on the (unit) sphere 42 coords_t r; 43 /// Weight 44 double w; 45 } angular_grid_t; 46 47 /// Form angular grid 48 std::vector<angular_grid_t> form_angular_grid(int lmax); 49 /// Compute values of spherical harmonics Ylm[ngrid][l,m] 50 std::vector< std::vector< std::complex<double> > > compute_spherical_harmonics(const std::vector<angular_grid_t> & grid, int lmax); 51 52 /// Expansion of orbitals 53 typedef struct { 54 /// Radial grid 55 std::vector<radial_grid_t> grid; 56 /// Expansion coefficients of orbitals clm[norbs][l,m][nrad] 57 std::vector< std::vector< std::vector< std::complex<double> > > > clm; 58 } expansion_t; 59 60 /// Compute expansion of orbitals around cen, return clm[norbs][l,m][nrad] up to lmax. Quadrature uses lquad:th order 61 expansion_t expand_orbitals(const arma::mat & C, const BasisSet & bas, const coords_t & cen, bool verbose=true, size_t Nrad=200, int lmax=5, int lquad=30); 62 63 /// Expansion of orbitals 64 typedef struct { 65 /// Radial grid 66 std::vector<radial_grid_t> grid; 67 /// Expansion coefficients of orbitals clm[norbs][l,m][nrad] 68 std::vector< std::vector< std::vector< double > > > clm; 69 } real_expansion_t; 70 71 /// Compute expansion of orbitals around cen, return clm[norbs][l,m][nrad] up to lmax. Quadrature uses lquad:th order 72 real_expansion_t expand_orbitals_real(const arma::mat & C, const BasisSet & bas, const coords_t & cen, bool verbose=true, size_t Nrad=200, int lmax=5, int lquad=30); 73 74 /// Compute weight decomposition. If total, last column contains total norm of orbital 75 arma::mat weight_decomposition(const real_expansion_t & exp, bool total=true); 76 /// Compute weight decomposition. If total, last column contains total norm of orbital 77 arma::mat weight_decomposition(const expansion_t & exp, bool total=true); 78 79 /// Compute angular decomposition - weights for (l,m) 80 arma::mat angular_decomposition(const real_expansion_t & exp, int l); 81 82 #endif 83