1 ////////////////////////////////////////////////////////////////////////////////////// 2 // This file is distributed under the University of Illinois/NCSA Open Source License. 3 // See LICENSE file in top directory for details. 4 // 5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. 6 // 7 // File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 8 // Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign 9 // Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory 10 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign 11 // Ye Luo, yeluo@anl.gov, Argonne National Laboratory 12 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory 13 // 14 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 15 ////////////////////////////////////////////////////////////////////////////////////// 16 17 18 #ifndef ATOMIC_ORBITAL_H 19 #define ATOMIC_ORBITAL_H 20 21 #include "config/stdlib/math.hpp" 22 #include "einspline/multi_bspline.h" 23 #include "QMCWaveFunctions/SPOSet.h" 24 #include "Lattice/CrystalLattice.h" 25 #include <Configuration.h> 26 #include "Utilities/TimerManager.h" 27 28 29 namespace qmcplusplus 30 { 31 /****************************************************************** 32 // This is just a template trick to avoid template specialization // 33 // in AtomicOrbital. // 34 ******************************************************************/ 35 36 template<typename StorageType> 37 struct AtomicOrbitalTraits 38 {}; 39 template<> 40 struct AtomicOrbitalTraits<double> 41 { 42 typedef multi_UBspline_1d_d SplineType; 43 }; 44 template<> 45 struct AtomicOrbitalTraits<std::complex<double>> 46 { 47 typedef multi_UBspline_1d_z SplineType; 48 }; 49 50 inline void EinsplineMultiEval(multi_UBspline_1d_d* spline, double x, double* val) 51 { 52 eval_multi_UBspline_1d_d(spline, x, val); 53 } 54 inline void EinsplineMultiEval(multi_UBspline_1d_z* spline, double x, std::complex<double>* val) 55 { 56 eval_multi_UBspline_1d_z(spline, x, val); 57 } 58 inline void EinsplineMultiEval(multi_UBspline_1d_d* spline, double x, double* val, double* grad, double* lapl) 59 { 60 eval_multi_UBspline_1d_d_vgl(spline, x, val, grad, lapl); 61 } 62 inline void EinsplineMultiEval(multi_UBspline_1d_z* spline, 63 double x, 64 std::complex<double>* val, 65 std::complex<double>* grad, 66 std::complex<double>* lapl) 67 { 68 eval_multi_UBspline_1d_z_vgl(spline, x, val, grad, lapl); 69 } 70 71 72 template<typename StorageType> 73 class AtomicOrbital 74 { 75 public: 76 typedef QMCTraits::PosType PosType; 77 typedef QMCTraits::RealType RealType; 78 typedef CrystalLattice<RealType, OHMMS_DIM> UnitCellType; 79 typedef Vector<double> RealValueVector_t; 80 typedef Vector<TinyVector<double, OHMMS_DIM>> RealGradVector_t; 81 typedef Vector<std::complex<double>> ComplexValueVector_t; 82 typedef Vector<TinyVector<std::complex<double>, OHMMS_DIM>> ComplexGradVector_t; 83 typedef Vector<Tensor<double, OHMMS_DIM>> RealHessVector_t; 84 typedef Vector<Tensor<std::complex<double>, OHMMS_DIM>> ComplexHessVector_t; 85 typedef typename AtomicOrbitalTraits<StorageType>::SplineType SplineType; 86 87 private: 88 // Store in order 89 // Index = l*(l+1) + m. There are (lMax+1)^2 Ylm's 90 std::vector<StorageType> YlmVec, dYlm_dthetaVec, dYlm_dphiVec, ulmVec, dulmVec, d2ulmVec; 91 92 SplineType* RadialSpline; 93 // The first index is n in r^n, the second is lm = l*(l+1)+m 94 Array<StorageType, 3> PolyCoefs; 95 NewTimer &YlmTimer, &SplineTimer, &SumTimer; 96 RealType rmagLast; 97 std::vector<PosType> TwistAngles; 98 99 public: 100 PosType Pos; 101 RealType CutoffRadius, SplineRadius, PolyRadius; 102 int SplinePoints; 103 int PolyOrder; 104 int lMax, Numlm, NumBands; 105 UnitCellType Lattice; 106 107 inline void set_pos(PosType pos) { Pos = pos; } 108 inline void set_lmax(int lmax) { lMax = lmax; } 109 inline void set_cutoff(RealType cutoff) { CutoffRadius = cutoff; } 110 inline void set_spline(RealType radius, int points) 111 { 112 SplineRadius = radius; 113 SplinePoints = points; 114 } 115 inline void set_polynomial(RealType radius, int order) 116 { 117 PolyRadius = radius; 118 PolyOrder = order; 119 } 120 inline void set_num_bands(int num_bands) { NumBands = num_bands; } 121 SplineType* get_radial_spline() { return RadialSpline; } 122 Array<StorageType, 3>& get_poly_coefs() { return PolyCoefs; } 123 124 inline void registerTimers() 125 { 126 YlmTimer.reset(); 127 SplineTimer.reset(); 128 } 129 130 void allocate(); 131 132 void set_band(int band, 133 Array<std::complex<double>, 2>& spline_data, 134 Array<std::complex<double>, 2>& poly_coefs, 135 PosType twist); 136 inline void CalcYlm(PosType rhat, 137 std::vector<std::complex<double>>& Ylm, 138 std::vector<std::complex<double>>& dYlm_dtheta, 139 std::vector<std::complex<double>>& dYlm_dphi); 140 141 inline void CalcYlm(PosType rhat, 142 std::vector<double>& Ylm, 143 std::vector<double>& dYlm_dtheta, 144 std::vector<double>& dYlm_dphi); 145 146 inline bool evaluate(PosType r, ComplexValueVector_t& vals); 147 inline bool evaluate(PosType r, ComplexValueVector_t& val, ComplexGradVector_t& grad, ComplexValueVector_t& lapl); 148 inline bool evaluate(PosType r, ComplexValueVector_t& val, ComplexGradVector_t& grad, ComplexHessVector_t& lapl); 149 inline bool evaluate(PosType r, RealValueVector_t& vals); 150 inline bool evaluate(PosType r, RealValueVector_t& val, RealGradVector_t& grad, RealValueVector_t& lapl); 151 inline bool evaluate(PosType r, RealValueVector_t& val, RealGradVector_t& grad, RealHessVector_t& lapl); 152 153 154 AtomicOrbital() 155 : RadialSpline(NULL), 156 YlmTimer(*timer_manager.createTimer("AtomicOrbital::CalcYlm")), 157 SplineTimer(*timer_manager.createTimer("AtomicOrbital::1D spline")), 158 SumTimer(*timer_manager.createTimer("AtomicOrbital::Summation")), 159 rmagLast(std::numeric_limits<RealType>::max()) 160 { 161 // Nothing else for now 162 } 163 }; 164 165 166 template<typename StorageType> 167 inline bool AtomicOrbital<StorageType>::evaluate(PosType r, ComplexValueVector_t& vals) 168 { 169 PosType dr = r - Pos; 170 PosType u = Lattice.toUnit(dr); 171 PosType img; 172 for (int i = 0; i < OHMMS_DIM; i++) 173 { 174 img[i] = round(u[i]); 175 u[i] -= img[i]; 176 } 177 dr = Lattice.toCart(u); 178 double r2 = dot(dr, dr); 179 if (r2 > CutoffRadius * CutoffRadius) 180 return false; 181 double rmag = std::sqrt(r2); 182 PosType rhat = (1.0 / rmag) * dr; 183 // Evaluate Ylm's 184 CalcYlm(rhat, YlmVec, dYlm_dthetaVec, dYlm_dphiVec); 185 if (std::abs(rmag - rmagLast) > 1.0e-6) 186 { 187 // Evaluate radial functions 188 if (rmag > PolyRadius) 189 EinsplineMultiEval(RadialSpline, rmag, &(ulmVec[0])); 190 else 191 { 192 for (int index = 0; index < ulmVec.size(); index++) 193 ulmVec[index] = StorageType(); 194 double r2n = 1.0; 195 for (int n = 0; n <= PolyOrder; n++) 196 { 197 int index = 0; 198 for (int i = 0; i < vals.size(); i++) 199 for (int lm = 0; lm < Numlm; lm++) 200 ulmVec[index++] += r2n * PolyCoefs(n, i, lm); 201 r2n *= rmag; 202 } 203 } 204 rmagLast = rmag; 205 } 206 SumTimer.start(); 207 int index = 0; 208 for (int i = 0; i < vals.size(); i++) 209 { 210 vals[i] = std::complex<double>(); 211 for (int lm = 0; lm < Numlm; lm++) 212 vals[i] += ulmVec[index++] * YlmVec[lm]; 213 double phase = -2.0 * M_PI * dot(TwistAngles[i], img); 214 // fprintf (stderr, "phase[%d] = %1.2f pi\n", i, phase/M_PI); 215 // fprintf (stderr, "img = [%f,%f,%f]\n", img[0], img[1], img[2]); 216 double s, c; 217 qmcplusplus::sincos(phase, &s, &c); 218 vals[i] *= std::complex<double>(c, s); 219 } 220 SumTimer.stop(); 221 return true; 222 } 223 224 225 template<typename StorageType> 226 inline bool AtomicOrbital<StorageType>::evaluate(PosType r, RealValueVector_t& vals) 227 { 228 PosType dr = r - Pos; 229 PosType u = Lattice.toUnit(dr); 230 PosType img; 231 for (int i = 0; i < OHMMS_DIM; i++) 232 { 233 img[i] = round(u[i]); 234 u[i] -= img[i]; 235 } 236 dr = Lattice.toCart(u); 237 double r2 = dot(dr, dr); 238 if (r2 > CutoffRadius * CutoffRadius) 239 return false; 240 double rmag = std::sqrt(r2); 241 PosType rhat = (1.0 / rmag) * dr; 242 // Evaluate Ylm's 243 CalcYlm(rhat, YlmVec, dYlm_dthetaVec, dYlm_dphiVec); 244 if (std::abs(rmag - rmagLast) > 1.0e-6) 245 { 246 // Evaluate radial functions 247 if (rmag > PolyRadius) 248 { 249 SplineTimer.start(); 250 EinsplineMultiEval(RadialSpline, rmag, &(ulmVec[0])); 251 SplineTimer.stop(); 252 } 253 else 254 { 255 for (int index = 0; index < ulmVec.size(); index++) 256 ulmVec[index] = StorageType(); 257 double r2n = 1.0; 258 for (int n = 0; n <= PolyOrder; n++) 259 { 260 int index = 0; 261 for (int i = 0; i < vals.size(); i++) 262 for (int lm = 0; lm < Numlm; lm++) 263 ulmVec[index++] += r2n * PolyCoefs(n, i, lm); 264 r2n *= rmag; 265 } 266 } 267 rmagLast = rmag; 268 } 269 SumTimer.start(); 270 int index = 0; 271 for (int i = 0; i < vals.size(); i++) 272 { 273 vals[i] = 0.0; 274 StorageType tmp = 0.0; 275 for (int lm = 0; lm < Numlm; lm++, index++) 276 tmp += ulmVec[index] * YlmVec[lm]; 277 //vals[i] += real(ulmVec[index++] * YlmVec[lm]); 278 // vals[i] += (ulmVec[index].real() * YlmVec[lm].real() - 279 // ulmVec[index].imag() * YlmVec[lm].imag()); 280 double phase = -2.0 * M_PI * dot(TwistAngles[i], img); 281 double s, c; 282 qmcplusplus::sincos(phase, &s, &c); 283 vals[i] = real(std::complex<double>(c, s) * tmp); 284 } 285 SumTimer.stop(); 286 return true; 287 } 288 289 template<typename StorageType> 290 inline bool AtomicOrbital<StorageType>::evaluate(PosType r, 291 RealValueVector_t& vals, 292 RealGradVector_t& grads, 293 RealHessVector_t& hess) 294 { 295 APP_ABORT(" AtomicOrbital<StorageType>::evaluate not implemented for Hess. \n"); 296 return true; 297 } 298 299 300 template<typename StorageType> 301 inline bool AtomicOrbital<StorageType>::evaluate(PosType r, 302 RealValueVector_t& vals, 303 RealGradVector_t& grads, 304 RealValueVector_t& lapl) 305 { 306 PosType dr = r - Pos; 307 PosType u = Lattice.toUnit(dr); 308 PosType img; 309 for (int i = 0; i < OHMMS_DIM; i++) 310 { 311 img[i] = round(u[i]); 312 u[i] -= img[i]; 313 } 314 dr = Lattice.toCart(u); 315 double r2 = dot(dr, dr); 316 if (r2 > CutoffRadius * CutoffRadius) 317 return false; 318 double rmag = std::sqrt(r2); 319 double rInv = 1.0 / rmag; 320 PosType rhat = rInv * dr; 321 double costheta = rhat[2]; 322 double sintheta = std::sqrt(1.0 - costheta * costheta); 323 double cosphi = rhat[0] / sintheta; 324 double sinphi = rhat[1] / sintheta; 325 PosType thetahat = PosType(costheta * cosphi, costheta * sinphi, -sintheta); 326 PosType phihat = PosType(-sinphi, cosphi, 0.0); 327 // Evaluate Ylm's 328 CalcYlm(rhat, YlmVec, dYlm_dthetaVec, dYlm_dphiVec); 329 // Evaluate radial functions 330 if (rmag > PolyRadius) 331 { 332 SplineTimer.start(); 333 EinsplineMultiEval(RadialSpline, rmag, &(ulmVec[0]), &(dulmVec[0]), &(d2ulmVec[0])); 334 SplineTimer.stop(); 335 } 336 else 337 { 338 for (int index = 0; index < ulmVec.size(); index++) 339 { 340 ulmVec[index] = StorageType(); 341 dulmVec[index] = StorageType(); 342 d2ulmVec[index] = StorageType(); 343 } 344 double r2n = 1.0, r2nm1 = 0.0, r2nm2 = 0.0; 345 double dn = 0.0; 346 double dnm1 = -1.0; 347 for (int n = 0; n <= PolyOrder; n++) 348 { 349 int index = 0; 350 for (int i = 0; i < vals.size(); i++) 351 for (int lm = 0; lm < Numlm; lm++, index++) 352 { 353 StorageType c = PolyCoefs(n, i, lm); 354 ulmVec[index] += r2n * c; 355 dulmVec[index] += dn * r2nm1 * c; 356 d2ulmVec[index] += dn * dnm1 * r2nm2 * c; 357 } 358 dn += 1.0; 359 dnm1 += 1.0; 360 r2nm2 = r2nm1; 361 r2nm1 = r2n; 362 r2n *= rmag; 363 } 364 } 365 SumTimer.start(); 366 int index = 0; 367 for (int i = 0; i < vals.size(); i++) 368 { 369 vals[i] = 0.0; 370 for (int j = 0; j < OHMMS_DIM; j++) 371 grads[i][j] = 0.0; 372 lapl[i] = 0.0; 373 // Compute e^{-i k.L} phase factor 374 double phase = -2.0 * M_PI * dot(TwistAngles[i], img); 375 double s, c; 376 qmcplusplus::sincos(phase, &s, &c); 377 std::complex<double> e2mikr(c, s); 378 StorageType tmp_val, tmp_lapl, grad_rhat, grad_thetahat, grad_phihat; 379 tmp_val = tmp_lapl = grad_rhat = grad_thetahat = grad_phihat = StorageType(); 380 int lm = 0; 381 for (int l = 0; l <= lMax; l++) 382 for (int m = -l; m <= l; m++, lm++, index++) 383 { 384 std::complex<double> im(0.0, (double)m); 385 tmp_val += ulmVec[index] * YlmVec[lm]; 386 grad_rhat += dulmVec[index] * YlmVec[lm]; 387 grad_thetahat += ulmVec[index] * rInv * dYlm_dthetaVec[lm]; 388 grad_phihat += (ulmVec[index] * dYlm_dphiVec[lm]) / (rmag * sintheta); 389 //grad_phihat += (ulmVec[index] * im *YlmVec[lm])/(rmag*sintheta); 390 tmp_lapl += YlmVec[lm] * 391 (-(double)(l * (l + 1)) * rInv * rInv * ulmVec[index] + d2ulmVec[index] + 2.0 * rInv * dulmVec[index]); 392 } 393 vals[i] = real(e2mikr * tmp_val); 394 lapl[i] = real(e2mikr * tmp_lapl); 395 grads[i] = (real(e2mikr * grad_rhat) * rhat + real(e2mikr * grad_thetahat) * thetahat + 396 real(e2mikr * grad_phihat) * phihat); 397 } 398 SumTimer.stop(); 399 rmagLast = rmag; 400 return true; 401 } 402 403 template<typename StorageType> 404 inline bool AtomicOrbital<StorageType>::evaluate(PosType r, 405 ComplexValueVector_t& vals, 406 ComplexGradVector_t& grads, 407 ComplexHessVector_t& hess) 408 { 409 APP_ABORT(" AtomicOrbital<StorageType>::evaluate not implemented for Hess. \n"); 410 return true; 411 } 412 413 template<typename StorageType> 414 inline bool AtomicOrbital<StorageType>::evaluate(PosType r, 415 ComplexValueVector_t& vals, 416 ComplexGradVector_t& grads, 417 ComplexValueVector_t& lapl) 418 { 419 PosType dr = r - Pos; 420 PosType u = Lattice.toUnit(dr); 421 PosType img; 422 for (int i = 0; i < OHMMS_DIM; i++) 423 { 424 img[i] = round(u[i]); 425 u[i] -= img[i]; 426 } 427 dr = Lattice.toCart(u); 428 double r2 = dot(dr, dr); 429 if (r2 > CutoffRadius * CutoffRadius) 430 return false; 431 double rmag = std::sqrt(r2); 432 double rInv = 1.0 / rmag; 433 PosType rhat = rInv * dr; 434 double costheta = rhat[2]; 435 double sintheta = std::sqrt(1.0 - costheta * costheta); 436 double cosphi = rhat[0] / sintheta; 437 double sinphi = rhat[1] / sintheta; 438 PosType thetahat = PosType(costheta * cosphi, costheta * sinphi, -sintheta); 439 PosType phihat = PosType(-sinphi, cosphi, 0.0); 440 // Evaluate Ylm's 441 CalcYlm(rhat, YlmVec, dYlm_dthetaVec, dYlm_dphiVec); 442 // Evaluate radial functions 443 if (rmag > PolyRadius) 444 { 445 SplineTimer.start(); 446 EinsplineMultiEval(RadialSpline, rmag, &(ulmVec[0]), &(dulmVec[0]), &(d2ulmVec[0])); 447 SplineTimer.stop(); 448 } 449 else 450 { 451 for (int index = 0; index < ulmVec.size(); index++) 452 { 453 ulmVec[index] = StorageType(); 454 dulmVec[index] = StorageType(); 455 d2ulmVec[index] = StorageType(); 456 } 457 double r2n = 1.0, r2nm1 = 0.0, r2nm2 = 0.0; 458 double dn = 0.0; 459 double dnm1 = -1.0; 460 for (int n = 0; n <= PolyOrder; n++) 461 { 462 int index = 0; 463 for (int i = 0; i < vals.size(); i++) 464 for (int lm = 0; lm < Numlm; lm++, index++) 465 { 466 StorageType c = PolyCoefs(n, i, lm); 467 ulmVec[index] += r2n * c; 468 dulmVec[index] += dn * r2nm1 * c; 469 d2ulmVec[index] += dn * dnm1 * r2nm2 * c; 470 } 471 dn += 1.0; 472 dnm1 += 1.0; 473 r2nm2 = r2nm1; 474 r2nm1 = r2n; 475 r2n *= rmag; 476 } 477 } 478 SumTimer.start(); 479 int index = 0; 480 for (int i = 0; i < vals.size(); i++) 481 { 482 vals[i] = 0.0; 483 for (int j = 0; j < OHMMS_DIM; j++) 484 grads[i][j] = 0.0; 485 lapl[i] = 0.0; 486 int lm = 0; 487 StorageType grad_rhat, grad_thetahat, grad_phihat; 488 // Compute e^{-i k.L} phase factor 489 double phase = -2.0 * M_PI * dot(TwistAngles[i], img); 490 double s, c; 491 qmcplusplus::sincos(phase, &s, &c); 492 std::complex<double> e2mikr(c, s); 493 for (int l = 0; l <= lMax; l++) 494 for (int m = -l; m <= l; m++, lm++, index++) 495 { 496 std::complex<double> im(0.0, (double)m); 497 vals[i] += ulmVec[index] * YlmVec[lm]; 498 grad_rhat += dulmVec[index] * YlmVec[lm]; 499 grad_thetahat += ulmVec[index] * rInv * dYlm_dthetaVec[lm]; 500 grad_phihat += (ulmVec[index] * im * YlmVec[lm]) / (rmag * sintheta); 501 lapl[i] += YlmVec[lm] * 502 (-(double)(l * (l + 1)) * rInv * rInv * ulmVec[index] + d2ulmVec[index] + 2.0 * rInv * dulmVec[index]); 503 } 504 vals[i] *= e2mikr; 505 lapl[i] *= e2mikr; 506 for (int j = 0; j < OHMMS_DIM; j++) 507 { 508 grads[i][j] = e2mikr * (grad_rhat * rhat[j] + grad_thetahat * thetahat[j] + grad_phihat * phihat[j]); 509 } 510 } 511 SumTimer.stop(); 512 rmagLast = rmag; 513 return true; 514 } 515 516 517 // Fast implementation 518 // See Geophys. J. Int. (1998) 135,pp.307-309 519 template<typename StorageType> 520 inline void AtomicOrbital<StorageType>::CalcYlm(PosType rhat, 521 std::vector<std::complex<double>>& Ylm, 522 std::vector<std::complex<double>>& dYlm_dtheta, 523 std::vector<std::complex<double>>& dYlm_dphi) 524 { 525 YlmTimer.start(); 526 const double fourPiInv = 0.0795774715459477; 527 double costheta = rhat[2]; 528 double sintheta = std::sqrt(1.0 - costheta * costheta); 529 double cottheta = costheta / sintheta; 530 double cosphi, sinphi; 531 cosphi = rhat[0] / sintheta; 532 sinphi = rhat[1] / sintheta; 533 std::complex<double> e2iphi(cosphi, sinphi); 534 double lsign = 1.0; 535 double dl = 0.0; 536 std::vector<double> XlmVec(2 * lMax + 1), dXlmVec(2 * lMax + 1); 537 for (int l = 0; l <= lMax; l++) 538 { 539 XlmVec[2 * l] = lsign; 540 dXlmVec[2 * l] = dl * cottheta * XlmVec[2 * l]; 541 XlmVec[0] = lsign * XlmVec[2 * l]; 542 dXlmVec[0] = lsign * dXlmVec[2 * l]; 543 double dm = dl; 544 double msign = lsign; 545 for (int m = l; m > 0; m--) 546 { 547 double tmp = std::sqrt((dl + dm) * (dl - dm + 1.0)); 548 XlmVec[l + m - 1] = -(dXlmVec[l + m] + dm * cottheta * XlmVec[l + m]) / tmp; 549 dXlmVec[l + m - 1] = (dm - 1.0) * cottheta * XlmVec[l + m - 1] + XlmVec[l + m] * tmp; 550 // Copy to negative m 551 XlmVec[l - (m - 1)] = -msign * XlmVec[l + m - 1]; 552 dXlmVec[l - (m - 1)] = -msign * dXlmVec[l + m - 1]; 553 msign *= -1.0; 554 dm -= 1.0; 555 } 556 double sum = 0.0; 557 for (int m = -l; m <= l; m++) 558 sum += XlmVec[l + m] * XlmVec[l + m]; 559 // Now, renormalize the Ylms for this l 560 double norm = std::sqrt((2.0 * dl + 1.0) * fourPiInv / sum); 561 for (int m = -l; m <= l; m++) 562 { 563 XlmVec[l + m] *= norm; 564 dXlmVec[l + m] *= norm; 565 } 566 // Multiply by azimuthal phase and store in Ylm 567 std::complex<double> e2imphi(1.0, 0.0); 568 std::complex<double> eye(0.0, 1.0); 569 for (int m = 0; m <= l; m++) 570 { 571 Ylm[l * (l + 1) + m] = XlmVec[l + m] * e2imphi; 572 Ylm[l * (l + 1) - m] = XlmVec[l - m] * qmcplusplus::conj(e2imphi); 573 dYlm_dphi[l * (l + 1) + m] = (double)m * eye * XlmVec[l + m] * e2imphi; 574 dYlm_dphi[l * (l + 1) - m] = -(double)m * eye * XlmVec[l - m] * qmcplusplus::conj(e2imphi); 575 dYlm_dtheta[l * (l + 1) + m] = dXlmVec[l + m] * e2imphi; 576 dYlm_dtheta[l * (l + 1) - m] = dXlmVec[l - m] * qmcplusplus::conj(e2imphi); 577 e2imphi *= e2iphi; 578 } 579 dl += 1.0; 580 lsign *= -1.0; 581 } 582 YlmTimer.stop(); 583 } 584 585 // Fast implementation 586 // See Geophys. J. Int. (1998) 135,pp.307-309 587 template<typename StorageType> 588 inline void AtomicOrbital<StorageType>::CalcYlm(PosType rhat, 589 std::vector<double>& Ylm, 590 std::vector<double>& dYlm_dtheta, 591 std::vector<double>& dYlm_dphi) 592 { 593 YlmTimer.start(); 594 const double fourPiInv = 0.0795774715459477; 595 double costheta = rhat[2]; 596 double sintheta = std::sqrt(1.0 - costheta * costheta); 597 double cottheta = costheta / sintheta; 598 double cosphi, sinphi; 599 cosphi = rhat[0] / sintheta; 600 sinphi = rhat[1] / sintheta; 601 std::complex<double> e2iphi(cosphi, sinphi); 602 double lsign = 1.0; 603 double dl = 0.0; 604 std::vector<double> XlmVec(2 * lMax + 1), dXlmVec(2 * lMax + 1); 605 for (int l = 0; l <= lMax; l++) 606 { 607 XlmVec[2 * l] = lsign; 608 dXlmVec[2 * l] = dl * cottheta * XlmVec[2 * l]; 609 XlmVec[0] = lsign * XlmVec[2 * l]; 610 dXlmVec[0] = lsign * dXlmVec[2 * l]; 611 double dm = dl; 612 double msign = lsign; 613 for (int m = l; m > 0; m--) 614 { 615 double tmp = std::sqrt((dl + dm) * (dl - dm + 1.0)); 616 XlmVec[l + m - 1] = -(dXlmVec[l + m] + dm * cottheta * XlmVec[l + m]) / tmp; 617 dXlmVec[l + m - 1] = (dm - 1.0) * cottheta * XlmVec[l + m - 1] + XlmVec[l + m] * tmp; 618 // Copy to negative m 619 XlmVec[l - (m - 1)] = -msign * XlmVec[l + m - 1]; 620 dXlmVec[l - (m - 1)] = -msign * dXlmVec[l + m - 1]; 621 msign *= -1.0; 622 dm -= 1.0; 623 } 624 double sum = 0.0; 625 for (int m = -l; m <= l; m++) 626 sum += XlmVec[l + m] * XlmVec[l + m]; 627 // Now, renormalize the Ylms for this l 628 double norm = std::sqrt((2.0 * dl + 1.0) * fourPiInv / sum); 629 for (int m = -l; m <= l; m++) 630 { 631 XlmVec[l + m] *= norm; 632 dXlmVec[l + m] *= norm; 633 } 634 // Multiply by azimuthal phase and store in Ylm 635 Ylm[l * (l + 1)] = XlmVec[l]; 636 dYlm_dphi[l * (l + 1)] = 0.0; 637 dYlm_dtheta[l * (l + 1)] = dXlmVec[l]; 638 std::complex<double> e2imphi = e2iphi; 639 for (int m = 1; m <= l; m++) 640 { 641 Ylm[l * (l + 1) + m] = XlmVec[l + m] * e2imphi.real(); 642 Ylm[l * (l + 1) - m] = XlmVec[l + m] * e2imphi.imag(); 643 dYlm_dphi[l * (l + 1) + m] = -(double)m * XlmVec[l + m] * e2imphi.imag(); 644 dYlm_dphi[l * (l + 1) - m] = (double)m * XlmVec[l + m] * e2imphi.real(); 645 dYlm_dtheta[l * (l + 1) + m] = dXlmVec[l + m] * e2imphi.real(); 646 dYlm_dtheta[l * (l + 1) - m] = dXlmVec[l + m] * e2imphi.imag(); 647 e2imphi *= e2iphi; 648 } 649 dl += 1.0; 650 lsign *= -1.0; 651 } 652 YlmTimer.stop(); 653 } 654 655 656 } // namespace qmcplusplus 657 #endif 658