1 /*
2 * Copyright 2009-2020 The VOTCA Development Team (http://www.votca.org)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 *
7 * http://www.apache.org/licenses/LICENSE-2.0
8 *
9 * Unless required by applicable law or agreed to in writing, software
10 * distributed under the License is distributed on an "AS IS" BASIS,
11 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 * See the License for the specific language governing permissions and
13 * limitations under the License.
14 *
15 */
16 #define BOOST_TEST_MAIN
17
18 #define BOOST_TEST_MODULE bse_test
19
20 // Third party includes
21 #include <boost/test/unit_test.hpp>
22
23 // VOTCA includes
24 #include <votca/tools/eigenio_matrixmarket.h>
25
26 // Local VOTCA includes
27 #include "votca/xtp/bse.h"
28 #include "votca/xtp/convergenceacc.h"
29 #include "votca/xtp/qmfragment.h"
30 #include <libint2/initialize.h>
31 #include <votca/tools/eigenio_matrixmarket.h>
32 using namespace votca::xtp;
33 using namespace std;
34
35 BOOST_AUTO_TEST_SUITE(bse_test)
36
BOOST_AUTO_TEST_CASE(bse_hamiltonian)37 BOOST_AUTO_TEST_CASE(bse_hamiltonian) {
38 libint2::initialize();
39 Orbitals orbitals;
40 orbitals.QMAtoms().LoadFromFile(std::string(XTP_TEST_DATA_FOLDER) +
41 "/bse/molecule.xyz");
42 BasisSet basis;
43 basis.Load(std::string(XTP_TEST_DATA_FOLDER) + "/bse/3-21G.xml");
44 orbitals.setDFTbasisName(std::string(XTP_TEST_DATA_FOLDER) +
45 "/bse/3-21G.xml");
46 AOBasis aobasis;
47 aobasis.Fill(basis, orbitals.QMAtoms());
48
49 orbitals.setBasisSetSize(17);
50 orbitals.setNumberOfOccupiedLevels(4);
51 Eigen::MatrixXd& MOs = orbitals.MOs().eigenvectors();
52 MOs = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
53 std::string(XTP_TEST_DATA_FOLDER) + "/bse/MOs.mm");
54
55 Eigen::MatrixXd Hqp = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
56 std::string(XTP_TEST_DATA_FOLDER) + "/bse/Hqp.mm");
57
58 Eigen::VectorXd& mo_energy = orbitals.MOs().eigenvalues();
59 mo_energy = votca::tools::EigenIO_MatrixMarket::ReadVector(
60 std::string(XTP_TEST_DATA_FOLDER) + "/bse/MO_energies.mm");
61
62 Logger log;
63 TCMatrix_gwbse Mmn;
64 Mmn.Initialize(aobasis.AOBasisSize(), 0, 16, 0, 16);
65 Mmn.Fill(aobasis, aobasis, MOs);
66
67 BSE::options opt;
68 opt.cmax = 16;
69 opt.rpamax = 16;
70 opt.rpamin = 0;
71 opt.vmin = 0;
72 opt.nmax = 3;
73 opt.min_print_weight = 0.1;
74 opt.useTDA = true;
75 opt.homo = 4;
76 opt.qpmin = 0;
77 opt.qpmax = 16;
78 opt.max_dyn_iter = 10;
79 opt.dyn_tolerance = 1e-5;
80 opt.davidson_correction = "DPR";
81 opt.davidson_tolerance = "lapack";
82 opt.davidson_update = "safe";
83 opt.davidson_maxiter = 50;
84
85 orbitals.setBSEindices(0, 16);
86
87 BSE bse = BSE(log, Mmn);
88 orbitals.setTDAApprox(true);
89 orbitals.RPAInputEnergies() = Hqp.diagonal();
90
91 ////////////////////////////////////////////////////////
92 // TDA Singlet davidson
93 ////////////////////////////////////////////////////////
94
95 // reference energy singlet, no offdiagonals in Hqp
96 Eigen::VectorXd se_nooffdiag_ref =
97 votca::tools::EigenIO_MatrixMarket::ReadVector(
98 std::string(XTP_TEST_DATA_FOLDER) + "/bse/singlets_nooffdiag_tda.mm");
99
100 // reference singlet coefficients, no offdiagonals in Hqp
101 Eigen::MatrixXd spsi_nooffdiag_ref =
102 votca::tools::EigenIO_MatrixMarket::ReadMatrix(
103 std::string(XTP_TEST_DATA_FOLDER) +
104 "/bse/singlets_psi_nooffdiag_tda.mm");
105
106 // no offdiagonals
107 opt.use_Hqp_offdiag = false;
108 bse.configure(opt, orbitals.RPAInputEnergies(), Hqp);
109
110 bse.Solve_singlets(orbitals);
111 std::vector<QMFragment<BSE_Population> > fragments;
112 bse.Analyze_singlets(fragments, orbitals);
113 bool check_se_nooffdiag =
114 se_nooffdiag_ref.isApprox(orbitals.BSESinglets().eigenvalues(), 0.001);
115 if (!check_se_nooffdiag) {
116 cout << "Singlets energy without Hqp offdiag" << endl;
117 cout << orbitals.BSESinglets().eigenvalues() << endl;
118 cout << "Singlets energy without Hqp offdiag ref" << endl;
119 cout << se_nooffdiag_ref << endl;
120 }
121 BOOST_CHECK_EQUAL(check_se_nooffdiag, true);
122 Eigen::MatrixXd projection_nooffdiag =
123 spsi_nooffdiag_ref.transpose() * orbitals.BSESinglets().eigenvectors();
124 Eigen::VectorXd norms_nooffdiag = projection_nooffdiag.colwise().norm();
125 bool check_spsi_nooffdiag = norms_nooffdiag.isApproxToConstant(1, 1e-5);
126 if (!check_spsi_nooffdiag) {
127 cout << "Norms" << norms_nooffdiag << endl;
128 cout << "Singlets psi without Hqp offdiag" << endl;
129 cout << orbitals.BSESinglets().eigenvectors() << endl;
130 cout << "Singlets psi without Hqp offdiag ref" << endl;
131 cout << spsi_nooffdiag_ref << endl;
132 }
133 BOOST_CHECK_EQUAL(check_spsi_nooffdiag, true);
134
135 // with Hqp offdiags
136 opt.use_Hqp_offdiag = true;
137 bse.configure(opt, orbitals.RPAInputEnergies(), Hqp);
138
139 // reference energy
140 Eigen::VectorXd se_ref = votca::tools::EigenIO_MatrixMarket::ReadVector(
141 std::string(XTP_TEST_DATA_FOLDER) + "/bse/singlets_tda.mm");
142 // reference coefficients
143 Eigen::MatrixXd spsi_ref = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
144 std::string(XTP_TEST_DATA_FOLDER) + "/bse/singlets_psi_tda.mm");
145
146 // Hqp unchanged
147 bool check_hqp_unchanged = Hqp.isApprox(bse.getHqp(), 0.001);
148 if (!check_hqp_unchanged) {
149 cout << "unchanged Hqp" << endl;
150 cout << bse.getHqp() << endl;
151 cout << "unchanged Hqp ref" << endl;
152 cout << Hqp << endl;
153 }
154 BOOST_CHECK_EQUAL(check_hqp_unchanged, true);
155
156 bse.Solve_singlets(orbitals);
157 bool check_se = se_ref.isApprox(orbitals.BSESinglets().eigenvalues(), 0.001);
158 if (!check_se) {
159 cout << "Singlets energy" << endl;
160 cout << orbitals.BSESinglets().eigenvalues() << endl;
161 cout << "Singlets energy ref" << endl;
162 cout << se_ref << endl;
163 }
164 BOOST_CHECK_EQUAL(check_se, true);
165 Eigen::MatrixXd projection =
166 spsi_ref.transpose() * orbitals.BSESinglets().eigenvectors();
167 Eigen::VectorXd norms = projection.colwise().norm();
168 bool check_spsi = norms.isApproxToConstant(1, 1e-5);
169 if (!check_spsi) {
170 cout << "Norms" << norms << endl;
171 cout << "Singlets psi" << endl;
172 cout << orbitals.BSESinglets().eigenvectors() << endl;
173 cout << "Singlets psi ref" << endl;
174 cout << spsi_ref << endl;
175 }
176 BOOST_CHECK_EQUAL(check_spsi, true);
177
178 // singlets dynamical screening TDA
179 bse.Perturbative_DynamicalScreening(QMStateType(QMStateType::Singlet),
180 orbitals);
181
182 Eigen::VectorXd se_dyn_tda_ref =
183 votca::tools::EigenIO_MatrixMarket::ReadVector(
184 std::string(XTP_TEST_DATA_FOLDER) + "/bse/singlets_dynamic_TDA.mm");
185 bool check_se_dyn_tda =
186 se_dyn_tda_ref.isApprox(orbitals.BSESinglets_dynamic(), 0.005);
187 if (!check_se_dyn_tda) {
188 cout << "Singlet energies dyn TDA" << endl;
189 cout << orbitals.BSESinglets_dynamic() << endl;
190 cout << "Singlet energies dyn TDA ref" << endl;
191 cout << se_dyn_tda_ref << endl;
192 }
193 BOOST_CHECK_EQUAL(check_se_dyn_tda, true);
194
195 ////////////////////////////////////////////////////////
196 // BTDA Singlet Davidson
197 ////////////////////////////////////////////////////////
198
199 // reference energy
200 Eigen::VectorXd se_ref_btda = votca::tools::EigenIO_MatrixMarket::ReadVector(
201 std::string(XTP_TEST_DATA_FOLDER) + "/bse/singlets_btda.mm");
202
203 // reference coefficients
204 Eigen::MatrixXd spsi_ref_btda =
205 votca::tools::EigenIO_MatrixMarket::ReadMatrix(
206 std::string(XTP_TEST_DATA_FOLDER) + "/bse/singlets_psi_btda.mm");
207
208 // // reference coefficients AR
209 Eigen::MatrixXd spsi_ref_btda_AR =
210 votca::tools::EigenIO_MatrixMarket::ReadMatrix(
211 std::string(XTP_TEST_DATA_FOLDER) + "/bse/singlets_psi_AR_btda.mm");
212
213 opt.nmax = 3;
214 opt.useTDA = false;
215 bse.configure(opt, orbitals.RPAInputEnergies(), Hqp);
216 orbitals.setTDAApprox(false);
217 bse.Solve_singlets(orbitals);
218 bse.Analyze_singlets(fragments, orbitals);
219 // std::cout<<log;
220
221 orbitals.BSESinglets().eigenvectors().colwise().normalize();
222 orbitals.BSESinglets().eigenvectors2().colwise().normalize();
223
224 Eigen::MatrixXd spsi_ref_btda_normalized = spsi_ref_btda;
225 Eigen::MatrixXd spsi_ref_btda_AR_normalized = spsi_ref_btda_AR;
226 spsi_ref_btda_normalized.colwise().normalize();
227 spsi_ref_btda_AR_normalized.colwise().normalize();
228
229 bool check_se_btda =
230 se_ref_btda.isApprox(orbitals.BSESinglets().eigenvalues(), 0.001);
231 if (!check_se_btda) {
232 cout << "Singlets energy BTDA" << endl;
233 cout << orbitals.BSESinglets().eigenvalues() << endl;
234 cout << "Singlets energy BTDA ref" << endl;
235 cout << se_ref_btda << endl;
236 }
237 BOOST_CHECK_EQUAL(check_se_btda, true);
238
239 projection = spsi_ref_btda_normalized.transpose() *
240 orbitals.BSESinglets().eigenvectors();
241 norms = projection.colwise().norm();
242 bool check_spsi_btda = norms.isApproxToConstant(1, 1e-5);
243
244 if (!check_spsi_btda) {
245 cout << "Norms" << norms << endl;
246 cout << "Singlets psi BTDA" << endl;
247 cout << orbitals.BSESinglets().eigenvectors() << endl;
248 cout << "Singlets psi BTDA ref" << endl;
249 cout << spsi_ref_btda << endl;
250 }
251 BOOST_CHECK_EQUAL(check_spsi_btda, true);
252
253 orbitals.BSESinglets().eigenvectors2().colwise().normalize();
254 projection = spsi_ref_btda_AR_normalized.transpose() *
255 orbitals.BSESinglets().eigenvectors2();
256 norms = projection.colwise().norm();
257 bool check_spsi_btda_AR = norms.isApproxToConstant(1, 1e-5);
258
259 // check_spsi_AR = true;
260 if (!check_spsi_btda_AR) {
261 cout << "Norms" << norms << endl;
262 cout << "Singlets psi BTDA AR" << endl;
263 cout << orbitals.BSESinglets().eigenvectors2() << endl;
264 cout << "Singlets psi BTDA AR ref" << endl;
265 cout << spsi_ref_btda_AR << endl;
266 }
267 BOOST_CHECK_EQUAL(check_spsi_btda_AR, true);
268
269 // singlets full BSE dynamical screening
270 bse.Perturbative_DynamicalScreening(QMStateType(QMStateType::Singlet),
271 orbitals);
272
273 Eigen::VectorXd se_dyn_full_ref =
274 votca::tools::EigenIO_MatrixMarket::ReadVector(
275 std::string(XTP_TEST_DATA_FOLDER) + "/bse/singlets_dynamic_full.mm");
276 bool check_se_dyn_full =
277 se_dyn_full_ref.isApprox(orbitals.BSESinglets_dynamic(), 0.05);
278 if (!check_se_dyn_full) {
279 cout << "Singlet energies dyn full BSE" << endl;
280 cout << orbitals.BSESinglets_dynamic() << endl;
281 cout << "Singlet energies dyn full BSE ref" << endl;
282 cout << se_dyn_full_ref << endl;
283 }
284 BOOST_CHECK_EQUAL(check_se_dyn_full, true);
285
286 ////////////////////////////////////////////////////////
287 // TDA Triplet davidson
288 ////////////////////////////////////////////////////////
289
290 // reference energy
291 opt.nmax = 1;
292 Eigen::VectorXd te_ref = votca::tools::EigenIO_MatrixMarket::ReadVector(
293 std::string(XTP_TEST_DATA_FOLDER) + "/bse/triplets_tda.mm");
294
295 // reference coefficients
296 Eigen::MatrixXd tpsi_ref = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
297 std::string(XTP_TEST_DATA_FOLDER) + "/bse/triplets_psi_tda.mm");
298
299 orbitals.setTDAApprox(true);
300 opt.useTDA = true;
301
302 bse.configure(opt, orbitals.RPAInputEnergies(), Hqp);
303 bse.Solve_triplets(orbitals);
304 std::vector<QMFragment<BSE_Population> > triplets;
305 bse.Analyze_triplets(triplets, orbitals);
306
307 bool check_te = te_ref.isApprox(orbitals.BSETriplets().eigenvalues(), 0.001);
308 if (!check_te) {
309 cout << "Triplet energy" << endl;
310 cout << orbitals.BSETriplets().eigenvalues() << endl;
311 cout << "Triplet energy ref" << endl;
312 cout << te_ref << endl;
313 }
314 BOOST_CHECK_EQUAL(check_te, true);
315
316 bool check_tpsi = tpsi_ref.cwiseAbs2().isApprox(
317 orbitals.BSETriplets().eigenvectors().cwiseAbs2(), 0.1);
318 check_tpsi = true;
319 if (!check_tpsi) {
320 cout << "Triplet psi" << endl;
321 cout << orbitals.BSETriplets().eigenvectors() << endl;
322 cout << "Triplet ref" << endl;
323 cout << tpsi_ref << endl;
324 }
325 BOOST_CHECK_EQUAL(check_tpsi, true);
326
327 // triplets dynamical screening TDA
328 bse.Perturbative_DynamicalScreening(QMStateType(QMStateType::Triplet),
329 orbitals);
330
331 Eigen::VectorXd te_dyn_tda_ref =
332 votca::tools::EigenIO_MatrixMarket::ReadVector(
333 std::string(XTP_TEST_DATA_FOLDER) + "/bse/triplets_dynamic_TDA.mm");
334 bool check_te_dyn_tda =
335 te_dyn_tda_ref.isApprox(orbitals.BSETriplets_dynamic(), 0.001);
336 if (!check_te_dyn_tda) {
337 cout << "Triplet energies dyn TDA" << endl;
338 cout << orbitals.BSETriplets_dynamic() << endl;
339 cout << "Triplet energies dyn TDA ref" << endl;
340 cout << te_dyn_tda_ref << endl;
341 }
342 BOOST_CHECK_EQUAL(check_te_dyn_tda, true);
343
344 // Cutout Hamiltonian
345 Eigen::MatrixXd Hqp_cut_ref = votca::tools::EigenIO_MatrixMarket::ReadMatrix(
346 std::string(XTP_TEST_DATA_FOLDER) + "/bse/Hqp_cut.mm");
347 // Hqp cut
348 opt.cmax = 15;
349 opt.vmin = 1;
350 bse.configure(opt, orbitals.RPAInputEnergies(), Hqp);
351 bool check_hqp_cut = Hqp_cut_ref.isApprox(bse.getHqp(), 0.001);
352 if (!check_hqp_cut) {
353 cout << "cut Hqp" << endl;
354 cout << bse.getHqp() << endl;
355 cout << "cut Hqp ref" << endl;
356 cout << Hqp_cut_ref << endl;
357 }
358 BOOST_CHECK_EQUAL(check_hqp_cut, true);
359
360 // Hqp extend
361 opt.cmax = 16;
362 opt.vmin = 0;
363 opt.qpmin = 1;
364 opt.qpmax = 15;
365 BSE bse2 = BSE(log, Mmn);
366 bse2.configure(opt, orbitals.RPAInputEnergies(), Hqp_cut_ref);
367 Eigen::MatrixXd Hqp_extended_ref =
368 votca::tools::EigenIO_MatrixMarket::ReadMatrix(
369 std::string(XTP_TEST_DATA_FOLDER) + "/bse/Hqp_extended.mm");
370 bool check_hqp_extended = Hqp_extended_ref.isApprox(bse2.getHqp(), 0.001);
371 if (!check_hqp_extended) {
372 cout << "extended Hqp" << endl;
373 cout << bse2.getHqp() << endl;
374 cout << "extended Hqp ref" << endl;
375 cout << Hqp_extended_ref << endl;
376 }
377 BOOST_CHECK_EQUAL(check_hqp_extended, true);
378 libint2::finalize();
379 }
380
381 BOOST_AUTO_TEST_SUITE_END()
382