1 /*
2  *            Copyright 2009-2020 The VOTCA Development Team
3  *                       (http://www.votca.org)
4  *
5  *      Licensed under the Apache License, Version 2.0 (the "License")
6  *
7  * You may not use this file except in compliance with the License.
8  * You may obtain a copy of the License at
9  *
10  *              http://www.apache.org/licenses/LICENSE-2.0
11  *
12  * Unless required by applicable law or agreed to in writing, software
13  * distributed under the License is distributed on an "AS IS" BASIS,
14  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15  * See the License for the specific language governing permissions and
16  * limitations under the License.
17  *
18  */
19 
20 // Local VOTCA includes
21 #include "votca/xtp/gridbox.h"
22 #include "votca/xtp/aobasis.h"
23 
24 namespace votca {
25 namespace xtp {
26 
27 void GridBox::FindSignificantShells(const AOBasis& basis) {
28 
29   for (const AOShell& store : basis) {
30     const double decay = store.getMinDecay();
31     const Eigen::Vector3d& shellpos = store.getPos();
32     for (const auto& point : grid_pos) {
33       Eigen::Vector3d dist = shellpos - point;
34       double distsq = dist.squaredNorm();
35       // if contribution is smaller than -ln(1e-10), add shell to list
36       if ((decay * distsq) < 20.7) {
37         addShell(&store);
38         break;
39       }
40     }
41   }
42 }
43 
44 AOShell::AOValues GridBox::CalcAOValues(const Eigen::Vector3d& point) const {
45   AOShell::AOValues result(Matrixsize());
46   for (Index j = 0; j < Shellsize(); ++j) {
47     const AOShell::AOValues val = significant_shells[j]->EvalAOspace(point);
48     result.derivatives.middleRows(aoranges[j].start, aoranges[j].size) =
49         val.derivatives;
50     result.values.segment(aoranges[j].start, aoranges[j].size) = val.values;
51   }
52   return result;
53 }
54 
55 void GridBox::AddtoBigMatrix(Eigen::MatrixXd& bigmatrix,
56                              const Eigen::MatrixXd& smallmatrix) const {
57   for (Index i = 0; i < Index(ranges.size()); i++) {
58     for (Index j = 0; j < Index(ranges.size()); j++) {
59       bigmatrix.block(ranges[i].start, ranges[j].start, ranges[i].size,
60                       ranges[j].size) +=
61           smallmatrix.block(inv_ranges[i].start, inv_ranges[j].start,
62                             inv_ranges[i].size, inv_ranges[j].size);
63     }
64   }
65   return;
66 }
67 
68 Eigen::MatrixXd GridBox::ReadFromBigMatrix(
69     const Eigen::MatrixXd& bigmatrix) const {
70   Eigen::MatrixXd matrix = Eigen::MatrixXd(matrix_size, matrix_size);
71   for (Index i = 0; i < Index(ranges.size()); i++) {
72     for (Index j = 0; j < Index(ranges.size()); j++) {
73       matrix.block(inv_ranges[i].start, inv_ranges[j].start, inv_ranges[i].size,
74                    inv_ranges[j].size) =
75           bigmatrix.block(ranges[i].start, ranges[j].start, ranges[i].size,
76                           ranges[j].size);
77     }
78   }
79   return matrix;
80 }
81 
82 Eigen::VectorXd GridBox::ReadFromBigVector(
83     const Eigen::VectorXd& bigvector) const {
84   Eigen::VectorXd vector = Eigen::VectorXd(matrix_size);
85   for (Index i = 0; i < Index(ranges.size()); i++) {
86     vector.segment(inv_ranges[i].start, inv_ranges[i].size) =
87         bigvector.segment(ranges[i].start, ranges[i].size);
88   }
89   return vector;
90 }
91 
92 void GridBox::PrepareForIntegration() {
93   Index index = 0;
94   aoranges = std::vector<GridboxRange>(0);
95   ranges = std::vector<GridboxRange>(0);
96   inv_ranges = std::vector<GridboxRange>(0);
97   std::vector<Index> start;
98   std::vector<Index> end;
99 
100   for (const auto shell : significant_shells) {
101     GridboxRange temp;
102     temp.size = shell->getNumFunc();
103     temp.start = index;
104     aoranges.push_back(temp);
105     index += shell->getNumFunc();
106     start.push_back(shell->getStartIndex());
107     end.push_back(shell->getStartIndex() + shell->getNumFunc());
108   }
109   std::vector<Index> startindex;
110   std::vector<Index> endindex;
111 
112   if (start.size() > 1) {
113     startindex.push_back(start[0]);
114 
115     for (Index i = 0; i < Index(start.size()) - 1; ++i) {
116 
117       if (end[i] != start[i + 1]) {
118         startindex.push_back(start[i + 1]);
119         endindex.push_back(end[i]);
120       }
121     }
122     endindex.push_back(end[end.size() - 1]);
123   } else {
124     startindex = start;
125     endindex = end;
126   }
127   Index shellstart = 0;
128   for (Index i = 0; i < Index(startindex.size()); ++i) {
129     Index size = endindex[i] - startindex[i];
130     GridboxRange temp;
131     temp.size = size;
132     temp.start = startindex[i];
133     ranges.push_back(temp);
134     GridboxRange temp2;
135     temp2.size = size;
136     temp2.start = shellstart;
137     inv_ranges.push_back(temp2);
138     shellstart += size;
139   }
140   return;
141 }
142 
143 }  // namespace xtp
144 }  // namespace votca
145