/*
* Copyright 2009-2018 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/orbitals.h>
#include <votca/xtp/aomatrix.h>
using namespace votca::xtp;
using namespace votca;
using namespace std;
BOOST_AUTO_TEST_SUITE(aomatrix_test)
BOOST_AUTO_TEST_CASE(aomatrices_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().LoadFromXYZ("molecule.xyz");
BasisSet basis;
basis.LoadBasisSet("3-21G.xml");
AOBasis aobasis;
aobasis.AOBasisFill(basis,orbitals.QMAtoms());
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);
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);
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);
if(!check_coulomb){
cout<<"ref"<<endl;
cout<<coulomb_ref<<endl;
cout<<"result"<<endl;
cout<<coulomb.Matrix()<<endl;
}
BOOST_CHECK_EQUAL(check_coulomb, 1);
coulomb.Pseudo_InvSqrt(1e-7);
Eigen::MatrixXd coulombinvsqrt_ref=Eigen::MatrixXd::Zero(17,17);
coulombinvsqrt_ref<<1.6659153595,4.0151975735,0,0,0,5.9651489891,0,0,0,1.8658410271,4.8359823061,1.8658410271,4.8359823061,1.8658410271,4.8359823061,1.8658410271,4.8359823061,
4.0151975735,18.378952637,0,0,0,31.374258628,0,0,0,10.303062989,26.609117042,10.303062989,26.609117042,10.303062989,26.609117042,10.303062989,26.609117042,
0,0,4.7594244166,0,0,0,7.0293902091,0,0,1.9658834037,2.4872154307,1.9658834037,2.4872154307,-1.9658834037,-2.4872154307,-1.9658834037,-2.4872154307,
0,0,0,4.7594244166,0,0,0,7.0293902091,0,1.9658834037,2.4872154307,-1.9658834037,-2.4872154307,-1.9658834037,-2.4872154307,1.9658834037,2.4872154307,
0,0,0,0,4.7594244166,0,0,0,7.0293902091,1.9658834037,2.4872154307,-1.9658834037,-2.4872154307,1.9658834037,2.4872154307,-1.9658834037,-2.4872154307,
5.9651489891,31.374258628,0,0,0,64.160947091,0,0,0,21.367082884,58.423867269,21.367082884,58.423867269,21.367082884,58.423867269,21.367082884,58.423867269,
0,0,7.0293902091,0,0,0,21.386982364,0,0,5.0791868054,8.8864083237,5.0791868054,8.8864083237,-5.0791868054,-8.8864083237,-5.0791868054,-8.8864083237,
0,0,0,7.0293902091,0,0,0,21.386982364,0,5.0791868054,8.8864083237,-5.0791868054,-8.8864083237,-5.0791868054,-8.8864083237,5.0791868054,8.8864083237,
0,0,0,0,7.0293902091,0,0,0,21.386982364,5.0791868054,8.8864083237,-5.0791868054,-8.8864083237,5.0791868054,8.8864083237,-5.0791868054,-8.8864083237,
1.8658410271,10.303062989,1.9658834037,1.9658834037,1.9658834037,21.367082884,5.0791868054,5.0791868054,5.0791868054,13.908456951,26.861509511,5.5466937828,17.041106037,5.5466937828,17.041106037,5.5466937828,17.041106037,
4.8359823061,26.609117042,2.4872154307,2.4872154307,2.4872154307,58.423867269,8.8864083237,8.8864083237,8.8864083237,26.861509511,68.596721551,17.041106037,50.770206191,17.041106037,50.770206191,17.041106037,50.770206191,
1.8658410271,10.303062989,1.9658834037,-1.9658834037,-1.9658834037,21.367082884,5.0791868054,-5.0791868054,-5.0791868054,5.5466937828,17.041106037,13.908456951,26.861509511,5.5466937828,17.041106037,5.5466937828,17.041106037,
4.8359823061,26.609117042,2.4872154307,-2.4872154307,-2.4872154307,58.423867269,8.8864083237,-8.8864083237,-8.8864083237,17.041106037,50.770206191,26.861509511,68.596721551,17.041106037,50.770206191,17.041106037,50.770206191,
1.8658410271,10.303062989,-1.9658834037,-1.9658834037,1.9658834037,21.367082884,-5.0791868054,-5.0791868054,5.0791868054,5.5466937828,17.041106037,5.5466937828,17.041106037,13.908456951,26.861509511,5.5466937828,17.041106037,
4.8359823061,26.609117042,-2.4872154307,-2.4872154307,2.4872154307,58.423867269,-8.8864083237,-8.8864083237,8.8864083237,17.041106037,50.770206191,17.041106037,50.770206191,26.861509511,68.596721551,17.041106037,50.770206191,
1.8658410271,10.303062989,-1.9658834037,1.9658834037,-1.9658834037,21.367082884,-5.0791868054,5.0791868054,-5.0791868054,5.5466937828,17.041106037,5.5466937828,17.041106037,5.5466937828,17.041106037,13.908456951,26.861509511,
4.8359823061,26.609117042,-2.4872154307,2.4872154307,-2.4872154307,58.423867269,-8.8864083237,8.8864083237,-8.8864083237,17.041106037,50.770206191,17.041106037,50.770206191,17.041106037,50.770206191,26.861509511,68.596721551;
bool check_coulombinvsqrt = coulomb.Matrix().isApprox(coulombinvsqrt_ref,0.00001);
if(!check_coulombinvsqrt){
cout<<"ref"<<endl;
cout<<coulombinvsqrt_ref<<endl;
cout<<"result"<<endl;
cout<<coulomb.Matrix()<<endl;
}
BOOST_CHECK_EQUAL(check_coulombinvsqrt, 1);
coulomb.Fill(aobasis);
coulomb.Pseudo_InvSqrt_GWBSE(overlap,1e-7);
Eigen::MatrixXd coulombinvsqrtgw_ref=Eigen::MatrixXd::Zero(17,17);
coulombinvsqrtgw_ref<<1.6659153595,4.0151975735,0,0,0,5.9651489891,0,0,0,1.8658410271,4.8359823061,1.8658410271,4.8359823061,1.8658410271,4.8359823061,1.8658410271,4.8359823061,
4.0151975735,18.378952637,0,0,0,31.374258628,0,0,0,10.303062989,26.609117042,10.303062989,26.609117042,10.303062989,26.609117042,10.303062989,26.609117042,
0,0,4.7594244166,0,0,0,7.0293902091,0,0,1.9658834037,2.4872154307,1.9658834037,2.4872154307,-1.9658834037,-2.4872154307,-1.9658834037,-2.4872154307,
0,0,0,4.7594244166,0,0,0,7.0293902091,0,1.9658834037,2.4872154307,-1.9658834037,-2.4872154307,-1.9658834037,-2.4872154307,1.9658834037,2.4872154307,
0,0,0,0,4.7594244166,0,0,0,7.0293902091,1.9658834037,2.4872154307,-1.9658834037,-2.4872154307,1.9658834037,2.4872154307,-1.9658834037,-2.4872154307,
5.9651489891,31.374258628,0,0,0,64.160947091,0,0,0,21.367082884,58.423867269,21.367082884,58.423867269,21.367082884,58.423867269,21.367082884,58.423867269,
0,0,7.0293902091,0,0,0,21.386982364,0,0,5.0791868054,8.8864083237,5.0791868054,8.8864083237,-5.0791868054,-8.8864083237,-5.0791868054,-8.8864083237,
0,0,0,7.0293902091,0,0,0,21.386982364,0,5.0791868054,8.8864083237,-5.0791868054,-8.8864083237,-5.0791868054,-8.8864083237,5.0791868054,8.8864083237,
0,0,0,0,7.0293902091,0,0,0,21.386982364,5.0791868054,8.8864083237,-5.0791868054,-8.8864083237,5.0791868054,8.8864083237,-5.0791868054,-8.8864083237,
1.8658410271,10.303062989,1.9658834037,1.9658834037,1.9658834037,21.367082884,5.0791868054,5.0791868054,5.0791868054,13.908456951,26.861509511,5.5466937828,17.041106037,5.5466937828,17.041106037,5.5466937828,17.041106037,
4.8359823061,26.609117042,2.4872154307,2.4872154307,2.4872154307,58.423867269,8.8864083237,8.8864083237,8.8864083237,26.861509511,68.596721551,17.041106037,50.770206191,17.041106037,50.770206191,17.041106037,50.770206191,
1.8658410271,10.303062989,1.9658834037,-1.9658834037,-1.9658834037,21.367082884,5.0791868054,-5.0791868054,-5.0791868054,5.5466937828,17.041106037,13.908456951,26.861509511,5.5466937828,17.041106037,5.5466937828,17.041106037,
4.8359823061,26.609117042,2.4872154307,-2.4872154307,-2.4872154307,58.423867269,8.8864083237,-8.8864083237,-8.8864083237,17.041106037,50.770206191,26.861509511,68.596721551,17.041106037,50.770206191,17.041106037,50.770206191,
1.8658410271,10.303062989,-1.9658834037,-1.9658834037,1.9658834037,21.367082884,-5.0791868054,-5.0791868054,5.0791868054,5.5466937828,17.041106037,5.5466937828,17.041106037,13.908456951,26.861509511,5.5466937828,17.041106037,
4.8359823061,26.609117042,-2.4872154307,-2.4872154307,2.4872154307,58.423867269,-8.8864083237,-8.8864083237,8.8864083237,17.041106037,50.770206191,17.041106037,50.770206191,26.861509511,68.596721551,17.041106037,50.770206191,
1.8658410271,10.303062989,-1.9658834037,1.9658834037,-1.9658834037,21.367082884,-5.0791868054,5.0791868054,-5.0791868054,5.5466937828,17.041106037,5.5466937828,17.041106037,5.5466937828,17.041106037,13.908456951,26.861509511,
4.8359823061,26.609117042,-2.4872154307,2.4872154307,-2.4872154307,58.423867269,-8.8864083237,8.8864083237,-8.8864083237,17.041106037,50.770206191,17.041106037,50.770206191,17.041106037,50.770206191,26.861509511,68.596721551;
bool check_coulombinvsqrtgw = coulomb.Matrix().isApprox(coulombinvsqrtgw_ref,0.00001);
if(!check_coulombinvsqrtgw){
cout<<"ref"<<endl;
cout<<coulombinvsqrtgw_ref<<endl;
cout<<"result"<<endl;
cout<<coulomb.Matrix()<<endl;
}
BOOST_CHECK_EQUAL(check_coulombinvsqrtgw, 1);
AOESP esp;
esp.Fillnucpotential(aobasis,orbitals.QMAtoms());
Eigen::MatrixXd esp_ref=Eigen::MatrixXd::Zero(17,17);
esp_ref<<0.485633,0.0929737,-0.0112854,0.0112854,-0.0112854,0.0875665,-0.00167437,0.00167437,-0.00167437,0.00881413,0.0388865,0.00881413,0.0388865,0.00881413,0.0388865,0.010454,0.0404184,
0.0929737,0.485422,-0.0825499,0.0825499,-0.0825499,0.367955,-0.0505522,0.0505522,-0.0505522,0.0731397,0.178738,0.0731397,0.178738,0.0731397,0.178738,0.194659,0.240204,
-0.0112854,-0.0825499,0.485023,-0.0299598,0.0299598,-0.0781392,0.252799,-0.0256949,0.0256949,0.0563721,0.0217104,0.0563721,0.0217104,-0.0677234,-0.0882558,-0.197548,-0.133855,
0.0112854,0.0825499,-0.0299598,0.485023,-0.0299598,0.0781392,-0.0256949,0.252799,-0.0256949,0.0677234,0.0882558,-0.0563721,-0.0217104,-0.0563721,-0.0217104,0.197548,0.133855,
-0.0112854,-0.0825499,0.0299598,-0.0299598,0.485023,-0.0781392,0.0256949,-0.0256949,0.252799,0.0563721,0.0217104,-0.0677234,-0.0882558,0.0563721,0.0217104,-0.197548,-0.133855,
0.0875665,0.367955,-0.0781392,0.0781392,-0.0781392,0.452433,-0.100822,0.100822,-0.100822,0.110741,0.253817,0.110741,0.253817,0.110741,0.253817,0.391022,0.405943,
-0.00167437,-0.0505522,0.252799,-0.0256949,0.0256949,-0.100822,0.407715,-0.0510951,0.0510951,0.084227,0.0617612,0.084227,0.0617612,-0.0960714,-0.152532,-0.362371,-0.271961,
0.00167437,0.0505522,-0.0256949,0.252799,-0.0256949,0.100822,-0.0510951,0.407715,-0.0510951,0.0960714,0.152532,-0.084227,-0.0617612,-0.084227,-0.0617612,0.362371,0.271961,
-0.00167437,-0.0505522,0.0256949,-0.0256949,0.252799,-0.100822,0.0510951,-0.0510951,0.407715,0.084227,0.0617612,-0.0960714,-0.152532,0.084227,0.0617612,-0.362371,-0.271961,
0.00881413,0.0731397,0.0563721,0.0677234,0.0563721,0.110741,0.084227,0.0960714,0.084227,0.297388,0.192083,0.0026714,0.0376198,0.0026714,0.0376198,0.00464633,0.042269,
0.0388865,0.178738,0.0217104,0.0882558,0.0217104,0.253817,0.0617612,0.152532,0.0617612,0.192083,0.2962,0.0376198,0.120354,0.0376198,0.120354,0.124795,0.179447,
0.00881413,0.0731397,0.0563721,-0.0563721,-0.0677234,0.110741,0.084227,-0.084227,-0.0960714,0.0026714,0.0376198,0.297388,0.192083,0.0026714,0.0376198,0.00464633,0.042269,
0.0388865,0.178738,0.0217104,-0.0217104,-0.0882558,0.253817,0.0617612,-0.0617612,-0.152532,0.0376198,0.120354,0.192083,0.2962,0.0376198,0.120354,0.124795,0.179447,
0.00881413,0.0731397,-0.0677234,-0.0563721,0.0563721,0.110741,-0.0960714,-0.084227,0.084227,0.0026714,0.0376198,0.0026714,0.0376198,0.297388,0.192083,0.00464633,0.042269,
0.0388865,0.178738,-0.0882558,-0.0217104,0.0217104,0.253817,-0.152532,-0.0617612,0.0617612,0.0376198,0.120354,0.0376198,0.120354,0.192083,0.2962,0.124795,0.179447,
0.010454,0.194659,-0.197548,0.197548,-0.197548,0.391022,-0.362371,0.362371,-0.362371,0.00464633,0.124795,0.00464633,0.124795,0.00464633,0.124795,1.72092,0.782663,
0.0404184,0.240204,-0.133855,0.133855,-0.133855,0.405943,-0.271961,0.271961,-0.271961,0.042269,0.179447,0.042269,0.179447,0.042269,0.179447,0.782663,0.683004;
bool check_esp = esp.Matrix().isApprox(esp_ref,0.00001);
BOOST_CHECK_EQUAL(check_esp, 1);
ofstream ecpfile("ecp.xml");
ecpfile << "<pseudopotential name=\"ECP_STUTTGART\">" << endl;
ecpfile << " <element name=\"C\" lmax=\"3\" ncore=\"2\">" << endl;
ecpfile << " <shell type=\"F\"><constant power=\"2\" decay=\"1.0\" contraction=\"0.0\"></constant></shell>" << endl;
ecpfile << " <shell type=\"S\"><constant power=\"2\" decay=\"6.40105200\" contraction=\"33.12163800\"></constant></shell>" << endl;
ecpfile << " <shell type=\"P\"><constant power=\"2\" decay=\"7.30774700\" contraction=\"-1.98625700\"></constant></shell>" << endl;
ecpfile << " <shell type=\"D\"><constant power=\"2\" decay=\"5.96179600\" contraction=\"-9.45431800\"></constant></shell>" << endl;
ecpfile << " </element>" << endl;
ecpfile << "</pseudopotential>" << endl;
BasisSet ecps;
ecps.LoadPseudopotentialSet("ecp.xml");
AOBasis ecpbasis;
ecpbasis.ECPFill(ecps,orbitals.QMAtoms());
AOECP ecp;
ecp.setECP(&ecpbasis);
ecp.Fill(aobasis);
Eigen::MatrixXd ecp_ref= Eigen::MatrixXd::Zero(17,17);
ecp_ref<<21.6188,1.34835,0,0,0,2.29744,0,0,0,0.209711,1.01592,0.209711,1.01592,0.209711,1.01592,0.209711,1.01592,
1.34835,0.702249,0,0,0,0.4993,0,0,0,0.0564639,0.225665,0.0564639,0.225665,0.0564639,0.225665,0.0564639,0.225665,
0,0,-0.0737545,0,0,0,-0.00882987,0,0,-0.00178626,-0.00193605,-0.00178626,-0.00193605,0.00178626,0.00193605,0.00178626,0.00193605,
0,0,0,-0.0737545,0,0,0,-0.00882987,0,-0.00178626,-0.00193605,0.00178626,0.00193605,0.00178626,0.00193605,-0.00178626,-0.00193605,
0,0,0,0,-0.0737545,0,0,0,-0.00882987,-0.00178626,-0.00193605,0.00178626,0.00193605,-0.00178626,-0.00193605,0.00178626,0.00193605,
2.29744,0.4993,0,0,0,0.458665,0,0,0,0.0477375,0.20545,0.0477375,0.20545,0.0477375,0.20545,0.0477375,0.20545,
0,0,-0.00882987,0,0,0,-0.0011596,0,0,-0.000240513,-0.000255319,-0.000240513,-0.000255319,0.000240513,0.000255319,0.000240513,0.000255319,
0,0,0,-0.00882987,0,0,0,-0.0011596,0,-0.000240513,-0.000255319,0.000240513,0.000255319,0.000240513,0.000255319,-0.000240513,-0.000255319,
0,0,0,0,-0.00882987,0,0,0,-0.0011596,-0.000240513,-0.000255319,0.000240513,0.000255319,-0.000240513,-0.000255319,0.000240513,0.000255319,
0.209711,0.0564639,-0.00178626,-0.00178626,-0.00178626,0.0477375,-0.000240513,-0.000240513,-0.000240513,0.00468574,0.0212243,0.0052935,0.0215396,0.0052935,0.0215396,0.0052935,0.0215396,
1.01592,0.225665,-0.00193605,-0.00193605,-0.00193605,0.20545,-0.000255319,-0.000255319,-0.000255319,0.0212243,0.0918741,0.0215396,0.0921252,0.0215396,0.0921252,0.0215396,0.0921252,
0.209711,0.0564639,-0.00178626,0.00178626,0.00178626,0.0477375,-0.000240513,0.000240513,0.000240513,0.0052935,0.0215396,0.00468574,0.0212243,0.0052935,0.0215396,0.0052935,0.0215396,
1.01592,0.225665,-0.00193605,0.00193605,0.00193605,0.20545,-0.000255319,0.000255319,0.000255319,0.0215396,0.0921252,0.0212243,0.0918741,0.0215396,0.0921252,0.0215396,0.0921252,
0.209711,0.0564639,0.00178626,0.00178626,-0.00178626,0.0477375,0.000240513,0.000240513,-0.000240513,0.0052935,0.0215396,0.0052935,0.0215396,0.00468574,0.0212243,0.0052935,0.0215396,
1.01592,0.225665,0.00193605,0.00193605,-0.00193605,0.20545,0.000255319,0.000255319,-0.000255319,0.0215396,0.0921252,0.0215396,0.0921252,0.0212243,0.0918741,0.0215396,0.0921252,
0.209711,0.0564639,0.00178626,-0.00178626,0.00178626,0.0477375,0.000240513,-0.000240513,0.000240513,0.0052935,0.0215396,0.0052935,0.0215396,0.0052935,0.0215396,0.00468574,0.0212243,
1.01592,0.225665,0.00193605,-0.00193605,0.00193605,0.20545,0.000255319,-0.000255319,0.000255319,0.0215396,0.0921252,0.0215396,0.0921252,0.0215396,0.0921252,0.0212243,0.0918741;
bool check_ecp = ecp.Matrix().isApprox(ecp_ref,0.00001);
BOOST_CHECK_EQUAL(check_ecp, 1);
}
BOOST_AUTO_TEST_CASE(externalmatrices_test){
Orbitals orbitals;
orbitals.QMAtoms().LoadFromXYZ("molecule.xyz");
BasisSet basis;
basis.LoadBasisSet("3-21G.xml");
AOBasis aobasis;
aobasis.AOBasisFill(basis,orbitals.QMAtoms());
ofstream mpsfile("polarsite.mps");
mpsfile<<"! One Site"<<endl;
mpsfile<<"! N=1 "<<endl;
mpsfile<<"Units angstrom"<<endl;
mpsfile<<" C +0 0 3 Rank 2"<<endl;
mpsfile<<"+0"<<endl;
mpsfile<<"10 0 0"<<endl;
mpsfile<<" 100 0 0 0 0"<<endl;
mpsfile<<"P +1.9445387 +0.0000000 +0.0000000 +1.9445387 +0.0000000 +1.9445387 "<<endl;
auto polar_segments = std::make_shared<MMRegion>();
PolarSegment seg=PolarSegment("",0);
seg.LoadFromMPS("polarsite.mps");
polar_segments->push_back(seg);
AODipole_Potential dip;
dip.Fillextpotential(aobasis,polar_segments);
Eigen::MatrixXd dip_ref= Eigen::MatrixXd::Zero(17,17);
dip_ref<<0.31114997753,0.059568868026,0.0090978711864,0,0,0.056104697636,0.0013498178976,0,0,0.0061933281198,0.025459181656,0.0061933281198,0.025459181656,0.0056130806569,0.024860733171,0.0056130806569,0.024860733171,
0.059568868026,0.31114997753,0.066842196408,0,0,0.2368963398,0.042609848798,0,0,0.073695437658,0.13588175121,0.073695437658,0.13588175121,0.047134059517,0.11392427425,0.047134059517,0.11392427425,
0.0090978711864,0.066842196408,0.32666220712,0,0,0.065599720802,0.17980265473,0,0,0.075224551189,0.083547839473,0.075224551189,0.083547839473,-0.035337083996,-0.0087046753266,-0.035337083996,-0.0087046753266,
0,0,0,0.30339386273,0,0,0,0.15697695635,0,0.061346561258,0.043126415371,-0.061346561258,-0.043126415371,-0.040820728861,-0.037257483146,0.040820728861,0.037257483146,
0,0,0,0,0.30339386273,0,0,0,0.15697695635,0.061346561258,0.043126415371,-0.061346561258,-0.043126415371,0.040820728861,0.037257483146,-0.040820728861,-0.037257483146,
0.056104697636,0.2368963398,0.065599720802,0,0,0.31114557005,0.12399363244,0,0,0.13566608278,0.24813270723,0.13566608278,0.24813270723,0.072222745084,0.16730389611,0.072222745084,0.16730389611,
0.0013498178976,0.042609848798,0.17980265473,0,0,0.12399363244,0.38517412728,0,0,0.13766652824,0.2355220863,0.13766652824,0.2355220863,-0.05335528609,-0.024078433641,-0.05335528609,-0.024078433641,
0,0,0,0.15697695635,0,0,0,0.27407784996,0,0.10959962132,0.10744652669,-0.10959962132,-0.10744652669,-0.060016250362,-0.076591108649,0.060016250362,0.076591108649,
0,0,0,0,0.15697695635,0,0,0,0.27407784996,0.10959962132,0.10744652669,-0.10959962132,-0.10744652669,0.060016250362,0.076591108649,-0.060016250362,-0.076591108649,
0.0061933281198,0.073695437658,0.075224551189,0.061346561258,0.061346561258,0.13566608278,0.13766652824,0.10959962132,0.10959962132,0.40885081105,0.26407630555,0.0038749978449,0.053453882722,0.002270661607,0.043208755867,0.002270661607,0.043208755867,
0.025459181656,0.13588175121,0.083547839473,0.043126415371,0.043126415371,0.24813270723,0.2355220863,0.10744652669,0.10744652669,0.26407630555,0.40853020343,0.053453882722,0.17647932404,0.026292260419,0.10354578683,0.026292260419,0.10354578683,
0.0061933281198,0.073695437658,0.075224551189,-0.061346561258,-0.061346561258,0.13566608278,0.13766652824,-0.10959962132,-0.10959962132,0.0038749978449,0.053453882722,0.40885081105,0.26407630555,0.002270661607,0.043208755867,0.002270661607,0.043208755867,
0.025459181656,0.13588175121,0.083547839473,-0.043126415371,-0.043126415371,0.24813270723,0.2355220863,-0.10744652669,-0.10744652669,0.053453882722,0.17647932404,0.26407630555,0.40853020343,0.026292260419,0.10354578683,0.026292260419,0.10354578683,
0.0056130806569,0.047134059517,-0.035337083996,-0.040820728861,0.040820728861,0.072222745084,-0.05335528609,-0.060016250362,0.060016250362,0.002270661607,0.026292260419,0.002270661607,0.026292260419,0.19479973371,0.12582094161,0.001654409346,0.023957628139,
0.024860733171,0.11392427425,-0.0087046753266,-0.037257483146,0.037257483146,0.16730389611,-0.024078433641,-0.076591108649,0.076591108649,0.043208755867,0.10354578683,0.043208755867,0.10354578683,0.12582094161,0.19479972247,0.023957628139,0.075477416664,
0.0056130806569,0.047134059517,-0.035337083996,0.040820728861,-0.040820728861,0.072222745084,-0.05335528609,0.060016250362,-0.060016250362,0.002270661607,0.026292260419,0.002270661607,0.026292260419,0.001654409346,0.023957628139,0.19479973371,0.12582094161,
0.024860733171,0.11392427425,-0.0087046753266,0.037257483146,-0.037257483146,0.16730389611,-0.024078433641,0.076591108649,-0.076591108649,0.043208755867,0.10354578683,0.043208755867,0.10354578683,0.023957628139,0.075477416664,0.12582094161,0.19479972247;
bool dip_check=dip_ref.isApprox(dip.Matrix(),1e-4);
BOOST_CHECK_EQUAL(dip_check, 1);
if(!dip_check){
std::cout<<"dip Ref"<<endl;
std::cout<<dip_ref<<endl;
std::cout<<"Dip"<<endl;
std::cout<<dip.Matrix()<<endl;
}
AOQuadrupole_Potential quad;
quad.Fillextpotential(aobasis,polar_segments);
Eigen::MatrixXd quad_ref= Eigen::MatrixXd::Zero(17,17);
quad_ref<<-0.54885754461,-0.10507737426,-0.024072484017,0,0,-0.098966700337,-0.0035715464751,0,0,-0.011179037808,-0.045172910673,-0.011179037808,-0.045172910673,-0.0096466958333,-0.043589632176,-0.0096466958333,-0.043589632176,
-0.10507737426,-0.54885754461,-0.17686090204,0,0,-0.41787675648,-0.11274339711,0,0,-0.13976724915,-0.24922034653,-0.13976724915,-0.24922034653,-0.072122590024,-0.19138301686,-0.072122590024,-0.19138301686,
-0.024072484017,-0.17686090204,-0.60358359226,0,0,-0.17357337756,-0.34400799982,0,0,-0.15272516987,-0.18718118563,-0.15272516987,-0.18718118563,0.049819674744,-0.010744186374,0.049819674744,-0.010744186374,
0,0,0,-0.52149452079,0,0,0,-0.26348055222,0,-0.11179427451,-0.074503486364,0.11179427451,0.074503486364,0.061741084451,0.060158153131,-0.061741084451,-0.060158153131,
0,0,0,0,-0.52149452079,0,0,0,-0.26348055222,-0.11179427451,-0.074503486364,0.11179427451,0.074503486364,-0.061741084451,-0.060158153131,0.061741084451,0.060158153131,
-0.098966700337,-0.41787675648,-0.17357337756,0,0,-0.54878692407,-0.32778400061,0,0,-0.25462737209,-0.46846679465,-0.25462737209,-0.46846679465,-0.10214301985,-0.26126296427,-0.10214301985,-0.26126296427,
-0.0035715464751,-0.11274339711,-0.34400799982,0,0,-0.32778400061,-0.80862341594,0,0,-0.27550032049,-0.54212154917,-0.27550032049,-0.54212154917,0.069746857669,-0.008774945581,0.069746857669,-0.008774945581,
0,0,0,-0.26348055222,0,0,0,-0.41807747815,0,-0.19230399467,-0.16717978258,0.19230399467,0.16717978258,0.082086844704,0.10623084333,-0.082086844704,-0.10623084333,
0,0,0,0,-0.26348055222,0,0,0,-0.41807747815,-0.19230399467,-0.16717978258,0.19230399467,0.16717978258,-0.082086844704,-0.10623084333,0.082086844704,0.10623084333,
-0.011179037808,-0.13976724915,-0.15272516987,-0.11179427451,-0.11179427451,-0.25462737209,-0.27550032049,-0.19230399467,-0.19230399467,-0.74359826632,-0.48028932852,-0.008644060733,-0.10930647188,-0.0037525383686,-0.077959498143,-0.0037525383686,-0.077959498143,
-0.045172910673,-0.24922034653,-0.18718118563,-0.074503486364,-0.074503486364,-0.46846679465,-0.54212154917,-0.16717978258,-0.16717978258,-0.48028932852,-0.73992811949,-0.10930647188,-0.39011966636,-0.038016753712,-0.1710846496,-0.038016753712,-0.1710846496,
-0.011179037808,-0.13976724915,-0.15272516987,0.11179427451,0.11179427451,-0.25462737209,-0.27550032049,0.19230399467,0.19230399467,-0.008644060733,-0.10930647188,-0.74359826632,-0.48028932852,-0.0037525383686,-0.077959498143,-0.0037525383686,-0.077959498143,
-0.045172910673,-0.24922034653,-0.18718118563,0.074503486364,0.074503486364,-0.46846679465,-0.54212154917,0.16717978258,0.16717978258,-0.10930647188,-0.39011966636,-0.48028932852,-0.73992811949,-0.038016753712,-0.1710846496,-0.038016753712,-0.1710846496,
-0.0096466958333,-0.072122590024,0.049819674744,0.061741084451,-0.061741084451,-0.10214301985,0.069746857669,0.082086844704,-0.082086844704,-0.0037525383686,-0.038016753712,-0.0037525383686,-0.038016753712,-0.25989585987,-0.16786646053,-0.0024117505537,-0.033638947738,
-0.043589632176,-0.19138301686,-0.010744186374,0.060158153131,-0.060158153131,-0.26126296427,-0.008774945581,0.10623084333,-0.10623084333,-0.077959498143,-0.1710846496,-0.077959498143,-0.1710846496,-0.16786646053,-0.259895667,-0.033638947738,-0.11005918435,
-0.0096466958333,-0.072122590024,0.049819674744,-0.061741084451,0.061741084451,-0.10214301985,0.069746857669,-0.082086844704,0.082086844704,-0.0037525383686,-0.038016753712,-0.0037525383686,-0.038016753712,-0.0024117505537,-0.033638947738,-0.25989585987,-0.16786646053,
-0.043589632176,-0.19138301686,-0.010744186374,-0.060158153131,0.060158153131,-0.26126296427,-0.008774945581,-0.10623084333,0.10623084333,-0.077959498143,-0.1710846496,-0.077959498143,-0.1710846496,-0.033638947738,-0.11005918435,-0.16786646053,-0.259895667;
bool quad_check=quad_ref.isApprox(quad.Matrix(),1e-4);
BOOST_CHECK_EQUAL(quad_check, 1);
if(!quad_check){
std::cout<<"Quad Ref"<<endl;
std::cout<<quad_ref<<endl;
std::cout<<"Quad"<<endl;
std::cout<<quad.Matrix()<<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();
Orbitals orbitals;
orbitals.QMAtoms().LoadFromXYZ("C.xyz");
BasisSet basis;
basis.LoadBasisSet("contracted.xml");
AOBasis aobasis;
aobasis.AOBasisFill(basis,orbitals.QMAtoms());
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) {
Orbitals orbitals;
orbitals.QMAtoms().LoadFromXYZ("molecule.xyz");
BasisSet basis;
basis.LoadBasisSet("3-21G.xml");
AOBasis aobasis;
aobasis.AOBasisFill(basis,orbitals.QMAtoms());
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_SUITE_END()