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;
}
|
