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 diis_restricted.cc
31
32 @brief DIISManagerRestricted class implementing direct inversion
33 in the iterative subspace (DIIS) for restricted SCF calculations.
34
35 @author: Elias Rudberg <em>responsible</em>.
36 */
37
38 #include <string.h>
39
40 #include "diis_restricted.h"
41
42 #include "output.h"
43 #include "solve_lin_eq_syst.h"
44 #include "utilities.h"
45
46
DIISManagerRestricted()47 DIISManagerRestricted::DIISManagerRestricted()
48 : DIISManager()
49 {
50 }
51
~DIISManagerRestricted()52 DIISManagerRestricted::~DIISManagerRestricted()
53 {
54 ClearList();
55 }
56
AddIterationToList(symmMatrix & F,normalMatrix & E)57 int DIISManagerRestricted::AddIterationToList(symmMatrix & F, normalMatrix & E)
58 {
59 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "entering DIISManagerRestricted::AddIterationToList, IterCount = %2i", IterCount);
60
61 Util::TimeMeter timeMeter;
62
63 if(IterCount > MaxNoOfIters)
64 {
65 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error:(IterCount > MaxNoOfIters)");
66 return -1;
67 }
68
69 // check if the list is full
70 if(IterCount == MaxNoOfIters)
71 {
72 // remove oldest iteration
73 delete F_list[0][0];
74 F_list[0][0] = NULL;
75 delete E_list[0][0];
76 E_list[0][0] = NULL;
77 int i;
78 for(i = 0; i < IterCount-1; i++)
79 {
80 F_list[0][i] = F_list[0][i+1];
81 F_list[0][i+1] = NULL;
82 E_list[0][i] = E_list[0][i+1];
83 E_list[0][i+1] = NULL;
84 }
85 RemoveOldestIteration(); /* note that this changes the value of IterCount */
86 }
87
88 F_list[0][IterCount] = new symmMatrix(F);
89 F_list[0][IterCount]->writeToFile();
90 E_list[0][IterCount] = new normalMatrix(E);
91 E_list[0][IterCount]->writeToFile();
92
93 // Create new B matrix
94 int dimB = IterCount + 1;
95 int dimBnew = IterCount + 2;
96 ergo_real* Bnew = new ergo_real[dimBnew*dimBnew];
97 memset(Bnew, 0, dimBnew*dimBnew*sizeof(ergo_real));
98 int i, j;
99 for(i = 0; i < dimB; i++)
100 for(j = 0; j < dimB; j++)
101 Bnew[i*dimBnew+j] = B[i*dimB+j];
102 // Set two matrix elements to -1
103 Bnew[0*dimBnew+dimBnew-1] = -1;
104 Bnew[(dimBnew-1)*dimBnew+0] = -1;
105 // Now it remains to complete B with scalar products of error matrices
106 for(i = 0; i < IterCount; i++)
107 {
108 // compute dot product of error matrix i and E
109 E_list[0][i]->readFromFile();
110 ergo_real scalarProd = DoScalarProductOfErrorMatrices(E, *E_list[0][i]);
111 E_list[0][i]->writeToFile();
112 Bnew[(dimBnew-1)*dimBnew+(1+i)] = scalarProd;
113 Bnew[(1+i)*dimBnew+(dimBnew-1)] = scalarProd;
114 }
115 // Do scalar product of the new E with itself
116 Bnew[(dimBnew-1)*dimBnew+(dimBnew-1)] = DoScalarProductOfErrorMatrices(E, E);
117
118 // Copy Bnew to B
119 memcpy(B, Bnew, dimBnew*dimBnew*sizeof(ergo_real));
120
121 delete [] Bnew;
122
123 IterCount++;
124
125 do_output(LOG_CAT_INFO, LOG_AREA_SCF, "DIISManagerRestricted::AddIterationToList ending OK.");
126 timeMeter.print(LOG_AREA_SCF, "DIISManagerRestricted::AddIterationToList");
127
128 return 0;
129 }
130
ClearList()131 int DIISManagerRestricted::ClearList()
132 {
133 int i;
134 for(i = 0; i < IterCount; i++)
135 {
136 delete F_list[0][i];
137 F_list[0][i] = NULL;
138 delete E_list[0][i];
139 E_list[0][i] = NULL;
140 }
141 IterCount = 0;
142 return 0;
143 }
144
GetCombinedFockMatrix(symmMatrix & result)145 int DIISManagerRestricted::GetCombinedFockMatrix(symmMatrix & result)
146 {
147 if(IterCount <= 0)
148 {
149 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in DIISManagerRestricted::GetCombinedFockMatrix: (IterCount <= 0)");
150 return -1;
151 }
152
153 int dimB = IterCount + 1;
154 ergo_real* RHS = new ergo_real[dimB];
155 ergo_real* cVector = new ergo_real[dimB];
156
157 // Construct vector RHS
158 RHS[0] = -1;
159 int i;
160 for(i = 0; i < IterCount; i++)
161 RHS[i+1] = 0;
162
163 // Solve equation system B*x = HL
164 if(solve_linear_equation_system(dimB, B, RHS, cVector) != 0)
165 {
166 do_output(LOG_CAT_ERROR, LOG_AREA_SCF, "error in solve_linear_equation_system");
167 return -1;
168 }
169
170 // Create linear combination of Fock matrices using coefficients from cVector.
171 F_list[0][0]->readFromFile();
172 result = *F_list[0][0];
173 F_list[0][0]->writeToFile();
174 result *= cVector[1];
175 for(i = 1; i < IterCount; i++)
176 {
177 F_list[0][i]->readFromFile();
178 result += cVector[1+i] * (*F_list[0][i]);
179 F_list[0][i]->writeToFile();
180 } // END FOR i
181
182
183 delete [] RHS;
184 delete [] cVector;
185
186 return 0;
187 }
188
189
190
191