1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
|
/* Ergo, version 3.8.2, a program for linear scaling electronic structure
* calculations.
* Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
* and Anastasia Kruchinina.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* Primary academic reference:
* Ergo: An open-source program for linear-scaling electronic structure
* calculations,
* Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
* Kruchinina,
* SoftwareX 7, 107 (2018),
* <http://dx.doi.org/10.1016/j.softx.2018.03.005>
*
* For further information about Ergo, see <http://www.ergoscf.org>.
*/
/** @file integrals_2el_J_mm_kernel.cc
\brief Code for multipole method computational kernel for
computing the Coulomb matrix J.
@author: Elias Rudberg <em>responsible</em>.
*/
#include "integrals_2el_J_mm_kernel.h"
int
do_multipole_interaction_between_2_boxes_branches(const distr_list_description_struct & distrDescription_1,
const multipole_struct_large & branchMultipole,
const multipole_struct_small* multipoleList_1,
ergo_real* result_J_list, // NULL if not used
ResultMatContrib* resultMatContrib, // NULL if not used
ergo_real threshold,
int* largest_L_used_so_far, // optional output, NULL if not used
MMInteractor & interactor,
const MMLimitTable & mmLimitTable)
{
int batchCount_1 = distrDescription_1.org.batchList.size();
if(batchCount_1 == 0)
return 0;
const batch_struct* batchList_1 = &distrDescription_1.org.batchList[0];
const cluster_struct* clusterList_1 = &distrDescription_1.org.clusterList[0];
const distr_group_struct* groupList_1 = &distrDescription_1.org.groupList[0];
const minimal_distr_struct* minimalDistrList_1 = &distrDescription_1.org.minimalDistrList[0];
const basis_func_pair_struct* basisFuncPairList = &distrDescription_1.org.basisFuncPairList[0];
// Prepare local result list, if needed
std::vector<ergo_real> result_J_list_local;
if(resultMatContrib) {
int result_J_list_local_size = distrDescription_1.org.basisFuncPairList.size();
result_J_list_local.resize(result_J_list_local_size);
memset(&result_J_list_local[0], 0x00, result_J_list_local_size*sizeof(ergo_real));
}
int distrCountTot = 0;
for(int batchIndex_1 = 0; batchIndex_1 < batchCount_1; batchIndex_1++) {
int clusterCount_1 = batchList_1[batchIndex_1].noOfClusters;
int cluster_start_1 = batchList_1[batchIndex_1].clusterStartIndex;
for(int clusterIndex_1 = cluster_start_1; clusterIndex_1 < cluster_start_1 + clusterCount_1; clusterIndex_1++) {
int group_start_1 = clusterList_1[clusterIndex_1].groupStartIndex;
int group_end_1 = group_start_1 + clusterList_1[clusterIndex_1].noOfGroups;
for(int groupIndex_1 = group_start_1; groupIndex_1 < group_end_1; groupIndex_1++) {
const distr_group_struct* currGroup_1 = &groupList_1[groupIndex_1];
const multipole_struct_small & multipoleCurrGroup = distrDescription_1.org_mm.multipoleListForGroups[groupIndex_1];
ergo_real dx = branchMultipole.centerCoords[0] - currGroup_1->centerCoords[0];
ergo_real dy = branchMultipole.centerCoords[1] - currGroup_1->centerCoords[1];
ergo_real dz = branchMultipole.centerCoords[2] - currGroup_1->centerCoords[2];
ergo_real r = template_blas_sqrt(dx*dx + dy*dy + dz*dz);
// loop over distrs of 1 (and at the same time over multipoles for those distrs)
// in order to find largest norm for each subvector.
int distr_start = currGroup_1->startIndex;
int distr_end = distr_start + currGroup_1->distrCount;
ergo_real maxMomentVectorNormForDistrsListCurrGroup[MAX_MULTIPOLE_DEGREE_BASIC+1];
for(int l = 0; l <= MAX_MULTIPOLE_DEGREE_BASIC; l++)
maxMomentVectorNormForDistrsListCurrGroup[l] = 0;
int maxDegreeForDistrs = 0;
for(int distrIndex = distr_start; distrIndex < distr_end; distrIndex++) {
const multipole_struct_small* distrMultipole = &multipoleList_1[distrCountTot + distrIndex - distr_start];
if(distrMultipole->degree > maxDegreeForDistrs)
maxDegreeForDistrs = distrMultipole->degree;
for(int l = 0; l <= distrMultipole->degree; l++) {
int startIndex = l*l;
int endIndex = (l+1)*(l+1);
ergo_real sum = 0;
for(int A = startIndex; A < endIndex; A++)
sum += distrMultipole->momentList[A]*distrMultipole->momentList[A];
ergo_real subNorm = template_blas_sqrt(sum);
if(subNorm > maxMomentVectorNormForDistrsListCurrGroup[l])
maxMomentVectorNormForDistrsListCurrGroup[l] = subNorm;
}
}
// check which degree is needed
int degreeNeeded = mmLimitTable.get_minimum_multipole_degree_needed(r, &branchMultipole, maxDegreeForDistrs,
maxMomentVectorNormForDistrsListCurrGroup, threshold);
if(degreeNeeded < 0)
return -1;
if(largest_L_used_so_far != NULL) {
if(degreeNeeded > *largest_L_used_so_far)
*largest_L_used_so_far = degreeNeeded;
}
int branchNoOfMoments = (degreeNeeded+1)*(degreeNeeded+1);
// create interaction matrix
ergo_real T[multipoleCurrGroup.noOfMoments * branchNoOfMoments];
interactor.getInteractionMatrix(dx, dy, dz, multipoleCurrGroup.degree, degreeNeeded, T);
ergo_real tempVector[MAX_NO_OF_MOMENTS_PER_MULTIPOLE];
for(int A = 0; A < multipoleCurrGroup.noOfMoments; A++) {
ergo_real sum = 0;
for(int B = 0; B < branchNoOfMoments; B++)
sum += branchMultipole.momentList[B] * T[A*branchNoOfMoments+B];
tempVector[A] = sum;
}
// loop over distrs of 1 (and at the same time over multipoles for those distrs)
for(int distrIndex = distr_start; distrIndex < distr_end; distrIndex++) {
const multipole_struct_small* distrMultipole = &multipoleList_1[distrCountTot];
distrCountTot++;
ergo_real sum = 0;
for(int A = 0; A < distrMultipole->noOfMoments; A++)
sum += tempVector[A] * distrMultipole->momentList[A];
int basisFuncPairIndex = minimalDistrList_1[distrIndex].basisFuncPairIndex;
int i1 = batchList_1[batchIndex_1].basisFuncPairListIndex+basisFuncPairIndex;
int pairIndex = basisFuncPairList[i1].pairIndex;
if(result_J_list)
result_J_list[pairIndex] += sum;
else
result_J_list_local[i1] += sum;
} // END FOR distrIndex
}
}
}
if(result_J_list == NULL) {
// Transfer results from result_J_list_local to resultMatContrib
assert(resultMatContrib != NULL);
for(int batchIndex_1 = 0; batchIndex_1 < batchCount_1; batchIndex_1++) {
int noOfBasisFuncPairs = batchList_1[batchIndex_1].noOfBasisFuncPairs;
for(int i = 0; i < noOfBasisFuncPairs; i++) {
int k = batchList_1[batchIndex_1].basisFuncPairListIndex+i;
int a = basisFuncPairList[k].index_1;
int b = basisFuncPairList[k].index_2;
ergo_real currContrib = result_J_list_local[k];
if(currContrib != 0)
resultMatContrib->addContrib(a, b, currContrib);
}
}
}
return 0;
}
|