/*
* Copyright 2009-2019 The VOTCA Development Team (http://www.votca.org)
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*/
#define BOOST_TEST_MAIN
#define BOOST_TEST_MODULE sigma_test
#include <boost/test/unit_test.hpp>
#include <fstream>
#include <votca/xtp/aobasis.h>
#include <votca/xtp/orbitals.h>
#include <votca/xtp/rpa.h>
#include <votca/xtp/sigma_exact.h>
#include <votca/xtp/threecenter.h>
using namespace votca::xtp;
using namespace std;
BOOST_AUTO_TEST_SUITE(sigma_test)
BOOST_AUTO_TEST_CASE(sigma_full) {
ofstream xyzfile("molecule.xyz");
xyzfile << " 5" << endl;
xyzfile << " methane" << endl;
xyzfile << " C .000000 .000000 .000000" << endl;
xyzfile << " H .629118 .629118 .629118" << endl;
xyzfile << " H -.629118 -.629118 .629118" << endl;
xyzfile << " H .629118 -.629118 -.629118" << endl;
xyzfile << " H -.629118 .629118 -.629118" << endl;
xyzfile.close();
ofstream basisfile("3-21G.xml");
basisfile << "<basis name=\"3-21G\">" << endl;
basisfile << " <element name=\"H\">" << endl;
basisfile << " <shell scale=\"1.0\" type=\"S\">" << endl;
basisfile << " <constant decay=\"5.447178e+00\">" << endl;
basisfile << " <contractions factor=\"1.562850e-01\" type=\"S\"/>"
<< endl;
basisfile << " </constant>" << endl;
basisfile << " <constant decay=\"8.245470e-01\">" << endl;
basisfile << " <contractions factor=\"9.046910e-01\" type=\"S\"/>"
<< endl;
basisfile << " </constant>" << endl;
basisfile << " </shell>" << endl;
basisfile << " <shell scale=\"1.0\" type=\"S\">" << endl;
basisfile << " <constant decay=\"1.831920e-01\">" << endl;
basisfile << " <contractions factor=\"1.000000e+00\" type=\"S\"/>"
<< endl;
basisfile << " </constant>" << endl;
basisfile << " </shell>" << endl;
basisfile << " </element>" << endl;
basisfile << " <element name=\"C\">" << endl;
basisfile << " <shell scale=\"1.0\" type=\"S\">" << endl;
basisfile << " <constant decay=\"1.722560e+02\">" << endl;
basisfile << " <contractions factor=\"6.176690e-02\" type=\"S\"/>"
<< endl;
basisfile << " </constant>" << endl;
basisfile << " <constant decay=\"2.591090e+01\">" << endl;
basisfile << " <contractions factor=\"3.587940e-01\" type=\"S\"/>"
<< endl;
basisfile << " </constant>" << endl;
basisfile << " <constant decay=\"5.533350e+00\">" << endl;
basisfile << " <contractions factor=\"7.007130e-01\" type=\"S\"/>"
<< endl;
basisfile << " </constant>" << endl;
basisfile << " </shell>" << endl;
basisfile << " <shell scale=\"1.0\" type=\"SP\">" << endl;
basisfile << " <constant decay=\"3.664980e+00\">" << endl;
basisfile << " <contractions factor=\"-3.958970e-01\" type=\"S\"/>"
<< endl;
basisfile << " <contractions factor=\"2.364600e-01\" type=\"P\"/>"
<< endl;
basisfile << " </constant>" << endl;
basisfile << " <constant decay=\"7.705450e-01\">" << endl;
basisfile << " <contractions factor=\"1.215840e+00\" type=\"S\"/>"
<< endl;
basisfile << " <contractions factor=\"8.606190e-01\" type=\"P\"/>"
<< endl;
basisfile << " </constant>" << endl;
basisfile << " </shell>" << endl;
basisfile << " <shell scale=\"1.0\" type=\"SP\">" << endl;
basisfile << " <constant decay=\"1.958570e-01\">" << endl;
basisfile << " <contractions factor=\"1.000000e+00\" type=\"S\"/>"
<< endl;
basisfile << " <contractions factor=\"1.000000e+00\" type=\"P\"/>"
<< endl;
basisfile << " </constant>" << endl;
basisfile << " </shell>" << endl;
basisfile << " </element>" << endl;
basisfile << "</basis>" << endl;
basisfile.close();
Orbitals orbitals;
orbitals.QMAtoms().LoadFromFile("molecule.xyz");
BasisSet basis;
basis.Load("3-21G.xml");
AOBasis aobasis;
aobasis.Fill(basis, orbitals.QMAtoms());
Eigen::VectorXd mo_energy = Eigen::VectorXd::Zero(17);
mo_energy << 0.0468207, 0.0907801, 0.0907801, 0.104563, 0.592491, 0.663355,
0.663355, 0.768373, 1.69292, 1.97724, 1.97724, 2.50877, 2.98732, 3.4418,
3.4418, 4.81084, 17.1838;
Eigen::MatrixXd MOs = Eigen::MatrixXd::Zero(17, 17);
MOs << 0.0185815, 2.9133e-17, 8.49354e-17, -0.00312916, 0.0420075,
1.11356e-16, 1.85886e-17, -0.0334732, 0.0485113, -8.71556e-18,
-3.79994e-17, -0.0346485, -0.0248392, -3.32286e-22, -4.62643e-17,
0.0144472, 0.996183, 0.166534, 2.80578e-16, 7.85515e-16, -0.0299557,
0.409156, 1.10045e-15, 1.67666e-16, -0.336895, 0.608646, -1.00343e-16,
-5.04622e-16, -0.437568, -0.299751, -2.77682e-17, -6.06812e-16, 0.176594,
-0.0866677, 0.0010572, -0.0210402, -0.0345975, 0.035778, 0.0611836,
0.0374747, -0.154443, 0.0892921, 0.0842611, -0.35309, -0.0759572,
0.278374, -0.409082, -0.64367, 0.308248, -0.261525, -0.000315534,
-0.0010572, -0.0404824, 0.000922593, -0.035778, -0.0611836, -0.115015,
-0.109676, -0.0892921, -0.0842611, -0.242326, 0.267806, -0.278374,
0.409082, -0.0548848, 0.711558, 0.261525, 0.000315534, 0.0010572,
-0.0194422, 0.0355201, 0.035778, 0.0611836, -0.152489, 0.0447677,
0.0892921, 0.0842611, 0.110764, 0.343764, 0.278374, -0.409082, 0.588785,
0.403311, -0.261525, -0.000315534, -0.823783, -9.8891e-16, -3.34692e-15,
0.103497, -0.277613, -5.51463e-16, 4.95594e-17, 0.163544, 0.121215,
-7.23985e-17, -2.63149e-16, -0.259891, -0.284396, -1.74149e-16,
-6.65818e-16, 0.208987, 0.00782842, -0.0333718, 0.22696, 0.373203,
-0.337332, -0.251625, -0.144131, 0.594004, -0.329076, -0.0456626, 0.18588,
0.0399869, -0.0631275, -0.0704844, -0.231899, 0.111054, -0.189161,
0.000129868, 0.0333718, 0.436683, -0.009952, 0.337332, 0.251625, 0.442357,
0.421824, 0.329076, 0.0456626, 0.12757, -0.140984, 0.0631275, 0.0704844,
-0.0197737, 0.256357, 0.189161, -0.000129868, -0.0333718, 0.209723,
-0.383155, -0.337332, -0.251625, 0.586489, -0.172181, -0.329076,
-0.0456626, -0.0583106, -0.180971, -0.0631275, -0.0704844, 0.212125,
0.145303, -0.189161, 0.000129868, -0.00177478, 0.0553645, -0.00126176,
-0.0164247, 0.23154, -0.262519, -0.250334, -0.0135392, -0.429472, 0.45567,
-0.503583, -0.223493, -0.211802, -0.020461, 0.265268, 0.0023362,
-0.00241145, 0.294363, -0.686239, 0.0156394, 0.204055, -0.360136,
0.267096, 0.254698, 0.074687, -0.0228668, 0.132236, -0.14614, -0.174986,
-0.185046, -0.0109958, 0.142556, 0.0661743, 0.0022999, -0.00177478,
-0.0265895, 0.0485779, -0.0164247, 0.23154, 0.348055, -0.102182,
-0.0135392, -0.429472, 0.208281, 0.646413, -0.223493, -0.211802,
-0.219498, -0.150354, 0.0023362, -0.00241145, 0.294363, 0.329576,
-0.60212, 0.204055, -0.360136, -0.354123, 0.103963, 0.074687, -0.0228668,
0.0604434, 0.18759, -0.174986, -0.185046, -0.117959, -0.0808008,
0.0661743, 0.0022999, -0.00177478, -0.028775, -0.0473162, -0.0164247,
0.23154, -0.0855356, 0.352515, -0.0135392, -0.429472, -0.663951, -0.14283,
-0.223493, -0.211802, 0.239959, -0.114914, 0.0023362, -0.00241145,
0.294363, 0.356664, 0.586481, 0.204055, -0.360136, 0.0870267, -0.358661,
0.074687, -0.0228668, -0.192679, -0.0414494, -0.174986, -0.185046,
0.128955, -0.0617554, 0.0661743, 0.0022999, 0.00741062, -3.87173e-16,
-4.31863e-16, 0.0468488, -0.0476991, 7.27357e-16, 1.23654e-15, -0.43422,
-0.159247, -4.34945e-17, 1.2743e-16, 0.503528, -0.228856, -7.97629e-17,
-2.53026e-16, 0.689669, -0.00301027, 0.173046, 7.91486e-15, 8.39419e-15,
-0.717804, -0.0195249, -1.10754e-15, -1.66789e-15, 0.551371, 0.0684292,
4.15572e-17, -1.84233e-16, 0.0105378, -0.148396, -1.63792e-16,
-4.6499e-16, 0.351571, 0.00210309;
Logger log;
TCMatrix_gwbse Mmn{log};
Mmn.Initialize(aobasis.AOBasisSize(), 0, 16, 0, 16);
Mmn.Fill(aobasis, aobasis, MOs);
RPA rpa(log, Mmn);
rpa.setRPAInputEnergies(mo_energy);
rpa.configure(4, 0, 16);
Sigma_Exact sigma = Sigma_Exact(Mmn, rpa);
Sigma_Exact::options opt;
opt.homo = 4;
opt.qpmin = 0;
opt.qpmax = 16;
opt.rpamin = 0;
opt.rpamax = 16;
sigma.configure(opt);
Eigen::MatrixXd x = sigma.CalcExchangeMatrix();
Eigen::MatrixXd x_ref = Eigen::MatrixXd::Zero(17, 17);
x_ref << -0.00370412, -1.59566e-08, -6.75073e-09, -0.000684112, 0.0188911,
-4.4624e-09, 1.22949e-08, -0.000795391, 0.00504307, 1.72119e-09,
-9.19106e-09, 0.0140303, 0.01811, -2.89711e-09, 3.40918e-09, -0.00266758,
2.00221e-05, -1.59566e-08, -0.00682224, -1.70145e-09, 2.72016e-09,
1.6257e-07, 0.00561074, 0.00559988, 3.98481e-09, 2.85634e-08, 0.00746612,
-0.00788248, 1.34435e-07, 1.84461e-07, -0.00064411, 0.0064346,
-2.38641e-08, -1.82782e-09, -6.75073e-09, -1.70145e-09, -0.00682222,
7.21988e-10, 6.40818e-08, -0.00559984, 0.00561072, 1.36378e-09,
1.16052e-08, 0.00788246, 0.00746608, 5.28253e-08, 7.22461e-08,
-0.00643455, -0.000644091, -9.27453e-09, -7.38479e-10, -0.000684112,
2.72016e-09, 7.21988e-10, -0.00595334, -0.00819783, 6.55678e-09,
-4.38787e-09, -0.0026928, 0.000460416, -8.03933e-10, 4.27309e-09,
-0.00228116, -0.0103568, -2.40859e-09, 1.68126e-09, -0.00102159,
-0.000111773, 0.0188911, 1.6257e-07, 6.40818e-08, -0.00819783, -0.209862,
8.00987e-08, -1.22485e-07, -0.0159641, -0.0364151, -1.40586e-08,
1.05127e-07, -0.154293, -0.21813, 5.43914e-09, -5.48907e-09, 0.00986884,
5.90199e-05, -4.4624e-09, 0.00561074, -0.00559984, 6.55678e-09,
8.00987e-08, -0.0291401, -5.4019e-09, 8.80569e-09, 1.14844e-08,
0.000642242, 0.024554, 5.87742e-08, 8.70793e-08, -0.0162291, -0.0198789,
-2.58402e-09, 4.9345e-10, 1.22949e-08, 0.00559988, 0.00561072,
-4.38787e-09, -1.22485e-07, -5.4019e-09, -0.0291401, -1.60476e-08,
-2.48907e-08, -0.024554, 0.000642261, -9.04449e-08, -1.27343e-07,
0.0198789, -0.0162291, -2.8592e-09, 4.71519e-10, -0.000795391,
3.98481e-09, 1.36378e-09, -0.0026928, -0.0159641, 8.80569e-09,
-1.60476e-08, -0.0115818, -0.00171152, -2.83121e-09, 7.39872e-09,
-0.00921298, -0.0142798, -2.05408e-09, -4.97901e-09, -0.0100267,
-9.13292e-05, 0.00504307, 2.85634e-08, 1.16052e-08, 0.000460416,
-0.0364151, 1.14844e-08, -2.48907e-08, -0.00171152, -0.0090557,
-4.06853e-09, 1.74842e-08, -0.0276487, -0.0363408, 4.8481e-09,
-6.02765e-09, 0.0021409, -2.62996e-05, 1.72119e-09, 0.00746612,
0.00788246, -8.03933e-10, -1.40586e-08, 0.000642242, -0.024554,
-2.83121e-09, -4.06853e-09, -0.0261176, 9.98743e-09, -2.24351e-08,
-2.98729e-08, 0.0176435, -0.0136511, 5.75818e-09, 1.88005e-09,
-9.19106e-09, -0.00788248, 0.00746608, 4.27309e-09, 1.05127e-07, 0.024554,
0.000642261, 7.39872e-09, 1.74842e-08, 9.98743e-09, -0.0261176,
8.24848e-08, 1.17485e-07, 0.0136511, 0.0176435, -8.51658e-09,
-8.53563e-10, 0.0140303, 1.34435e-07, 5.28253e-08, -0.00228116, -0.154293,
5.87742e-08, -9.04449e-08, -0.00921298, -0.0276487, -2.24351e-08,
8.24848e-08, -0.117991, -0.162036, 6.73442e-09, -6.66614e-10, 0.0100442,
0.000214376, 0.01811, 1.84461e-07, 7.22461e-08, -0.0103568, -0.21813,
8.70793e-08, -1.27343e-07, -0.0142798, -0.0363408, -2.98729e-08,
1.17485e-07, -0.162036, -0.234311, 6.81854e-09, -4.59511e-11, 0.0131863,
-9.13286e-05, -2.89711e-09, -0.00064411, -0.00643455, -2.40859e-09,
5.43914e-09, -0.0162291, 0.0198789, -2.05408e-09, 4.8481e-09, 0.0176435,
0.0136511, 6.73442e-09, 6.81854e-09, -0.0242219, 6.69066e-09, -1.5142e-09,
-5.43482e-10, 3.40918e-09, 0.0064346, -0.000644091, 1.68126e-09,
-5.48907e-09, -0.0198789, -0.0162291, -4.97901e-09, -6.02765e-09,
-0.0136511, 0.0176435, -6.66614e-10, -4.59511e-11, 6.69066e-09,
-0.0242219, -9.38585e-09, 1.1836e-09, -0.00266758, -2.38641e-08,
-9.27453e-09, -0.00102159, 0.00986884, -2.58402e-09, -2.8592e-09,
-0.0100267, 0.0021409, 5.75818e-09, -8.51658e-09, 0.0100442, 0.0131863,
-1.5142e-09, -9.38585e-09, -0.0129496, -0.000241474, 2.00221e-05,
-1.82782e-09, -7.38479e-10, -0.000111773, 5.90199e-05, 4.9345e-10,
4.71519e-10, -9.13292e-05, -2.62996e-05, 1.88005e-09, -8.53563e-10,
0.000214376, -9.13286e-05, -5.43482e-10, 1.1836e-09, -0.000241474,
-0.000593496;
bool check_x = x_ref.isApprox(x, 1e-5);
if (!check_x) {
cout << "Sigma X" << endl;
cout << x << endl;
cout << "Sigma X ref" << endl;
cout << x_ref << endl;
}
BOOST_CHECK_EQUAL(check_x, true);
sigma.PrepareScreening();
Eigen::MatrixXd c = sigma.CalcCorrelationOffDiag(mo_energy);
c.diagonal() = sigma.CalcCorrelationDiag(mo_energy);
Eigen::MatrixXd c_ref = Eigen::MatrixXd::Zero(17, 17);
c_ref << 0.00049504, 5.25671e-10, 1.14604e-09, -0.000237725, -0.00170732,
2.17028e-09, 2.21063e-09, 0.0287528, 0.00236979, 7.93769e-10,
-3.48238e-09, 0.00389123, 0.00385479, -1.30823e-08, 7.63278e-09,
-0.00924927, 0.000467862, 5.25671e-10, -0.000254971, -5.93549e-10,
-1.80089e-09, 2.50201e-08, 0.001731, 0.00172765, 2.55389e-07, 7.89582e-09,
0.00836737, -0.00883398, 1.54842e-07, 2.25869e-07, -0.00225727, 0.0225497,
-9.57052e-08, -4.6909e-09, 1.14604e-09, -5.93549e-10, -0.000254963,
-9.17251e-11, 6.17814e-09, -0.00172762, 0.00173098, 1.47149e-07,
4.59616e-09, 0.00883394, 0.00836729, 5.84135e-08, 8.17386e-08, -0.0225494,
-0.00225715, -9.92156e-09, -1.8038e-09, -0.000237725, -1.80089e-09,
-9.17251e-11, -1.9489e-05, 0.00250374, 3.8445e-09, 2.74965e-09, 0.0305465,
0.00162667, -1.60766e-09, -1.79559e-09, 0.00369192, -0.000977068,
1.07777e-08, 2.14276e-08, -0.0063554, -3.46277e-05, -0.00170732,
2.50201e-08, 6.17814e-09, 0.00250374, -0.0077376, 1.36329e-08,
-3.98221e-08, -0.178039, -0.00688451, -5.01666e-09, 5.81575e-08,
-0.0987428, -0.154697, 2.07288e-07, -1.48835e-08, 0.116353, 0.00188979,
2.17028e-09, 0.001731, -0.00172762, 3.8445e-09, 1.36329e-08, -0.032511,
-9.0345e-09, 6.78846e-08, 2.07284e-08, 0.0014161, 0.0541402, 7.79436e-08,
2.04687e-07, -0.176453, -0.216136, -6.15182e-08, -1.73517e-10,
2.21063e-09, 0.00172765, 0.00173098, 2.74965e-09, -3.98221e-08,
-9.0345e-09, -0.032511, -2.49668e-07, -5.74159e-08, -0.0541402,
0.00141615, -2.29796e-07, -2.76846e-07, 0.216136, -0.176453, -7.59654e-08,
1.01268e-09, 0.0287528, 2.55389e-07, 1.47149e-07, 0.0305465, -0.178039,
6.78846e-08, -2.49668e-07, 0.220689, -0.0179985, -2.10184e-07,
6.56681e-09, -0.16273, -0.196948, 3.59899e-07, -8.49758e-09, 0.109104,
-0.00412267, 0.00236979, 7.89582e-09, 4.59616e-09, 0.00162667,
-0.00688451, 2.07284e-08, -5.74159e-08, -0.0179985, -0.0509681,
-4.97808e-08, -6.40726e-09, -0.0242414, -0.0200948, 2.30011e-07,
-6.0392e-08, 0.0568319, -0.00858586, 7.93769e-10, 0.00836737, 0.00883394,
-1.60766e-09, -5.01666e-09, 0.0014161, -0.0541402, -2.10184e-07,
-4.97808e-08, -0.161675, -5.67728e-09, -6.7963e-08, -4.25365e-08,
0.219819, -0.170079, -1.03912e-06, -2.0088e-09, -3.48238e-09, -0.00883398,
0.00836729, -1.79559e-09, 5.81575e-08, 0.0541402, 0.00141615, 6.56681e-09,
-6.40726e-09, -5.67728e-09, -0.161674, 1.31128e-07, 2.3699e-07, 0.170079,
0.219819, 1.53369e-07, -6.59807e-09, 0.00389123, 1.54842e-07, 5.84135e-08,
0.00369192, -0.0987428, 7.79436e-08, -2.29796e-07, -0.16273, -0.0242414,
-6.7963e-08, 1.31128e-07, -0.342331, -0.470923, 4.83395e-07, -1.98519e-07,
0.333295, 0.0168582, 0.00385479, 2.25869e-07, 8.17386e-08, -0.000977068,
-0.154697, 2.04687e-07, -2.76846e-07, -0.196948, -0.0200948, -4.25365e-08,
2.3699e-07, -0.470923, -0.750008, 6.00089e-07, 1.39026e-07, 0.386832,
0.0253594, -1.30823e-08, -0.00225727, -0.0225494, 1.07777e-08,
2.07288e-07, -0.176453, 0.216136, 3.59899e-07, 2.30011e-07, 0.219819,
0.170079, 4.83395e-07, 6.00089e-07, -0.653838, 1.58273e-07, 4.20355e-07,
7.2481e-09, 7.63278e-09, 0.0225497, -0.00225715, 2.14276e-08,
-1.48835e-08, -0.216136, -0.176453, -8.49758e-09, -6.0392e-08, -0.170079,
0.219819, -1.98519e-07, 1.39026e-07, 1.58273e-07, -0.653838, -7.90656e-08,
-7.14016e-11, -0.00924927, -9.57052e-08, -9.92156e-09, -0.0063554,
0.116353, -6.15182e-08, -7.59654e-08, 0.109104, 0.0568319, -1.03912e-06,
1.53369e-07, 0.333295, 0.386832, 4.20355e-07, -7.90656e-08, -0.389354,
-0.0138153, 0.000467862, -4.6909e-09, -1.8038e-09, -3.46277e-05,
0.00188979, -1.73517e-10, 1.01268e-09, -0.00412267, -0.00858586,
-2.0088e-09, -6.59807e-09, 0.0168582, 0.0253594, 7.2481e-09, -7.14016e-11,
-0.0138153, -0.0743819;
bool check_c_diag = c.diagonal().isApprox(c_ref.diagonal(), 1e-5);
if (!check_c_diag) {
cout << "Sigma C" << endl;
cout << c.diagonal() << endl;
cout << "Sigma C ref" << endl;
cout << c_ref.diagonal() << endl;
}
BOOST_CHECK_EQUAL(check_c_diag, true);
bool check_c = c.isApprox(c_ref, 1e-5);
if (!check_c) {
cout << "Sigma C" << endl;
cout << c << endl;
cout << "Sigma C ref" << endl;
cout << c_ref << endl;
}
BOOST_CHECK_EQUAL(check_c, true);
}
BOOST_AUTO_TEST_SUITE_END()