/*
* 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 aomatrix_test
#include <boost/test/unit_test.hpp>
#include <votca/xtp/aomatrix.h>
#include <votca/xtp/orbitals.h>
using namespace votca::xtp;
using namespace votca;
using namespace std;
BOOST_AUTO_TEST_SUITE(aomatrix_test)
QMMolecule Methane() {
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();
QMMolecule mol(" ", 0);
mol.LoadFromFile("molecule.xyz");
return mol;
}
BOOST_AUTO_TEST_CASE(aomatrices_test) {
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();
QMMolecule mol = Methane();
BasisSet basis;
basis.Load("3-21G.xml");
AOBasis aobasis;
aobasis.Fill(basis, mol);
AOOverlap overlap;
overlap.Fill(aobasis);
Eigen::MatrixXd overlap_ref = Eigen::MatrixXd::Zero(17, 17);
overlap_ref << 1, 0.191448, 0, 0, 0, 0.180314, 0, 0, 0, 0.0189724, 0.0808612,
0.0189724, 0.0808612, 0.0189724, 0.0808612, 0.0189724, 0.0808612,
0.191448, 1.00001, 0, 0, 0, 0.761361, 0, 0, 0, 0.194748, 0.401447,
0.194748, 0.401447, 0.194748, 0.401447, 0.194748, 0.401447, 0, 0, 1, 0, 0,
0, 0.528959, 0, 0, 0.169584, 0.135615, 0.169584, 0.135615, -0.169584,
-0.135615, -0.169584, -0.135615, 0, 0, 0, 1, 0, 0, 0, 0.528959, 0,
0.169584, 0.135615, -0.169584, -0.135615, -0.169584, -0.135615, 0.169584,
0.135615, 0, 0, 0, 0, 1, 0, 0, 0, 0.528959, 0.169584, 0.135615, -0.169584,
-0.135615, 0.169584, 0.135615, -0.169584, -0.135615, 0.180314, 0.761361,
0, 0, 0, 1, 0, 0, 0, 0.338796, 0.668849, 0.338796, 0.668849, 0.338796,
0.668849, 0.338796, 0.668849, 0, 0, 0.528959, 0, 0, 0, 1, 0, 0, 0.290649,
0.340149, 0.290649, 0.340149, -0.290649, -0.340149, -0.290649, -0.340149,
0, 0, 0, 0.528959, 0, 0, 0, 1, 0, 0.290649, 0.340149, -0.290649,
-0.340149, -0.290649, -0.340149, 0.290649, 0.340149, 0, 0, 0, 0, 0.528959,
0, 0, 0, 1, 0.290649, 0.340149, -0.290649, -0.340149, 0.290649, 0.340149,
-0.290649, -0.340149, 0.0189724, 0.194748, 0.169584, 0.169584, 0.169584,
0.338796, 0.290649, 0.290649, 0.290649, 1, 0.645899, 0.00778321, 0.116994,
0.00778321, 0.116994, 0.00778321, 0.116994, 0.0808612, 0.401447, 0.135615,
0.135615, 0.135615, 0.668849, 0.340149, 0.340149, 0.340149, 0.645899, 1,
0.116994, 0.354983, 0.116994, 0.354983, 0.116994, 0.354983, 0.0189724,
0.194748, 0.169584, -0.169584, -0.169584, 0.338796, 0.290649, -0.290649,
-0.290649, 0.00778321, 0.116994, 1, 0.645899, 0.00778321, 0.116994,
0.00778321, 0.116994, 0.0808612, 0.401447, 0.135615, -0.135615, -0.135615,
0.668849, 0.340149, -0.340149, -0.340149, 0.116994, 0.354983, 0.645899, 1,
0.116994, 0.354983, 0.116994, 0.354983, 0.0189724, 0.194748, -0.169584,
-0.169584, 0.169584, 0.338796, -0.290649, -0.290649, 0.290649, 0.00778321,
0.116994, 0.00778321, 0.116994, 1, 0.645899, 0.00778321, 0.116994,
0.0808612, 0.401447, -0.135615, -0.135615, 0.135615, 0.668849, -0.340149,
-0.340149, 0.340149, 0.116994, 0.354983, 0.116994, 0.354983, 0.645899, 1,
0.116994, 0.354983, 0.0189724, 0.194748, -0.169584, 0.169584, -0.169584,
0.338796, -0.290649, 0.290649, -0.290649, 0.00778321, 0.116994,
0.00778321, 0.116994, 0.00778321, 0.116994, 1, 0.645899, 0.0808612,
0.401447, -0.135615, 0.135615, -0.135615, 0.668849, -0.340149, 0.340149,
-0.340149, 0.116994, 0.354983, 0.116994, 0.354983, 0.116994, 0.354983,
0.645899, 1;
bool check_overlap = overlap.Matrix().isApprox(overlap_ref, 0.0001);
BOOST_CHECK_EQUAL(check_overlap, 1);
if (!check_overlap) {
cout << "ref" << endl;
cout << overlap_ref << endl;
cout << "result" << endl;
cout << overlap.Matrix() << endl;
}
AOKinetic kinetic;
kinetic.Fill(aobasis);
Eigen::MatrixXd kinetic_ref = Eigen::MatrixXd::Zero(17, 17);
kinetic_ref << 16.579, -1.43503, 0, 0, 0, 0.10275, 0, 0, 0, -0.0439437,
0.0214514, -0.0439437, 0.0214514, -0.0439437, 0.0214514, -0.0439437,
0.0214514, -1.43503, 1.35738, 0, 0, 0, 0.346414, 0, 0, 0, -0.0154143,
0.103305, -0.0154143, 0.103305, -0.0154143, 0.103305, -0.0154143,
0.103305, 0, 0, 2.58667, 0, 0, 0, 0.41751, 0, 0, 0.0928809, 0.0755731,
0.0928809, 0.0755731, -0.0928809, -0.0755731, -0.0928809, -0.0755731, 0,
0, 0, 2.58667, 0, 0, 0, 0.41751, 0, 0.0928809, 0.0755731, -0.0928809,
-0.0755731, -0.0928809, -0.0755731, 0.0928809, 0.0755731, 0, 0, 0, 0,
2.58667, 0, 0, 0, 0.41751, 0.0928809, 0.0755731, -0.0928809, -0.0755731,
0.0928809, 0.0755731, -0.0928809, -0.0755731, 0.10275, 0.346414, 0, 0, 0,
0.293786, 0, 0, 0, 0.0889197, 0.139112, 0.0889197, 0.139112, 0.0889197,
0.139112, 0.0889197, 0.139112, 0, 0, 0.41751, 0, 0, 0, 0.489642, 0, 0,
0.169257, 0.135141, 0.169257, 0.135141, -0.169257, -0.135141, -0.169257,
-0.135141, 0, 0, 0, 0.41751, 0, 0, 0, 0.489642, 0, 0.169257, 0.135141,
-0.169257, -0.135141, -0.169257, -0.135141, 0.169257, 0.135141, 0, 0, 0,
0, 0.41751, 0, 0, 0, 0.489642, 0.169257, 0.135141, -0.169257, -0.135141,
0.169257, 0.135141, -0.169257, -0.135141, -0.0439437, -0.0154143,
0.0928809, 0.0928809, 0.0928809, 0.0889197, 0.169257, 0.169257, 0.169257,
1.5494, 0.293152, -0.0206172, -0.00736852, -0.0206172, -0.00736852,
-0.0206172, -0.00736852, 0.0214514, 0.103305, 0.0755731, 0.0755731,
0.0755731, 0.139112, 0.135141, 0.135141, 0.135141, 0.293152, 0.274788,
-0.00736852, 0.0301943, -0.00736852, 0.0301943, -0.00736852, 0.0301943,
-0.0439437, -0.0154143, 0.0928809, -0.0928809, -0.0928809, 0.0889197,
0.169257, -0.169257, -0.169257, -0.0206172, -0.00736852, 1.5494, 0.293152,
-0.0206172, -0.00736852, -0.0206172, -0.00736852, 0.0214514, 0.103305,
0.0755731, -0.0755731, -0.0755731, 0.139112, 0.135141, -0.135141,
-0.135141, -0.00736852, 0.0301943, 0.293152, 0.274788, -0.00736852,
0.0301943, -0.00736852, 0.0301943, -0.0439437, -0.0154143, -0.0928809,
-0.0928809, 0.0928809, 0.0889197, -0.169257, -0.169257, 0.169257,
-0.0206172, -0.00736852, -0.0206172, -0.00736852, 1.5494, 0.293152,
-0.0206172, -0.00736852, 0.0214514, 0.103305, -0.0755731, -0.0755731,
0.0755731, 0.139112, -0.135141, -0.135141, 0.135141, -0.00736852,
0.0301943, -0.00736852, 0.0301943, 0.293152, 0.274788, -0.00736852,
0.0301943, -0.0439437, -0.0154143, -0.0928809, 0.0928809, -0.0928809,
0.0889197, -0.169257, 0.169257, -0.169257, -0.0206172, -0.00736852,
-0.0206172, -0.00736852, -0.0206172, -0.00736852, 1.5494, 0.293152,
0.0214514, 0.103305, -0.0755731, 0.0755731, -0.0755731, 0.139112,
-0.135141, 0.135141, -0.135141, -0.00736852, 0.0301943, -0.00736852,
0.0301943, -0.00736852, 0.0301943, 0.293152, 0.274788;
bool check_kinetic = kinetic.Matrix().isApprox(kinetic_ref, 0.00001);
BOOST_CHECK_EQUAL(check_kinetic, 1);
if (!check_kinetic) {
cout << "ref" << endl;
cout << kinetic_ref << endl;
cout << "result" << endl;
cout << kinetic.Matrix() << endl;
}
AOCoulomb coulomb;
coulomb.Fill(aobasis);
Eigen::MatrixXd coulomb_ref = Eigen::MatrixXd::Zero(17, 17);
coulomb_ref << 1.66592, 4.0152, 0, 0, 0, 5.96515, 0, 0, 0, 1.86584, 4.83598,
1.86584, 4.83598, 1.86584, 4.83598, 1.86584, 4.83598, 4.0152, 18.379, 0,
0, 0, 31.3743, 0, 0, 0, 10.3031, 26.6091, 10.3031, 26.6091, 10.3031,
26.6091, 10.3031, 26.6091, 0, 0, 4.75942, 0, 0, 0, 7.02939, 0, 0, 1.96588,
2.48722, 1.96588, 2.48722, -1.96588, -2.48722, -1.96588, -2.48722, 0, 0,
0, 4.75942, 0, 0, 0, 7.02939, 0, 1.96588, 2.48722, -1.96588, -2.48722,
-1.96588, -2.48722, 1.96588, 2.48722, 0, 0, 0, 0, 4.75942, 0, 0, 0,
7.02939, 1.96588, 2.48722, -1.96588, -2.48722, 1.96588, 2.48722, -1.96588,
-2.48722, 5.96515, 31.3743, 0, 0, 0, 64.1609, 0, 0, 0, 21.3671, 58.4239,
21.3671, 58.4239, 21.3671, 58.4239, 21.3671, 58.4239, 0, 0, 7.02939, 0, 0,
0, 21.387, 0, 0, 5.07919, 8.88641, 5.07919, 8.88641, -5.07919, -8.88641,
-5.07919, -8.88641, 0, 0, 0, 7.02939, 0, 0, 0, 21.387, 0, 5.07919,
8.88641, -5.07919, -8.88641, -5.07919, -8.88641, 5.07919, 8.88641, 0, 0,
0, 0, 7.02939, 0, 0, 0, 21.387, 5.07919, 8.88641, -5.07919, -8.88641,
5.07919, 8.88641, -5.07919, -8.88641, 1.86584, 10.3031, 1.96588, 1.96588,
1.96588, 21.3671, 5.07919, 5.07919, 5.07919, 13.9085, 26.8615, 5.54669,
17.0411, 5.54669, 17.0411, 5.54669, 17.0411, 4.83598, 26.6091, 2.48722,
2.48722, 2.48722, 58.4239, 8.88641, 8.88641, 8.88641, 26.8615, 68.5967,
17.0411, 50.7702, 17.0411, 50.7702, 17.0411, 50.7702, 1.86584, 10.3031,
1.96588, -1.96588, -1.96588, 21.3671, 5.07919, -5.07919, -5.07919,
5.54669, 17.0411, 13.9085, 26.8615, 5.54669, 17.0411, 5.54669, 17.0411,
4.83598, 26.6091, 2.48722, -2.48722, -2.48722, 58.4239, 8.88641, -8.88641,
-8.88641, 17.0411, 50.7702, 26.8615, 68.5967, 17.0411, 50.7702, 17.0411,
50.7702, 1.86584, 10.3031, -1.96588, -1.96588, 1.96588, 21.3671, -5.07919,
-5.07919, 5.07919, 5.54669, 17.0411, 5.54669, 17.0411, 13.9085, 26.8615,
5.54669, 17.0411, 4.83598, 26.6091, -2.48722, -2.48722, 2.48722, 58.4239,
-8.88641, -8.88641, 8.88641, 17.0411, 50.7702, 17.0411, 50.7702, 26.8615,
68.5967, 17.0411, 50.7702, 1.86584, 10.3031, -1.96588, 1.96588, -1.96588,
21.3671, -5.07919, 5.07919, -5.07919, 5.54669, 17.0411, 5.54669, 17.0411,
5.54669, 17.0411, 13.9085, 26.8615, 4.83598, 26.6091, -2.48722, 2.48722,
-2.48722, 58.4239, -8.88641, 8.88641, -8.88641, 17.0411, 50.7702, 17.0411,
50.7702, 17.0411, 50.7702, 26.8615, 68.5967;
bool check_coulomb = coulomb.Matrix().isApprox(coulomb_ref, 0.00001);
BOOST_CHECK_EQUAL(check_coulomb, 1);
if (!check_coulomb) {
cout << "ref" << endl;
cout << coulomb_ref << endl;
cout << "result" << endl;
cout << coulomb.Matrix() << endl;
}
Eigen::MatrixXd ps_invSqrt = coulomb.Pseudo_InvSqrt(1e-7);
Eigen::MatrixXd coulombinvsqrt_ref = Eigen::MatrixXd::Zero(17, 17);
coulombinvsqrt_ref << 1.16163, -0.307153, 9.97835e-16, -1.10548e-16,
1.00508e-16, 0.0492417, -1.93564e-15, 1.85948e-16, -1.33024e-16,
0.0226533, -0.00856968, 0.0226533, -0.00856968, 0.0226533, -0.00856968,
0.0226533, -0.00856968, -0.307153, 0.829138, -1.96961e-15, 3.70613e-16,
-3.72916e-16, -0.608979, 3.60427e-15, -5.74697e-16, 5.22424e-16,
-0.0218795, 0.0808469, -0.0218795, 0.0808469, -0.0218795, 0.0808469,
-0.0218795, 0.0808469, 9.97835e-16, -1.96961e-15, 0.679881, 2.75893e-16,
3.26907e-17, 2.17372e-15, -0.207947, -3.5126e-16, -4.68835e-16,
-0.0646361, 0.0554385, -0.0646361, 0.0554385, 0.0646361, -0.0554385,
0.0646361, -0.0554385, -1.10548e-16, 3.70613e-16, 2.76025e-16, 0.679881,
-9.91215e-17, -6.96261e-16, 8.15636e-17, -0.207947, 3.66591e-16,
-0.0646361, 0.0554385, 0.0646361, -0.0554385, 0.0646361, -0.0554385,
-0.0646361, 0.0554385, 1.00508e-16, -3.72916e-16, 3.28582e-17,
-9.70689e-17, 0.679881, 5.26717e-16, 1.97761e-16, 5.74318e-16, -0.207947,
-0.0646361, 0.0554385, 0.0646361, -0.0554385, -0.0646361, 0.0554385,
0.0646361, -0.0554385, 0.0492417, -0.608979, 2.17372e-15, -6.96261e-16,
5.26717e-16, 1.36693, -3.70734e-15, 9.17285e-16, -8.05922e-16, -0.119143,
-0.233467, -0.119143, -0.233467, -0.119143, -0.233467, -0.119143,
-0.233467, -1.93564e-15, 3.60427e-15, -0.207947, 8.14087e-17, 1.96659e-16,
-3.70734e-15, 0.526054, 1.22817e-16, 2.41716e-16, -0.012717, -0.162726,
-0.012717, -0.162726, 0.012717, 0.162726, 0.012717, 0.162726, 1.85948e-16,
-5.74697e-16, -3.50753e-16, -0.207947, 5.70762e-16, 9.17285e-16,
1.1573e-16, 0.526054, -7.70072e-16, -0.012717, -0.162726, 0.012717,
0.162726, 0.012717, 0.162726, -0.012717, -0.162726, -1.33024e-16,
5.22424e-16, -4.6736e-16, 3.63935e-16, -0.207947, -8.05922e-16,
2.42516e-16, -7.66731e-16, 0.526054, -0.012717, -0.162726, 0.012717,
0.162726, -0.012717, -0.162726, 0.012717, 0.162726, 0.0226533, -0.0218795,
-0.0646361, -0.0646361, -0.0646361, -0.119143, -0.012717, -0.012717,
-0.012717, 0.606742, -0.164568, 0.0157438, 0.0286141, 0.0157438,
0.0286141, 0.0157438, 0.0286141, -0.00856968, 0.0808469, 0.0554385,
0.0554385, 0.0554385, -0.233467, -0.162726, -0.162726, -0.162726,
-0.164568, 0.512794, 0.0286141, -0.0727664, 0.0286141, -0.0727664,
0.0286141, -0.0727664, 0.0226533, -0.0218795, -0.0646361, 0.0646361,
0.0646361, -0.119143, -0.012717, 0.012717, 0.012717, 0.0157438, 0.0286141,
0.606742, -0.164568, 0.0157438, 0.0286141, 0.0157438, 0.0286141,
-0.00856968, 0.0808469, 0.0554385, -0.0554385, -0.0554385, -0.233467,
-0.162726, 0.162726, 0.162726, 0.0286141, -0.0727664, -0.164568, 0.512794,
0.0286141, -0.0727664, 0.0286141, -0.0727664, 0.0226533, -0.0218795,
0.0646361, 0.0646361, -0.0646361, -0.119143, 0.012717, 0.012717,
-0.012717, 0.0157438, 0.0286141, 0.0157438, 0.0286141, 0.606742,
-0.164568, 0.0157438, 0.0286141, -0.00856968, 0.0808469, -0.0554385,
-0.0554385, 0.0554385, -0.233467, 0.162726, 0.162726, -0.162726,
0.0286141, -0.0727664, 0.0286141, -0.0727664, -0.164568, 0.512794,
0.0286141, -0.0727664, 0.0226533, -0.0218795, 0.0646361, -0.0646361,
0.0646361, -0.119143, 0.012717, -0.012717, 0.012717, 0.0157438, 0.0286141,
0.0157438, 0.0286141, 0.0157438, 0.0286141, 0.606742, -0.164568,
-0.00856968, 0.0808469, -0.0554385, 0.0554385, -0.0554385, -0.233467,
0.162726, -0.162726, 0.162726, 0.0286141, -0.0727664, 0.0286141,
-0.0727664, 0.0286141, -0.0727664, -0.164568, 0.512794;
bool check_coulombinvsqrt = ps_invSqrt.isApprox(coulombinvsqrt_ref, 0.00001);
BOOST_CHECK_EQUAL(check_coulombinvsqrt, 1);
if (!check_coulombinvsqrt) {
cout << "ref" << endl;
cout << coulombinvsqrt_ref << endl;
cout << "result" << endl;
cout << ps_invSqrt << endl;
}
Eigen::MatrixXd ps_invSqrtgw = coulomb.Pseudo_InvSqrt_GWBSE(overlap, 1e-7);
Eigen::MatrixXd coulombinvsqrtgw_ref = Eigen::MatrixXd::Zero(17, 17);
coulombinvsqrtgw_ref << 1.18216, -0.221181, 3.31007e-16, -2.77017e-17,
5.4648e-17, -0.024762, -4.67451e-17, 2.42791e-18, -7.39467e-17, 0.0191022,
0.00100653, 0.0191022, 0.00100653, 0.0191022, 0.00100653, 0.0191022,
0.00100653, -0.387336, 0.905821, -9.70665e-16, 1.12468e-16, -7.20284e-17,
-0.458094, -2.51973e-16, 2.48351e-16, 4.71999e-16, 0.00887867, 0.00238888,
0.00887867, 0.00238888, 0.00887867, 0.00238888, 0.00887867, 0.00238888,
4.23087e-16, -6.09524e-16, 0.705798, 4.81316e-16, 1.3367e-16, 3.47315e-17,
-0.13188, -1.27084e-16, 3.48185e-17, -0.0613791, 0.0311124, -0.0613791,
0.0311124, 0.0613791, -0.0311124, 0.0613791, -0.0311124, 1.0619e-17,
4.2506e-17, 4.78159e-16, 0.705798, -1.81479e-16, 7.98158e-17,
-7.45323e-17, -0.13188, 1.77653e-16, -0.0613791, 0.0311124, 0.0613791,
-0.0311124, 0.0613791, -0.0311124, -0.0613791, 0.0311124, 8.8482e-17,
-1.10927e-16, 7.25119e-17, -1.27572e-16, 0.705798, 9.94042e-17,
1.74247e-16, 9.21421e-17, -0.13188, -0.0613791, 0.0311124, 0.0613791,
-0.0311124, -0.0613791, 0.0311124, 0.0613791, -0.0311124, 0.140772,
-0.845289, 1.08307e-15, -1.38689e-16, 2.19957e-17, 1.29026, 6.23914e-16,
-5.71606e-16, -7.72672e-16, -0.162867, -0.0533785, -0.162867, -0.0533785,
-0.162867, -0.0533785, -0.162867, -0.0533785, -2.36379e-16, 1.27991e-16,
-0.281562, -1.88063e-16, 1.68519e-16, 4.23859e-16, 0.519419, 2.15099e-16,
-2.20189e-16, -0.00395714, -0.139109, -0.00395714, -0.139109, 0.00395714,
0.139109, 0.00395714, 0.139109, -2.23714e-16, 4.63617e-16, -3.5789e-16,
-0.281562, 1.691e-16, -6.9076e-16, 4.32912e-16, 0.519419, -1.50764e-16,
-0.00395714, -0.139109, 0.00395714, 0.139109, 0.00395714, 0.139109,
-0.00395714, -0.139109, -1.29103e-16, 2.23031e-16, -1.79135e-16,
1.09785e-16, -0.281562, -3.18918e-16, -2.12058e-16, 3.09506e-17, 0.519419,
-0.00395714, -0.139109, 0.00395714, 0.139109, -0.00395714, -0.139109,
0.00395714, 0.139109, 0.0210071, -0.0124102, -0.0704645, -0.0704645,
-0.0704645, -0.0878611, -0.00143061, -0.00143061, -0.00143061, 0.628099,
-0.0864583, 0.0165705, 0.0141384, 0.0165705, 0.0141384, 0.0165705,
0.0141384, -0.0234403, 0.126287, 0.0856204, 0.0856204, 0.0856204,
-0.249203, -0.182259, -0.182259, -0.182259, -0.210269, 0.423453,
0.0481652, -0.0956113, 0.0481652, -0.0956113, 0.0481652, -0.0956113,
0.0210071, -0.0124102, -0.0704645, 0.0704645, 0.0704645, -0.0878611,
-0.00143061, 0.00143061, 0.00143061, 0.0165705, 0.0141384, 0.628099,
-0.0864583, 0.0165705, 0.0141384, 0.0165705, 0.0141384, -0.0234403,
0.126287, 0.0856204, -0.0856204, -0.0856204, -0.249203, -0.182259,
0.182259, 0.182259, 0.0481652, -0.0956113, -0.210269, 0.423453, 0.0481652,
-0.0956113, 0.0481652, -0.0956113, 0.0210071, -0.0124102, 0.0704645,
0.0704645, -0.0704645, -0.0878611, 0.00143061, 0.00143061, -0.00143061,
0.0165705, 0.0141384, 0.0165705, 0.0141384, 0.628099, -0.0864583,
0.0165705, 0.0141384, -0.0234403, 0.126287, -0.0856204, -0.0856204,
0.0856204, -0.249203, 0.182259, 0.182259, -0.182259, 0.0481652,
-0.0956113, 0.0481652, -0.0956113, -0.210269, 0.423453, 0.0481652,
-0.0956113, 0.0210071, -0.0124102, 0.0704645, -0.0704645, 0.0704645,
-0.0878611, 0.00143061, -0.00143061, 0.00143061, 0.0165705, 0.0141384,
0.0165705, 0.0141384, 0.0165705, 0.0141384, 0.628099, -0.0864583,
-0.0234403, 0.126287, -0.0856204, 0.0856204, -0.0856204, -0.249203,
0.182259, -0.182259, 0.182259, 0.0481652, -0.0956113, 0.0481652,
-0.0956113, 0.0481652, -0.0956113, -0.210269, 0.423453;
bool check_coulombinvsqrtgw =
ps_invSqrtgw.isApprox(coulombinvsqrtgw_ref, 0.00001);
BOOST_CHECK_EQUAL(check_coulombinvsqrtgw, 1);
if (!check_coulombinvsqrtgw) {
cout << "ref" << endl;
cout << coulombinvsqrtgw_ref << endl;
cout << "result" << endl;
cout << ps_invSqrtgw << endl;
}
}
BOOST_AUTO_TEST_CASE(aomatrices_contracted_test) {
std::ofstream basisfile("contracted.xml");
basisfile << "<basis name=\"cc-pVTZ\">" << std::endl;
basisfile << "<element name=\"C\">" << std::endl;
basisfile << " <shell scale=\"1.0\" type=\"S\">" << std::endl;
basisfile << " <constant decay=\"8236.0\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"0.0005424302\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"1235.0\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"0.0041964279\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"280.8\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"0.0215409141\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"79.27\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"0.0836149496\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"25.59\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"0.2398716189\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"8.997\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"0.4437518201\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"3.319\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"0.3535796965\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"0.3643\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"-0.0091763661\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " </shell>" << std::endl;
basisfile << " <shell scale=\"1.0\" type=\"S\">" << std::endl;
basisfile << " <constant decay=\"8236.0\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"-0.0001963922\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"1235.0\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"-0.0015259503\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"280.8\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"-0.007890449\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"79.27\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"-0.0315148705\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"25.59\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"-0.0969100083\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"8.997\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"-0.2205415263\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"3.319\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"-0.2960691129\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"0.3643\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"1.0405034329\" type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " </shell>" << std::endl;
basisfile << " <shell scale=\"1.0\" type=\"S\">" << std::endl;
basisfile << " <constant decay=\"0.9059\">"
<< std::endl;
basisfile << " <contractions factor=\"1.0\" "
"type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " </shell>" << std::endl;
basisfile << " <shell scale=\"1.0\" type=\"S\">" << std::endl;
basisfile << " <constant decay=\"0.1285\">"
<< std::endl;
basisfile << " <contractions factor=\"1.0\" "
"type=\"S\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " </shell>" << std::endl;
basisfile << " <shell scale=\"1.0\" type=\"P\">" << std::endl;
basisfile << " <constant decay=\"18.71\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"0.0394263872\" type=\"P\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"4.133\">"
<< std::endl;
basisfile << " <contractions "
"factor=\"0.2440889849\" type=\"P\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " <constant decay=\"1.2\">" << std::endl;
basisfile << " <contractions "
"factor=\"0.8154920089\" type=\"P\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " </shell>" << std::endl;
basisfile << " <shell scale=\"1.0\" type=\"P\">" << std::endl;
basisfile << " <constant decay=\"0.3827\">"
<< std::endl;
basisfile << " <contractions factor=\"1.0\" "
"type=\"P\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " </shell>" << std::endl;
basisfile << " <shell scale=\"1.0\" type=\"P\">" << std::endl;
basisfile << " <constant decay=\"0.1209\">"
<< std::endl;
basisfile << " <contractions factor=\"1.0\" "
"type=\"P\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " </shell>" << std::endl;
basisfile << " <shell scale=\"1.0\" type=\"D\">" << std::endl;
basisfile << " <constant decay=\"1.097\">"
<< std::endl;
basisfile << " <contractions factor=\"1.0\" "
"type=\"D\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " </shell>" << std::endl;
basisfile << " <shell scale=\"1.0\" type=\"D\">" << std::endl;
basisfile << " <constant decay=\"0.318\">"
<< std::endl;
basisfile << " <contractions factor=\"1.0\" "
"type=\"D\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " </shell>" << std::endl;
basisfile << " <shell scale=\"1.0\" type=\"F\">" << std::endl;
basisfile << " <constant decay=\"0.761\">"
<< std::endl;
basisfile << " <contractions factor=\"1.0\" "
"type=\"F\"/>"
<< std::endl;
basisfile << " </constant>" << std::endl;
basisfile << " </shell>" << std::endl;
basisfile << " </element>" << std::endl;
basisfile << "</basis>" << std::endl;
basisfile.close();
std::ofstream xyzfile("C.xyz");
xyzfile << " 1" << std::endl;
xyzfile << " C" << std::endl;
xyzfile << " C .000000 .000000 .000000" << std::endl;
xyzfile.close();
QMMolecule mol("C", 0);
mol.LoadFromFile("C.xyz");
BasisSet basis;
basis.Load("contracted.xml");
AOBasis aobasis;
aobasis.Fill(basis, mol);
AOOverlap overlap;
overlap.Fill(aobasis);
Eigen::MatrixXd overlap_ref =
Eigen::MatrixXd::Zero(aobasis.AOBasisSize(), aobasis.AOBasisSize());
overlap_ref << 1, -0.271918, 0.510892, 0.140478, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.271918, 1, 0.553776,
0.755963, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0.510892, 0.553776, 1, 0.535796, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.140478, 0.755963,
0.535796, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0.611573, 0, 0, 0.221927, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, -2.7742e-17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0.611573, 0, 0, 0.221927, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
2.95039e-17, 0, 0, 0, 5.55112e-17, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0.611573,
0, 0, 0.221927, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-5.55112e-17, 0, 0, 0, 0, 0.611573, 0, 0, 1, 0, 0, 0.674476, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1.11022e-16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.611573, 0, 0, 1, 0, 0, 0.674476, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-5.55112e-17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.611573, 0, 0, 1, 0, 0,
0.674476, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.11022e-16, 0, 0, 0, 0, 0,
0, 0, 0, 0.221927, 0, 0, 0.674476, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1.11022e-16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.221927, 0, 0,
0.674476, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4.16334e-17, 0, 0,
0, -5.55112e-17, 0, 0, 0, 0, 0, 0, 0, 0.221927, 0, 0, 0.674476, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -5.55112e-17, 0, 0, 0, 5.55112e-17, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0.531577, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0.531577, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0.531577, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0.531577, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0.531577, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0.531577, 0, 0, 0, 0, 1, 0, 0, 0, -5.55112e-17, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.531577, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.531577, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.531577, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.531577,
-5.55112e-17, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.7742e-17, 0,
0, 1.11022e-16, 0, 0, 1.11022e-16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, -5.93751e-18, 0, 0, 0, 0, 0, 0, 0, 2.95039e-17, 0, 0,
-5.55112e-17, 0, 0, -4.16334e-17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 8.37926e-18, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.11022e-16, 0, 0,
-5.55112e-17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 8.37926e-18,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, -5.93751e-18, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 5.55112e-17, 0, 0,
0, 0, 0, -5.55112e-17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8.37926e-18, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -5.55112e-17, 0, 0, 0, 0, 0, 5.55112e-17, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8.37926e-18, 0, 0, 0, 1;
bool check_overlap = overlap.Matrix().isApprox(overlap_ref, 0.0001);
if (!check_overlap) {
std::cout << std::endl;
std::cout << "Ref" << std::endl;
std::cout << overlap_ref << std::endl;
std::cout << "Result" << std::endl;
std::cout << overlap.Matrix() << std::endl;
}
BOOST_CHECK_EQUAL(check_overlap, 1);
}
BOOST_AUTO_TEST_CASE(aocoulomb_inv_test) {
QMMolecule mol = Methane();
BasisSet basis;
basis.Load("3-21G.xml");
AOBasis aobasis;
aobasis.Fill(basis, mol);
AOCoulomb cou;
cou.Fill(aobasis);
Eigen::MatrixXd PseudoInvSqrt = cou.Pseudo_InvSqrt(1e-7);
Eigen::MatrixXd Reformed = PseudoInvSqrt * PseudoInvSqrt * cou.Matrix();
bool check_inv = Reformed.isApprox(Eigen::MatrixXd::Identity(17, 17), 0.0001);
if (!check_inv) {
std::cout << "reformed" << endl;
std::cout << Reformed << endl;
}
BOOST_CHECK_EQUAL(check_inv, 1);
}
BOOST_AUTO_TEST_CASE(large_l_test) {
std::ofstream xyzfile("C2.xyz");
xyzfile << " 2" << std::endl;
xyzfile << " C2" << std::endl;
xyzfile << " C .000000 .000000 .000000" << std::endl;
xyzfile << " C 1.000000 .000000 .000000" << std::endl;
xyzfile.close();
QMMolecule mol("C", 0);
mol.LoadFromFile("C2.xyz");
ofstream basisfile("G.xml");
basisfile << "<basis name=\"G\">" << endl;
basisfile << " <element name=\"C\">" << endl;
basisfile << " <shell scale=\"1.0\" type=\"G\">" << endl;
basisfile << " <constant decay=\"5.447178e+00\">" << endl;
basisfile << " <contractions factor=\"1.562850e-01\" type=\"G\"/>"
<< endl;
basisfile << " </constant>" << endl;
basisfile << " </shell>" << endl;
basisfile << " </element>" << endl;
basisfile << "</basis>" << std::endl;
basisfile.close();
ofstream auxbasisfile("I.xml");
auxbasisfile << "<basis name=\"I\">" << endl;
auxbasisfile << " <element name=\"C\">" << endl;
auxbasisfile << " <shell scale=\"1.0\" type=\"I\">" << endl;
auxbasisfile << " <constant decay=\"5.447178e+00\">" << endl;
auxbasisfile << " <contractions factor=\"1.562850e-01\" type=\"I\"/>"
<< endl;
auxbasisfile << " </constant>" << endl;
auxbasisfile << " </shell>" << endl;
auxbasisfile << " </element>" << endl;
auxbasisfile << "</basis>" << std::endl;
auxbasisfile.close();
BasisSet basisset;
basisset.Load("G.xml");
BasisSet auxbasisset;
auxbasisset.Load("I.xml");
AOBasis dftbasis;
dftbasis.Fill(basisset, mol);
AOBasis auxbasis;
auxbasis.Fill(auxbasisset, mol);
AOOverlap overlap;
overlap.Fill(auxbasis);
Index auxbasissize = 26;
Index dftbasissize = 18;
Eigen::MatrixXd overlap_ref =
Eigen::MatrixXd::Zero(auxbasissize, auxbasissize);
overlap_ref << 1, 0, 0, 0, 0.00980207, 0, 0, 0, 4.12241e-16, 0, 0, 0,
1.249e-16, 0.0237511, 0, 0, 0, -0.0176617, 0, 0, 0, 0.0148613, 0, 0, 0,
0.00194066, 0, 1, 0, 0, 0, 2.88549e-16, 0, 0, 0, 3.43734e-16, 0, 0, 0, 0,
0.00415193, 0, 0, 0, -0.00971203, 0, 0, 0, 0.00833903, 0, 0, 0, 0, 0, 1,
0, 0, 0, 9.21638e-18, 0, 0, 0, 4.41139e-16, 0, 0, 0, 0, -0.0433016, 0, 0,
0, 0.0310482, 0, 0, 0, -0.0111976, 0, 0, 0, 0, 0, 1, 0.00473485, 0, 0,
-8.01766e-17, 0, 0, 0, -2.61143e-16, 0, 0, 0, 0, -0.0138466, -0.00256163,
0, 0, 0.0185842, 0, 0, 0, -0.00613319, 0, 0.00980207, 0, 0, 0.00473485,
1.19531, 0, 0, -0.0103735, 0.0155603, 0, 0, 0.0210687, -0.0210687,
-0.0176617, 0, 0, -0.00256163, 0.05203, 0, 0, 0.00561224, -0.00671312, 0,
0, -0.0113985, -0.0213299, 0, 2.88549e-16, 0, 0, 0, 1, 0, 0, 0,
1.2767e-16, 0, 0, 0, 0, -0.00971203, 0, 0, 0, 0.0246417, 0, 0, 0,
-0.0263528, 0, 0, 0, 0, 0, 9.21638e-18, 0, 0, 0, 1, 0, 0, 0, 4.28107e-17,
0, 0, 0, 0, 0.0310482, 0, 0, 0, -0.0296094, 0, 0, 0, 0.0227903, 0, 0, 0,
0, 0, -8.01766e-17, -0.0103735, 0, 0, 1, 0, 0, 0, -6.61878e-17, 0, 0, 0,
0, 0.0185842, 0.00561224, 0, 0, -0.030451, 0, 0, 0, 0.0211324, 0,
4.12241e-16, 0, 0, 0, 0.0155603, 0, 0, 0, 1, 0, 0, 0, 1.11022e-16,
0.0148613, 0, 0, 0, -0.00671312, 0, 0, 0, 0.0281051, 0, 0, 0, -0.0185897,
0, 3.43734e-16, 0, 0, 0, 1.2767e-16, 0, 0, 0, 1, 0, 0, 0, 0, 0.00833903,
0, 0, 0, -0.0263528, 0, 0, 0, 0.0451492, 0, 0, 0, 0, 0, 4.41139e-16, 0, 0,
0, 4.28107e-17, 0, 0, 0, 1, 0, 0, 0, 0, -0.0111976, 0, 0, 0, 0.0227903, 0,
0, 0, -0.0414569, 0, 0, 0, 0, 0, -2.61143e-16, 0.0210687, 0, 0,
-6.61878e-17, 0, 0, 0, 1, 0, 0, 0, 0, -0.00613319, -0.0113985, 0, 0,
0.0211324, 0, 0, 0, -0.0700703, 0, 1.249e-16, 0, 0, 0, -0.0210687, 0, 0,
0, 1.11022e-16, 0, 0, 0, 1, 0.00194066, 0, 0, 0, -0.0213299, 0, 0, 0,
-0.0185897, 0, 0, 0, 0.0703197, 0.0237511, 0, 0, 0, -0.0176617, 0, 0, 0,
0.0148613, 0, 0, 0, 0.00194066, 1, 0, 0, 0, 0.00980207, 0, 0, 0,
4.12241e-16, 0, 0, 0, 1.249e-16, 0, 0.00415193, 0, 0, 0, -0.00971203, 0,
0, 0, 0.00833903, 0, 0, 0, 0, 1, 0, 0, 0, 2.88549e-16, 0, 0, 0,
3.43734e-16, 0, 0, 0, 0, 0, -0.0433016, 0, 0, 0, 0.0310482, 0, 0, 0,
-0.0111976, 0, 0, 0, 0, 1, 0, 0, 0, 9.21638e-18, 0, 0, 0, 4.41139e-16, 0,
0, 0, 0, 0, -0.0138466, -0.00256163, 0, 0, 0.0185842, 0, 0, 0,
-0.00613319, 0, 0, 0, 0, 1, 0.00473485, 0, 0, -8.01766e-17, 0, 0, 0,
-2.61143e-16, 0, -0.0176617, 0, 0, -0.00256163, 0.05203, 0, 0, 0.00561224,
-0.00671312, 0, 0, -0.0113985, -0.0213299, 0.00980207, 0, 0, 0.00473485,
1.19531, 0, 0, -0.0103735, 0.0155603, 0, 0, 0.0210687, -0.0210687, 0,
-0.00971203, 0, 0, 0, 0.0246417, 0, 0, 0, -0.0263528, 0, 0, 0, 0,
2.88549e-16, 0, 0, 0, 1, 0, 0, 0, 1.2767e-16, 0, 0, 0, 0, 0, 0.0310482, 0,
0, 0, -0.0296094, 0, 0, 0, 0.0227903, 0, 0, 0, 0, 9.21638e-18, 0, 0, 0, 1,
0, 0, 0, 4.28107e-17, 0, 0, 0, 0, 0, 0.0185842, 0.00561224, 0, 0,
-0.030451, 0, 0, 0, 0.0211324, 0, 0, 0, 0, -8.01766e-17, -0.0103735, 0, 0,
1, 0, 0, 0, -6.61878e-17, 0, 0.0148613, 0, 0, 0, -0.00671312, 0, 0, 0,
0.0281051, 0, 0, 0, -0.0185897, 4.12241e-16, 0, 0, 0, 0.0155603, 0, 0, 0,
1, 0, 0, 0, 1.11022e-16, 0, 0.00833903, 0, 0, 0, -0.0263528, 0, 0, 0,
0.0451492, 0, 0, 0, 0, 3.43734e-16, 0, 0, 0, 1.2767e-16, 0, 0, 0, 1, 0, 0,
0, 0, 0, -0.0111976, 0, 0, 0, 0.0227903, 0, 0, 0, -0.0414569, 0, 0, 0, 0,
4.41139e-16, 0, 0, 0, 4.28107e-17, 0, 0, 0, 1, 0, 0, 0, 0, 0, -0.00613319,
-0.0113985, 0, 0, 0.0211324, 0, 0, 0, -0.0700703, 0, 0, 0, 0,
-2.61143e-16, 0.0210687, 0, 0, -6.61878e-17, 0, 0, 0, 1, 0, 0.00194066, 0,
0, 0, -0.0213299, 0, 0, 0, -0.0185897, 0, 0, 0, 0.0703197, 1.249e-16, 0,
0, 0, -0.0210687, 0, 0, 0, 1.11022e-16, 0, 0, 0, 1;
bool check_overlap = overlap.Matrix().isApprox(overlap_ref, 0.00001);
BOOST_CHECK_EQUAL(check_overlap, 1);
if (!check_overlap) {
cout << "ref" << endl;
cout << overlap_ref << endl;
cout << "result" << endl;
cout << overlap.Matrix() << endl;
}
AOCoulomb coulomb;
coulomb.Fill(auxbasis);
Eigen::MatrixXd coulomb_ref =
Eigen::MatrixXd::Zero(auxbasissize, auxbasissize);
coulomb_ref << 0.177458, 0, 0, 0, 0.00173945, 0, 0, 0, -5.0899e-16, 0, 0, 0,
-6.59195e-16, 0.00737579, 0, 0, 0, -0.0046258, 0, 0, 0, 0.00251989, 0, 0,
0, 0.00101389, 0, 0.177458, 0, 0, 0, 5.72107e-16, 0, 0, 0, 2.60091e-16, 0,
0, 0, 0, 0.00186685, 0, 0, 0, -0.00372032, 0, 0, 0, 0.00167281, 0, 0, 0,
0, 0, 0.177458, 0, 0, 0, -4.81256e-16, 0, 0, 0, 4.5071e-16, 0, 0, 0, 0,
-0.0128239, 0, 0, 0, 0.00741287, 0, 0, 0, -0.000167197, 0, 0, 0, 0, 0,
0.177458, 0.000840236, 0, 0, 1.46075e-16, 0, 0, 0, 5.76574e-17, 0, 0, 0,
0, -0.00579144, -0.000946348, 0, 0, 0.00629574, 0, 0, 0, -9.15777e-05, 0,
0.00173945, 0, 0, 0.000840236, 0.659262, 0, 0, -0.00184086, 0.0027613, 0,
0, 0.00373881, -0.00373881, -0.0046258, 0, 0, -0.000946348, 0.267231, 0,
0, 0.00207334, -0.000303246, 0, 0, -0.00421098, -0.00776951, 0,
5.72107e-16, 0, 0, 0, 0.177458, 0, 0, 0, 1.67775e-16, 0, 0, 0, 0,
-0.00372032, 0, 0, 0, 0.00868561, 0, 0, 0, -0.00648608, 0, 0, 0, 0, 0,
-4.81256e-16, 0, 0, 0, 0.177458, 0, 0, 0, -1.18905e-16, 0, 0, 0, 0,
0.00741287, 0, 0, 0, -0.0078581, 0, 0, 0, 0.00504498, 0, 0, 0, 0, 0,
1.46075e-16, -0.00184086, 0, 0, 0.177458, 0, 0, 0, -1.11731e-17, 0, 0, 0,
0, 0.00629574, 0.00207334, 0, 0, -0.0104345, 0, 0, 0, 0.00329104, 0,
-5.0899e-16, 0, 0, 0, 0.0027613, 0, 0, 0, 0.177458, 0, 0, 0, 1.74513e-15,
0.00251989, 0, 0, 0, -0.000303246, 0, 0, 0, 0.00930186, 0, 0, 0,
-0.00349634, 0, 2.60091e-16, 0, 0, 0, 1.67775e-16, 0, 0, 0, 0.177458, 0,
0, 0, 0, 0.00167281, 0, 0, 0, -0.00648608, 0, 0, 0, 0.0142938, 0, 0, 0, 0,
0, 4.5071e-16, 0, 0, 0, -1.18905e-16, 0, 0, 0, 0.177458, 0, 0, 0, 0,
-0.000167197, 0, 0, 0, 0.00504498, 0, 0, 0, -0.0138951, 0, 0, 0, 0, 0,
5.76574e-17, 0.00373881, 0, 0, -1.11731e-17, 0, 0, 0, 0.177458, 0, 0, 0,
0, -9.15777e-05, -0.00421098, 0, 0, 0.00329104, 0, 0, 0, -0.0183512, 0,
-6.59195e-16, 0, 0, 0, -0.00373881, 0, 0, 0, 1.74513e-15, 0, 0, 0,
0.177458, 0.00101389, 0, 0, 0, -0.00776951, 0, 0, 0, -0.00349634, 0, 0, 0,
0.0190278, 0.00737579, 0, 0, 0, -0.0046258, 0, 0, 0, 0.00251989, 0, 0, 0,
0.00101389, 0.177458, 0, 0, 0, 0.00173945, 0, 0, 0, -5.0899e-16, 0, 0, 0,
-6.59195e-16, 0, 0.00186685, 0, 0, 0, -0.00372032, 0, 0, 0, 0.00167281, 0,
0, 0, 0, 0.177458, 0, 0, 0, 5.72107e-16, 0, 0, 0, 2.60091e-16, 0, 0, 0, 0,
0, -0.0128239, 0, 0, 0, 0.00741287, 0, 0, 0, -0.000167197, 0, 0, 0, 0,
0.177458, 0, 0, 0, -4.81256e-16, 0, 0, 0, 4.5071e-16, 0, 0, 0, 0, 0,
-0.00579144, -0.000946348, 0, 0, 0.00629574, 0, 0, 0, -9.15777e-05, 0, 0,
0, 0, 0.177458, 0.000840236, 0, 0, 1.46075e-16, 0, 0, 0, 5.76574e-17, 0,
-0.0046258, 0, 0, -0.000946348, 0.267231, 0, 0, 0.00207334, -0.000303246,
0, 0, -0.00421098, -0.00776951, 0.00173945, 0, 0, 0.000840236, 0.659262,
0, 0, -0.00184086, 0.0027613, 0, 0, 0.00373881, -0.00373881, 0,
-0.00372032, 0, 0, 0, 0.00868561, 0, 0, 0, -0.00648608, 0, 0, 0, 0,
5.72107e-16, 0, 0, 0, 0.177458, 0, 0, 0, 1.67775e-16, 0, 0, 0, 0, 0,
0.00741287, 0, 0, 0, -0.0078581, 0, 0, 0, 0.00504498, 0, 0, 0, 0,
-4.81256e-16, 0, 0, 0, 0.177458, 0, 0, 0, -1.18905e-16, 0, 0, 0, 0, 0,
0.00629574, 0.00207334, 0, 0, -0.0104345, 0, 0, 0, 0.00329104, 0, 0, 0, 0,
1.46075e-16, -0.00184086, 0, 0, 0.177458, 0, 0, 0, -1.11731e-17, 0,
0.00251989, 0, 0, 0, -0.000303246, 0, 0, 0, 0.00930186, 0, 0, 0,
-0.00349634, -5.0899e-16, 0, 0, 0, 0.0027613, 0, 0, 0, 0.177458, 0, 0, 0,
1.74513e-15, 0, 0.00167281, 0, 0, 0, -0.00648608, 0, 0, 0, 0.0142938, 0,
0, 0, 0, 2.60091e-16, 0, 0, 0, 1.67775e-16, 0, 0, 0, 0.177458, 0, 0, 0, 0,
0, -0.000167197, 0, 0, 0, 0.00504498, 0, 0, 0, -0.0138951, 0, 0, 0, 0,
4.5071e-16, 0, 0, 0, -1.18905e-16, 0, 0, 0, 0.177458, 0, 0, 0, 0, 0,
-9.15777e-05, -0.00421098, 0, 0, 0.00329104, 0, 0, 0, -0.0183512, 0, 0, 0,
0, 5.76574e-17, 0.00373881, 0, 0, -1.11731e-17, 0, 0, 0, 0.177458, 0,
0.00101389, 0, 0, 0, -0.00776951, 0, 0, 0, -0.00349634, 0, 0, 0,
0.0190278, -6.59195e-16, 0, 0, 0, -0.00373881, 0, 0, 0, 1.74513e-15, 0, 0,
0, 0.177458;
bool check_coulomb = coulomb.Matrix().isApprox(coulomb_ref, 0.00001);
BOOST_CHECK_EQUAL(check_coulomb, 1);
if (!check_coulomb) {
cout << "ref" << endl;
cout << coulomb_ref << endl;
cout << "result" << endl;
cout << coulomb.Matrix() << endl;
}
AOKinetic kinetic;
kinetic.Fill(dftbasis);
Eigen::MatrixXd kinetic_ref =
Eigen::MatrixXd::Zero(dftbasissize, dftbasissize);
kinetic_ref << 29.9595, 0, 0, 0, 1.03944e-15, 0, 0, 0, -1.98878e-16,
-0.011401, 0, 0, 0, -0.000410714, 0, 0, 0, 0.0614917, 0, 29.9595, 0, 0, 0,
1.19856e-15, 0, 0, 0, 0, -0.0116766, 0, 0, 0, 0.0272569, 0, 0, 0, 0, 0,
29.9595, 0, 0, 0, -1.87629e-15, 0, 0, 0, 0, 0.0573089, 0, 0, 0,
-0.0337547, 0, 0, 0, 0, 0, 29.9595, 0, 0, 0, 3.30788e-16, 0, 0, 0, 0,
0.0275401, 0, 0, 0, -0.0225031, 0, 1.03944e-15, 0, 0, 0, 29.9595, 0, 0, 0,
-7.11277e-16, -0.000410714, 0, 0, 0, 0.012974, 0, 0, 0, -0.0734951, 0,
1.19856e-15, 0, 0, 0, 29.9595, 0, 0, 0, 0, 0.0272569, 0, 0, 0, -0.0734894,
0, 0, 0, 0, 0, -1.87629e-15, 0, 0, 0, 29.9595, 0, 0, 0, 0, -0.0337547, 0,
0, 0, 0.0488035, 0, 0, 0, 0, 0, 3.30788e-16, 0, 0, 0, 29.9595, 0, 0, 0, 0,
-0.0225031, 0, 0, 0, 0.0785723, 0, -1.98878e-16, 0, 0, 0, -7.11277e-16, 0,
0, 0, 29.9595, 0.0614917, 0, 0, 0, -0.0734951, 0, 0, 0, 0.0237354,
-0.011401, 0, 0, 0, -0.000410714, 0, 0, 0, 0.0614917, 29.9595, 0, 0, 0,
1.03944e-15, 0, 0, 0, -1.98878e-16, 0, -0.0116766, 0, 0, 0, 0.0272569, 0,
0, 0, 0, 29.9595, 0, 0, 0, 1.19856e-15, 0, 0, 0, 0, 0, 0.0573089, 0, 0, 0,
-0.0337547, 0, 0, 0, 0, 29.9595, 0, 0, 0, -1.87629e-15, 0, 0, 0, 0, 0,
0.0275401, 0, 0, 0, -0.0225031, 0, 0, 0, 0, 29.9595, 0, 0, 0, 3.30788e-16,
0, -0.000410714, 0, 0, 0, 0.012974, 0, 0, 0, -0.0734951, 1.03944e-15, 0,
0, 0, 29.9595, 0, 0, 0, -7.11277e-16, 0, 0.0272569, 0, 0, 0, -0.0734894,
0, 0, 0, 0, 1.19856e-15, 0, 0, 0, 29.9595, 0, 0, 0, 0, 0, -0.0337547, 0,
0, 0, 0.0488035, 0, 0, 0, 0, -1.87629e-15, 0, 0, 0, 29.9595, 0, 0, 0, 0,
0, -0.0225031, 0, 0, 0, 0.0785723, 0, 0, 0, 0, 3.30788e-16, 0, 0, 0,
29.9595, 0, 0.0614917, 0, 0, 0, -0.0734951, 0, 0, 0, 0.0237354,
-1.98878e-16, 0, 0, 0, -7.11277e-16, 0, 0, 0, 29.9595;
bool check_kinetic = kinetic.Matrix().isApprox(kinetic_ref, 0.00001);
BOOST_CHECK_EQUAL(check_kinetic, 1);
if (!check_kinetic) {
cout << "ref" << endl;
cout << kinetic_ref << endl;
cout << "result" << endl;
cout << kinetic.Matrix() << endl;
}
}
BOOST_AUTO_TEST_SUITE_END()