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