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) 2019 QMCPACK developers. 6 // 7 // File developed by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab 8 // ChangMo Yang, nichthierwohne@gmail.com, University of Illinois at Urbana-Champaign 9 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 10 // Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory 11 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign 12 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory 13 // 14 // File created by: ChangMo Yang, nichthierwohne@gmail.com, University of Illinois at Urbana-Champaign 15 ////////////////////////////////////////////////////////////////////////////////////// 16 17 18 /**@file 19 */ 20 #ifndef QMCPLUSPLUS_MULTIDIRACDETERMINANTCALCULATOR_H 21 #define QMCPLUSPLUS_MULTIDIRACDETERMINANTCALCULATOR_H 22 23 #include "OhmmsPETE/OhmmsMatrix.h" 24 25 namespace qmcplusplus 26 { 27 /** Function class calculates multi determinant ratio from matrix elements 28 * needs to be the size of your result matrix or larger 29 * includes manual expansions for smaller evaluations 30 */ 31 template<typename T> 32 struct MultiDiracDeterminantCalculator 33 { 34 std::vector<T> M; 35 std::vector<int> Pivot; 36 resizeMultiDiracDeterminantCalculator37 void resize(int n) 38 { 39 M.resize(n * n); 40 Pivot.resize(n); 41 } 42 evaluateMultiDiracDeterminantCalculator43 inline T evaluate(T a11, T a12, T a21, T a22) { return a11 * a22 - a21 * a12; } 44 evaluateMultiDiracDeterminantCalculator45 inline T evaluate(T a11, T a12, T a13, T a21, T a22, T a23, T a31, T a32, T a33) 46 { 47 return (a11 * (a22 * a33 - a32 * a23) - a21 * (a12 * a33 - a32 * a13) + a31 * (a12 * a23 - a22 * a13)); 48 } 49 evaluateMultiDiracDeterminantCalculator50 inline T evaluate(T a11, 51 T a12, 52 T a13, 53 T a14, 54 T a21, 55 T a22, 56 T a23, 57 T a24, 58 T a31, 59 T a32, 60 T a33, 61 T a34, 62 T a41, 63 T a42, 64 T a43, 65 T a44) 66 { 67 return (a11 * (a22 * (a33 * a44 - a43 * a34) - a32 * (a23 * a44 - a43 * a24) + a42 * (a23 * a34 - a33 * a24)) - 68 a21 * (a12 * (a33 * a44 - a43 * a34) - a32 * (a13 * a44 - a43 * a14) + a42 * (a13 * a34 - a33 * a14)) + 69 a31 * (a12 * (a23 * a44 - a43 * a24) - a22 * (a13 * a44 - a43 * a14) + a42 * (a13 * a24 - a23 * a14)) - 70 a41 * (a12 * (a23 * a34 - a33 * a24) - a22 * (a13 * a34 - a33 * a14) + a32 * (a13 * a24 - a23 * a14))); 71 } 72 evaluateMultiDiracDeterminantCalculator73 inline T evaluate(T a11, 74 T a12, 75 T a13, 76 T a14, 77 T a15, 78 T a21, 79 T a22, 80 T a23, 81 T a24, 82 T a25, 83 T a31, 84 T a32, 85 T a33, 86 T a34, 87 T a35, 88 T a41, 89 T a42, 90 T a43, 91 T a44, 92 T a45, 93 T a51, 94 T a52, 95 T a53, 96 T a54, 97 T a55) 98 { 99 return ( 100 a11 * 101 (a22 * (a33 * (a44 * a55 - a54 * a45) - a43 * (a34 * a55 - a54 * a35) + a53 * (a34 * a45 - a44 * a35)) - 102 a32 * (a23 * (a44 * a55 - a54 * a45) - a43 * (a24 * a55 - a54 * a25) + a53 * (a24 * a45 - a44 * a25)) + 103 a42 * (a23 * (a34 * a55 - a54 * a35) - a33 * (a24 * a55 - a54 * a25) + a53 * (a24 * a35 - a34 * a25)) - 104 a52 * (a23 * (a34 * a45 - a44 * a35) - a33 * (a24 * a45 - a44 * a25) + a43 * (a24 * a35 - a34 * a25))) - 105 a21 * 106 (a12 * (a33 * (a44 * a55 - a54 * a45) - a43 * (a34 * a55 - a54 * a35) + a53 * (a34 * a45 - a44 * a35)) - 107 a32 * (a13 * (a44 * a55 - a54 * a45) - a43 * (a14 * a55 - a54 * a15) + a53 * (a14 * a45 - a44 * a15)) + 108 a42 * (a13 * (a34 * a55 - a54 * a35) - a33 * (a14 * a55 - a54 * a15) + a53 * (a14 * a35 - a34 * a15)) - 109 a52 * (a13 * (a34 * a45 - a44 * a35) - a33 * (a14 * a45 - a44 * a15) + a43 * (a14 * a35 - a34 * a15))) + 110 a31 * 111 (a12 * (a23 * (a44 * a55 - a54 * a45) - a43 * (a24 * a55 - a54 * a25) + a53 * (a24 * a45 - a44 * a25)) - 112 a22 * (a13 * (a44 * a55 - a54 * a45) - a43 * (a14 * a55 - a54 * a15) + a53 * (a14 * a45 - a44 * a15)) + 113 a42 * (a13 * (a24 * a55 - a54 * a25) - a23 * (a14 * a55 - a54 * a15) + a53 * (a14 * a25 - a24 * a15)) - 114 a52 * (a13 * (a24 * a45 - a44 * a25) - a23 * (a14 * a45 - a44 * a15) + a43 * (a14 * a25 - a24 * a15))) - 115 a41 * 116 (a12 * (a23 * (a34 * a55 - a54 * a35) - a33 * (a24 * a55 - a54 * a25) + a53 * (a24 * a35 - a34 * a25)) - 117 a22 * (a13 * (a34 * a55 - a54 * a35) - a33 * (a14 * a55 - a54 * a15) + a53 * (a14 * a35 - a34 * a15)) + 118 a32 * (a13 * (a24 * a55 - a54 * a25) - a23 * (a14 * a55 - a54 * a15) + a53 * (a14 * a25 - a24 * a15)) - 119 a52 * (a13 * (a24 * a35 - a34 * a25) - a23 * (a14 * a35 - a34 * a15) + a33 * (a14 * a25 - a24 * a15))) + 120 a51 * 121 (a12 * (a23 * (a34 * a45 - a44 * a35) - a33 * (a24 * a45 - a44 * a25) + a43 * (a24 * a35 - a34 * a25)) - 122 a22 * (a13 * (a34 * a45 - a44 * a35) - a33 * (a14 * a45 - a44 * a15) + a43 * (a14 * a35 - a34 * a15)) + 123 a32 * (a13 * (a24 * a45 - a44 * a25) - a23 * (a14 * a45 - a44 * a15) + a43 * (a14 * a25 - a24 * a15)) - 124 a42 * (a13 * (a24 * a35 - a34 * a25) - a23 * (a14 * a35 - a34 * a15) + a33 * (a14 * a25 - a24 * a15)))); 125 } 126 127 /** default implementation of MultiDiracDeterminant::CalculateRatioFromMatrixElements 128 * If you don't have a perfect square below 25 of dots this is what you hit 129 * dots must be a 2n by 2n Matrix 130 * iter must be size n^2 131 * this object must have been resized to n 132 */ 133 template<typename ITER> evaluateMultiDiracDeterminantCalculator134 inline T evaluate(Matrix<T>& dots, ITER it, int n) 135 { 136 typename std::vector<T>::iterator d = M.begin(); 137 for (int i = 0; i < n; i++) 138 for (int j = 0; j < n; j++) 139 { 140 //performance through proper iterator indistiquishable from data pointer 141 *(d++) = dots(*(it + i), *(it + n + j)); 142 } 143 return Determinant(M.data(), n, n, Pivot.data()); 144 } 145 }; 146 147 } // namespace qmcplusplus 148 #endif 149