1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: test_molden.cc
4 // Copyright (C) 2012 Toru Shiozaki
5 //
6 // Author: Shane Parker < shane.parker@u.northwestern.edu >
7 // Maintainer: Shiozaki group
8 //
9 // This file is part of the BAGEL package.
10 //
11 // This program is free software: you can redistribute it and/or modify
12 // it under the terms of the GNU General Public License as published by
13 // the Free Software Foundation, either version 3 of the License, or
14 // (at your option) any later version.
15 //
16 // This program is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU General Public License
22 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
23 //
24 
25 
26 #include <src/util/io/moldenout.h>
27 #include <src/util/io/moldenin.h>
28 #include <src/scf/hf/fock.h>
29 
molden_out_energy(std::string inp1,std::string inp2)30 double molden_out_energy(std::string inp1, std::string inp2) {
31 
32   auto ofs = std::make_shared<std::ofstream>(inp1 + ".testout", std::ios::trunc);
33   std::streambuf* backup_stream = std::cout.rdbuf(ofs->rdbuf());
34 
35   {
36     std::stringstream ss; ss << location__ << inp1 << ".json";
37     auto idata = std::make_shared<const PTree>(ss.str());
38     auto keys = idata->get_child("bagel");
39     std::shared_ptr<Geometry> geom;
40 
41     std::shared_ptr<const Reference> ref;
42 
43     for (auto& itree : *keys) {
44       const std::string method = to_lower(itree->get<std::string>("title", ""));
45 
46       if (method == "molecule") {
47         geom = std::make_shared<Geometry>(itree);
48 
49       } else if (method == "hf") {
50         auto scf = std::make_shared<RHF>(itree, geom);
51         scf->compute();
52         ref = scf->conv_to_ref();
53       }
54       else if (method == "print") {
55         std::string out_file = itree->get<std::string>("file", inp1 + ".molden");
56 
57         MoldenOut mfs(out_file);
58         mfs << geom;
59         mfs << ref;
60 
61       }
62     }
63   }
64 
65   double energy = 0.0;
66   {
67     std::stringstream ss; ss << location__ << inp2 << ".json";
68     auto idata = std::make_shared<const PTree>(ss.str());
69     std::shared_ptr<const PTree> keys = idata->get_child("bagel");
70     std::shared_ptr<const PTree> mol = *keys->begin();
71 
72     const std::string method = to_lower(mol->get<std::string>("title", ""));
73     if (method != "molecule") throw std::logic_error("broken test case");
74     auto geom = std::make_shared<Geometry>(mol);
75 
76     auto coeff = std::make_shared<Coeff>(geom);
77     auto tmp = std::make_shared<VectorB>(geom->nbasis());
78     auto tmp2 = std::make_shared<VectorB>(geom->nbasis());
79 
80     std::string filename = inp1 + ".molden";
81     MoldenIn mfs(filename, geom->spherical());
82     mfs.read();
83     mfs >> tie(coeff, geom, tmp, tmp2);
84 
85     std::shared_ptr<const Matrix> ao_density = coeff->form_density_rhf(geom->nele()/2);
86     auto hcore = std::make_shared<const Hcore>(geom);
87     auto fock = std::make_shared<const Fock<1>>(geom, hcore, ao_density, geom->schwarz());
88 
89     auto hcore_fock = std::make_shared<const Matrix>(*hcore + *fock);
90     energy = ((*ao_density)*(*hcore_fock)).trace();
91     energy = 0.5*energy + geom->nuclear_repulsion();
92 
93   }
94   std::cout.rdbuf(backup_stream);
95   return energy;
96 }
97 
98 BOOST_AUTO_TEST_SUITE(TEST_MOLDEN)
99 
BOOST_AUTO_TEST_CASE(MOLDEN)100 BOOST_AUTO_TEST_CASE(MOLDEN) {
101     BOOST_CHECK(compare(molden_out_energy("hf_write_mol_sph", "hf_read_mol_sph"),   -99.84772354 ));
102     BOOST_CHECK(compare(molden_out_energy("hf_write_mol_cart", "hf_read_mol_cart"), -99.84911270 ));
103 }
104 
105 BOOST_AUTO_TEST_SUITE_END()
106