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