1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /** @file SCF_unrestricted.cc
31
32 @brief Class for self-consistent field (SCF) procedure;
33 spin-unrestricted case.
34
35 @author: Elias Rudberg <em>responsible</em>.
36 */
37
38 #include <sstream>
39 #include "SCF_unrestricted.h"
40 #include "output.h"
41 #include "scf_utils.h"
42 #include "utilities.h"
43 #include "diis_unrestricted.h"
44 #include "density_projection.h"
45 #include "densfromf_full.h"
46 #include "density_description_file.h"
47 #include "densitymanager.h"
48 #include "matrix_utilities.h"
49 #include "units.h"
50 #include "atom_labels.h"
51 #include "dipole_moment.h"
52
53
SCF_unrestricted(const Molecule & molecule_,const Molecule & extraCharges_,const BasisInfoStruct & basisInfo_,const IntegralInfo & integralInfo_,const char * guessDmatFileName_,const JK::Params & J_K_params_,const Dft::GridParams & gridParams_,const SCF::Options & scfopts,const SCF::MatOptions & matOpts,ergo_real threshold_integrals_1el_input,int alpha_beta_diff_input)54 SCF_unrestricted::SCF_unrestricted(const Molecule& molecule_,
55 const Molecule& extraCharges_,
56 const BasisInfoStruct & basisInfo_,
57 const IntegralInfo& integralInfo_,
58 const char* guessDmatFileName_,
59 const JK::Params& J_K_params_,
60 const Dft::GridParams& gridParams_,
61 const SCF::Options& scfopts,
62 const SCF::MatOptions& matOpts,
63 ergo_real threshold_integrals_1el_input,
64 int alpha_beta_diff_input)
65 : SCF_general(molecule_,
66 extraCharges_,
67 basisInfo_,
68 integralInfo_,
69 guessDmatFileName_,
70 J_K_params_,
71 gridParams_,
72 scfopts,
73 matOpts,
74 threshold_integrals_1el_input),
75 alpha_beta_diff(alpha_beta_diff_input)
76 {
77 DIIS = new DIISManagerUnrestricted;
78 if(determine_number_of_electrons_unrestricted(noOfElectrons, alpha_beta_diff, &noOfElectrons_alpha, &noOfElectrons_beta) != 0)
79 {
80 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in determine_number_of_electrons_unrestricted");
81 throw "error in determine_number_of_electrons_unrestricted";
82 }
83
84 DensFromFock.do_unrestricted_calculations(); // set factor = 1
85
86 }
87
~SCF_unrestricted()88 SCF_unrestricted::~SCF_unrestricted()
89 {
90 delete ((DIISManagerUnrestricted*)DIIS);
91 }
92
93
get_Fock_matrices(symmMatrix & FockMatrix_a,symmMatrix & FockMatrix_b)94 void SCF_unrestricted::get_Fock_matrices(symmMatrix & FockMatrix_a, symmMatrix & FockMatrix_b)
95 {
96 FockMatrix_alpha.readFromFile();
97 FockMatrix_a = FockMatrix_alpha;
98 FockMatrix_alpha.writeToFile();
99 FockMatrix_beta.readFromFile();
100 FockMatrix_b = FockMatrix_beta;
101 FockMatrix_beta.writeToFile();
102 }
103
104
get_no_of_electrons(int & noOfElectrons_a,int & noOfElectrons_b)105 void SCF_unrestricted::get_no_of_electrons(int & noOfElectrons_a, int & noOfElectrons_b)
106 {
107 noOfElectrons_a = noOfElectrons_alpha;
108 noOfElectrons_b = noOfElectrons_beta;
109 }
110
111
initialize_matrices()112 void SCF_unrestricted::initialize_matrices()
113 {
114 densityMatrix_alpha.resetSizesAndBlocks(matOpts.size_block_info,
115 matOpts.size_block_info);
116 densityMatrix_beta.resetSizesAndBlocks(matOpts.size_block_info,
117 matOpts.size_block_info);
118 FockMatrix_alpha.resetSizesAndBlocks(matOpts.size_block_info,
119 matOpts.size_block_info);
120 FockMatrix_beta.resetSizesAndBlocks(matOpts.size_block_info,
121 matOpts.size_block_info);
122 Fprev_alpha.resetSizesAndBlocks(matOpts.size_block_info,
123 matOpts.size_block_info);
124 Fprev_beta.resetSizesAndBlocks(matOpts.size_block_info,
125 matOpts.size_block_info);
126 // Dprev_alpha.resetSizesAndBlocks(matOpts.size_block_info,
127 // matOpts.size_block_info);
128 // Dprev_beta.resetSizesAndBlocks(matOpts.size_block_info,
129 // matOpts.size_block_info);
130 F_ort_prev_alpha.resetSizesAndBlocks(matOpts.size_block_info,
131 matOpts.size_block_info);
132 F_ort_prev_beta.resetSizesAndBlocks(matOpts.size_block_info,
133 matOpts.size_block_info);
134 bestFockMatrixSoFar_alpha.resetSizesAndBlocks(matOpts.size_block_info,
135 matOpts.size_block_info);
136 bestFockMatrixSoFar_beta.resetSizesAndBlocks(matOpts.size_block_info,
137 matOpts.size_block_info);
138 bestFockMatrixSoFar2_alpha.resetSizesAndBlocks(matOpts.size_block_info,
139 matOpts.size_block_info);
140 bestFockMatrixSoFar2_beta.resetSizesAndBlocks(matOpts.size_block_info,
141 matOpts.size_block_info);
142 ErrorMatrix_alpha.resetSizesAndBlocks(matOpts.size_block_info,
143 matOpts.size_block_info);
144 ErrorMatrix_beta.resetSizesAndBlocks(matOpts.size_block_info,
145 matOpts.size_block_info);
146 G_alpha.resetSizesAndBlocks(matOpts.size_block_info,
147 matOpts.size_block_info);
148 G_beta.resetSizesAndBlocks(matOpts.size_block_info,
149 matOpts.size_block_info);
150 }
151
152
153
check_params()154 void SCF_unrestricted::check_params()
155 {
156 }
157
158
159
160
get_starting_guess_density()161 void SCF_unrestricted::get_starting_guess_density()
162 {
163 // set up starting guess
164
165 int n = basisInfo.noOfBasisFuncs;
166 DensFromFock.set_SCF_step(SCF_step);
167
168 if(guessDmatFileName != NULL)
169 {
170 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "getting starting guess density from file '%s'", guessDmatFileName);
171 int noOfDensityMatrices = 2;
172
173
174 symmMatrix* matrixList[2];
175 matrixList[0] = &densityMatrix_alpha;
176 matrixList[1] = &densityMatrix_beta;
177
178 int noOfElectronsList[2];
179 noOfElectronsList[0] = noOfElectrons_alpha + scfopts.starting_guess_spin_diff;
180 noOfElectronsList[1] = noOfElectrons_beta - scfopts.starting_guess_spin_diff;
181
182 if(load_density_and_project_sparse(DensFromFock,
183 guessDmatFileName,
184 noOfDensityMatrices,
185 &integralInfo,
186 basisInfo,
187 S_symm,
188 matrixList,
189 noOfElectronsList,
190 matOpts.size_block_info,
191 matOpts.permutationHML,
192 matOpts.sparse_threshold) != 0)
193 {
194 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in load_density_and_project_sparse");
195 throw "error in load_density_and_project_sparse";
196 }
197 }
198 else
199 {
200 if(scfopts.use_simple_starting_guess == 1)
201 {
202 if(get_simple_starting_guess_sparse(n, noOfElectrons_alpha, densityMatrix_alpha) != 0)
203 throw "error in get_simple_starting_guess_sparse";
204 densityMatrix_alpha.writeToFile();
205 if(get_simple_starting_guess_sparse(n, noOfElectrons_beta, densityMatrix_beta) != 0)
206 throw "error in get_simple_starting_guess_sparse";
207 densityMatrix_beta.writeToFile();
208 }
209 else if(scfopts.use_diag_guess_from_file == 1)
210 {
211 if(get_diag_matrix_from_file(n, densityMatrix_alpha, "diagdens_alpha.txt",matOpts.permutationHML) != 0)
212 throw "error in get_diag_matrix_from_file";
213 if(get_diag_matrix_from_file(n, densityMatrix_beta , "diagdens_beta.txt",matOpts.permutationHML) != 0)
214 throw "error in get_diag_matrix_from_file";
215 }
216 else
217 {
218 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
219 "calling get_dens_from_fock for alpha and beta to diagonalize H_core for starting guesses, n = %i, matOpts.sparse_threshold = %g",
220 n, (double)matOpts.sparse_threshold);
221
222
223
224 // save flags
225 int use_diag_on_error = DensFromFock.get_use_diag_on_error();
226
227 if(DensFromFock.get_use_diag_on_error_guess() == 1)
228 DensFromFock.set_use_diag_on_error();
229
230
231 // ALPHA
232 symmMatrix F_ort_prev_dummy_alpha;
233 F_ort_prev_dummy_alpha.resetSizesAndBlocks(matOpts.size_block_info,
234 matOpts.size_block_info);
235 F_ort_prev_dummy_alpha.writeToFile();
236 densityMatrix_alpha.writeToFile();
237
238 DensFromFock.set_no_occupied_orbs(noOfElectrons_alpha);
239 DensFromFock.clean_eigs_intervals();
240
241 if(DensFromFock.get_dens_from_fock(H_core_Matrix,
242 densityMatrix_alpha,
243 F_ort_prev_dummy_alpha) != 0)
244 {
245 throw "SCF_unrestricted::get_starting_guess_density: Error in get_dens_from_fock for alpha.";
246 }
247
248
249 // BETA
250 symmMatrix F_ort_prev_dummy_beta;
251 F_ort_prev_dummy_beta.resetSizesAndBlocks(matOpts.size_block_info,
252 matOpts.size_block_info);
253 F_ort_prev_dummy_beta.writeToFile();
254 densityMatrix_beta.writeToFile();
255
256 DensFromFock.set_no_occupied_orbs(noOfElectrons_beta);
257 DensFromFock.clean_eigs_intervals();
258
259 if(DensFromFock.get_dens_from_fock(H_core_Matrix,
260 densityMatrix_beta,
261 F_ort_prev_dummy_beta) != 0)
262 {
263 throw "SCF_unrestricted::get_starting_guess_density: Error in get_dens_from_fock for beta.";
264 }
265
266 // return values of flags
267 if( use_diag_on_error != 1 )
268 DensFromFock.unset_use_diag_on_error();
269
270
271
272 } // END ELSE use H_core to get starting guesses
273 } // END ELSE no dmat given
274
275 densityMatrix_alpha.readFromFile();
276 densityMatrix_beta.readFromFile();
277 output_sparsity_symm(n, densityMatrix_alpha, "starting guess density matrix (alpha)");
278 output_sparsity_symm(n, densityMatrix_beta , "starting guess density matrix (beta )");
279 densityMatrix_alpha.writeToFile();
280 densityMatrix_beta.writeToFile();
281
282 }
283
284
285
add_random_disturbance_to_starting_guess()286 void SCF_unrestricted::add_random_disturbance_to_starting_guess()
287 {
288 if(scfopts.sg_disturb_specific_elements > SCF::DISTURB_ELEMENT_MAX_COUNT)
289 throw "Error in SCF_unrestricted::add_random_disturbance_to_starting_guess: (scfopts.sg_disturb_specific_elements > SCF::DISTURB_ELEMENT_MAX_COUNT)";
290 int n = basisInfo.noOfBasisFuncs;
291 densityMatrix_alpha.readFromFile();
292 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
293 "SCF_unrestricted::add_random_disturbance_to_starting_guess, scfopts.starting_guess_disturbance = %7.3f",
294 scfopts.starting_guess_disturbance);
295 add_disturbance_to_matrix(n,
296 densityMatrix_alpha,
297 scfopts.starting_guess_disturbance,
298 scfopts.sg_disturb_specific_elements,
299 scfopts.disturbedElementIndexVector,
300 matOpts.permutationHML);
301 densityMatrix_alpha.writeToFile();
302 densityMatrix_beta.readFromFile();
303 add_disturbance_to_matrix(n,
304 densityMatrix_beta ,
305 -1 * scfopts.starting_guess_disturbance,
306 scfopts.sg_disturb_specific_elements,
307 scfopts.disturbedElementIndexVector,
308 matOpts.permutationHML);
309 densityMatrix_beta.writeToFile();
310 }
311
312
313
initialize_homo_lumo_limits()314 void SCF_unrestricted::initialize_homo_lumo_limits()
315 {
316 intervalType hugeInterval(-1e22, 1e22);
317 homoInterval_F_ort_prev_alpha = hugeInterval;
318 lumoInterval_F_ort_prev_alpha = hugeInterval;
319 homoInterval_F_ort_prev_beta = hugeInterval;
320 lumoInterval_F_ort_prev_beta = hugeInterval;
321 homoInterval_Fprev_alpha = hugeInterval;
322 lumoInterval_Fprev_alpha = hugeInterval;
323 homoInterval_Fprev_beta = hugeInterval;
324 lumoInterval_Fprev_beta = hugeInterval;
325 }
326
327
328
write_matrices_to_file()329 void SCF_unrestricted::write_matrices_to_file()
330 {
331 FockMatrix_alpha.writeToFile();
332 FockMatrix_beta.writeToFile();
333
334 Fprev_alpha.writeToFile();
335 Fprev_beta.writeToFile();
336
337 Dprev_alpha.writeToFile();
338 Dprev_beta.writeToFile();
339
340 bestFockMatrixSoFar_alpha.writeToFile();
341 bestFockMatrixSoFar_beta.writeToFile();
342
343 bestFockMatrixSoFar2_alpha.writeToFile();
344 bestFockMatrixSoFar2_beta.writeToFile();
345
346 F_ort_prev_alpha.writeToFile();
347 F_ort_prev_beta.writeToFile();
348 }
349
350
351
352 template<typename T>
353 static void
printMat(const T & m,const char * msg,int nbast)354 printMat(const T& m, const char *msg, int nbast)
355 {
356 #if 0
357 ergo_real* m_full = new ergo_real[nbast*nbast];
358 m.fullmatrix(m_full, nbast, nbast);
359 puts(msg);
360 for(int row=0; row<nbast; row++) {
361 for(int col=0; col<nbast; col++)
362 printf(" %10.6f", m_full[row + col*nbast]);
363 puts("");
364 }
365 delete []m_full;
366 #endif
367 }
368
get_2e_part_and_energy()369 void SCF_unrestricted::get_2e_part_and_energy()
370 {
371 // calculate Fock matrices
372 if(!scfopts.force_restricted) {
373 // F_alpha = H_core + G_alpha
374 // F_beta = H_core + G_beta
375 densityMatrix_alpha.readFromFile();
376 densityMatrix_beta.readFromFile();
377 if(get_2e_matrices_and_energy_sparse_unrestricted(basisInfo,
378 molecule,
379 integralInfo,
380 CAM_params,
381 G_alpha,
382 G_beta,
383 densityMatrix_alpha,
384 densityMatrix_beta,
385 J_K_params,
386 gridParams,
387 scfopts.use_dft,
388 &energy_2el,
389 noOfElectrons,
390 matOpts.size_block_info,
391 matOpts.permutationHML,
392 matOpts.inversePermutationHML) != 0)
393 {
394 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in get_2e_matrices_and_energy_sparse_unrestricted");
395 throw "error in get_2e_matrices_and_energy_sparse_unrestricted";
396 }
397 densityMatrix_alpha.writeToFile();
398 densityMatrix_beta.writeToFile();
399
400 H_core_Matrix.readFromFile();
401 FockMatrix_alpha.readFromFile();
402 FockMatrix_alpha = H_core_Matrix + G_alpha;
403 FockMatrix_alpha.frob_thresh(matOpts.sparse_threshold);
404 printMat(FockMatrix_alpha, "FockMatrix_alpha", basisInfo.noOfBasisFuncs);
405 FockMatrix_alpha.writeToFile();
406 FockMatrix_beta.readFromFile();
407 FockMatrix_beta = H_core_Matrix + G_beta;
408 FockMatrix_beta.frob_thresh(matOpts.sparse_threshold);
409 FockMatrix_beta.writeToFile();
410 H_core_Matrix.writeToFile();
411 } else {
412 // F_alpha = F_beta = Frestricted
413 densityMatrix_alpha.readFromFile();
414 densityMatrix_beta.readFromFile();
415 if(get_2e_matrices_and_energy_restricted_open(basisInfo,
416 molecule,
417 integralInfo,
418 CAM_params,
419 G_alpha, /*J_a+J_b+X_a+X_b*/
420 G_beta, /*J_a+J_b+X_a*/
421 densityMatrix_alpha,
422 densityMatrix_beta,
423 J_K_params,
424 gridParams,
425 scfopts.use_dft,
426 &energy_2el,
427 noOfElectrons,
428 matOpts.size_block_info,
429 matOpts.permutationHML,
430 matOpts.inversePermutationHML) != 0)
431 {
432 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in get_2e_matrices_and_energy_restricted_open");
433 throw "error in get_2e_matrices_and_energy_restricted_open";
434 }
435
436 printMat(densityMatrix_alpha, "densityMatrix_alpha",
437 basisInfo.noOfBasisFuncs);
438 printMat(densityMatrix_beta, "densityMatrix_beta",
439 basisInfo.noOfBasisFuncs);
440
441 symmMatrix S, Fc, Fo, FcMinusFo, fockConstrained, syTmp;
442 normalMatrix Sc, So, tmp, tmp1, tmp2;
443 Sc.resetSizesAndBlocks(matOpts.size_block_info,
444 matOpts.size_block_info);
445 So.resetSizesAndBlocks(matOpts.size_block_info,
446 matOpts.size_block_info);
447 tmp.resetSizesAndBlocks(matOpts.size_block_info,
448 matOpts.size_block_info);
449 tmp1.resetSizesAndBlocks(matOpts.size_block_info,
450 matOpts.size_block_info);
451 tmp2.resetSizesAndBlocks(matOpts.size_block_info,
452 matOpts.size_block_info);
453 fockConstrained.resetSizesAndBlocks(matOpts.size_block_info,
454 matOpts.size_block_info);
455
456 H_core_Matrix.readFromFile();
457 FockMatrix_alpha.readFromFile();
458 FockMatrix_beta.readFromFile();
459
460 get_overlap_matrix(S);
461 printMat(S, "overlap", basisInfo.noOfBasisFuncs);
462
463 syTmp = densityMatrix_alpha - densityMatrix_beta;
464 So = syTmp*S;
465
466 Sc = densityMatrix_beta*S;
467 printMat(Sc, "Sc", basisInfo.noOfBasisFuncs);
468
469 Fc = H_core_Matrix + G_alpha;
470 Fo = H_core_Matrix + G_beta; Fo *= 0.5;
471
472 printMat(Fc, "Fc", basisInfo.noOfBasisFuncs);
473 printMat(Fo, "Fo", basisInfo.noOfBasisFuncs);
474 tmp = 1;
475 tmp -= So;
476 // tmp = (-1.0)*So + one;
477 tmp1 = Fc*tmp;
478 tmp2 = transpose(tmp1)*tmp;
479 fockConstrained = tmp2;
480
481 tmp = 1;
482 tmp -= Sc;
483 tmp1 = Fo * tmp;
484 tmp2 = transpose(tmp1) * tmp;
485 syTmp = tmp2;
486 printMat(tmp2, "a2", basisInfo.noOfBasisFuncs);
487 fockConstrained += syTmp;
488
489 printMat(fockConstrained, "a1+a2", basisInfo.noOfBasisFuncs);
490
491 FcMinusFo = Fc - Fo;
492
493 tmp2 = FcMinusFo*Sc;
494 tmp1 = transpose(tmp2)*So;
495 tmp = transpose(tmp1);
496 tmp1 += tmp;
497 syTmp = tmp1;
498 fockConstrained += syTmp;
499 printMat(syTmp, "a3", basisInfo.noOfBasisFuncs);
500
501 printMat(fockConstrained, "a1+a2+a3", basisInfo.noOfBasisFuncs);
502
503 FockMatrix_alpha = fockConstrained;
504 FockMatrix_beta = fockConstrained;
505 FockMatrix_alpha.writeToFile();
506 FockMatrix_beta.writeToFile();
507 H_core_Matrix.writeToFile();
508 densityMatrix_alpha.writeToFile();
509 densityMatrix_beta.writeToFile();
510 }
511 }
512
513
514
515
output_sparsity_S_F_D(SCF_statistics & stats)516 void SCF_unrestricted::output_sparsity_S_F_D(SCF_statistics & stats)
517 {
518 int n = basisInfo.noOfBasisFuncs;
519 S_symm.readFromFile();
520 output_sparsity_symm(n, S_symm, "S");
521 S_symm.writeToFile();
522 FockMatrix_alpha.readFromFile();
523 output_sparsity_symm(n, FockMatrix_alpha, "F_alpha");
524 FockMatrix_alpha.writeToFile();
525 FockMatrix_beta.readFromFile();
526 output_sparsity_symm(n, FockMatrix_beta, "F_beta ");
527 FockMatrix_beta.writeToFile();
528 densityMatrix_alpha.readFromFile();
529 output_sparsity_symm(n, densityMatrix_alpha, "D_alpha");
530 densityMatrix_alpha.writeToFile();
531 densityMatrix_beta.readFromFile();
532 output_sparsity_symm(n, densityMatrix_beta, "D_beta ");
533 densityMatrix_beta.writeToFile();
534 }
535
536
537
538
calculate_energy()539 void SCF_unrestricted::calculate_energy()
540 {
541 // calculate energy
542 H_core_Matrix.readFromFile();
543 densityMatrix_alpha.readFromFile();
544 densityMatrix_beta.readFromFile();
545 energy = 0;
546 energy += symmMatrix::trace_ab(densityMatrix_alpha, H_core_Matrix);
547 energy += symmMatrix::trace_ab(densityMatrix_beta , H_core_Matrix);
548 energy += energy_2el;
549 densityMatrix_alpha.writeToFile();
550 densityMatrix_beta.writeToFile();
551 H_core_Matrix.writeToFile();
552 energy += nuclearEnergy;
553 }
554
555
556
557
get_FDSminusSDF()558 void SCF_unrestricted::get_FDSminusSDF()
559 {
560 int n = basisInfo.noOfBasisFuncs;
561
562 // ALPHA
563 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "calling compute_FDSminusSDF_sparse for ALPHA, n = %i", n);
564 densityMatrix_alpha.readFromFile();
565 S_symm.readFromFile();
566 FockMatrix_alpha.readFromFile();
567 compute_FDSminusSDF_sparse(n, FockMatrix_alpha, densityMatrix_alpha, S_symm,
568 ErrorMatrix_alpha, matOpts.sparse_threshold);
569 S_symm.writeToFile();
570 FockMatrix_alpha.writeToFile();
571 densityMatrix_alpha.writeToFile();
572 // write to file and read back again to reduce memory fragmentation.
573 ErrorMatrix_alpha.writeToFile();
574 ErrorMatrix_alpha.readFromFile();
575 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "compute_FDSminusSDF_sparse for ALPHA finished.");
576
577 // BETA
578 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "calling compute_FDSminusSDF_sparse for BETA, n = %i", n);
579 densityMatrix_beta.readFromFile();
580 S_symm.readFromFile();
581 FockMatrix_beta.readFromFile();
582 compute_FDSminusSDF_sparse(n, FockMatrix_beta, densityMatrix_beta, S_symm,
583 ErrorMatrix_beta, matOpts.sparse_threshold);
584 S_symm.writeToFile();
585 FockMatrix_beta.writeToFile();
586 densityMatrix_beta.writeToFile();
587 // write to file and read back again to reduce memory fragmentation.
588 ErrorMatrix_beta.writeToFile();
589 ErrorMatrix_beta.readFromFile();
590 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "compute_FDSminusSDF_sparse for BETA finished.");
591
592 output_sparsity(n, ErrorMatrix_alpha, "FDS-SDF (alpha)");
593 output_sparsity(n, ErrorMatrix_beta , "FDS-SDF (beta )");
594 }
595
596
597
598
get_error_measure()599 void SCF_unrestricted::get_error_measure()
600 {
601 ergo_real error_maxabs_alpha = compute_maxabs_sparse(ErrorMatrix_alpha);
602 ergo_real error_maxabs_beta = compute_maxabs_sparse(ErrorMatrix_beta );
603
604 ergo_real error_frob_alpha = ErrorMatrix_alpha.frob();
605 ergo_real error_frob_beta = ErrorMatrix_beta .frob();
606
607 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "maxabs FDS-SDF (alpha) is %8.3g", (double)error_maxabs_alpha);
608 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "maxabs FDS-SDF (beta ) is %8.3g", (double)error_maxabs_beta);
609 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "frob FDS-SDF (alpha) is %8.3g", (double)error_frob_alpha);
610 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "frob FDS-SDF (beta ) is %8.3g", (double)error_frob_beta);
611
612 errorMeasure = 0.5 * (error_maxabs_alpha + error_maxabs_beta);
613 }
614
615
616
617
618
add_to_DIIS_list()619 void SCF_unrestricted::add_to_DIIS_list()
620 {
621 FockMatrix_alpha.readFromFile();
622 FockMatrix_beta.readFromFile();
623 if(((DIISManagerUnrestricted*)DIIS)->AddIterationToList(FockMatrix_alpha,
624 FockMatrix_beta,
625 ErrorMatrix_alpha,
626 ErrorMatrix_beta) != 0)
627 {
628 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in DIIS AddIterationToList");
629 throw "error in DIIS AddIterationToList";
630 }
631 FockMatrix_alpha.writeToFile();
632 FockMatrix_beta.writeToFile();
633 }
634
635
636
637
638
update_best_fock_so_far()639 void SCF_unrestricted::update_best_fock_so_far()
640 {
641 // ALPHA
642 Fprev_alpha.readFromFile();
643 bestFockMatrixSoFar_alpha.readFromFile();
644 bestFockMatrixSoFar_alpha = Fprev_alpha;
645 Fprev_alpha.writeToFile();
646 bestFockMatrixSoFar_alpha.writeToFile();
647 FockMatrix_alpha.readFromFile();
648 bestFockMatrixSoFar2_alpha.readFromFile();
649 bestFockMatrixSoFar2_alpha = FockMatrix_alpha;
650 FockMatrix_alpha.writeToFile();
651 bestFockMatrixSoFar2_alpha.writeToFile();
652 // BETA
653 Fprev_beta.readFromFile();
654 bestFockMatrixSoFar_beta.readFromFile();
655 bestFockMatrixSoFar_beta = Fprev_beta;
656 Fprev_beta.writeToFile();
657 bestFockMatrixSoFar_beta.writeToFile();
658 FockMatrix_beta.readFromFile();
659 bestFockMatrixSoFar2_beta.readFromFile();
660 bestFockMatrixSoFar2_beta = FockMatrix_beta;
661 FockMatrix_beta.writeToFile();
662 bestFockMatrixSoFar2_beta.writeToFile();
663 }
664
665
666
667
combine_old_fock_matrices(ergo_real stepLength)668 void SCF_unrestricted::combine_old_fock_matrices(ergo_real stepLength)
669 {
670 // ALPHA
671 bestFockMatrixSoFar_alpha.readFromFile();
672 bestFockMatrixSoFar2_alpha.readFromFile();
673 FockMatrix_alpha.readFromFile();
674 FockMatrix_alpha = 0;
675 FockMatrix_alpha += stepLength * bestFockMatrixSoFar2_alpha;
676 FockMatrix_alpha += (1 - stepLength) * bestFockMatrixSoFar_alpha;
677 FockMatrix_alpha.writeToFile();
678 bestFockMatrixSoFar_alpha.writeToFile();
679 bestFockMatrixSoFar2_alpha.writeToFile();
680 // BETA
681 bestFockMatrixSoFar_beta.readFromFile();
682 bestFockMatrixSoFar2_beta.readFromFile();
683 FockMatrix_beta.readFromFile();
684 FockMatrix_beta = 0;
685 FockMatrix_beta += stepLength * bestFockMatrixSoFar2_beta;
686 FockMatrix_beta += (1 - stepLength) * bestFockMatrixSoFar_beta;
687 FockMatrix_beta.writeToFile();
688 bestFockMatrixSoFar_beta.writeToFile();
689 bestFockMatrixSoFar2_beta.writeToFile();
690 }
691
692
693
694
use_diis_to_get_new_fock_matrix()695 void SCF_unrestricted::use_diis_to_get_new_fock_matrix()
696 {
697 symmMatrix newFsymm_alpha;
698 symmMatrix newFsymm_beta;
699 if(((DIISManagerUnrestricted*)DIIS)->GetCombinedFockMatrices(newFsymm_alpha, newFsymm_beta) != 0)
700 {
701 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in DIIS GetCombinedFockMatrices");
702 throw "error in DIIS GetCombinedFockMatrices";
703 }
704 // ALPHA
705 FockMatrix_alpha.readFromFile();
706 FockMatrix_alpha = newFsymm_alpha;
707 FockMatrix_alpha.frob_thresh(matOpts.sparse_threshold);
708 FockMatrix_alpha.writeToFile();
709 // BETA
710 FockMatrix_beta.readFromFile();
711 FockMatrix_beta = newFsymm_beta;
712 FockMatrix_beta.frob_thresh(matOpts.sparse_threshold);
713 FockMatrix_beta.writeToFile();
714 }
715
716
717
718
719
clear_diis_list()720 void SCF_unrestricted::clear_diis_list()
721 {
722 ((DIISManagerUnrestricted*)DIIS)->ClearList();
723 }
724
725
726
727
clear_error_matrices()728 void SCF_unrestricted::clear_error_matrices()
729 {
730 ErrorMatrix_alpha.clear();
731 ErrorMatrix_beta.clear();
732 }
733
734
735
736
save_current_fock_as_fprev()737 void SCF_unrestricted::save_current_fock_as_fprev()
738 {
739 // ALPHA
740 FockMatrix_alpha.readFromFile();
741 Fprev_alpha.readFromFile();
742 Fprev_alpha = FockMatrix_alpha;
743 FockMatrix_alpha.writeToFile();
744 Fprev_alpha.writeToFile();
745 // BETA
746 FockMatrix_beta.readFromFile();
747 Fprev_beta.readFromFile();
748 Fprev_beta = FockMatrix_beta;
749 FockMatrix_beta.writeToFile();
750 Fprev_beta.writeToFile();
751 }
752
753
754
755
get_new_density_matrix()756 void SCF_unrestricted::get_new_density_matrix()
757 {
758 do_output(LOG_CAT_INFO, LOG_AREA_SCF,
759 "SCF_unrestricted::get_new_density_matrix, noOfElectrons_alpha noOfElectrons_beta = %5i %5i",
760 noOfElectrons_alpha, noOfElectrons_beta);
761
762 DensFromFock.set_SCF_step(SCF_step);
763
764 // ALPHA
765 DensFromFock.set_no_occupied_orbs( noOfElectrons_alpha);
766 DensFromFock.set_eigs_F_ort_prev(homoInterval_F_ort_prev_alpha, lumoInterval_F_ort_prev_alpha);
767 DensFromFock.set_generate_figures("_alpha");
768
769 // if eigenvectors are needed, set params
770 int use_init_guess = 0;
771 if(scfopts.use_prev_vector_as_initial_guess == 1 &&
772 scfopts.min_number_of_iterations-1 <= SCF_step &&
773 SCF_step > 1 &&
774 !eigVecUNOCC_alpha.empty() && !eigVecOCC_alpha.empty()) // ensure that we computed vectors in previous cycle
775 use_init_guess = 1;
776
777 // if we use new purification
778 if(SCF_step > 1 &&
779 (scfopts.min_number_of_iterations-1 <= SCF_step) &&
780 DensFromFock.get_use_purification() == 1 &&
781 DensFromFock.get_output_homo_and_lumo_eigenvectors() == 1)
782 {
783 DensFromFock.compute_eigenvectors(scfopts.eigenvectors_method,
784 scfopts.eigenvectors_iterative_method,
785 scfopts.eigensolver_accuracy,
786 scfopts.eigensolver_maxiter,
787 use_init_guess,
788 scfopts.try_eigv_on_next_iteration_if_fail);
789
790 DensFromFock.compute_eigenvectors_extra(scfopts.puri_compute_eigv_in_each_iteration, scfopts.run_shift_and_square_method_on_F);
791 }
792
793
794 if(DensFromFock.get_dens_from_fock(FockMatrix_alpha,
795 densityMatrix_alpha,
796 F_ort_prev_alpha) != 0)
797 {
798 throw "SCF_unrestricted::get_new_density_matrix: Error in get_dens_from_fock for alpha";
799 }
800
801 // return empty vectors in case nothing is computed
802 DensFromFock.get_computed_eigenpairs(eigVecUNOCC_alpha, eigVecOCC_alpha, eigValUNOCC_alpha, eigValOCC_alpha);
803
804 DensFromFock.get_eigs_F_ort_prev(homoInterval_F_ort_prev_alpha, lumoInterval_F_ort_prev_alpha);
805 DensFromFock.get_eigs_Fprev(homoInterval_Fprev_alpha, lumoInterval_Fprev_alpha);
806
807 DensFromFock.unset_generate_figures();
808 ergo_real electronicEntropyTerm_alpha = DensFromFock.get_result_entropy_term();
809
810
811 // BETA
812 DensFromFock.set_no_occupied_orbs( noOfElectrons_beta);
813 DensFromFock.set_eigs_F_ort_prev(homoInterval_F_ort_prev_beta, lumoInterval_F_ort_prev_beta);
814 DensFromFock.set_generate_figures("_beta");
815
816 // if eigenvectors are needed, set params
817 use_init_guess = 0;
818 if(scfopts.use_prev_vector_as_initial_guess == 1 &&
819 scfopts.min_number_of_iterations-1 <= SCF_step &&
820 SCF_step > 1 &&
821 !eigVecUNOCC_beta.empty() && !eigVecOCC_beta.empty()) // ensure that we computed vectors in previous cycle
822 use_init_guess = 1;
823
824 // if we use new purification
825 if(SCF_step > 1 &&
826 (scfopts.min_number_of_iterations-1 <= SCF_step) &&
827 DensFromFock.get_use_purification() == 1 &&
828 DensFromFock.get_output_homo_and_lumo_eigenvectors() == 1)
829 {
830 DensFromFock.compute_eigenvectors(scfopts.eigenvectors_method,
831 scfopts.eigenvectors_iterative_method,
832 scfopts.eigensolver_accuracy,
833 scfopts.eigensolver_maxiter,
834 use_init_guess,
835 scfopts.try_eigv_on_next_iteration_if_fail);
836
837 DensFromFock.compute_eigenvectors_extra(scfopts.puri_compute_eigv_in_each_iteration, scfopts.run_shift_and_square_method_on_F);
838 }
839
840
841 if(DensFromFock.get_dens_from_fock(FockMatrix_beta,
842 densityMatrix_beta,
843 F_ort_prev_beta) != 0)
844 {
845 throw "SCF_unrestricted::get_new_density_matrix: Error in get_dens_from_fock for beta";
846 }
847 // return empty vectors in case nothing is computed
848 DensFromFock.get_computed_eigenpairs(eigVecUNOCC_beta, eigVecOCC_beta, eigValUNOCC_beta, eigValOCC_beta);
849
850 DensFromFock.get_eigs_F_ort_prev(homoInterval_F_ort_prev_beta, lumoInterval_F_ort_prev_beta);
851 DensFromFock.get_eigs_Fprev(homoInterval_Fprev_beta, lumoInterval_Fprev_beta);
852
853 DensFromFock.unset_generate_figures();
854 ergo_real electronicEntropyTerm_beta = DensFromFock.get_result_entropy_term();
855
856
857 electronicEntropyTerm = electronicEntropyTerm_alpha + electronicEntropyTerm_beta;
858
859 // Report trace(DS) and spin info.
860 S_symm.readFromFile();
861 densityMatrix_alpha.readFromFile();
862 densityMatrix_beta.readFromFile();
863 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Tr( D_alpha * S ) = %22.11f", (double)symmMatrix::trace_ab(densityMatrix_alpha, S_symm));
864 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Tr( D_beta * S ) = %22.11f", (double)symmMatrix::trace_ab(densityMatrix_beta , S_symm));
865 symmMatrix spinDensityMatrix(densityMatrix_alpha);
866 spinDensityMatrix += (ergo_real)(-1.0) * densityMatrix_beta;
867 densityMatrix_alpha.writeToFile();
868 densityMatrix_beta.writeToFile();
869 S_symm.writeToFile();
870 normalMatrix spinDensityMatrix_normal(spinDensityMatrix);
871 spinDensityMatrix.clear();
872 ergo_real maxabs = compute_maxabs_sparse(spinDensityMatrix_normal);
873 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "maxabs element in spin density matrix : %22.11f", (double)maxabs);
874
875 ergo_real S2_exact, S2;
876 get_S2(S2_exact, S2);
877 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "<S2>_exact = %7.3f, <S2> = %22.11f", (double)S2_exact, (double)S2);
878 }
879
880
881
prep_matrix_description_struct(matrix_description_struct & mat,int nvalues,std::vector<int> & rowind,std::vector<int> & colind,std::vector<ergo_real> & values)882 static void prep_matrix_description_struct(matrix_description_struct & mat,
883 int nvalues,
884 std::vector<int> & rowind,
885 std::vector<int> & colind,
886 std::vector<ergo_real> & values) {
887 mat.nvalues = nvalues;
888 if(nvalues > 0) {
889 mat.rowind = &rowind[0];
890 mat.colind = &colind[0];
891 mat.values = &values[0];
892 }
893 else {
894 mat.rowind = NULL;
895 mat.colind = NULL;
896 mat.values = NULL;
897 }
898 }
899
write_density_to_file()900 void SCF_unrestricted::write_density_to_file()
901 {
902 matrix_description_struct matrixList[2];
903 // ALPHA
904 densityMatrix_alpha.readFromFile();
905 int nvalues_alpha = densityMatrix_alpha.nvalues();
906 std::vector<int> rowind_alpha;
907 rowind_alpha.reserve(nvalues_alpha);
908 std::vector<int> colind_alpha;
909 colind_alpha.reserve(nvalues_alpha);
910 std::vector<ergo_real> values_alpha;
911 values_alpha.reserve(nvalues_alpha);
912 densityMatrix_alpha.get_all_values(rowind_alpha,
913 colind_alpha,
914 values_alpha,
915 matOpts.inversePermutationHML,
916 matOpts.inversePermutationHML);
917 densityMatrix_alpha.writeToFile();
918 prep_matrix_description_struct(matrixList[0], nvalues_alpha, rowind_alpha, colind_alpha, values_alpha);
919 // BETA
920 densityMatrix_beta.readFromFile();
921 int nvalues_beta = densityMatrix_beta.nvalues();
922 std::vector<int> rowind_beta;
923 rowind_beta.reserve(nvalues_beta);
924 std::vector<int> colind_beta;
925 colind_beta.reserve(nvalues_beta);
926 std::vector<ergo_real> values_beta;
927 values_beta.reserve(nvalues_beta);
928 densityMatrix_beta.get_all_values(rowind_beta,
929 colind_beta,
930 values_beta,
931 matOpts.inversePermutationHML,
932 matOpts.inversePermutationHML);
933 densityMatrix_beta.writeToFile();
934 prep_matrix_description_struct(matrixList[1], nvalues_beta, rowind_beta, colind_beta, values_beta);
935
936 if(ddf_writeShellListAndDensityMatricesToFile_sparse(basisInfo,
937 2,
938 matrixList,
939 "density.bin") != 0)
940 {
941 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in ddf_writeShellListAndDensityMatricesToFile_sparse");
942 throw "error in ddf_writeShellListAndDensityMatricesToFile_sparse";
943 }
944
945 }
946
947
save_final_potential()948 void SCF_unrestricted::save_final_potential()
949 {
950 do_output(LOG_CAT_ERROR, LOG_AREA_SCF,
951 "error: save_final_potential not implemented for unrestricted case.");
952 throw "error: save_final_potential not implemented for unrestricted case.";
953 }
954
955
output_expected_values_pos_operator()956 void SCF_unrestricted::output_expected_values_pos_operator()
957 {
958 do_output(LOG_CAT_ERROR, LOG_AREA_SCF,
959 "error: output_expected_values_pos_operator not implemented for unrestricted case.");
960 throw "error: output_expected_values_pos_operator not implemented for unrestricted case.";
961 }
962
963
output_density_images_orbital(generalVector & eigVec,const std::string & filename_id)964 void SCF_unrestricted::output_density_images_orbital(generalVector &eigVec, const std::string &filename_id)
965 {
966 do_output(LOG_CAT_ERROR, LOG_AREA_SCF,
967 "error: output_density_images_orbital not implemented for unrestricted case.");
968 throw "error: output_density_images_orbital not implemented for unrestricted case.";
969 }
970
output_density_images()971 void SCF_unrestricted::output_density_images()
972 {
973 Util::TimeMeter timeMeter;
974
975 int n = basisInfo.noOfBasisFuncs;
976
977 ergo_real* densityMatrixFull_tot = new ergo_real[n*n];
978 ergo_real* densityMatrixFull_spin = new ergo_real[n*n];
979
980 // Get full matrix versions of density matrices
981 {
982 std::vector<ergo_real> densityMatrixFull_a(n*n);
983 std::vector<ergo_real> densityMatrixFull_b(n*n);
984
985 densityMatrix_alpha.readFromFile();
986 densityMatrix_alpha.fullMatrix(densityMatrixFull_a,
987 matOpts.inversePermutationHML,
988 matOpts.inversePermutationHML);
989 densityMatrix_alpha.writeToFile();
990
991 densityMatrix_beta.readFromFile();
992 densityMatrix_beta.fullMatrix(densityMatrixFull_b,
993 matOpts.inversePermutationHML,
994 matOpts.inversePermutationHML);
995 densityMatrix_beta.writeToFile();
996
997 for(int i = 0; i < n*n; i++)
998 {
999 densityMatrixFull_tot [i] = densityMatrixFull_a[i] + densityMatrixFull_b[i];
1000 densityMatrixFull_spin[i] = densityMatrixFull_a[i] - densityMatrixFull_b[i];
1001 }
1002 }
1003
1004 do_density_images(basisInfo,
1005 molecule,
1006 densityMatrixFull_tot,
1007 densityMatrixFull_spin,
1008 scfopts.output_density_images_boxwidth);
1009
1010 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_unrestricted::output_density_images finished OK.");
1011 timeMeter.print(LOG_AREA_SCF, "SCF_unrestricted::output_density_images");
1012 }
1013
1014
do_spin_flip(int atomCount)1015 void SCF_unrestricted::do_spin_flip(int atomCount)
1016 {
1017 int n = basisInfo.noOfBasisFuncs;
1018
1019 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_unrestricted::do_spin_flip, n = %6i", n);
1020
1021 ergo_real* densityMatrixFull_tot = new ergo_real[n*n];
1022 ergo_real* densityMatrixFull_spin = new ergo_real[n*n];
1023
1024 // Get full matrix versions of density matrices
1025 std::vector<ergo_real> densityMatrixFull_a(n*n);
1026 std::vector<ergo_real> densityMatrixFull_b(n*n);
1027
1028 densityMatrix_alpha.readFromFile();
1029 densityMatrix_beta.readFromFile();
1030 S_symm.readFromFile();
1031
1032 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Before spinflip: Tr( D_alpha * S ) = %22.11f", (double)symmMatrix::trace_ab(densityMatrix_alpha, S_symm));
1033 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Before spinflip: Tr( D_beta * S ) = %22.11f", (double)symmMatrix::trace_ab(densityMatrix_beta , S_symm));
1034
1035 densityMatrix_alpha.fullMatrix(densityMatrixFull_a,
1036 matOpts.inversePermutationHML,
1037 matOpts.inversePermutationHML);
1038 densityMatrix_beta.fullMatrix(densityMatrixFull_b,
1039 matOpts.inversePermutationHML,
1040 matOpts.inversePermutationHML);
1041
1042 for(int i = 0; i < n*n; i++)
1043 {
1044 densityMatrixFull_tot [i] = densityMatrixFull_a[i] + densityMatrixFull_b[i];
1045 densityMatrixFull_spin[i] = densityMatrixFull_a[i] - densityMatrixFull_b[i];
1046 }
1047
1048 // Now modify spin density matrix
1049 for(int i = 0; i < n; i++)
1050 for(int j = 0; j < n; j++)
1051 {
1052 // get atom index for basis funcs i and j
1053 int atomIndex_i = -1;
1054 int atomIndex_j = -1;
1055 for(int k = 0; k < molecule.getNoOfAtoms(); k++)
1056 {
1057 int iok = 1;
1058 int jok = 1;
1059 for(int coord = 0; coord < 3; coord++)
1060 {
1061 if(template_blas_fabs(molecule.getAtom(k).coords[coord] - basisInfo.basisFuncList[i].centerCoords[coord]) > 1e-5)
1062 iok = 0;
1063 if(template_blas_fabs(molecule.getAtom(k).coords[coord] - basisInfo.basisFuncList[j].centerCoords[coord]) > 1e-5)
1064 jok = 0;
1065 }
1066 if(iok == 1)
1067 atomIndex_i = k;
1068 if(jok == 1)
1069 atomIndex_j = k;
1070 } // END FOR k
1071
1072 // Change spin density matrix element if both basis functions belong to a spin-flip atom.
1073 if(atomIndex_i < atomCount && atomIndex_j < atomCount)
1074 densityMatrixFull_spin[i*n+j] *= -1;
1075
1076 } // END FOR i j
1077
1078 // Now the spin density matrix has been modified. Recreate new alpha- and beta- density matrices.
1079 for(int i = 0; i < n*n; i++)
1080 {
1081 densityMatrixFull_a[i] = 0.5 * (densityMatrixFull_tot[i] + densityMatrixFull_spin[i]);
1082 densityMatrixFull_b[i] = 0.5 * (densityMatrixFull_tot[i] - densityMatrixFull_spin[i]);
1083 }
1084
1085 densityMatrix_alpha.assignFromFull(densityMatrixFull_a,
1086 matOpts.permutationHML,
1087 matOpts.permutationHML);
1088 densityMatrix_beta.assignFromFull(densityMatrixFull_b,
1089 matOpts.permutationHML,
1090 matOpts.permutationHML);
1091
1092 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "After spinflip: Tr( D_alpha * S ) = %22.11f", (double)symmMatrix::trace_ab(densityMatrix_alpha, S_symm));
1093 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "After spinflip: Tr( D_beta * S ) = %22.11f", (double)symmMatrix::trace_ab(densityMatrix_beta , S_symm));
1094
1095 densityMatrix_alpha.writeToFile();
1096 densityMatrix_beta.writeToFile();
1097 S_symm.writeToFile();
1098
1099 delete [] densityMatrixFull_tot;
1100 delete [] densityMatrixFull_spin;
1101 }
1102
1103
1104
write_diag_dens_to_file()1105 void SCF_unrestricted::write_diag_dens_to_file()
1106 {
1107 int n = basisInfo.noOfBasisFuncs;
1108
1109 densityMatrix_alpha.readFromFile();
1110 write_diag_elements_to_file(n, densityMatrix_alpha, "diagdens_alpha.txt",
1111 matOpts.permutationHML);
1112 densityMatrix_alpha.writeToFile();
1113
1114 densityMatrix_beta.readFromFile();
1115 write_diag_elements_to_file(n, densityMatrix_beta , "diagdens_beta.txt",
1116 matOpts.permutationHML);
1117 densityMatrix_beta.writeToFile();
1118 }
1119
1120
save_full_matrices_for_matlab()1121 void SCF_unrestricted::save_full_matrices_for_matlab()
1122 {
1123 throw "error: SCF_unrestricted::save_full_matrices_for_matlab not implemented.";
1124 }
1125
1126
report_final_results()1127 void SCF_unrestricted::report_final_results()
1128 {
1129 ergo_real S2_exact, S2;
1130 get_S2(S2_exact, S2);
1131 do_output(LOG_CAT_RESULTS, LOG_AREA_SCF, "FINAL <S2> = %22.11f", (double)S2);
1132 }
1133
1134
get_S2(ergo_real & S2_exact,ergo_real & S2)1135 void SCF_unrestricted::get_S2(ergo_real & S2_exact, ergo_real & S2)
1136 {
1137 // FIXME: the <S2> stuff should work also if N_alpha < N_beta.
1138 if(noOfElectrons_alpha >= noOfElectrons_beta)
1139 {
1140 // Compute <S2>
1141 S_symm.readFromFile();
1142 densityMatrix_alpha.readFromFile();
1143 densityMatrix_beta.readFromFile();
1144
1145 normalMatrix B_alpha;
1146 B_alpha.resetSizesAndBlocks(matOpts.size_block_info,
1147 matOpts.size_block_info);
1148 B_alpha = (ergo_real)1.0 * densityMatrix_alpha * S_symm;
1149
1150 normalMatrix B_beta;
1151 B_beta.resetSizesAndBlocks(matOpts.size_block_info,
1152 matOpts.size_block_info);
1153 B_beta = (ergo_real)1.0 * densityMatrix_beta * S_symm;
1154
1155 ergo_real x = normalMatrix::trace_ab(B_alpha, B_beta);
1156 S2_exact = ( (ergo_real)(noOfElectrons_alpha - noOfElectrons_beta) / 2) * ( ( (ergo_real)(noOfElectrons_alpha - noOfElectrons_beta) / 2) + 1);
1157 S2 = S2_exact + noOfElectrons_beta - x;
1158
1159 densityMatrix_alpha.writeToFile();
1160 densityMatrix_beta.writeToFile();
1161 S_symm.writeToFile();
1162 }
1163 else
1164 {
1165 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "Error: <S2> computation not impl for na < nb");
1166 S2_exact = 0;
1167 S2 = 0;
1168 }
1169 }
1170
1171
save_density_as_prevdens()1172 void SCF_unrestricted::save_density_as_prevdens()
1173 {
1174 // alpha
1175 densityMatrix_alpha.readFromFile();
1176 Dprev_alpha.readFromFile();
1177 Dprev_alpha = densityMatrix_alpha;
1178 densityMatrix_alpha.writeToFile();
1179 Dprev_alpha.writeToFile();
1180 // beta
1181 densityMatrix_beta.readFromFile();
1182 Dprev_beta.readFromFile();
1183 Dprev_beta = densityMatrix_beta;
1184 densityMatrix_beta.writeToFile();
1185 Dprev_beta.writeToFile();
1186 }
1187
1188
1189
report_density_difference()1190 void SCF_unrestricted::report_density_difference()
1191 {
1192 // alpha
1193 densityMatrix_alpha.readFromFile();
1194 Dprev_alpha.readFromFile();
1195 symmMatrix diff_alpha(densityMatrix_alpha);
1196 diff_alpha += (ergo_real)-1.0 * Dprev_alpha;
1197 ergo_real diff_eucl_alpha = GetEuclideanNormOfMatrix(diff_alpha);
1198 densityMatrix_alpha.writeToFile();
1199 Dprev_alpha.writeToFile();
1200
1201 // beta
1202 densityMatrix_beta.readFromFile();
1203 Dprev_beta.readFromFile();
1204 symmMatrix diff_beta(densityMatrix_beta);
1205 diff_beta += (ergo_real)-1.0 * Dprev_beta;
1206 ergo_real diff_eucl_beta = GetEuclideanNormOfMatrix(diff_beta);
1207 densityMatrix_beta.writeToFile();
1208 Dprev_beta.writeToFile();
1209
1210 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_unrestricted::report_density_difference, diff_eucl_alpha = %22.11f", (double)diff_eucl_alpha);
1211 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_unrestricted::report_density_difference, diff_eucl_beta = %22.11f", (double)diff_eucl_beta );
1212 }
1213
compute_dipole_moment()1214 void SCF_unrestricted::compute_dipole_moment()
1215 {
1216 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_unrestricted::compute_dipole_moment");
1217 densityMatrix_alpha.readFromFile();
1218 densityMatrix_beta.readFromFile();
1219 symmMatrix densityMatrix_tot(densityMatrix_alpha);
1220 densityMatrix_tot += (ergo_real)1.0 * densityMatrix_beta;
1221 get_dipole_moment(densityMatrix_tot, basisInfo, matOpts.size_block_info, matOpts.permutationHML, molecule, LOG_AREA_SCF, "SCF");
1222 densityMatrix_alpha.writeToFile();
1223 densityMatrix_beta.writeToFile();
1224 }
1225
do_mulliken_pop_stuff()1226 void SCF_unrestricted::do_mulliken_pop_stuff()
1227 {
1228 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "SCF_unrestricted::do_mulliken_pop_stuff");
1229 densityMatrix_alpha.readFromFile();
1230 densityMatrix_beta.readFromFile();
1231 S_symm.readFromFile();
1232 symmMatrix densityMatrix_tot(densityMatrix_alpha);
1233 densityMatrix_tot += (ergo_real)1.0 * densityMatrix_beta;
1234 do_mulliken_atomic_charges(densityMatrix_tot,
1235 S_symm,
1236 basisInfo,
1237 matOpts.size_block_info,
1238 matOpts.permutationHML,
1239 matOpts.inversePermutationHML,
1240 molecule);
1241 densityMatrix_tot.clear(); /* This is to reduce memory usage. */
1242 symmMatrix spinDensityMatrix(densityMatrix_alpha);
1243 spinDensityMatrix += (ergo_real)-1.0 * densityMatrix_beta;
1244 do_mulliken_spin_densities(spinDensityMatrix,
1245 S_symm,
1246 basisInfo,
1247 matOpts.size_block_info,
1248 matOpts.permutationHML,
1249 matOpts.inversePermutationHML,
1250 molecule);
1251 densityMatrix_alpha.writeToFile();
1252 densityMatrix_beta.writeToFile();
1253 S_symm.writeToFile();
1254 }
1255
create_mtx_files_F(int const scfIter)1256 void SCF_unrestricted::create_mtx_files_F(int const scfIter) {
1257 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Creating mtx files for Fock matrices");
1258 { // alpha
1259 std::stringstream ss_fileName;
1260 ss_fileName << "F_matrix_alpha_" << scfIter;
1261 std::stringstream ss_id;
1262 ss_id << scfopts.calculation_identifier << " - effective Hamiltonian matrix (alpha), SCF cycle " << scfIter;
1263 FockMatrix_alpha.readFromFile();
1264 write_matrix_in_matrix_market_format( FockMatrix_alpha, matOpts.inversePermutationHML, ss_fileName.str(),
1265 ss_id.str(), scfopts.method_and_basis_set );
1266 FockMatrix_alpha.writeToFile();
1267 }
1268 { // beta
1269 std::stringstream ss_fileName;
1270 ss_fileName << "F_matrix_beta_" << scfIter;
1271 std::stringstream ss_id;
1272 ss_id << scfopts.calculation_identifier << " - effective Hamiltonian matrix (beta), SCF cycle " << scfIter;
1273 FockMatrix_beta.readFromFile();
1274 write_matrix_in_matrix_market_format( FockMatrix_beta, matOpts.inversePermutationHML, ss_fileName.str(),
1275 ss_id.str(), scfopts.method_and_basis_set );
1276 FockMatrix_beta.writeToFile();
1277 }
1278 }
create_mtx_files_D(int const scfIter)1279 void SCF_unrestricted::create_mtx_files_D(int const scfIter) {
1280 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Creating mtx files for density matrices");
1281 { // alpha
1282 std::stringstream ss_fileName;
1283 ss_fileName << "D_matrix_alpha_" << scfIter;
1284 std::stringstream ss_id;
1285 ss_id << scfopts.calculation_identifier << " - density matrix (alpha), SCF cycle " << scfIter;
1286 densityMatrix_alpha.readFromFile();
1287 write_matrix_in_matrix_market_format( densityMatrix_alpha, matOpts.inversePermutationHML, ss_fileName.str(),
1288 ss_id.str(), scfopts.method_and_basis_set );
1289 densityMatrix_alpha.writeToFile();
1290 }
1291 { // beta
1292 std::stringstream ss_fileName;
1293 ss_fileName << "D_matrix_beta_" << scfIter;
1294 std::stringstream ss_id;
1295 ss_id << scfopts.calculation_identifier << " - density matrix (beta), SCF cycle " << scfIter;
1296 densityMatrix_beta.readFromFile();
1297 write_matrix_in_matrix_market_format( densityMatrix_beta, matOpts.inversePermutationHML, ss_fileName.str(),
1298 ss_id.str(), scfopts.method_and_basis_set );
1299 densityMatrix_beta.writeToFile();
1300 }
1301 }
1302
1303
create_eigenvalues_files() const1304 void SCF_unrestricted::create_eigenvalues_files() const
1305 {
1306 if(! eigValOCC_alpha.empty() )
1307 {
1308 std::ofstream ff("occupied_spectrum_alpha.txt");
1309 for (size_t ind = 0; ind < eigValOCC_alpha.size(); ind++)
1310 {
1311 ff << (double)eigValOCC_alpha[ind] << std::endl;
1312 }
1313 ff.close();
1314 }
1315
1316 if(! eigValOCC_beta.empty() )
1317 {
1318 std::ofstream ff("occupied_spectrum_beta.txt");
1319 for (size_t ind = 0; ind < eigValOCC_beta.size(); ind++)
1320 {
1321 ff << (double)eigValOCC_beta[ind] << std::endl;
1322 }
1323 ff.close();
1324 }
1325
1326 if(! eigValUNOCC_alpha.empty() )
1327 {
1328 std::ofstream ff("unoccupied_spectrum_alpha.txt");
1329 for (size_t ind = 0; ind < eigValUNOCC_alpha.size(); ind++)
1330 {
1331 ff << (double)eigValUNOCC_alpha[ind] << std::endl;
1332 }
1333 ff.close();
1334 }
1335 if(! eigValUNOCC_beta.empty() )
1336 {
1337 std::ofstream ff("unoccupied_spectrum_beta.txt");
1338 for (size_t ind = 0; ind < eigValUNOCC_beta.size(); ind++)
1339 {
1340 ff << (double)eigValUNOCC_beta[ind] << std::endl;
1341 }
1342 ff.close();
1343 }
1344 }
1345
1346
create_eigenvectors_files() const1347 void SCF_unrestricted::create_eigenvectors_files() const
1348 {
1349 create_eigenvalues_files();
1350
1351 if(! eigVecOCC_alpha.empty() && ! eigVecOCC_beta.empty() )
1352 {
1353 create_eigvec_file(eigVecOCC_alpha[0], eigVecOCC_beta[0],
1354 "HOMO",
1355 "homo_coefficient_vec");
1356 assert(eigVecOCC_alpha.size() == eigVecOCC_beta.size());
1357 for (size_t i = 1; i < eigVecOCC_alpha.size(); i++)
1358 {
1359 std::stringstream ss_name, ss_filename;
1360 ss_name << "HOMO-" << i;
1361 ss_filename << "occ_" << i << "_coefficient_vec";
1362 create_eigvec_file(eigVecOCC_alpha[i], eigVecOCC_beta[i],
1363 ss_name.str().c_str(),
1364 ss_filename.str().c_str());
1365 }
1366 }
1367
1368 if(! eigVecUNOCC_alpha.empty() && ! eigVecUNOCC_beta.empty() )
1369 {
1370 create_eigvec_file(eigVecUNOCC_alpha[0], eigVecUNOCC_beta[0],
1371 "LUMO",
1372 "lumo_coefficient_vec");
1373 assert(eigVecUNOCC_alpha.size() == eigVecUNOCC_beta.size());
1374 for (size_t i = 1; i < eigVecUNOCC_alpha.size(); i++)
1375 {
1376 std::stringstream ss_name, ss_filename;
1377 ss_name << "LUMO-" << i;
1378 ss_filename << "unocc_" << i << "_coefficient_vec";
1379 create_eigvec_file(eigVecUNOCC_alpha[i], eigVecUNOCC_beta[i],
1380 ss_name.str().c_str(),
1381 ss_filename.str().c_str());
1382 }
1383 }
1384
1385 }
1386
create_eigvec_file(const generalVector & eigVec_alpha,const generalVector & eigVec_beta,const char * vector_name,const char * filename_id) const1387 void SCF_unrestricted::create_eigvec_file(const generalVector &eigVec_alpha,
1388 const generalVector &eigVec_beta,
1389 const char *vector_name,
1390 const char *filename_id) const
1391 {
1392 if (eigVec_alpha.is_empty()) {
1393 do_output(LOG_CAT_WARNING, LOG_AREA_SCF, "Failed to output %s (alpha) eigenvector to file. "
1394 "No %s (alpha) eigenvector stored.", vector_name, vector_name);
1395 }
1396 else{
1397 char ffname[888];
1398 sprintf(ffname, "%s_alpha.txt", filename_id);
1399 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Storing %s (alpha) eigenvector to file "
1400 "%s.", vector_name, ffname);
1401 int n = basisInfo.noOfBasisFuncs;
1402 std::vector<ergo_real> vec_perm(n);
1403 eigVec_alpha.fullvector(vec_perm);
1404 // now we have the permuted vector
1405 std::vector<ergo_real> vec(n);
1406 for (int ind = 0; ind < n; ind++)
1407 vec[ind] = vec_perm[matOpts.inversePermutationHML[ind]];
1408 std::ofstream ff(ffname);
1409 for (int ind = 0; ind < n; ind++)
1410 ff << (double)vec[ind] << std::endl;
1411 ff.close();
1412 }
1413
1414
1415 if (eigVec_beta.is_empty()) {
1416 do_output(LOG_CAT_WARNING, LOG_AREA_SCF, "Failed to output %s (beta) eigenvector to file. "
1417 "No %s (beta) eigenvector stored.", vector_name, vector_name);
1418 }
1419 else{
1420 char ffname[888];
1421 sprintf(ffname, "%s_beta.txt", filename_id);
1422 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Storing %s (beta) eigenvector to file "
1423 "%s.", vector_name, ffname);
1424 int n = basisInfo.noOfBasisFuncs;
1425 std::vector<ergo_real> vec_perm(n);
1426 eigVec_beta.fullvector(vec_perm);
1427 // now we have the permuted vector
1428 std::vector<ergo_real> vec(n);
1429 for (int ind = 0; ind < n; ind++)
1430 vec[ind] = vec_perm[matOpts.inversePermutationHML[ind]];
1431 std::ofstream ff(ffname);
1432 for (int ind = 0; ind < n; ind++)
1433 ff << (double)vec[ind] << std::endl;
1434 ff.close();
1435 }
1436
1437 }
1438
1439
1440
1441
1442 static void
output_orbital_coeffs_in_gabedit_order(const BasisInfoStruct & basisInfo,std::vector<int> const & shellIdxList,std::ofstream & ff,std::vector<ergo_real> const & orbital_vec)1443 output_orbital_coeffs_in_gabedit_order(const BasisInfoStruct & basisInfo,
1444 std::vector<int> const & shellIdxList,
1445 std::ofstream & ff,
1446 std::vector<ergo_real> const & orbital_vec) {
1447 int count = 0;
1448 for(int i = 0; i < basisInfo.noOfShells; i++) {
1449 int k = shellIdxList[i];
1450 int startIdx = basisInfo.shellList[k].startIndexInMatrix;
1451 switch(basisInfo.shellList[k].shellType) {
1452 case 0: // s-type shell
1453 ff << count+1 << " " << (double)orbital_vec[startIdx] << std::endl; count++;
1454 break;
1455 case 1: // p-type shell
1456 ff << count+1 << " " << (double)orbital_vec[startIdx+2] << std::endl; count++;
1457 ff << count+1 << " " << (double)orbital_vec[startIdx+0] << std::endl; count++;
1458 ff << count+1 << " " << (double)orbital_vec[startIdx+1] << std::endl; count++;
1459 break;
1460 case 2: // d-type shell
1461 ff << count+1 << " " << (double)orbital_vec[startIdx+2] << std::endl; count++;
1462 ff << count+1 << " " << (double)orbital_vec[startIdx+3] << std::endl; count++;
1463 ff << count+1 << " " << (double)orbital_vec[startIdx+1] << std::endl; count++;
1464 ff << count+1 << " " << (double)orbital_vec[startIdx+4] << std::endl; count++;
1465 ff << count+1 << " " << (double)orbital_vec[startIdx+0] << std::endl; count++;
1466 break;
1467 case 3: // f-type shell
1468 ff << count+1 << " " << (double)orbital_vec[startIdx+3] << std::endl; count++;
1469 ff << count+1 << " " << (double)orbital_vec[startIdx+4] << std::endl; count++;
1470 ff << count+1 << " " << (double)orbital_vec[startIdx+2] << std::endl; count++;
1471 ff << count+1 << " " << (double)orbital_vec[startIdx+5] << std::endl; count++;
1472 ff << count+1 << " " << (double)orbital_vec[startIdx+1] << std::endl; count++;
1473 ff << count+1 << " " << (double)orbital_vec[startIdx+6] << std::endl; count++;
1474 ff << count+1 << " " << (double)orbital_vec[startIdx+0] << std::endl; count++;
1475 break;
1476 case 4: // g-type shell
1477 ff << count+1 << " " << (double)orbital_vec[startIdx+4] << std::endl; count++;
1478 ff << count+1 << " " << (double)orbital_vec[startIdx+5] << std::endl; count++;
1479 ff << count+1 << " " << (double)orbital_vec[startIdx+3] << std::endl; count++;
1480 ff << count+1 << " " << (double)orbital_vec[startIdx+6] << std::endl; count++;
1481 ff << count+1 << " " << (double)orbital_vec[startIdx+2] << std::endl; count++;
1482 ff << count+1 << " " << (double)orbital_vec[startIdx+7] << std::endl; count++;
1483 ff << count+1 << " " << (double)orbital_vec[startIdx+1] << std::endl; count++;
1484 ff << count+1 << " " << (double)orbital_vec[startIdx+8] << std::endl; count++;
1485 ff << count+1 << " " << (double)orbital_vec[startIdx+0] << std::endl; count++;
1486 break;
1487 default: throw "error in output_orbital_coeffs_in_gabedit_order: shell types beyond g not implemented!";
1488 }
1489 }
1490 if(count != basisInfo.noOfBasisFuncs)
1491 throw "error in output_orbital_coeffs_in_gabedit_order: (count != basisInfo.noOfBasisFuncs)";
1492 }
1493
1494
1495
create_gabedit_file() const1496 void SCF_unrestricted::create_gabedit_file() const
1497 {
1498 if (eigVecOCC_alpha.empty() || eigVecUNOCC_alpha.empty() ||
1499 eigVecOCC_beta.empty() || eigVecUNOCC_beta.empty()) {
1500 do_output(LOG_CAT_WARNING, LOG_AREA_SCF, "Failed to output eigenvectors to gabedit file; no eigenvectors info available.");
1501 return;
1502 }
1503 if(basisInfo.use_6_d_funcs == 1) {
1504 do_output(LOG_CAT_WARNING, LOG_AREA_SCF, "Failed to output eigenvectors to gabedit file; not implemented for use_6_d_funcs case.");
1505 return;
1506 }
1507
1508 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Creating Gabedit file with eigenvector info.");
1509 int n = basisInfo.noOfBasisFuncs;
1510
1511 // Create Gabedit file.
1512 const char fileName [] = "gabeditfile.gab";
1513 std::ofstream ff(fileName);
1514 /* FIXME: check if we should use "Cart" or "Sphe" here. That is,
1515 should we use use_6_d_funcs? */
1516 int use_6_d_funcs = 0;
1517 ff << "[Gabedit Format] Sphe" << std::endl;
1518 ff << "[Atoms] Angs" << std::endl;
1519 for(int i = 0; i < molecule.getNoOfAtoms(); i++) {
1520 char atomLabelString[4];
1521 get_atom_label_from_charge_int(molecule.getAtom(i).charge, atomLabelString, 4);
1522 ff << atomLabelString << " " << i+1 << " " << (double)molecule.getAtom(i).charge
1523 << " " << (double)(molecule.getAtom(i).coords[0] / UNIT_one_Angstrom)
1524 << " " << (double)(molecule.getAtom(i).coords[1] / UNIT_one_Angstrom)
1525 << " " << (double)(molecule.getAtom(i).coords[2] / UNIT_one_Angstrom)
1526 << std::endl;
1527 }
1528
1529 std::vector<int> shellIdxList(basisInfo.noOfShells);
1530 int shellIdxCounter = 0;
1531 SquareFuncIntegrator sfi;
1532 ff << "[Basis]" << std::endl;
1533 for(int i = 0; i < molecule.getNoOfAtoms(); i++) {
1534 ff << i+1 << " 0" << std::endl;
1535 // Now output info about shells for this atom.
1536 for(int k = 0; k < basisInfo.noOfShells; k++) {
1537 // Check if this shell belongs to the current atom.
1538 ergo_real absdx = template_blas_fabs(basisInfo.shellList[k].centerCoords[0] - molecule.getAtom(i).coords[0]);
1539 ergo_real absdy = template_blas_fabs(basisInfo.shellList[k].centerCoords[1] - molecule.getAtom(i).coords[1]);
1540 ergo_real absdz = template_blas_fabs(basisInfo.shellList[k].centerCoords[2] - molecule.getAtom(i).coords[2]);
1541 ergo_real distlimit = 0.01;
1542 if(absdx > distlimit || absdy > distlimit || absdz > distlimit)
1543 continue;
1544 // OK, now we know this shell is at least very near the current atom.
1545 shellIdxList[shellIdxCounter] = k;
1546 shellIdxCounter++;
1547 char shellChar = 'x';
1548 int shellType = basisInfo.shellList[k].shellType;
1549 switch(shellType) {
1550 case 0: shellChar = 's'; break;
1551 case 1: shellChar = 'p'; break;
1552 case 2: shellChar = 'd'; break;
1553 case 3: shellChar = 'f'; break;
1554 case 4: shellChar = 'g'; break;
1555 default: throw "SCF_restricted::create_gabedit_file error: shell types beyond g not implemented!";
1556 }
1557 ff << shellChar << " " << basisInfo.shellList[k].noOfContr << " 1.00" << std::endl;
1558 for(int contridx = 0; contridx < basisInfo.shellList[k].noOfContr; contridx++) {
1559 ergo_real exponent = basisInfo.shellList[k].exponentList[contridx];
1560 ergo_real shellFactor = sfi.getShellFactor(integralInfo, exponent, shellType, use_6_d_funcs);
1561 ergo_real coeff = basisInfo.shellList[k].coeffList[contridx] / shellFactor;
1562 ff << (double)exponent << " " << (double)coeff << std::endl;
1563 }
1564 }
1565 // Blank line before shells for next atom.
1566 ff << std::endl;
1567 }
1568
1569 if(shellIdxCounter != basisInfo.noOfShells)
1570 throw "Error: (shellIdxCounter != basisInfo.noOfShells)";
1571
1572
1573 // MO section.
1574 ff << "[MO]" << std::endl;
1575
1576 // ALPHA
1577 // Occupied orbitals
1578 for (int orb = (int)eigVecOCC_alpha.size()-1; orb >= 0; orb--) // note that orb can be negative
1579 {
1580 std::vector<ergo_real> orb_vec_perm(n);
1581 eigVecOCC_alpha[orb].fullvector(orb_vec_perm);
1582 // now we have the permuted vector
1583 std::vector<ergo_real> orb_vec(n);
1584 for (int ind = 0; ind < n; ind++)
1585 {
1586 orb_vec[matOpts.inversePermutationHML[ind]] = orb_vec_perm[ind];
1587 }
1588 ff << "Spin=Alpha" << std::endl;
1589 if((int)eigValOCC_alpha.size() > orb) // in case we do not have eigenvalues
1590 ff << "Energy= " << (double)eigValOCC_alpha[orb] << std::endl;
1591 ff << "Occup= 1.000000" << std::endl;
1592 output_orbital_coeffs_in_gabedit_order(basisInfo, shellIdxList, ff, orb_vec);
1593 }
1594
1595 // Unoccupied orbitals
1596 for (size_t orb = 0; orb < eigVecUNOCC_alpha.size(); orb++)
1597 {
1598 std::vector<ergo_real> orb_vec_perm(n);
1599 eigVecUNOCC_alpha[orb].fullvector(orb_vec_perm);
1600 // now we have the permuted vector
1601 std::vector<ergo_real> orb_vec(n);
1602 for (int ind = 0; ind < n; ind++)
1603 {
1604 orb_vec[matOpts.inversePermutationHML[ind]] = orb_vec_perm[ind];
1605 }
1606 ff << "Spin=Alpha" << std::endl;
1607 if(eigValUNOCC_alpha.size() > orb) // in case we do not have eigenvalues
1608 ff << "Energy= " << (double)eigValUNOCC_alpha[orb] << std::endl;
1609 ff << "Occup= 0.000000" << std::endl;
1610 output_orbital_coeffs_in_gabedit_order(basisInfo, shellIdxList, ff, orb_vec);
1611 }
1612
1613
1614
1615
1616
1617
1618
1619 // BETA
1620 // Occupied orbitals
1621 for (int orb = (int)eigVecOCC_beta.size()-1; orb >= 0; orb--) // note that orb can be negative
1622 {
1623 std::vector<ergo_real> orb_vec_perm(n);
1624 eigVecOCC_beta[orb].fullvector(orb_vec_perm);
1625 // now we have the permuted vector
1626 std::vector<ergo_real> orb_vec(n);
1627 for (int ind = 0; ind < n; ind++)
1628 {
1629 orb_vec[matOpts.inversePermutationHML[ind]] = orb_vec_perm[ind];
1630 }
1631 ff << "Spin=Beta" << std::endl;
1632 if((int)eigValOCC_beta.size() > orb) // in case we do not have eigenvalues
1633 ff << "Energy= " << (double)eigValOCC_beta[orb] << std::endl;
1634 ff << "Occup= 1.000000" << std::endl;
1635 output_orbital_coeffs_in_gabedit_order(basisInfo, shellIdxList, ff, orb_vec);
1636 }
1637
1638 // Unoccupied orbitals
1639 for (size_t orb = 0; orb < eigVecUNOCC_beta.size(); orb++)
1640 {
1641 std::vector<ergo_real> orb_vec_perm(n);
1642 eigVecUNOCC_beta[orb].fullvector(orb_vec_perm);
1643 // now we have the permuted vector
1644 std::vector<ergo_real> orb_vec(n);
1645 for (int ind = 0; ind < n; ind++)
1646 {
1647 orb_vec[matOpts.inversePermutationHML[ind]] = orb_vec_perm[ind];
1648 }
1649 ff << "Spin=Beta" << std::endl;
1650 if(eigValUNOCC_beta.size() > orb) // in case we do not have eigenvalues
1651 ff << "Energy= " << (double)eigValUNOCC_beta[orb] << std::endl;
1652 ff << "Occup= 0.000000" << std::endl;
1653 output_orbital_coeffs_in_gabedit_order(basisInfo, shellIdxList, ff, orb_vec);
1654 }
1655
1656
1657 // Blank line before end of file.
1658 ff << std::endl;
1659 // Close file.
1660 ff.close();
1661
1662 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "Gabedit file '%s' with eigenvectors info created OK.", fileName);
1663 }
1664
1665
update_subspace_diff()1666 void SCF_unrestricted::update_subspace_diff()
1667 {
1668 throw "Error: SCF_unrestricted::update_subspace_diff() not implemented";
1669 }
1670
1671
disturb_fock_matrix(ergo_real subspaceError)1672 void SCF_unrestricted::disturb_fock_matrix(ergo_real subspaceError)
1673 {
1674 throw "SCF_unrestricted::disturb_fock_matrix() not implemented";
1675 }
1676
disturb_dens_matrix(ergo_real subspaceError)1677 void SCF_unrestricted::disturb_dens_matrix(ergo_real subspaceError)
1678 {
1679 throw "SCF_unrestricted::disturb_dens_matrix() not implemented";
1680 }
1681
disturb_dens_matrix_exact(ergo_real subspaceError)1682 void SCF_unrestricted::disturb_dens_matrix_exact(ergo_real subspaceError)
1683 {
1684 throw "SCF_unrestricted::disturb_dens_matrix_exact() not implemented";
1685 }
1686
compute_gradient_fixeddens()1687 void SCF_unrestricted::compute_gradient_fixeddens()
1688 {
1689 throw "SCF_unrestricted::compute_gradient_fixeddens() not implemented";
1690 }
1691
1692