Kindfield
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Properties Friends Pages
Functions
development.cpp File Reference
#include "electronsystems/gaussian/gaussiancore.h"
#include "electronsystems/gaussian/gaussiansystem.h"
#include "solvers/unrestrictedhartreefocksolver.h"
#include "solvers/restrictedhartreefocksolver.h"
#include "basisfunctions/gaussian/integrals/gaussianelectroninteractionintegral.h"
#include "basisfunctions/gaussian/integrals/gaussianoverlapintegral.h"
#include "math/boysfunction.h"
#include <iostream>
#include <iomanip>
#include <unittest++/UnitTest++.h>
Include dependency graph for development.cpp:

Go to the source code of this file.

Functions

 SUITE (Development)
 

Function Documentation

SUITE ( Development  )

Definition at line 17 of file development.cpp.

17  {
18  TEST(Dummy) {
19  }
20  TEST(LinearWater) {
21  vector<GaussianCore> cores;
22  string basis = "6-311++Gdsds";
23  cores.push_back(GaussianCore({0,0,0}, "atom_8_basis_" + basis + ".tm"));
24  cores.push_back(GaussianCore({1.8, 0, 0}, "atom_1_basis_" + basis + ".tm"));
25  cores.push_back(GaussianCore({-0.44916, 1.743056, 0.0}, "atom_1_basis_" + basis + ".tm"));
26  GaussianSystem system;
27  for(const GaussianCore &core : cores) {
28  system.addCore(core);
29  }
30  UnrestrictedHartreeFockSolver solver(&system);
31 // solver.setDiisEnabled(true);
32 // solver.setDiisSampleCount(10);
33 // solver.setDiisStartingIteration(20);
34  mat coefficientsUp = randn(system.nBasisFunctions(), system.nParticlesUp());
35  mat coefficientsDown = randn(system.nBasisFunctions(), system.nParticlesDown());
36  solver.setInitialCoefficientMatrices(coefficientsUp, coefficientsDown);
37  solver.setConvergenceTreshold(1e-9);
38  solver.setNIterationsMax(1e4);
39  solver.setDensityMixFactor(0.95);
40  solver.solve();
41 
42  cout << "Ground state energy: " << solver.energy() << endl;
43 
44  vector<GaussianCore> cores2;
45  cores2.push_back(GaussianCore({0,0,0}, "atom_8_basis_" + basis + ".tm"));
46  cores2.push_back(GaussianCore({8,0,0}, "atom_1_basis_" + basis + ".tm"));
47  cores2.push_back(GaussianCore({-7.999999999971834,2.122871834682184e-5,0}, "atom_1_basis_" + basis + ".tm"));
48  GaussianSystem system2;
49  for(const GaussianCore &core : cores2) {
50  system2.addCore(core);
51  }
52  UnrestrictedHartreeFockSolver solver2(&system2);
53  solver2.setDiisEnabled(true);
54  solver2.setDiisSampleCount(10);
55  solver2.setDiisStartingIteration(20);
56  solver2.setInitialCoefficientMatrices(solver.coeffcientMatrixUp(), solver.coeffcientMatrixDown());
57  solver2.setConvergenceTreshold(1e-9);
58  solver2.setNIterationsMax(1e3);
59  solver2.setDensityMixFactor(0.0);
60  solver2.solve();
61  cout << std::setprecision(20);
62  cout << "Energy: " << solver2.energy() << endl;
63  cout << "Iterations: " << solver2.iterationsUsed() << endl;
64  // CHECK_CLOSE(-149.5117583638509, solver.energy(), 1e-5);
65  }
66 // TEST(LinearHydrogen) {
67 // vector<GaussianCore> cores;
68 // cores.push_back(GaussianCore({0,0,0}, "atom_1_basis_6-31Gdsds.tm"));
69 // cores.push_back(GaussianCore({12,0,0}, "atom_1_basis_6-31Gdsds.tm"));
70 // cores.push_back(GaussianCore({24,0,0}, "atom_1_basis_6-31Gdsds.tm"));
71 // GaussianSystem system;
72 // for(const GaussianCore &core : cores) {
73 // system.addCore(core);
74 // }
75 // mat C;
76 // UnrestrictedHartreeFockSolver solver(&system);
77 // mat coefficientsUp = randn(system.nBasisFunctions(), system.nParticlesUp());
78 // mat coefficientsDown = randn(system.nBasisFunctions(), system.nParticlesDown());
79 // solver.setInitialCoefficientMatrices(coefficientsUp, coefficientsDown);
80 // solver.setConvergenceTreshold(1e-9);
81 // solver.setNIterationsMax(1e3);
82 // solver.setDensityMixFactor(0.95);
83 // solver.solve();
84 // cout << solver.overlapMatrix() << endl;
85 // cout << std::setprecision(20);
86 // cout << "Energy: " << solver.energy() << endl;
87 // cout << "Iterations: " << solver.iterationsUsed() << endl;
88 // // CHECK_CLOSE(-149.5117583638509, solver.energy(), 1e-5);
89 // }
90 // TEST(OxygenSix) {
91 // vector<GaussianCore> cores;
92 // cores.push_back(GaussianCore({0,0,0}, "atom_6_basis_6-311++Gdsds.tm"));
93 // GaussianSystem system;
94 // for(const GaussianCore &core : cores) {
95 // system.addCore(core);
96 // }
97 // mat C;
98 // UnrestrictedHartreeFockSolver solver(&system);
99 // solver.setConvergenceTreshold(1e-12);
100 // solver.setNIterationsMax(1e4);
101 // solver.setDensityMixFactor(0.95);
102 // solver.solve();
103 // cout << std::setprecision(20);
104 // cout << solver.energy() << endl;
105 // cout << solver.iterationsUsed() << endl;
106 // // CHECK_CLOSE(-149.5117583638509, solver.energy(), 1e-5);
107 // }
108 // TEST(Water431G) {
109 // vector<GaussianCore> cores;
112 // cores.push_back(GaussianCore({0,0,0}, "atom_8_basis_4-31G.tm"));
113 // cores.push_back(GaussianCore({1.7970, 0, 0}, "atom_1_basis_4-31G.tm"));
114 // cores.push_back(GaussianCore({-0.448, 1.740, 0.0}, "atom_1_basis_4-31G.tm"));
115 // GaussianSystem system;
116 // for(const GaussianCore &core : cores) {
117 // system.addCore(core);
118 // }
119 // mat C;
120 // UnrestrictedHartreeFockSolver solver(&system);
121 // solver.setConvergenceTreshold(1e-12);
122 // solver.setNIterationsMax(1e4);
123 // solver.setDensityMixFactor(0.95);
124 // solver.solve();
125 // cout << std::setprecision(20);
126 // cout << solver.energy() << endl;
127 // cout << solver.iterationsUsed() << endl;
128 // // CHECK_CLOSE(-149.5117583638509, solver.energy(), 1e-5);
129 // }
130 // TEST(WaterSixPlusPlus) {
131 // vector<GaussianCore> cores;
132 // cores.push_back(GaussianCore({0,0,0}, "atom_8_basis_6-311++Gdsds.tm"));
133 // cores.push_back(GaussianCore({1.797, 0.0, 0.0}, "atom_1_basis_6-311++Gdsds.tm"));
134 // cores.push_back(GaussianCore({-0.448, 1.740, 0.0}, "atom_1_basis_6-311++Gdsds.tm"));
135 // GaussianSystem system;
136 // for(const GaussianCore &core : cores) {
137 // system.addCore(core);
138 // }
139 // mat C;
140 // UnrestrictedHartreeFockSolver solver(&system);
141 // solver.setConvergenceTreshold(1e-12);
142 // solver.setNIterationsMax(1e4);
143 // solver.setDensityMixFactor(0.95);
144 // solver.solve();
145 // cout << std::setprecision(20);
146 
147 // cout << solver.energy() << endl;
148 // cout << solver.iterationsUsed() << endl;
149 // // CHECK_CLOSE(-149.5117583638509, solver.energy(), 1e-5);
150 // }
151 // TEST(OxygenSix) {
152 // vector<GaussianCore> cores;
153 // cores.push_back(GaussianCore({0,0,0}, "atom_8_basis_6-31Gds.tm"));
154 // cores.push_back(GaussianCore({2.281,0,0}, "atom_8_basis_6-31Gds.tm"));
155 // GaussianSystem system;
156 // for(const GaussianCore &core : cores) {
157 // system.addCore(core);
158 // }
159 // system.setNParticlesDown(7);
160 // mat C;
161 // UnrestrictedHartreeFockSolver solver(&system);
162 // solver.setConvergenceTreshold(1e-12);
163 // solver.setNIterationsMax(1e4);
164 // solver.setDensityMixFactor(0.95);
165 // solver.solve();
166 // cout << std::setprecision(20);
167 // cout << solver.energy() << endl;
168 // cout << solver.iterationsUsed() << endl;
169 // // CHECK_CLOSE(-149.5117583638509, solver.energy(), 1e-5);
170 // }
171 // TEST(BoronHydrogen) {
172 // vector<GaussianCore> cores;
173 // cores.push_back(GaussianCore({0,0,0}, "atom_5_basis_STO-3G.tm"));
174 // cores.push_back(GaussianCore({1.2,0,0}, "atom_1_basis_STO-3G.tm"));
175 // GaussianSystem system;
176 // for(const GaussianCore &core : cores) {
177 // system.addCore(core);
178 // }
179 // mat C;
180 // system.setNParticlesDown(4);
181 // UnrestrictedHartreeFockSolver solver(&system);
182 // solver.setConvergenceTreshold(1e-12);
183 // solver.setNIterationsMax(1e4);
184 // solver.setDensityMixFactor(0.95);
185 // solver.solve();
186 // cout << std::setprecision(20);
187 // cout << solver.energy() << endl;
188 // cout << solver.iterationsUsed() << endl;
189 // // CHECK_CLOSE(-149.5117583638509, solver.energy(), 1e-5);
190 // }
191  // TEST(OxygenSix) {
192  // vector<GaussianCore> cores;
193  // cores.push_back(GaussianCore({0,0,0}, "atom_8_basis_4-31G.tm"));
194  // cores.push_back(GaussianCore({2.282,0,0}, "atom_8_basis_4-31G.tm"));
195  // GaussianSystem system;
196  // for(const GaussianCore &core : cores) {
197  // system.addCore(core);
198  // }
199  // mat C;
200  // UnrestrictedHartreeFockSolver solver(&system);
201  // solver.setConvergenceTreshold(1e-12);
202  // solver.setNIterationsMax(1e4);
203  // solver.setDensityMixFactor(0.5);
204  // solver.solve();
205  // cout << solver.energy() << endl;
206  // cout << solver.iterationsUsed() << endl;
208  // }
209  // TEST(Water) {
210  // vector<GaussianCore> cores;
211  // cores.push_back(GaussianCore({0,0,0}, "atom_8_basis_4-31G.tm"));
212  // cores.push_back(GaussianCore({-1.43,1.108,0}, "atom_1_basis_4-31G.tm"));
213  // cores.push_back(GaussianCore({1.43,1.108,0}, "atom_1_basis_4-31G.tm"));
214  // GaussianSystem system;
215  // for(const GaussianCore &core : cores) {
216  // system.addCore(core);
217  // }
218  // mat C;
219  // HartreeFockSolver solver(&system);
220  // solver.setConvergenceTreshold(1e-12);
221  // solver.setNIterationsMax(1e3);
222  // solver.setDensityMixFactor(0.5);
223  // solver.solve();
224  // cout << solver.energy() << endl;
225  // cout << solver.iterationsUsed() << endl;
226  // CHECK_CLOSE(-75.90736859918989, solver.energy(), 1e-6);
227  // }
228 
229  // TEST(OxygenSixAsterisk) {
230  // vector<GaussianCore> cores;
231  // cores.push_back(GaussianCore({0,0,0}, "atom_8_basis_6-31Gs.tm"));
232  // cores.push_back(GaussianCore({2.282,0,0}, "atom_8_basis_6-31Gs.tm"));
233  // GaussianSystem system;
234  // for(const GaussianCore &core : cores) {
235  // system.addCore(core);
236  // }
237  // mat C;
238  // UnrestrictedHartreeFockSolver solver(&system);
239  // solver.setConvergenceTreshold(1e-8);
240  // solver.setNIterationsMax(1e4);
241  // solver.setDensityMixFactor(0.5);
242  // solver.solve();
243  // cout << solver.iterationsUsed() << endl;
244  // CHECK_CLOSE(-149.5876095851103, solver.energy(), 1e-5);
245  // }
246 
247  // TEST(CheapOxygen) {
248  // vector<GaussianCore> cores;
249  // cores.push_back(GaussianCore({0,0,0}, "atom_8_basis_cheap.tm"));
250  // cores.push_back(GaussianCore({2.282,0,0}, "atom_8_basis_cheap.tm"));
251  // GaussianSystem system;
252  // for(const GaussianCore &core : cores) {
253  // system.addCore(core);
254  // }
255  // mat C;
256  // HartreeFockSolver solver(&system);
257  // solver.setConvergenceTreshold(1e-12);
258  // solver.setNIterationsMax(1e4);
259  // solver.setDensityMixFactor(0.0);
260  // solver.advance();
261  // // solver.solve();
262  // cout << "Used " << solver.iterationsUsed() << endl;
263  // cout << "Oxygen cheap: " << solver.energy() << endl;
264  // }
265 
266  // TEST(Variations) {
267  // for(int i = 0; i < 100; i++) {
268  // vector<GaussianCore> cores;
269  // cores.push_back(GaussianCore({0,0,0}, "atom_1_basis_4-31G.tm"));
270  // cores.push_back(GaussianCore({0.5 + i * 0.1,0,0}, "atom_1_basis_4-31G.tm"));
271  // GaussianSystem system;
272  // for(const GaussianCore &core : cores) {
273  // system.addCore(core);
274  // }
275  // UnrestrictedHartreeFockSolver solver(&system);
276  // solver.setConvergenceTreshold(1e-9);
277  // solver.setNIterationsMax(1e3);
278  // solver.setDensityMixFactor(0.5);
279  // solver.solve();
280  // cout << solver.energy() << endl;
281  // }
282  // }
283  // TEST(Water) {
284  // vector<GaussianCore> cores;
285  // cores.push_back(GaussianCore({0,0,0}, "atom_8_basis_3-21G.tm"));
286  // cores.push_back(GaussianCore({-1.43,1.108,0}, "atom_1_basis_3-21G.tm"));
287  // cores.push_back(GaussianCore({1.43,1.108,0}, "atom_1_basis_3-21G.tm"));
288  // GaussianSystem system;
289  // for(const GaussianCore &core : cores) {
290  // system.addCore(core);
291  // }
292  // mat C;
293  // HartreeFockSolver solver(&system);
294  // solver.setConvergenceTreshold(1e-12);
295  // solver.setNIterationsMax(1e3);
296  // solver.setDensityMixFactor(0.0);
297  // solver.solve();
298  // CHECK_CLOSE(solver.energy(), -75.90736859918989, 1e-6);
299  // }
300  // TEST(HydrogenMolecule) {
301  // vector<GaussianCore> cores;
302  // cores.push_back(GaussianCore({0,0,0}, "atom_1_basis_6-31Gds.tm"));
303  // cores.push_back(GaussianCore({1.4,0.0,0}, "atom_1_basis_6-31Gds.tm"));
304  // GaussianSystem system;
305  // for(const GaussianCore &core : cores) {
306  // system.addCore(core);
307  // }
308  // HartreeFockSolver solver(&system);
309  // solver.setConvergenceTreshold(1e-12);
310  // solver.setNIterationsMax(1e3);
311  // solver.setDensityMixFactor(0.5);
312  // solver.solve();
313  // CHECK_CLOSE(-1.122933363617109, solver.energy(), 1e-6);
314  // }
315  // TEST(WaterAsterix) {
316  // vector<GaussianCore> cores;
317  // cores.push_back(GaussianCore({0,0,0}, "atom_8_basis_6-31Gs.tm"));
318  // cores.push_back(GaussianCore({-1.43,1.108,0}, "atom_1_basis_3-21G.tm"));
319  // cores.push_back(GaussianCore({1.43,1.108,0}, "atom_1_basis_3-21G.tm"));
320  // GaussianSystem system;
321  // for(const GaussianCore &core : cores) {
322  // system.addCore(core);
323  // }
324  // mat C;
325  // HartreeFockSolver solver(&system);
326  // solver.setConvergenceTreshold(1e-12);
327  // solver.setNIterationsMax(1e3);
328  // solver.setDensityMixFactor(0.0);
329  // solver.solve();
330  // CHECK_CLOSE(-75.90736859918989, solver.energy() , 1e-6);
331  // }
332 
333 }