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 density_description_file.cc
31
32 @brief An interface file for writing and reading density matrices
33 to/from a file, including basis set information.
34
35 @author: Elias Rudberg <em>responsible</em>
36 */
37
38 #define _LARGEFILE_SOURCE 1
39 #include <stdlib.h>
40 #include <stdio.h>
41 #include <sys/types.h>
42 #include <sys/stat.h>
43 #include <unistd.h>
44 #include <math.h>
45 #include <errno.h>
46 #include <string>
47
48 #include "density_description_file.h"
49 #include "densfromf_full.h"
50 #include "integrals_general.h"
51 #include "matrix_algebra.h"
52 #include "memorymanag.h"
53 #include "output.h"
54 #include "utilities.h"
55
56 // Version number.
57 // This is written as part of the header.
58 // Should be changed if the implementation is changed in such a way
59 // that old density files will no longer work.
60 const int DENSITY_FILE_VERSION_NUMBER = 10002;
61
62
63 // Matrix storage types
64 #define MATRIX_STORAGE_TYPE_FULL 1
65 #define MATRIX_STORAGE_TYPE_TRIANGLE 2
66 #define MATRIX_STORAGE_TYPE_VECTORS 3
67
68 const ergo_real THRESHOLD_FOR_VECTOR_STORAGE = 1e-12;
69
70
71 // File header.
72 // We want this to contain an even number of integers,
73 // to avoid compilers inserting padding.
74 // We use double instead of size_t for sizes to ensure compatibility of density files.
75 typedef struct
76 {
77 int densityFileVersion;
78 int typeOfMatrixStorage;
79 int noOfShells;
80 int noOfBasisFuncs;
81 int noOfDensityMatrices;
82 int padding_int;
83 double matrixSize_1;
84 double matrixSize_2;
85 double fileSizeInBytes;
86 } densityFileHeaderStruct;
87
88
89 /** stores the upper triangle of a square matrix in the specified
90 memory block.
91 @param p memory block
92 @param n matrix dimension
93 @param matrix square matrix
94 */
95 static void
ddf_store_triangular_matrix(char * p,int n,const ergo_real * matrix)96 ddf_store_triangular_matrix(char* p, int n, const ergo_real* matrix)
97 {
98 ergo_real* resultPtr = (ergo_real*)p;
99 long count = 0;
100 for(long i = 0; i < n; i++)
101 for(long j = i; j < n; j++)
102 {
103 resultPtr[count] = matrix[i*n+j];
104 count++;
105 }
106 }
107
108
109 /** stores the upper triangle of the matrix given by the
110 matrix_description_struct in the specified memory block.
111 @param p memory block
112 @param n matrix dimension
113 @param matrix matrix
114 */
115 static void
ddf_store_triangular_matrix_sparse(char * p,long n,matrix_description_struct * matrix)116 ddf_store_triangular_matrix_sparse(char* p, long n, matrix_description_struct* matrix)
117 {
118 ergo_real* resultPtr = (ergo_real*)p;
119 long totCount = n*(n+1) / 2;
120 memset(resultPtr, 0, totCount*sizeof(ergo_real));
121 for(long i = 0; i < matrix->nvalues; i++)
122 {
123 long row = matrix->rowind[i];
124 long col = matrix->colind[i];
125 ergo_real value = matrix->values[i];
126 // Make sure (row,col) refers to upper triangle.
127 if(col < row)
128 {
129 long tmp = row;
130 row = col;
131 col = tmp;
132 }
133 long index = row*n+col-row*(row+1)/2;
134 resultPtr[index] = value;
135 }
136 }
137
138
139 static void
ddf_get_triangular_matrix_from_storage(const char * p,long n,ergo_real * resultMatrix)140 ddf_get_triangular_matrix_from_storage(const char* p, long n, ergo_real* resultMatrix)
141 {
142 const ergo_real* sourcePtr = (ergo_real*)p;
143 long count = 0;
144 for(long i = 0; i < n; i++)
145 for(long j = i; j < n; j++)
146 {
147 ergo_real a = sourcePtr[count];
148 resultMatrix[i*n+j] = a;
149 resultMatrix[j*n+i] = a;
150 count++;
151 }
152 }
153
154 static void
ddf_get_triangular_matrix_from_storage_sparse(const char * p,int n,int * rowind,int * colind,ergo_real * values)155 ddf_get_triangular_matrix_from_storage_sparse(const char* p,
156 int n,
157 int* rowind,
158 int* colind,
159 ergo_real* values)
160 {
161 const ergo_real* sourcePtr = (ergo_real*)p;
162 long count = 0;
163 for(long i = 0; i < n; i++)
164 for(long j = i; j < n; j++)
165 {
166 ergo_real a = sourcePtr[count];
167 rowind[count] = i;
168 colind[count] = j;
169 values[count] = a;
170 count++;
171 }
172 }
173
174 static int
ddf_get_nvalues_symm_matrix(int n,const ergo_real * matrix)175 ddf_get_nvalues_symm_matrix(int n, const ergo_real* matrix)
176 {
177 long nvalues = 0;
178 for(long i = 0; i < n; i++)
179 for(long j = i; j < n; j++)
180 {
181 if(template_blas_fabs(matrix[i*n+j]) > THRESHOLD_FOR_VECTOR_STORAGE)
182 nvalues++;
183 }
184 return nvalues;
185 }
186
187 static int
ddf_store_matrix_as_vectors(char * p,int n,const ergo_real * matrix,size_t sizeInBytes)188 ddf_store_matrix_as_vectors(char* p, int n, const ergo_real* matrix, size_t sizeInBytes)
189 {
190 long nvalues = ddf_get_nvalues_symm_matrix(n, matrix);
191 if(sizeInBytes != nvalues*(sizeof(ergo_real)+2*sizeof(int)))
192 return -1;
193 int* rowind = (int*)p;
194 int* colind = &rowind[nvalues];
195 // The values follow as a sequence of ergo_real after rowind and colind, but we use a char ptr and memcpy to avoid memory alignment problems e.g. in quad precision case.
196 char* valuesPtr = (char*)(&rowind[2*nvalues]);
197 long count = 0;
198 for(long i = 0; i < n; i++)
199 for(long j = i; j < n; j++)
200 {
201 if(template_blas_fabs(matrix[i*n+j]) > THRESHOLD_FOR_VECTOR_STORAGE)
202 {
203 if(count >= nvalues)
204 return -1;
205 rowind[count] = i;
206 colind[count] = j;
207 memcpy(valuesPtr+count*sizeof(ergo_real), &matrix[i*n+j], sizeof(ergo_real)); // Use memcpy to avoid memory alignment problems e.g. in quad precision case.
208 count++;
209 }
210 } // END FOR i j
211 return 0;
212 }
213
214
215 static int
ddf_store_matrix_as_vectors_sparse(char * p,int n,matrix_description_struct * matrix,size_t sizeInBytes)216 ddf_store_matrix_as_vectors_sparse(char* p, int n, matrix_description_struct* matrix, size_t sizeInBytes)
217 {
218 long nvalues = matrix->nvalues;
219 if(sizeInBytes != nvalues*(sizeof(ergo_real)+2*sizeof(int)))
220 return -1;
221 int* rowind = (int*)p;
222 int* colind = &rowind[nvalues];
223 ergo_real* values = (ergo_real*)(&rowind[2*nvalues]);
224 for(long i = 0; i < nvalues; i++)
225 {
226 rowind[i] = matrix->rowind[i];
227 colind[i] = matrix->colind[i];
228 values[i] = matrix->values[i];
229 }
230 return 0;
231 }
232
233
234 static int
ddf_get_matrix_from_vector_storage(const char * p,size_t sizeInBytes,int n,ergo_real * resultMatrix)235 ddf_get_matrix_from_vector_storage(const char* p, size_t sizeInBytes, int n, ergo_real* resultMatrix)
236 {
237 int nBytesPerElement = sizeof(ergo_real)+2*sizeof(int);
238 if(sizeInBytes % nBytesPerElement != 0)
239 return -1;
240 long nvalues = sizeInBytes / nBytesPerElement;
241 int* rowind = (int*)p;
242 int* colind = &rowind[nvalues];
243 // The values follow as a sequence of ergo_real after rowind and colind, but we use a char ptr and memcpy to avoid memory alignment problems e.g. in quad precision case.
244 char* valuesPtr = (char*)(&rowind[2*nvalues]);
245 memset(resultMatrix, 0, n*n*sizeof(ergo_real));
246 for(long i = 0; i < nvalues; i++)
247 {
248 long row = rowind[i];
249 long col = colind[i];
250 ergo_real value = 0;
251 memcpy(&value, &valuesPtr[i*sizeof(ergo_real)], sizeof(ergo_real)); // Use memcpy to avoid memory alignment problems e.g. in quad precision case.
252 resultMatrix[row*n+col] = value;
253 resultMatrix[col*n+row] = value;
254 }
255 return 0;
256 }
257
258 static int
ddf_get_matrix_from_vector_storage_sparse(const char * p,size_t sizeInBytes,int n,int * rowind2,int * colind2,ergo_real * values2)259 ddf_get_matrix_from_vector_storage_sparse(const char* p, size_t sizeInBytes, int n, int* rowind2, int* colind2, ergo_real* values2)
260 {
261 int nBytesPerElement = sizeof(ergo_real)+2*sizeof(int);
262 if(sizeInBytes % nBytesPerElement != 0)
263 return -1;
264 long nvalues = sizeInBytes / nBytesPerElement;
265 int* rowind = (int*)p;
266 int* colind = &rowind[nvalues];
267 ergo_real* values = (ergo_real*)(&rowind[2*nvalues]);
268 for(long i = 0; i < nvalues; i++)
269 {
270 rowind2[i] = rowind[i];
271 colind2[i] = colind[i];
272 values2[i] = values[i];
273 }
274 return 0;
275 }
276
277 static size_t
ddf_get_matrix_storage_size(int storageType,int n,const ergo_real * matrix)278 ddf_get_matrix_storage_size(int storageType, int n, const ergo_real* matrix)
279 {
280 size_t storageSize = 0;
281 switch(storageType)
282 {
283 case MATRIX_STORAGE_TYPE_FULL:
284 storageSize = n*n*sizeof(ergo_real);
285 break;
286 case MATRIX_STORAGE_TYPE_TRIANGLE:
287 storageSize = (n * (n+1) / 2) * sizeof(ergo_real);
288 break;
289 case MATRIX_STORAGE_TYPE_VECTORS:
290 storageSize = ddf_get_nvalues_symm_matrix(n, matrix) * (sizeof(ergo_real) + 2*sizeof(int));
291 break;
292 default:
293 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_storage_size: unknown storage type %i", storageType);
294 return 0;
295 }
296 return storageSize;
297 }
298
299
300 static size_t
ddf_get_matrix_storage_size_sparse(int storageType,int n,matrix_description_struct * matrix)301 ddf_get_matrix_storage_size_sparse(int storageType, int n, matrix_description_struct* matrix)
302 {
303 size_t storageSize = 0;
304 switch(storageType)
305 {
306 case MATRIX_STORAGE_TYPE_FULL:
307 storageSize = n*n*sizeof(ergo_real);
308 break;
309 case MATRIX_STORAGE_TYPE_TRIANGLE:
310 storageSize = (n * (n+1) / 2) * sizeof(ergo_real);
311 break;
312 case MATRIX_STORAGE_TYPE_VECTORS:
313 storageSize = matrix->nvalues * (sizeof(ergo_real) + 2*sizeof(int));
314 break;
315 default:
316 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_storage_size_sparse: unknown storage type %i", storageType);
317 return 0;
318 }
319 return storageSize;
320 }
321
322
323 static int
ddf_store_matrix(char * p,size_t sizeInBytes,int n,const ergo_real * matrix,int storageType)324 ddf_store_matrix(char* p, size_t sizeInBytes, int n, const ergo_real* matrix, int storageType)
325 {
326 size_t sizeNeeded = 0;
327 switch(storageType)
328 {
329 case MATRIX_STORAGE_TYPE_FULL:
330 sizeNeeded = n*n*sizeof(ergo_real);
331 if(sizeNeeded != sizeInBytes)
332 {
333 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix: (sizeNeeded != sizeInBytes)");
334 return -1;
335 }
336 memcpy(p, matrix, sizeNeeded);
337 break;
338 case MATRIX_STORAGE_TYPE_TRIANGLE:
339 sizeNeeded = (n * (n+1) / 2) * sizeof(ergo_real);
340 if(sizeNeeded != sizeInBytes)
341 {
342 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix: (sizeNeeded != sizeInBytes)");
343 return -1;
344 }
345 ddf_store_triangular_matrix(p, n, matrix);
346 break;
347 case MATRIX_STORAGE_TYPE_VECTORS:
348 if(ddf_store_matrix_as_vectors(p, n, matrix, sizeInBytes) != 0)
349 {
350 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix_as_vectors.");
351 return -1;
352 }
353 break;
354 default:
355 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix: unknown storage type %i", storageType);
356 return -1;
357 }
358 return 0;
359 }
360
361
362 static int
ddf_store_matrix_sparse(char * p,size_t sizeInBytes,int n,matrix_description_struct * matrix,int storageType)363 ddf_store_matrix_sparse(char* p,
364 size_t sizeInBytes,
365 int n,
366 matrix_description_struct* matrix,
367 int storageType)
368 {
369 size_t sizeNeeded = 0;
370 switch(storageType)
371 {
372 case MATRIX_STORAGE_TYPE_FULL:
373 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix_sparse: MATRIX_STORAGE_TYPE_FULL not supported.");
374 return -1;
375 break;
376 case MATRIX_STORAGE_TYPE_TRIANGLE:
377 sizeNeeded = (n * (n+1) / 2) * sizeof(ergo_real);
378 if(sizeNeeded != sizeInBytes)
379 {
380 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix: (sizeNeeded != sizeInBytes)");
381 return -1;
382 }
383 ddf_store_triangular_matrix_sparse(p, n, matrix);
384 break;
385 case MATRIX_STORAGE_TYPE_VECTORS:
386 if(ddf_store_matrix_as_vectors_sparse(p, n, matrix, sizeInBytes) != 0)
387 {
388 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix_as_vectors.");
389 return -1;
390 }
391 break;
392 default:
393 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix: unknown storage type %i", storageType);
394 return -1;
395 }
396 return 0;
397 }
398
399
400 static int
ddf_get_matrix_nvalues_from_storage(const char * p,size_t sizeInBytes,int n,long * result_nvalues,int storageType)401 ddf_get_matrix_nvalues_from_storage(const char* p, size_t sizeInBytes, int n,
402 long* result_nvalues, int storageType)
403 {
404 switch(storageType)
405 {
406 case MATRIX_STORAGE_TYPE_FULL:
407 // In this case the needed nvalues is n*n
408 *result_nvalues = n*n;
409 break;
410 case MATRIX_STORAGE_TYPE_TRIANGLE:
411 *result_nvalues = (n * (n+1) / 2);
412 break;
413 case MATRIX_STORAGE_TYPE_VECTORS:
414 {
415 int nBytesPerElement = sizeof(ergo_real)+2*sizeof(int);
416 if(sizeInBytes % nBytesPerElement != 0)
417 return -1;
418 *result_nvalues = sizeInBytes / nBytesPerElement;
419 }
420 break;
421 default:
422 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_nvalues_from_storage: unknown storage type %i", storageType);
423 return -1;
424 }
425 return 0;
426 }
427
428
429 static int
ddf_get_matrix_from_storage(const char * p,size_t sizeInBytes,int n,ergo_real * resultMatrix,int storageType)430 ddf_get_matrix_from_storage(const char* p, size_t sizeInBytes, int n,
431 ergo_real* resultMatrix, int storageType)
432 {
433 size_t sizeNeeded = 0;
434 switch(storageType)
435 {
436 case MATRIX_STORAGE_TYPE_FULL:
437 sizeNeeded = n*n*sizeof(ergo_real);
438 if(sizeNeeded != sizeInBytes)
439 {
440 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_from_storage: (sizeNeeded != sizeInBytes)");
441 return -1;
442 }
443 memcpy(resultMatrix, p, sizeNeeded);
444 break;
445 case MATRIX_STORAGE_TYPE_TRIANGLE:
446 sizeNeeded = (n * (n+1) / 2) * sizeof(ergo_real);
447 if(sizeNeeded != sizeInBytes)
448 {
449 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix: (sizeNeeded != sizeInBytes)");
450 return -1;
451 }
452 ddf_get_triangular_matrix_from_storage(p, n, resultMatrix);
453 break;
454 case MATRIX_STORAGE_TYPE_VECTORS:
455 if(ddf_get_matrix_from_vector_storage(p, sizeInBytes, n, resultMatrix) != 0)
456 {
457 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_from_vector_storage.");
458 return -1;
459 }
460 break;
461 default:
462 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_from_storage: unknown storage type %i", storageType);
463 return -1;
464 }
465 return 0;
466 }
467
468
469 static int
ddf_get_matrix_from_storage_sparse(const char * p,size_t sizeInBytes,int n,int * rowind,int * colind,ergo_real * values,int storageType)470 ddf_get_matrix_from_storage_sparse(const char* p,
471 size_t sizeInBytes,
472 int n,
473 int* rowind,
474 int* colind,
475 ergo_real* values,
476 int storageType)
477 {
478 size_t sizeNeeded = 0;
479 switch(storageType)
480 {
481 case MATRIX_STORAGE_TYPE_FULL:
482 {
483 sizeNeeded = n*n*sizeof(ergo_real);
484 if(sizeNeeded != sizeInBytes)
485 {
486 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_from_storage: (sizeNeeded != sizeInBytes)");
487 return -1;
488 }
489 ergo_real* A = (ergo_real*)p;
490 long count = 0;
491 for(long i = 0; i < n; i++)
492 for(long j = 0; j < n; j++)
493 {
494 rowind[count] = i;
495 colind[count] = j;
496 values[count] = A[i*n+j];
497 count++;
498 }
499 }
500 break;
501 case MATRIX_STORAGE_TYPE_TRIANGLE:
502 sizeNeeded = (n * (n+1) / 2) * sizeof(ergo_real);
503 if(sizeNeeded != sizeInBytes)
504 {
505 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix: (sizeNeeded != sizeInBytes)");
506 return -1;
507 }
508 ddf_get_triangular_matrix_from_storage_sparse(p, n, rowind, colind, values);
509 break;
510 case MATRIX_STORAGE_TYPE_VECTORS:
511 if(ddf_get_matrix_from_vector_storage_sparse(p, sizeInBytes, n, rowind, colind, values) != 0)
512 {
513 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_from_vector_storage.");
514 return -1;
515 }
516 break;
517 default:
518 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_from_storage: unknown storage type %i", storageType);
519 return -1;
520 }
521 return 0;
522 }
523
524
525 int
ddf_writeShellListAndDensityMatricesToFile(const BasisInfoStruct & basisInfo,int noOfDensityMatrices,ergo_real ** densityMatrixList,const char * fileName)526 ddf_writeShellListAndDensityMatricesToFile(const BasisInfoStruct & basisInfo,
527 int noOfDensityMatrices,
528 ergo_real** densityMatrixList,
529 const char* fileName) {
530 int n = basisInfo.noOfBasisFuncs;
531 int nShells = basisInfo.noOfShells;
532 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "ddf_writeShellListAndDensityMatricesToFile, n = %d, nShells = %d", n, nShells);
533 if(n <= 0 || nShells <= 0) {
534 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in writeShellListAndDensityMatricesToFile: "
535 "(n <= 0 || nShells <= 0)");
536 return -1;
537 }
538 if(noOfDensityMatrices <= 0 || noOfDensityMatrices > 2) {
539 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in writeShellListAndDensityMatricesToFile: "
540 "(noOfDensityMatrices <= 0 || noOfDensityMatrices > 2)");
541 return -1;
542 }
543 FILE* f = fopen(fileName, "wb");
544 if(f == NULL) {
545 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error opening file '%s' for writing", fileName);
546 return -1;
547 }
548
549 // decide which storage format to use
550 int matrixStorageType;
551 if(n < 22000) {
552 if(ddf_get_matrix_storage_size(MATRIX_STORAGE_TYPE_TRIANGLE, n, densityMatrixList[0]) < ddf_get_matrix_storage_size(MATRIX_STORAGE_TYPE_VECTORS, n, densityMatrixList[0])) {
553 matrixStorageType = MATRIX_STORAGE_TYPE_TRIANGLE;
554 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "Choosing matrixStorageType = MATRIX_STORAGE_TYPE_TRIANGLE.");
555 }
556 else
557 {
558 matrixStorageType = MATRIX_STORAGE_TYPE_VECTORS;
559 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "Choosing matrixStorageType = MATRIX_STORAGE_TYPE_VECTORS.");
560 }
561 }
562 else {
563 // triangle storage not possible, use vector storage.
564 matrixStorageType = MATRIX_STORAGE_TYPE_VECTORS;
565 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "Choosing matrixStorageType = MATRIX_STORAGE_TYPE_VECTORS.");
566 }
567
568 // Get sizes needed for matrix storage.
569 size_t matrixStorageSizeList[2];
570 matrixStorageSizeList[0] = 0;
571 matrixStorageSizeList[1] = 0;
572 for(int i = 0; i < noOfDensityMatrices; i++) {
573 size_t sizeInBytes = ddf_get_matrix_storage_size(matrixStorageType, n, densityMatrixList[i]);
574 if(sizeInBytes == 0) {
575 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_storage_size");
576 return -1;
577 }
578 matrixStorageSizeList[i] = sizeInBytes;
579 }
580
581 size_t matrixStorageSizeTotal = 0;
582 for(int i = 0; i < noOfDensityMatrices; i++)
583 matrixStorageSizeTotal += matrixStorageSizeList[i];
584
585 // Compute file size
586 size_t fileSize =
587 sizeof(densityFileHeaderStruct) +
588 nShells * sizeof(ShellSpecStruct) +
589 matrixStorageSizeTotal;
590 densityFileHeaderStruct fileHeader;
591
592 // Create header
593 fileHeader.densityFileVersion = DENSITY_FILE_VERSION_NUMBER;
594 fileHeader.typeOfMatrixStorage = matrixStorageType;
595 fileHeader.noOfShells = nShells;
596 fileHeader.noOfBasisFuncs = n;
597 fileHeader.noOfDensityMatrices = noOfDensityMatrices;
598 fileHeader.padding_int = 0; /* Set padding_int also to make valgrind happy. */
599 fileHeader.matrixSize_1 = matrixStorageSizeList[0];
600 fileHeader.matrixSize_2 = matrixStorageSizeList[1];
601 fileHeader.fileSizeInBytes = fileSize;
602
603 // prepare data to write
604 char* buffer = (char*)ergo_malloc(fileSize);
605 char* p = buffer;
606 memcpy(p, &fileHeader, sizeof(densityFileHeaderStruct));
607 p += sizeof(densityFileHeaderStruct);
608 memcpy(p, basisInfo.shellList, nShells * sizeof(ShellSpecStruct));
609 p += nShells * sizeof(ShellSpecStruct);
610 for(int i = 0; i < noOfDensityMatrices; i++) {
611 if(ddf_store_matrix(p, matrixStorageSizeList[i], n, densityMatrixList[i], matrixStorageType) != 0) {
612 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix");
613 return -1;
614 }
615 p += matrixStorageSizeList[i];
616 } // END FOR i
617 if(size_t(p - buffer) != fileSize) {
618 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error: ((p - buffer) != fileSize)");
619 return -1;
620 }
621
622 // write to file
623 size_t noOfBytesWritten = fwrite(buffer, 1, fileSize, f);
624 fclose(f);
625 if(noOfBytesWritten != fileSize) {
626 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in fwrite, fileSize = %i, noOfBytesWritten = %i",
627 fileSize, noOfBytesWritten);
628 return -1;
629 }
630
631 ergo_free(buffer);
632
633 // return 0 to indicate success
634 return 0;
635 }
636
637
638
639 int
ddf_writeShellListAndDensityMatricesToFile_sparse(const BasisInfoStruct & basisInfo,int noOfDensityMatrices,matrix_description_struct * densityMatrixList,const char * fileName)640 ddf_writeShellListAndDensityMatricesToFile_sparse(const BasisInfoStruct & basisInfo,
641 int noOfDensityMatrices,
642 matrix_description_struct* densityMatrixList,
643 const char* fileName)
644 {
645 int n = basisInfo.noOfBasisFuncs;
646
647 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "ddf_writeShellListAndDensityMatricesToFile_sparse, n = %6i, densityMatrixList[0].nvalues = %11i", n, densityMatrixList[0].nvalues);
648
649 int nShells = basisInfo.noOfShells;
650 if(n <= 0 || nShells <= 0) {
651 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in writeShellListAndDensityMatricesToFile: "
652 "(n <= 0 || nShells <= 0)");
653 return -1;
654 }
655 if(noOfDensityMatrices <= 0 || noOfDensityMatrices > 2) {
656 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in writeShellListAndDensityMatricesToFile: "
657 "(noOfDensityMatrices <= 0 || noOfDensityMatrices > 2)");
658 return -1;
659 }
660
661 /* We create a tmp file and rename it on success. This is useful
662 when a file of identical name exists: this method guarantees that
663 the created density file does not get corrupted when the process
664 is interrupted mid-stream. */
665 std::string tmpFName(fileName);
666 tmpFName += ".XXXXXX";
667 int tmpFD = mkstemp(&tmpFName[0]);
668 FILE* f;
669
670 if(tmpFD == -1 || (f = fdopen(tmpFD, "wb")) == NULL)
671 {
672 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN,
673 "error opening temporary file '%s' for writing",
674 tmpFName.c_str());
675 return -1;
676 }
677 do_output(LOG_CAT_INFO, LOG_AREA_MAIN,
678 "opened temporary file '%s' for writing",
679 tmpFName.c_str());
680
681 // decide which storage format to use
682 int matrixStorageType;
683 if(n < 22000)
684 {
685 if(ddf_get_matrix_storage_size_sparse(MATRIX_STORAGE_TYPE_TRIANGLE, n, &densityMatrixList[0]) < ddf_get_matrix_storage_size_sparse(MATRIX_STORAGE_TYPE_VECTORS, n, &densityMatrixList[0]))
686 {
687 matrixStorageType = MATRIX_STORAGE_TYPE_TRIANGLE;
688 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "Choosing matrixStorageType = MATRIX_STORAGE_TYPE_TRIANGLE.");
689 }
690 else
691 {
692 matrixStorageType = MATRIX_STORAGE_TYPE_VECTORS;
693 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "Choosing matrixStorageType = MATRIX_STORAGE_TYPE_VECTORS.");
694 }
695 }
696 else
697 {
698 // triangle storage not possible, use vector storage.
699 matrixStorageType = MATRIX_STORAGE_TYPE_VECTORS;
700 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "Choosing matrixStorageType = MATRIX_STORAGE_TYPE_VECTORS.");
701 }
702
703 // Get sizes needed for matrix storage.
704 size_t matrixStorageSizeList[2];
705 matrixStorageSizeList[0] = 0;
706 matrixStorageSizeList[1] = 0;
707 for(int i = 0; i < noOfDensityMatrices; i++)
708 {
709 size_t sizeInBytes = ddf_get_matrix_storage_size_sparse(matrixStorageType, n, &densityMatrixList[i]);
710 if(sizeInBytes <= 0)
711 {
712 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_storage_size_sparse");
713 return -1;
714 }
715 matrixStorageSizeList[i] = sizeInBytes;
716 }
717
718 size_t matrixStorageSizeTotal = 0;
719 for(int i = 0; i < noOfDensityMatrices; i++)
720 matrixStorageSizeTotal += matrixStorageSizeList[i];
721
722 // Compute file size
723 size_t fileSize =
724 sizeof(densityFileHeaderStruct) +
725 nShells * sizeof(ShellSpecStruct) +
726 matrixStorageSizeTotal;
727
728 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "fileSize = %15.0f bytes = %10.3f MegaBytes", (double)fileSize, (double)fileSize / 1000000);
729
730 densityFileHeaderStruct fileHeader;
731
732 // Create header
733 fileHeader.densityFileVersion = DENSITY_FILE_VERSION_NUMBER;
734 fileHeader.typeOfMatrixStorage = matrixStorageType;
735 fileHeader.noOfShells = nShells;
736 fileHeader.noOfBasisFuncs = n;
737 fileHeader.noOfDensityMatrices = noOfDensityMatrices;
738 fileHeader.padding_int = 0; /* Set padding_int also to make valgrind happy. */
739 fileHeader.matrixSize_1 = matrixStorageSizeList[0];
740 fileHeader.matrixSize_2 = matrixStorageSizeList[1];
741 fileHeader.fileSizeInBytes = fileSize;
742
743 // prepare data to write
744 char* buffer = (char*)ergo_malloc(fileSize);
745 char* p = buffer;
746 memcpy(p, &fileHeader, sizeof(densityFileHeaderStruct));
747 p += sizeof(densityFileHeaderStruct);
748 memcpy(p, basisInfo.shellList, nShells * sizeof(ShellSpecStruct));
749 p += nShells * sizeof(ShellSpecStruct);
750 for(int i = 0; i < noOfDensityMatrices; i++) {
751 if(ddf_store_matrix_sparse(p, matrixStorageSizeList[i], n, &densityMatrixList[i], matrixStorageType) != 0) {
752 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_store_matrix");
753 return -1;
754 }
755 p += matrixStorageSizeList[i];
756 } // END FOR i
757 if(size_t(p - buffer) != fileSize) {
758 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error: ((p - buffer) != fileSize)");
759 return -1;
760 }
761
762 // write to file
763 size_t noOfBytesWritten = fwrite(buffer, 1, fileSize, f);
764 if(fsync(fileno(f)) != 0) {
765 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in fsync on %s: %s",
766 fileName, strerror(errno));
767 return -1;
768 }
769 if(fclose(f) != 0) {
770 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in fclose on %s: %s",
771 fileName, strerror(errno));
772 return -1;
773 }
774
775 if(noOfBytesWritten != fileSize) {
776 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in fwrite, fileSize = %i, noOfBytesWritten = %i",
777 fileSize, noOfBytesWritten);
778 return -1;
779 }
780 if(rename(tmpFName.c_str(), fileName) != 0) {
781 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error renaming %s to %s: %s\n",
782 tmpFName.c_str(), fileName, strerror(errno));
783 return -1;
784 }
785
786 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "renamed file '%s' to '%s'",
787 tmpFName.c_str(), fileName);
788
789 /** Data loss was observed with large files on AFS, we do extra verification to detect it early on... */
790 struct stat fileStats;
791 if(stat(fileName, &fileStats) != 0) {
792 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "Cannot stat saved file %s: %s\n",
793 fileName, strerror(errno));
794 return -1;
795 }
796 if( size_t(fileStats.st_size) != fileSize) {
797 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "File %s size got altered from %lu to %lu\n",
798 fileName, (unsigned long)fileSize, (unsigned long)fileStats.st_size);
799 return -1;
800 }
801
802 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "ddf_writeShellListAndDensityMatricesToFile_sparse freeing buffer.");
803 ergo_free(buffer);
804
805 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "ddf_writeShellListAndDensityMatricesToFile_sparse returning OK.");
806 // return 0 to indicate success
807 return 0;
808 }
809
810
811
812 static int
ddf_load_density_getSizes(const char * fileName,int * result_noOfShells,int * result_noOfBasisFuncs,int * result_noOfDensitiesOnFile,long * result_noOfValuesList)813 ddf_load_density_getSizes(const char* fileName,
814 int* result_noOfShells,
815 int* result_noOfBasisFuncs,
816 int* result_noOfDensitiesOnFile,
817 long* result_noOfValuesList)
818 {
819 FILE* f = fopen(fileName, "rb");
820 if(f == NULL)
821 {
822 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error opening file '%s'", fileName);
823 return -1;
824 }
825
826 // check file size
827 fseeko(f, 0L, SEEK_END);
828 size_t fileSize = size_t(ftello(f));
829 if(fileSize < sizeof(densityFileHeaderStruct))
830 {
831 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN,
832 "error in ddf_load_density_getSizes: (fileSize < sizeof(densityFileHeaderStruct)) %lu < %lu",
833 (unsigned long)fileSize, (unsigned long)sizeof(densityFileHeaderStruct));
834 return -1;
835 }
836 fseeko(f, 0L, SEEK_SET);
837 char* buffer = (char*)ergo_malloc(fileSize);
838 if(fread(buffer, 1, fileSize, f) != fileSize)
839 {
840 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error reading file '%s'", fileName);
841 return -1;
842 }
843 fclose(f);
844
845 // OK, file contents are now in buffer.
846 densityFileHeaderStruct* headerPtr = (densityFileHeaderStruct*)buffer;
847 if(headerPtr->densityFileVersion != DENSITY_FILE_VERSION_NUMBER)
848 {
849 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_load_density_getSizes: (headerPtr->densityFileVersion != DENSITY_FILE_VERSION_NUMBER)");
850 return -1;
851 }
852 if(headerPtr->fileSizeInBytes != fileSize)
853 {
854 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN,
855 "error in ddf_load_density_getSizes: (headerPtr->fileSizeInBytes != fileSize) %lu %lu",
856 (unsigned long)headerPtr->fileSizeInBytes, (unsigned long)fileSize);
857 return -1;
858 }
859 int n = headerPtr->noOfBasisFuncs;
860 int nShells = headerPtr->noOfShells;
861
862 size_t fileSize_expected =
863 sizeof(densityFileHeaderStruct) +
864 nShells * sizeof(ShellSpecStruct) +
865 (size_t)headerPtr->matrixSize_1 + (size_t)headerPtr->matrixSize_2;
866 if(fileSize_expected != fileSize)
867 {
868 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_load_density_getSizes: (fileSize_expected != fileSize) %lu != %lu",
869 (unsigned long)fileSize_expected, (unsigned long)fileSize);
870 return -1;
871 }
872
873 *result_noOfShells = nShells;
874 *result_noOfBasisFuncs = n;
875 *result_noOfDensitiesOnFile = headerPtr->noOfDensityMatrices;
876 if(result_noOfValuesList != NULL)
877 {
878 // Get result_noOfValuesList.
879 char* p = buffer + sizeof(densityFileHeaderStruct);
880
881 // skip shell list
882 p += nShells*sizeof(ShellSpecStruct);
883
884 // get nvalues for stored matrices
885 size_t matrixSizeList[2];
886 matrixSizeList[0] = (size_t)headerPtr->matrixSize_1;
887 matrixSizeList[1] = (size_t)headerPtr->matrixSize_2;
888 for(int i = 0; i < headerPtr->noOfDensityMatrices; i++)
889 {
890 size_t currSize = matrixSizeList[i];
891 if(ddf_get_matrix_nvalues_from_storage(p, currSize, n, &result_noOfValuesList[i],
892 headerPtr->typeOfMatrixStorage) != 0)
893 {
894 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_nvalues_from_storage");
895 return -1;
896 }
897 p += currSize;
898 } // END FOR i
899
900 if(size_t(p - buffer) != fileSize)
901 {
902 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error: ((p - buffer) != fileSize)");
903 return -1;
904 }
905 }
906
907 ergo_free(buffer);
908 return 0;
909 }
910
911
912 /** ddf_read_shells_and_density_matrices() reads the basis set
913 information and requested number of density matrices from a
914 specified file. basisInfo needs to be allocated and zeroed in
915 advance. densityMatrixList must be properly allocated as well. */
916 static int
ddf_read_shells_and_density_matrices(BasisInfoStruct * basisInfo,int noOfDensityMatrices,ergo_real ** densityMatrixList,const char * fileName)917 ddf_read_shells_and_density_matrices(BasisInfoStruct* basisInfo,
918 int noOfDensityMatrices,
919 ergo_real** densityMatrixList,
920 const char* fileName)
921 {
922 FILE* f = fopen(fileName, "rb");
923 if(f == NULL)
924 {
925 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error opening file '%s'", fileName);
926 return -1;
927 }
928
929 // check file size
930 fseeko(f, 0L, SEEK_END);
931 size_t fileSize = size_t(ftello(f));
932 if(fileSize < sizeof(densityFileHeaderStruct))
933 {
934 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_read_shells_and_density_matrices: (fileSize < sizeof(densityFileHeaderStruct))");
935 return -1;
936 }
937 fseeko(f, 0L, SEEK_SET);
938 char* buffer = (char*)ergo_malloc(fileSize);
939 size_t readsz = fread(buffer, 1, fileSize, f);
940 fclose(f);
941 if(readsz != fileSize)
942 {
943 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error reading file '%s'", fileName);
944 ergo_free(buffer);
945 return -1;
946 }
947
948 // OK, file contents are now in buffer.
949 densityFileHeaderStruct* headerPtr = (densityFileHeaderStruct*)buffer;
950 if(headerPtr->densityFileVersion != DENSITY_FILE_VERSION_NUMBER)
951 {
952 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_read_shells_and_density_matrices: (headerPtr->densityFileVersion != DENSITY_FILE_VERSION_NUMBER)");
953 return -1;
954 }
955 if(headerPtr->fileSizeInBytes != fileSize)
956 {
957 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_read_shells_and_density_matrices: (headerPtr->fileSizeInBytes != fileSize)");
958 return -1;
959 }
960 if(headerPtr->noOfDensityMatrices != noOfDensityMatrices)
961 {
962 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_read_shells_and_density_matrices: (headerPtr->noOfDensityMatrices != noOfDensityMatrices)");
963 return -1;
964 }
965 int n = headerPtr->noOfBasisFuncs;
966 int nShells = headerPtr->noOfShells;
967
968 size_t fileSize_expected =
969 sizeof(densityFileHeaderStruct) +
970 nShells * sizeof(ShellSpecStruct) +
971 (size_t)headerPtr->matrixSize_1 + (size_t)headerPtr->matrixSize_2;
972 if(fileSize_expected != fileSize)
973 {
974 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_read_shells_and_density_matrices: (fileSize_expected != fileSize)");
975 return -1;
976 }
977
978 char* p = buffer + sizeof(densityFileHeaderStruct);
979
980 // get shell list
981 basisInfo->noOfShells = nShells;
982 memcpy(basisInfo->shellList, p, nShells*sizeof(ShellSpecStruct));
983 p += nShells*sizeof(ShellSpecStruct);
984
985 // get density matrices
986 size_t matrixSizeList[2];
987 matrixSizeList[0] = (size_t)headerPtr->matrixSize_1;
988 matrixSizeList[1] = (size_t)headerPtr->matrixSize_2;
989 for(int i = 0; i < noOfDensityMatrices; i++)
990 {
991 long currSize = matrixSizeList[i];
992 if(ddf_get_matrix_from_storage(p, currSize, n, densityMatrixList[i],
993 headerPtr->typeOfMatrixStorage) != 0)
994 {
995 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_from_storage");
996 return -1;
997 }
998 p += currSize;
999 } // END FOR i
1000
1001 if(size_t(p - buffer) != fileSize)
1002 {
1003 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error: ((p - buffer) != fileSize)");
1004 return -1;
1005 }
1006
1007 ergo_free(buffer);
1008 return 0;
1009 }
1010
1011
1012
1013 static int
ddf_read_shells_and_density_matrices_sparse(BasisInfoStruct ** basisInfo,int noOfDensityMatrices,int ** rowindList,int ** colindList,ergo_real ** valuesList,const char * fileName)1014 ddf_read_shells_and_density_matrices_sparse(BasisInfoStruct** basisInfo,
1015 int noOfDensityMatrices,
1016 int** rowindList,
1017 int** colindList,
1018 ergo_real** valuesList,
1019 const char* fileName)
1020 {
1021 FILE* f = fopen(fileName, "rb");
1022 if(f == NULL)
1023 {
1024 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error opening file '%s'", fileName);
1025 return -1;
1026 }
1027
1028 // check file size
1029 fseeko(f, 0L, SEEK_END);
1030 size_t fileSize = size_t(ftello(f));
1031 if(fileSize < sizeof(densityFileHeaderStruct))
1032 {
1033 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_read_shells_and_density_matrices: (fileSize < sizeof(densityFileHeaderStruct))");
1034 return -1;
1035 }
1036 fseeko(f, 0L, SEEK_SET);
1037 char* buffer = (char*)ergo_malloc(fileSize);
1038 size_t readsz = fread(buffer, 1, fileSize, f);
1039 fclose(f);
1040 if(readsz != fileSize)
1041 {
1042 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error reading file '%s'", fileName);
1043 ergo_free(buffer);
1044 return -1;
1045 }
1046
1047 // OK, file contents are now in buffer.
1048 densityFileHeaderStruct* headerPtr = (densityFileHeaderStruct*)buffer;
1049 if(headerPtr->densityFileVersion != DENSITY_FILE_VERSION_NUMBER)
1050 {
1051 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_read_shells_and_density_matrices: (headerPtr->densityFileVersion != DENSITY_FILE_VERSION_NUMBER)");
1052 return -1;
1053 }
1054 if(headerPtr->fileSizeInBytes != fileSize)
1055 {
1056 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_read_shells_and_density_matrices: (headerPtr->fileSizeInBytes != fileSize)");
1057 return -1;
1058 }
1059 if(headerPtr->noOfDensityMatrices != noOfDensityMatrices)
1060 {
1061 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_read_shells_and_density_matrices: (headerPtr->noOfDensityMatrices != noOfDensityMatrices)");
1062 return -1;
1063 }
1064 int n = headerPtr->noOfBasisFuncs;
1065 int nShells = headerPtr->noOfShells;
1066
1067 size_t fileSize_expected =
1068 sizeof(densityFileHeaderStruct) +
1069 nShells * sizeof(ShellSpecStruct) +
1070 (size_t)headerPtr->matrixSize_1 + (size_t)headerPtr->matrixSize_2;
1071 if(fileSize_expected != fileSize)
1072 {
1073 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_read_shells_and_density_matrices: (fileSize_expected != fileSize)");
1074 return -1;
1075 }
1076
1077 char* p = buffer + sizeof(densityFileHeaderStruct);
1078
1079 // get shell list
1080 ShellSpecStruct* shellList = new ShellSpecStruct[nShells];
1081 memcpy(shellList, p, nShells*sizeof(ShellSpecStruct));
1082 // Now check if use_6_d_funcs was used for the basis set.
1083 int use_6_d_funcs = 0;
1084 for(int i = 0; i < nShells; i++) {
1085 if(shellList[i].noOfBasisFuncs == 6)
1086 use_6_d_funcs = 1;
1087 }
1088 // Now create basisInfo object.
1089 *basisInfo = new BasisInfoStruct(use_6_d_funcs);
1090 (*basisInfo)->noOfShells = nShells;
1091 (*basisInfo)->shellList = shellList;
1092 // Move pointer past shell list
1093 p += nShells*sizeof(ShellSpecStruct);
1094
1095 // get density matrices
1096 size_t matrixSizeList[2];
1097 matrixSizeList[0] = (size_t)headerPtr->matrixSize_1;
1098 matrixSizeList[1] = (size_t)headerPtr->matrixSize_2;
1099 for(int i = 0; i < noOfDensityMatrices; i++)
1100 {
1101 size_t currSize = matrixSizeList[i];
1102 if(ddf_get_matrix_from_storage_sparse(p, currSize, n, rowindList[i], colindList[i], valuesList[i],
1103 headerPtr->typeOfMatrixStorage) != 0)
1104 {
1105 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_get_matrix_from_storage_sparse");
1106 return -1;
1107 }
1108 p += currSize;
1109 } // END FOR i
1110
1111 if(size_t(p - buffer) != fileSize)
1112 {
1113 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error: ((p - buffer) != fileSize)");
1114 return -1;
1115 }
1116
1117 ergo_free(buffer);
1118 return 0;
1119 }
1120
1121
1122
1123
1124 int
ddf_load_density(const char * densityFileName,int noOfDensityMatrices,const IntegralInfo & integralInfo,BasisInfoStruct ** basisInfo,ergo_real ** densityMatrix)1125 ddf_load_density(const char *densityFileName,
1126 int noOfDensityMatrices,
1127 const IntegralInfo& integralInfo,
1128 BasisInfoStruct **basisInfo,
1129 ergo_real **densityMatrix)
1130 {
1131 do_output(LOG_CAT_INFO, LOG_AREA_MAIN,
1132 "getting starting guess from file, noOfDensityMatrices = %i",
1133 noOfDensityMatrices);
1134
1135 // get basisInfo and density matrix (or matrices) for starting guess
1136 // FIXME: here it is assumed that use_6_d_funcs was not used for the saved density.
1137 int use_6_d_funcs = 0;
1138 *basisInfo = new BasisInfoStruct(use_6_d_funcs);
1139
1140 int noOfShells_sg = 0;
1141 int noOfBasisFuncs_sg = 0;
1142 int noOfDensitiesOnFile = 0;
1143 if(ddf_load_density_getSizes(densityFileName,
1144 &noOfShells_sg,
1145 &noOfBasisFuncs_sg,
1146 &noOfDensitiesOnFile,
1147 NULL) != 0)
1148 {
1149 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_load_density_getSizes");
1150 return -1;
1151 }
1152 if(noOfDensitiesOnFile != noOfDensityMatrices)
1153 {
1154 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN,
1155 "error: (headerPtr->noOfDensityMatrices != noOfDensityMatrices)");
1156 return -1;
1157 }
1158
1159
1160 (*basisInfo)->shellList = new ShellSpecStruct[noOfShells_sg];
1161 for(int i = 0; i < noOfDensityMatrices; i++)
1162 densityMatrix[i] =
1163 ergo_new(noOfBasisFuncs_sg*noOfBasisFuncs_sg,ergo_real);
1164
1165 if(ddf_read_shells_and_density_matrices(*basisInfo,
1166 noOfDensityMatrices,
1167 densityMatrix,
1168 densityFileName) != 0)
1169 {
1170 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_load_density");
1171 return -1;
1172 }
1173 if((*basisInfo)->get_basis_funcs() != 0)
1174 {
1175 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in get_basis_funcs");
1176 return -1;
1177 }
1178 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "get_basis_funcs for starting guess basis set returned OK,"
1179 " number of basis funcs: %i",
1180 (*basisInfo)->noOfBasisFuncs);
1181 if((*basisInfo)->getSimplePrimitivesAll(integralInfo) != 0)
1182 {
1183 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in getSimplePrimitivesAll");
1184 return -1;
1185 }
1186 do_output(LOG_CAT_INFO, LOG_AREA_MAIN,
1187 "getSimplePrimitivesAll for starting guess basis set "
1188 "returned OK, n = %i",
1189 (*basisInfo)->noOfSimplePrimitives);
1190 return 0;
1191 }
1192
1193
1194
1195 int
ddf_load_density_sparse(const char * densityFileName,const IntegralInfo & integralInfo,BasisInfoStruct ** basisInfo,int * noOfDensitiesRead,int ** rowindList,int ** colindList,ergo_real ** valuesList,long * nvaluesList)1196 ddf_load_density_sparse(const char *densityFileName,
1197 const IntegralInfo& integralInfo,
1198 BasisInfoStruct **basisInfo,
1199 int* noOfDensitiesRead,
1200 int** rowindList,
1201 int** colindList,
1202 ergo_real** valuesList,
1203 long* nvaluesList)
1204 {
1205 do_output(LOG_CAT_INFO, LOG_AREA_MAIN,
1206 "ddf_load_density_sparse(\"%s\")", densityFileName);
1207
1208 // get basisInfo and density matrix (or matrices) for starting guess
1209 // FIXME: here it is assumed that use_6_d_funcs was not used for the saved density.
1210 int use_6_d_funcs = 0;
1211 *basisInfo = new BasisInfoStruct(use_6_d_funcs);
1212
1213 int noOfShells_sg = 0;
1214 int noOfBasisFuncs_sg = 0;
1215 if(ddf_load_density_getSizes(densityFileName,
1216 &noOfShells_sg,
1217 &noOfBasisFuncs_sg,
1218 noOfDensitiesRead,
1219 nvaluesList) != 0)
1220 {
1221 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_load_density_getSizes");
1222 return -1;
1223 }
1224
1225 for(int i = 0; i < *noOfDensitiesRead; i++)
1226 {
1227 long nvalues = nvaluesList[i];
1228 rowindList[i] = new int[nvalues];
1229 colindList[i] = new int[nvalues];
1230 valuesList[i] = new ergo_real[nvalues];
1231 }
1232
1233 if(ddf_read_shells_and_density_matrices_sparse(basisInfo,
1234 *noOfDensitiesRead,
1235 rowindList,
1236 colindList,
1237 valuesList,
1238 densityFileName) != 0)
1239 {
1240 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in ddf_read_shells_and_density_matrices_sparse");
1241 return -1;
1242 }
1243
1244 if((*basisInfo)->get_basis_funcs() != 0)
1245 {
1246 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in get_basis_funcs");
1247 return -1;
1248 }
1249 do_output(LOG_CAT_INFO, LOG_AREA_MAIN, "get_basis_funcs for starting guess basis set returned OK,"
1250 " number of basis funcs: %i",
1251 (*basisInfo)->noOfBasisFuncs);
1252 if((*basisInfo)->getSimplePrimitivesAll(integralInfo) != 0)
1253 {
1254 do_output(LOG_CAT_ERROR, LOG_AREA_MAIN, "error in getSimplePrimitivesAll");
1255 return -1;
1256 }
1257 do_output(LOG_CAT_INFO, LOG_AREA_MAIN,
1258 "getSimplePrimitivesAll for starting guess basis set "
1259 "returned OK, n = %i",
1260 (*basisInfo)->noOfSimplePrimitives);
1261 return 0;
1262 }
1263