/*
* 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 aomatrix3d_test
#include <boost/test/unit_test.hpp>
#include <votca/xtp/aomatrix3d.h>
#include <votca/xtp/orbitals.h>
using namespace votca::xtp;
using namespace votca;
using namespace std;
BOOST_AUTO_TEST_SUITE(aomatrix3d_test)
BOOST_AUTO_TEST_CASE(aomatrices3d_test) {
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());
AODipole dip;
dip.Fill(aobasis);
Eigen::MatrixXd dip0_ref = Eigen::MatrixXd::Zero(17, 17);
dip0_ref << 0, 0, 0, 0, 0.082882, 0, 0, 0, 0.0122969, 0.00265068, 0.00272642,
-0.00265068, -0.00272642, 0.00265068, 0.00272642, -0.00265068,
-0.00272642, 0, 0, 0, 0, 0.608935, 0, 0, 0, 0.388177, 0.127733, 0.10069,
-0.127733, -0.10069, 0.127733, 0.10069, -0.127733, -0.10069, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0.104648, 0.029801, -0.104648, -0.029801, -0.104648,
-0.029801, 0.104648, 0.029801, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.104648,
0.029801, 0.104648, 0.029801, -0.104648, -0.029801, -0.104648, -0.029801,
0.082882, 0.608935, 0, 0, 0, 0.597616, 0, 0, 0, 0.187889, 0.341144,
0.187889, 0.341144, 0.187889, 0.341144, 0.187889, 0.341144, 0, 0, 0, 0,
0.597616, 0, 0, 0, 1.1298, 0.328374, 0.3843, -0.328374, -0.3843, 0.328374,
0.3843, -0.328374, -0.3843, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.282166, 0.195439,
-0.282166, -0.195439, -0.282166, -0.195439, 0.282166, 0.195439, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0.282166, 0.195439, 0.282166, 0.195439, -0.282166,
-0.195439, -0.282166, -0.195439, 0.0122969, 0.388177, 0, 0, 0, 1.1298, 0,
0, 0, 0.423587, 0.976351, 0.423587, 0.976351, 0.423587, 0.976351,
0.423587, 0.976351, 0.00265068, 0.127733, 0.104648, 0.104648, 0.187889,
0.328374, 0.282166, 0.282166, 0.423587, 1.18886, 0.767884, -3.38813e-20,
0.0901008, 0.00925314, 0.139089, -3.38813e-20, 0.0901008, 0.00272642,
0.10069, 0.029801, 0.029801, 0.341144, 0.3843, 0.195439, 0.195439,
0.976351, 0.767884, 1.18886, -0.0901008, 0, 0.139089, 0.422025,
-0.0901008, 0, -0.00265068, -0.127733, -0.104648, 0.104648, 0.187889,
-0.328374, -0.282166, 0.282166, 0.423587, 3.38813e-20, -0.0901008,
-1.18886, -0.767884, 3.38813e-20, -0.0901008, -0.00925314, -0.139089,
-0.00272642, -0.10069, -0.029801, 0.029801, 0.341144, -0.3843, -0.195439,
0.195439, 0.976351, 0.0901008, 0, -0.767884, -1.18886, 0.0901008, 0,
-0.139089, -0.422025, 0.00265068, 0.127733, -0.104648, -0.104648,
0.187889, 0.328374, -0.282166, -0.282166, 0.423587, 0.00925314, 0.139089,
-3.38813e-20, 0.0901008, 1.18886, 0.767884, -3.38813e-20, 0.0901008,
0.00272642, 0.10069, -0.029801, -0.029801, 0.341144, 0.3843, -0.195439,
-0.195439, 0.976351, 0.139089, 0.422025, -0.0901008, 0, 0.767884, 1.18886,
-0.0901008, 0, -0.00265068, -0.127733, 0.104648, -0.104648, 0.187889,
-0.328374, 0.282166, -0.282166, 0.423587, 3.38813e-20, -0.0901008,
-0.00925314, -0.139089, 3.38813e-20, -0.0901008, -1.18886, -0.767884,
-0.00272642, -0.10069, 0.029801, -0.029801, 0.341144, -0.3843, 0.195439,
-0.195439, 0.976351, 0.0901008, 0, -0.139089, -0.422025, 0.0901008, 0,
-0.767884, -1.18886;
Eigen::MatrixXd dip1_ref = Eigen::MatrixXd::Zero(17, 17);
dip1_ref << 0, 0, 0, 0.082882, 0, 0, 0, 0.0122969, 0, 0.00265068, 0.00272642,
-0.00265068, -0.00272642, -0.00265068, -0.00272642, 0.00265068,
0.00272642, 0, 0, 0, 0.608935, 0, 0, 0, 0.388177, 0, 0.127733, 0.10069,
-0.127733, -0.10069, -0.127733, -0.10069, 0.127733, 0.10069, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0.104648, 0.029801, -0.104648, -0.029801, 0.104648,
0.029801, -0.104648, -0.029801, 0.082882, 0.608935, 0, 0, 0, 0.597616, 0,
0, 0, 0.187889, 0.341144, 0.187889, 0.341144, 0.187889, 0.341144,
0.187889, 0.341144, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.104648, 0.029801,
0.104648, 0.029801, -0.104648, -0.029801, -0.104648, -0.029801, 0, 0, 0,
0.597616, 0, 0, 0, 1.1298, 0, 0.328374, 0.3843, -0.328374, -0.3843,
-0.328374, -0.3843, 0.328374, 0.3843, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.282166,
0.195439, -0.282166, -0.195439, 0.282166, 0.195439, -0.282166, -0.195439,
0.0122969, 0.388177, 0, 0, 0, 1.1298, 0, 0, 0, 0.423587, 0.976351,
0.423587, 0.976351, 0.423587, 0.976351, 0.423587, 0.976351, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0.282166, 0.195439, 0.282166, 0.195439, -0.282166, -0.195439,
-0.282166, -0.195439, 0.00265068, 0.127733, 0.104648, 0.187889, 0.104648,
0.328374, 0.282166, 0.423587, 0.282166, 1.18886, 0.767884, -3.38813e-20,
0.0901008, -3.38813e-20, 0.0901008, 0.00925314, 0.139089, 0.00272642,
0.10069, 0.029801, 0.341144, 0.029801, 0.3843, 0.195439, 0.976351,
0.195439, 0.767884, 1.18886, -0.0901008, 0, -0.0901008, 0, 0.139089,
0.422025, -0.00265068, -0.127733, -0.104648, 0.187889, 0.104648,
-0.328374, -0.282166, 0.423587, 0.282166, 3.38813e-20, -0.0901008,
-1.18886, -0.767884, -0.00925314, -0.139089, 3.38813e-20, -0.0901008,
-0.00272642, -0.10069, -0.029801, 0.341144, 0.029801, -0.3843, -0.195439,
0.976351, 0.195439, 0.0901008, 0, -0.767884, -1.18886, -0.139089,
-0.422025, 0.0901008, 0, -0.00265068, -0.127733, 0.104648, 0.187889,
-0.104648, -0.328374, 0.282166, 0.423587, -0.282166, 3.38813e-20,
-0.0901008, -0.00925314, -0.139089, -1.18886, -0.767884, 3.38813e-20,
-0.0901008, -0.00272642, -0.10069, 0.029801, 0.341144, -0.029801, -0.3843,
0.195439, 0.976351, -0.195439, 0.0901008, 0, -0.139089, -0.422025,
-0.767884, -1.18886, 0.0901008, 0, 0.00265068, 0.127733, -0.104648,
0.187889, -0.104648, 0.328374, -0.282166, 0.423587, -0.282166, 0.00925314,
0.139089, -3.38813e-20, 0.0901008, -3.38813e-20, 0.0901008, 1.18886,
0.767884, 0.00272642, 0.10069, -0.029801, 0.341144, -0.029801, 0.3843,
-0.195439, 0.976351, -0.195439, 0.139089, 0.422025, -0.0901008, 0,
-0.0901008, 0, 0.767884, 1.18886;
Eigen::MatrixXd dip2_ref = Eigen::MatrixXd::Zero(17, 17);
dip2_ref << 0, 0, 0.082882, 0, 0, 0, 0.0122969, 0, 0, 0.00265068, 0.00272642,
0.00265068, 0.00272642, -0.00265068, -0.00272642, -0.00265068,
-0.00272642, 0, 0, 0.608935, 0, 0, 0, 0.388177, 0, 0, 0.127733, 0.10069,
0.127733, 0.10069, -0.127733, -0.10069, -0.127733, -0.10069, 0.082882,
0.608935, 0, 0, 0, 0.597616, 0, 0, 0, 0.187889, 0.341144, 0.187889,
0.341144, 0.187889, 0.341144, 0.187889, 0.341144, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0.104648, 0.029801, -0.104648, -0.029801, 0.104648, 0.029801,
-0.104648, -0.029801, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.104648, 0.029801,
-0.104648, -0.029801, -0.104648, -0.029801, 0.104648, 0.029801, 0, 0,
0.597616, 0, 0, 0, 1.1298, 0, 0, 0.328374, 0.3843, 0.328374, 0.3843,
-0.328374, -0.3843, -0.328374, -0.3843, 0.0122969, 0.388177, 0, 0, 0,
1.1298, 0, 0, 0, 0.423587, 0.976351, 0.423587, 0.976351, 0.423587,
0.976351, 0.423587, 0.976351, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.282166,
0.195439, -0.282166, -0.195439, 0.282166, 0.195439, -0.282166, -0.195439,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.282166, 0.195439, -0.282166, -0.195439,
-0.282166, -0.195439, 0.282166, 0.195439, 0.00265068, 0.127733, 0.187889,
0.104648, 0.104648, 0.328374, 0.423587, 0.282166, 0.282166, 1.18886,
0.767884, 0.00925314, 0.139089, -3.38813e-20, 0.0901008, -3.38813e-20,
0.0901008, 0.00272642, 0.10069, 0.341144, 0.029801, 0.029801, 0.3843,
0.976351, 0.195439, 0.195439, 0.767884, 1.18886, 0.139089, 0.422025,
-0.0901008, 0, -0.0901008, 0, 0.00265068, 0.127733, 0.187889, -0.104648,
-0.104648, 0.328374, 0.423587, -0.282166, -0.282166, 0.00925314, 0.139089,
1.18886, 0.767884, -3.38813e-20, 0.0901008, -3.38813e-20, 0.0901008,
0.00272642, 0.10069, 0.341144, -0.029801, -0.029801, 0.3843, 0.976351,
-0.195439, -0.195439, 0.139089, 0.422025, 0.767884, 1.18886, -0.0901008,
0, -0.0901008, 0, -0.00265068, -0.127733, 0.187889, 0.104648, -0.104648,
-0.328374, 0.423587, 0.282166, -0.282166, 3.38813e-20, -0.0901008,
3.38813e-20, -0.0901008, -1.18886, -0.767884, -0.00925314, -0.139089,
-0.00272642, -0.10069, 0.341144, 0.029801, -0.029801, -0.3843, 0.976351,
0.195439, -0.195439, 0.0901008, 0, 0.0901008, 0, -0.767884, -1.18886,
-0.139089, -0.422025, -0.00265068, -0.127733, 0.187889, -0.104648,
0.104648, -0.328374, 0.423587, -0.282166, 0.282166, 3.38813e-20,
-0.0901008, 3.38813e-20, -0.0901008, -0.00925314, -0.139089, -1.18886,
-0.767884, -0.00272642, -0.10069, 0.341144, -0.029801, 0.029801, -0.3843,
0.976351, -0.195439, 0.195439, 0.0901008, 0, 0.0901008, 0, -0.139089,
-0.422025, -0.767884, -1.18886;
bool check_dip0 = dip.Matrix()[0].isApprox(dip0_ref, 0.0001);
BOOST_CHECK_EQUAL(check_dip0, 1);
if (!check_dip0) {
cout << "ref" << endl;
cout << dip0_ref << endl;
cout << "result" << endl;
cout << dip.Matrix()[0] << endl;
}
bool check_dip1 = dip.Matrix()[1].isApprox(dip1_ref, 0.0001);
BOOST_CHECK_EQUAL(check_dip1, 1);
if (!check_dip1) {
cout << "ref" << endl;
cout << dip1_ref << endl;
cout << "result" << endl;
cout << dip.Matrix()[1] << endl;
}
bool check_dip2 = dip.Matrix()[2].isApprox(dip2_ref, 0.0001);
BOOST_CHECK_EQUAL(check_dip2, 1);
if (!check_dip2) {
cout << "ref" << endl;
cout << dip2_ref << endl;
cout << "result" << endl;
cout << dip.Matrix()[2] << endl;
}
AOMomentum momentum;
momentum.Fill(aobasis);
Eigen::MatrixXd momentum0_ref = Eigen::MatrixXd::Zero(17, 17);
momentum0_ref << 0, 0, 0, 0, 1.12967, 0, 0, 0, 0.154782, 0.0328303, 0.0342226,
-0.0328303, -0.0342226, 0.0328303, 0.0342226, -0.0328303, -0.0342226, 0,
0, 0, 0, 0.539423, 0, 0, 0, 0.521834, 0.179053, 0.13797, -0.179053,
-0.13797, 0.179053, 0.13797, -0.179053, -0.13797, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0.170183, 0.0481524, -0.170183, -0.0481524, -0.170183, -0.0481524,
0.170183, 0.0481524, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.170183, 0.0481524,
0.170183, 0.0481524, -0.170183, -0.0481524, -0.170183, -0.0481524,
-1.12967, -0.539423, 0, 0, 0, -0.234095, 0, 0, 0, 0.0275385, -0.0659188,
0.0275385, -0.0659188, 0.0275385, -0.0659188, 0.0275385, -0.0659188, 0, 0,
0, 0, 0.234095, 0, 0, 0, 0.442557, 0.128629, 0.150536, -0.128629,
-0.150536, 0.128629, 0.150536, -0.128629, -0.150536, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0.110528, 0.0765563, -0.110528, -0.0765563, -0.110528, -0.0765563,
0.110528, 0.0765563, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.110528, 0.0765563,
0.110528, 0.0765563, -0.110528, -0.0765563, -0.110528, -0.0765563,
-0.154782, -0.521834, 0, 0, 0, -0.442557, 0, 0, 0, -0.133948, -0.209557,
-0.133948, -0.209557, -0.133948, -0.209557, -0.133948, -0.209557,
-0.0328303, -0.179053, -0.170183, -0.170183, -0.0275385, -0.128629,
-0.110528, -0.110528, 0.133948, -2.18791e-16, -4.98587e-17, -0.0153284,
-0.0839716, -6.65767e-20, -9.1556e-18, -0.0153284, -0.0839716, -0.0342226,
-0.13797, -0.0481524, -0.0481524, 0.0659188, -0.150536, -0.0765563,
-0.0765563, 0.209557, -2.24414e-16, 0, -0.0839716, -0.154623,
-4.12093e-17, 0, -0.0839716, -0.154623, 0.0328303, 0.179053, 0.170183,
-0.170183, -0.0275385, 0.128629, 0.110528, -0.110528, 0.133948, 0.0153284,
0.0839716, 2.18791e-16, 4.98587e-17, 0.0153284, 0.0839716, 6.65767e-20,
9.1556e-18, 0.0342226, 0.13797, 0.0481524, -0.0481524, 0.0659188,
0.150536, 0.0765563, -0.0765563, 0.209557, 0.0839716, 0.154623,
2.24414e-16, 0, 0.0839716, 0.154623, 4.12093e-17, 0, -0.0328303,
-0.179053, 0.170183, 0.170183, -0.0275385, -0.128629, 0.110528, 0.110528,
0.133948, -6.65767e-20, -9.1556e-18, -0.0153284, -0.0839716, -2.18791e-16,
-4.98587e-17, -0.0153284, -0.0839716, -0.0342226, -0.13797, 0.0481524,
0.0481524, 0.0659188, -0.150536, 0.0765563, 0.0765563, 0.209557,
-4.12093e-17, 0, -0.0839716, -0.154623, -2.24414e-16, 0, -0.0839716,
-0.154623, 0.0328303, 0.179053, -0.170183, 0.170183, -0.0275385, 0.128629,
-0.110528, 0.110528, 0.133948, 0.0153284, 0.0839716, 6.65767e-20,
9.1556e-18, 0.0153284, 0.0839716, 2.18791e-16, 4.98587e-17, 0.0342226,
0.13797, -0.0481524, 0.0481524, 0.0659188, 0.150536, -0.0765563,
0.0765563, 0.209557, 0.0839716, 0.154623, 4.12093e-17, 0, 0.0839716,
0.154623, 2.24414e-16, 0;
Eigen::MatrixXd momentum1_ref = Eigen::MatrixXd::Zero(17, 17);
momentum1_ref << 0, 0, 0, 1.12967, 0, 0, 0, 0.154782, 0, 0.0328303, 0.0342226,
-0.0328303, -0.0342226, -0.0328303, -0.0342226, 0.0328303, 0.0342226, 0,
0, 0, 0.539423, 0, 0, 0, 0.521834, 0, 0.179053, 0.13797, -0.179053,
-0.13797, -0.179053, -0.13797, 0.179053, 0.13797, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0.170183, 0.0481524, -0.170183, -0.0481524, 0.170183, 0.0481524,
-0.170183, -0.0481524, -1.12967, -0.539423, 0, 0, 0, -0.234095, 0, 0, 0,
0.0275385, -0.0659188, 0.0275385, -0.0659188, 0.0275385, -0.0659188,
0.0275385, -0.0659188, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.170183, 0.0481524,
0.170183, 0.0481524, -0.170183, -0.0481524, -0.170183, -0.0481524, 0, 0,
0, 0.234095, 0, 0, 0, 0.442557, 0, 0.128629, 0.150536, -0.128629,
-0.150536, -0.128629, -0.150536, 0.128629, 0.150536, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0.110528, 0.0765563, -0.110528, -0.0765563, 0.110528, 0.0765563,
-0.110528, -0.0765563, -0.154782, -0.521834, 0, 0, 0, -0.442557, 0, 0, 0,
-0.133948, -0.209557, -0.133948, -0.209557, -0.133948, -0.209557,
-0.133948, -0.209557, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.110528, 0.0765563,
0.110528, 0.0765563, -0.110528, -0.0765563, -0.110528, -0.0765563,
-0.0328303, -0.179053, -0.170183, -0.0275385, -0.170183, -0.128629,
-0.110528, 0.133948, -0.110528, -2.18791e-16, -4.98587e-17, -0.0153284,
-0.0839716, -0.0153284, -0.0839716, -6.65767e-20, -9.1556e-18, -0.0342226,
-0.13797, -0.0481524, 0.0659188, -0.0481524, -0.150536, -0.0765563,
0.209557, -0.0765563, -2.24414e-16, 0, -0.0839716, -0.154623, -0.0839716,
-0.154623, -4.12093e-17, 0, 0.0328303, 0.179053, 0.170183, -0.0275385,
-0.170183, 0.128629, 0.110528, 0.133948, -0.110528, 0.0153284, 0.0839716,
2.18791e-16, 4.98587e-17, 6.65767e-20, 9.1556e-18, 0.0153284, 0.0839716,
0.0342226, 0.13797, 0.0481524, 0.0659188, -0.0481524, 0.150536, 0.0765563,
0.209557, -0.0765563, 0.0839716, 0.154623, 2.24414e-16, 0, 4.12093e-17, 0,
0.0839716, 0.154623, 0.0328303, 0.179053, -0.170183, -0.0275385, 0.170183,
0.128629, -0.110528, 0.133948, 0.110528, 0.0153284, 0.0839716,
6.65767e-20, 9.1556e-18, 2.18791e-16, 4.98587e-17, 0.0153284, 0.0839716,
0.0342226, 0.13797, -0.0481524, 0.0659188, 0.0481524, 0.150536,
-0.0765563, 0.209557, 0.0765563, 0.0839716, 0.154623, 4.12093e-17, 0,
2.24414e-16, 0, 0.0839716, 0.154623, -0.0328303, -0.179053, 0.170183,
-0.0275385, 0.170183, -0.128629, 0.110528, 0.133948, 0.110528,
-6.65767e-20, -9.1556e-18, -0.0153284, -0.0839716, -0.0153284, -0.0839716,
-2.18791e-16, -4.98587e-17, -0.0342226, -0.13797, 0.0481524, 0.0659188,
0.0481524, -0.150536, 0.0765563, 0.209557, 0.0765563, -4.12093e-17, 0,
-0.0839716, -0.154623, -0.0839716, -0.154623, -2.24414e-16, 0;
Eigen::MatrixXd momentum2_ref = Eigen::MatrixXd::Zero(17, 17);
momentum2_ref << 0, 0, 1.12967, 0, 0, 0, 0.154782, 0, 0, 0.0328303, 0.0342226,
0.0328303, 0.0342226, -0.0328303, -0.0342226, -0.0328303, -0.0342226, 0,
0, 0.539423, 0, 0, 0, 0.521834, 0, 0, 0.179053, 0.13797, 0.179053,
0.13797, -0.179053, -0.13797, -0.179053, -0.13797, -1.12967, -0.539423, 0,
0, 0, -0.234095, 0, 0, 0, 0.0275385, -0.0659188, 0.0275385, -0.0659188,
0.0275385, -0.0659188, 0.0275385, -0.0659188, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.170183, 0.0481524, -0.170183, -0.0481524, 0.170183, 0.0481524,
-0.170183, -0.0481524, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.170183, 0.0481524,
-0.170183, -0.0481524, -0.170183, -0.0481524, 0.170183, 0.0481524, 0, 0,
0.234095, 0, 0, 0, 0.442557, 0, 0, 0.128629, 0.150536, 0.128629, 0.150536,
-0.128629, -0.150536, -0.128629, -0.150536, -0.154782, -0.521834, 0, 0, 0,
-0.442557, 0, 0, 0, -0.133948, -0.209557, -0.133948, -0.209557, -0.133948,
-0.209557, -0.133948, -0.209557, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.110528,
0.0765563, -0.110528, -0.0765563, 0.110528, 0.0765563, -0.110528,
-0.0765563, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.110528, 0.0765563, -0.110528,
-0.0765563, -0.110528, -0.0765563, 0.110528, 0.0765563, -0.0328303,
-0.179053, -0.0275385, -0.170183, -0.170183, -0.128629, 0.133948,
-0.110528, -0.110528, -2.18791e-16, -4.98587e-17, -6.65767e-20,
-9.1556e-18, -0.0153284, -0.0839716, -0.0153284, -0.0839716, -0.0342226,
-0.13797, 0.0659188, -0.0481524, -0.0481524, -0.150536, 0.209557,
-0.0765563, -0.0765563, -2.24414e-16, 0, -4.12093e-17, 0, -0.0839716,
-0.154623, -0.0839716, -0.154623, -0.0328303, -0.179053, -0.0275385,
0.170183, 0.170183, -0.128629, 0.133948, 0.110528, 0.110528, -6.65767e-20,
-9.1556e-18, -2.18791e-16, -4.98587e-17, -0.0153284, -0.0839716,
-0.0153284, -0.0839716, -0.0342226, -0.13797, 0.0659188, 0.0481524,
0.0481524, -0.150536, 0.209557, 0.0765563, 0.0765563, -4.12093e-17, 0,
-2.24414e-16, 0, -0.0839716, -0.154623, -0.0839716, -0.154623, 0.0328303,
0.179053, -0.0275385, -0.170183, 0.170183, 0.128629, 0.133948, -0.110528,
0.110528, 0.0153284, 0.0839716, 0.0153284, 0.0839716, 2.18791e-16,
4.98587e-17, 6.65767e-20, 9.1556e-18, 0.0342226, 0.13797, 0.0659188,
-0.0481524, 0.0481524, 0.150536, 0.209557, -0.0765563, 0.0765563,
0.0839716, 0.154623, 0.0839716, 0.154623, 2.24414e-16, 0, 4.12093e-17, 0,
0.0328303, 0.179053, -0.0275385, 0.170183, -0.170183, 0.128629, 0.133948,
0.110528, -0.110528, 0.0153284, 0.0839716, 0.0153284, 0.0839716,
6.65767e-20, 9.1556e-18, 2.18791e-16, 4.98587e-17, 0.0342226, 0.13797,
0.0659188, 0.0481524, -0.0481524, 0.150536, 0.209557, 0.0765563,
-0.0765563, 0.0839716, 0.154623, 0.0839716, 0.154623, 4.12093e-17, 0,
2.24414e-16, 0;
bool check_momentum0 = momentum.Matrix()[0].isApprox(momentum0_ref, 0.0001);
BOOST_CHECK_EQUAL(check_momentum0, 1);
if (!check_momentum0) {
cout << "ref" << endl;
cout << momentum0_ref << endl;
cout << "result" << endl;
cout << momentum.Matrix()[0] << endl;
}
bool check_momentum1 = momentum.Matrix()[1].isApprox(momentum1_ref, 0.0001);
BOOST_CHECK_EQUAL(check_momentum1, 1);
if (!check_momentum1) {
cout << "ref" << endl;
cout << momentum1_ref << endl;
cout << "result" << endl;
cout << momentum.Matrix()[1] << endl;
}
bool check_momentum2 = momentum.Matrix()[2].isApprox(momentum2_ref, 0.0001);
BOOST_CHECK_EQUAL(check_momentum2, 1);
if (!check_momentum2) {
cout << "ref" << endl;
cout << momentum2_ref << endl;
cout << "result" << endl;
cout << momentum.Matrix()[2] << endl;
}
}
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();
BasisSet basisset;
basisset.Load("G.xml");
AOBasis dftbasis;
dftbasis.Fill(basisset, mol);
AODipole dip;
dip.Fill(dftbasis);
Index dftbasissize = 18;
Eigen::MatrixXd dip0_ref = Eigen::MatrixXd::Zero(dftbasissize, dftbasissize);
dip0_ref << 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00698192, 0, 0, 0, -0.00896982, 0, 0,
0, 0.00649193, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.000964788, 0, 0, 0,
-0.00240334, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0136987, 0, 0,
0, 0.0111631, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.00385377, 0, 0,
0, 0.00744207, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.00896982, 0, 0, 0,
0.0122167, 0, 0, 0, -0.0112798, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.00240334,
0, 0, 0, 0.00641505, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0111631,
0, 0, 0, -0.0108859, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00744207,
0, 0, 0, -0.0207308, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00649193, 0, 0, 0,
-0.0112798, 0, 0, 0, 0.0213368, 0.00698192, 0, 0, 0, -0.00896982, 0, 0, 0,
0.00649193, 1.88973, 0, 0, 0, 3.46492e-17, 0, 0, 0, 1.94289e-16, 0,
0.000964788, 0, 0, 0, -0.00240334, 0, 0, 0, 0, 1.88973, 0, 0, 0,
1.92546e-16, 0, 0, 0, 0, 0, -0.0136987, 0, 0, 0, 0.0111631, 0, 0, 0, 0,
1.88973, 0, 0, 0, -2.44805e-16, 0, 0, 0, 0, 0, -0.00385377, 0, 0, 0,
0.00744207, 0, 0, 0, 0, 1.88973, 0, 0, 0, 1.2951e-17, 0, -0.00896982, 0,
0, 0, 0.0122167, 0, 0, 0, -0.0112798, 1.48984e-16, 0, 0, 0, 1.88973, 0, 0,
0, -1.36647e-16, 0, -0.00240334, 0, 0, 0, 0.00641505, 0, 0, 0, 0,
1.44582e-16, 0, 0, 0, 1.88973, 0, 0, 0, 0, 0, 0.0111631, 0, 0, 0,
-0.0108859, 0, 0, 0, 0, -2.15709e-16, 0, 0, 0, 1.88973, 0, 0, 0, 0, 0,
0.00744207, 0, 0, 0, -0.0207308, 0, 0, 0, 0, 3.50989e-17, 0, 0, 0,
1.88973, 0, 0.00649193, 0, 0, 0, -0.0112798, 0, 0, 0, 0.0213368,
4.31659e-16, 0, 0, 0, -3.98495e-18, 0, 0, 0, 1.88973;
Eigen::MatrixXd dip1_ref = Eigen::MatrixXd::Zero(dftbasissize, dftbasissize);
dip1_ref << 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00150191, 0, 0, 0,
-0.00258254, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.000964788, 0, 0, 0,
0.000386357, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.000964788, 0, 0, 0,
-0.00240334, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.00150191, 0, 0, 0,
0.00153947, 0, 0, 0, 0.000142264, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0.00153947, 0, 0, 0, 0.0031273, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.00240334, 0, 0, 0, -0.00119144, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0.000386357, 0, 0, 0, 0.00119144, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.00258254, 0, 0, 0, -0.0031273, 0, 0, 0, 0.00104397, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, -0.000142264, 0, 0, 0, -0.00104397, 0, 0, 0, 0,
-0.00150191, 0, 0, 0, 0.00258254, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.000964788, 0, 0, 0, -0.000386357, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0.000964788, 0, 0, 0, 0.00240334, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.00150191, 0, 0, 0, -0.00153947, 0, 0, 0, -0.000142264, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0.00153947, 0, 0, 0, -0.0031273, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, -0.00240334, 0, 0, 0, 0.00119144, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0.000386357, 0, 0, 0, -0.00119144, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, -0.00258254, 0, 0, 0, 0.0031273, 0, 0, 0, -0.00104397, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0.000142264, 0, 0, 0, 0.00104397, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0;
Eigen::MatrixXd dip2_ref = Eigen::MatrixXd::Zero(dftbasissize, dftbasissize);
dip2_ref << 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.000135096, 0, 0, 0,
0.000379317, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.000159912, 0, 0,
0, 0.000986306, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.000135096, 0, 0, 0,
0.000816107, 0, 0, 0, -0.00309704, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.000159912, 0, 0, 0, -0.000212013, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, -0.000816107, 0, 0, 0, 0.000133185, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0.000212013, 0, 0, 0, -0.0026893, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0.000379317, 0, 0, 0, -0.000133185, 0, 0, 0, 0.00342749, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, -0.000986306, 0, 0, 0, 0.0026893, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0.00309704, 0, 0, 0, -0.00342749, 0, 0, 0, 0,
-0.000135096, 0, 0, 0, -0.000379317, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0.000159912, 0, 0, 0, -0.000986306, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.000135096, 0, 0, 0, -0.000816107, 0, 0, 0, 0.00309704, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, -0.000159912, 0, 0, 0, 0.000212013, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0.000816107, 0, 0, 0, -0.000133185, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, -0.000212013, 0, 0, 0, 0.0026893, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0.000379317, 0, 0, 0, 0.000133185, 0, 0, 0, -0.00342749, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0.000986306, 0, 0, 0, -0.0026893, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, -0.00309704, 0, 0, 0, 0.00342749, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0;
bool check_dip0 = dip.Matrix()[0].isApprox(dip0_ref, 0.0001);
BOOST_CHECK_EQUAL(check_dip0, 1);
if (!check_dip0) {
cout << "ref" << endl;
cout << dip0_ref << endl;
cout << "result" << endl;
cout << dip.Matrix()[0] << endl;
}
bool check_dip1 = dip.Matrix()[1].isApprox(dip1_ref, 0.0001);
BOOST_CHECK_EQUAL(check_dip1, 1);
if (!check_dip1) {
cout << "ref" << endl;
cout << dip1_ref << endl;
cout << "result" << endl;
cout << dip.Matrix()[1] << endl;
}
bool check_dip2 = dip.Matrix()[2].isApprox(dip2_ref, 0.0001);
BOOST_CHECK_EQUAL(check_dip2, 1);
if (!check_dip2) {
cout << "ref" << endl;
cout << dip2_ref << endl;
cout << "result" << endl;
cout << dip.Matrix()[2] << endl;
}
AOMomentum momentum;
momentum.Fill(dftbasis);
Eigen::MatrixXd momentum0_ref =
Eigen::MatrixXd::Zero(dftbasissize, dftbasissize);
momentum0_ref << 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0452282, 0, 0, 0, -0.0551799, 0,
0, 0, 0.0280469, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00821239, 0, 0, 0,
-0.0201019, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0955273, 0, 0, 0,
0.07533, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0290925, 0, 0, 0,
0.05022, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0551799, 0, 0, 0, 0.0727415, 0,
0, 0, -0.0552111, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0201019, 0, 0, 0,
0.0537993, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.07533, 0, 0, 0,
-0.076546, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.05022, 0, 0, 0,
-0.142981, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0280469, 0, 0, 0, -0.0552111,
0, 0, 0, 0.126837, -0.0452282, 0, 0, 0, 0.0551799, 0, 0, 0, -0.0280469, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, -0.00821239, 0, 0, 0, 0.0201019, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0955273, 0, 0, 0, -0.07533, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0.0290925, 0, 0, 0, -0.05022, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0.0551799, 0, 0, 0, -0.0727415, 0, 0, 0, 0.0552111, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0.0201019, 0, 0, 0, -0.0537993, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, -0.07533, 0, 0, 0, 0.076546, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, -0.05022, 0, 0, 0, 0.142981, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0.0280469, 0, 0, 0, 0.0552111, 0, 0, 0, -0.126837, 0, 0, 0, 0, 0, 0, 0,
0, 0;
Eigen::MatrixXd momentum1_ref =
Eigen::MatrixXd::Zero(dftbasissize, dftbasissize);
momentum1_ref << 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0100472, 0, 0, 0,
-0.0145434, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00821239, 0, 0, 0,
-0.01029, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00821239, 0, 0, 0,
-0.0201019, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0100472, 0, 0, 0,
-0.0180007, 0, 0, 0, 0.0293066, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0.0180007, 0, 0, 0, 0.0336053, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0.0201019, 0, 0, 0, 0.0290685, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0.01029, 0, 0, 0, 0.0290685, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-0.0145434, 0, 0, 0, 0.0336053, 0, 0, 0, -0.0942411, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0.0293066, 0, 0, 0, -0.0942411, 0, 0, 0, 0, -0.0100472, 0,
0, 0, 0.0145434, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.00821239, 0, 0, 0,
0.01029, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.00821239, 0, 0, 0,
0.0201019, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0100472, 0, 0, 0,
0.0180007, 0, 0, 0, -0.0293066, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.0180007, 0, 0, 0, -0.0336053, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.0201019, 0, 0, 0, -0.0290685, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.01029, 0, 0, 0, -0.0290685, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.0145434, 0, 0, 0, -0.0336053, 0, 0, 0, 0.0942411, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, -0.0293066, 0, 0, 0, 0.0942411, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0;
Eigen::MatrixXd momentum2_ref =
Eigen::MatrixXd::Zero(dftbasissize, dftbasissize);
momentum2_ref << 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0472832, 0, 0, 0,
-0.0343107, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00824288, 0, 0, 0,
-0.0107452, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0472832, 0, 0, 0, -0.0583714,
0, 0, 0, 0.0337403, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00824288, 0, 0, 0,
-0.0195091, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0583714, 0, 0, 0,
0.0464027, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0195091, 0, 0, 0,
0.0292982, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0343107, 0, 0, 0, 0.0464027,
0, 0, 0, -0.0373403, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0107452, 0, 0, 0,
0.0292982, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0337403, 0, 0, 0,
-0.0373403, 0, 0, 0, 0, -0.0472832, 0, 0, 0, 0.0343107, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, -0.00824288, 0, 0, 0, 0.0107452, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, -0.0472832, 0, 0, 0, 0.0583714, 0, 0, 0, -0.0337403, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, -0.00824288, 0, 0, 0, 0.0195091, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0.0583714, 0, 0, 0, -0.0464027, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0.0195091, 0, 0, 0, -0.0292982, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0.0343107, 0, 0, 0, -0.0464027, 0, 0, 0, 0.0373403, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0.0107452, 0, 0, 0, -0.0292982, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, -0.0337403, 0, 0, 0, 0.0373403, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0;
bool check_momentum0 = momentum.Matrix()[0].isApprox(momentum0_ref, 0.0001);
BOOST_CHECK_EQUAL(check_momentum0, 1);
if (!check_momentum0) {
cout << "ref" << endl;
cout << momentum0_ref << endl;
cout << "result" << endl;
cout << momentum.Matrix()[0] << endl;
}
bool check_momentum1 = momentum.Matrix()[1].isApprox(momentum1_ref, 0.0001);
BOOST_CHECK_EQUAL(check_momentum1, 1);
if (!check_momentum1) {
cout << "ref" << endl;
cout << momentum1_ref << endl;
cout << "result" << endl;
cout << momentum.Matrix()[1] << endl;
}
bool check_momentum2 = momentum.Matrix()[2].isApprox(momentum2_ref, 0.0001);
BOOST_CHECK_EQUAL(check_momentum2, 1);
if (!check_momentum2) {
cout << "ref" << endl;
cout << momentum2_ref << endl;
cout << "result" << endl;
cout << momentum.Matrix()[2] << endl;
}
}
BOOST_AUTO_TEST_SUITE_END()