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 private VOTCA includes
21 #include "Density_filter.h"
22
23 namespace votca {
24 namespace xtp {
25
Initialize(const tools::Property & options)26 void Density_filter::Initialize(const tools::Property& options) {
27 threshold_ = options.get(".").as<double>();
28 }
29
Info(Logger & log) const30 void Density_filter::Info(Logger& log) const {
31 if (threshold_ == 0.0) {
32 XTP_LOG(Log::error, log)
33 << "Using density filter with no threshold " << std::flush;
34 } else {
35 XTP_LOG(Log::error, log)
36 << "Using density filter with threshold " << threshold_ << std::flush;
37 }
38 }
39
UpdateHist(const Orbitals & orb,QMState state)40 void Density_filter::UpdateHist(const Orbitals& orb, QMState state) {
41 laststate_dmat_ = orb.DensityMatrixFull(state);
42 }
43
CalculateDNorm(const Orbitals & orb,QMStateType type) const44 Eigen::VectorXd Density_filter::CalculateDNorm(const Orbitals& orb,
45 QMStateType type) const {
46 Index nostates = orb.NumberofStates(type);
47 Eigen::VectorXd norm = Eigen::VectorXd::Zero(nostates);
48 for (Index i = 0; i < nostates; i++) {
49 QMState state(type, i, false);
50 norm(i) = (orb.DensityMatrixFull(state) - laststate_dmat_).norm();
51 }
52 double lastnorm = laststate_dmat_.norm();
53 return norm / lastnorm;
54 }
55
CalcIndeces(const Orbitals & orb,QMStateType type) const56 std::vector<Index> Density_filter::CalcIndeces(const Orbitals& orb,
57 QMStateType type) const {
58 Eigen::VectorXd Overlap = CalculateDNorm(orb, type);
59 Index offset = 0;
60 if (type == QMStateType::DQPstate) {
61 offset = orb.getGWAmin();
62 }
63 return ReduceAndSortIndecesDown(Overlap, offset, threshold_);
64 }
65
WriteToCpt(CheckpointWriter & w)66 void Density_filter::WriteToCpt(CheckpointWriter& w) {
67 w(laststate_dmat_, "laststatedmat");
68 w(threshold_, "threshold");
69 }
70
ReadFromCpt(CheckpointReader & r)71 void Density_filter::ReadFromCpt(CheckpointReader& r) {
72 r(laststate_dmat_, "laststatedmat");
73 r(threshold_, "threshold");
74 }
75
76 } // namespace xtp
77 } // namespace votca
78