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/threecenter.h"
22 #include "votca/xtp/aomatrix.h"
23 #include "votca/xtp/openmp_cuda.h"
24 #include "votca/xtp/symmetric_matrix.h"
25 
26 namespace votca {
27 namespace xtp {
28 
29 void TCMatrix_gwbse::Initialize(Index basissize, Index mmin, Index mmax,
30                                 Index nmin, Index nmax) {
31 
32   // here as storage indices starting from zero
33   nmin_ = nmin;
34   nmax_ = nmax;
35   ntotal_ = nmax - nmin + 1;
36   mmin_ = mmin;
37   mmax_ = mmax;
38   mtotal_ = mmax - mmin + 1;
39   auxbasissize_ = basissize;
40 
41   // vector has mtotal elements
42   // largest object should be allocated in multithread fashion
43   matrix_ = std::vector<Eigen::MatrixXd>(mtotal_);
44 #pragma omp parallel for schedule(dynamic, 4)
45   for (Index i = 0; i < mtotal_; i++) {
46     matrix_[i] = Eigen::MatrixXd::Zero(ntotal_, auxbasissize_);
47   }
48 }
49 
50 /*
51  * Modify 3-center matrix elements consistent with use of symmetrized
52  * Coulomb interaction using either CUDA or Openmp.
53  */
54 void TCMatrix_gwbse::MultiplyRightWithAuxMatrix(const Eigen::MatrixXd& matrix) {
55   OpenMP_CUDA gemm;
56   gemm.setOperators(matrix_, matrix);
57 #pragma omp parallel
58   {
59     Index threadid = OPENMP::getThreadId();
60 #pragma omp for schedule(dynamic)
61     for (Index i = 0; i < msize(); i++) {
62       gemm.MultiplyRight(matrix_[i], threadid);
63     }
64   }
65 }
66 /*
67  * Fill the 3-center object by looping over shells of GW basis set and
68  * calling FillBlock, which calculates all 3-center overlap integrals
69  * associated to a particular shell, convoluted with the DFT orbital
70  * coefficients
71  */
72 void TCMatrix_gwbse::Fill(const AOBasis& auxbasis, const AOBasis& dftbasis,
73                           const Eigen::MatrixXd& dft_orbitals) {
74   // needed for Rebuild())
75   auxbasis_ = &auxbasis;
76   dftbasis_ = &dftbasis;
77   dft_orbitals_ = &dft_orbitals;
78 
79   Fill3cMO(auxbasis, dftbasis, dft_orbitals);
80 
81   AOOverlap auxoverlap;
82   auxoverlap.Fill(auxbasis);
83   AOCoulomb auxcoulomb;
84   auxcoulomb.Fill(auxbasis);
85   Eigen::MatrixXd inv_sqrt = auxcoulomb.Pseudo_InvSqrt_GWBSE(auxoverlap, 5e-7);
86   removedfunctions_ = auxcoulomb.Removedfunctions();
87   MultiplyRightWithAuxMatrix(inv_sqrt);
88 
89   return;
90 }
91 
92 }  // namespace xtp
93 }  // namespace votca
94