CheMPS2
FourIndex.cpp
1 /*
2  CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3  Copyright (C) 2013-2016 Sebastian Wouters
4 
5  This program is free software; you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation; either version 2 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License along
16  with this program; if not, write to the Free Software Foundation, Inc.,
17  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19 
20 #include <stdlib.h>
21 #include <assert.h>
22 #include <iostream>
23 #include <sstream>
24 #include <string>
25 #include <climits>
26 
27 #include "FourIndex.h"
28 #include "Lapack.h"
29 #include "MyHDF5.h"
30 #include "MPIchemps2.h"
31 
32 using namespace std;
33 
34 CheMPS2::FourIndex::FourIndex(const int nGroup, const int * IrrepSizes){
35 
36  SymmInfo.setGroup(nGroup);
37 
38  Isizes = new int[SymmInfo.getNumberOfIrreps()];
39  for (int Icenter=0; Icenter<SymmInfo.getNumberOfIrreps(); Icenter++){
40  Isizes[Icenter] = IrrepSizes[Icenter];
41  }
42 
43  arrayLength = calcNumberOfUniqueElements(true); //true means allocate the storage!
44  theElements = new double[arrayLength];
45 
46  Clear();
47 
48 }
49 
50 long long CheMPS2::FourIndex::calcNumberOfUniqueElements(const bool allocate){
51 
52  //The object size: see text above storage in Fourindex.h
53  long long theTotalSize = 0;
54 
55  if (allocate){ storage = new long long****[SymmInfo.getNumberOfIrreps()]; }
56  for (int Icenter=0; Icenter<SymmInfo.getNumberOfIrreps(); Icenter++){
57  if (allocate){ storage[Icenter] = new long long***[SymmInfo.getNumberOfIrreps()]; }
58  for (int I_i=0; I_i<SymmInfo.getNumberOfIrreps(); I_i++){
59  int I_j = Irreps::directProd(Icenter,I_i);
60  if ((Isizes[I_i]>0)&&(Isizes[I_j]>0)){
61  if (allocate){ storage[Icenter][I_i] = new long long**[SymmInfo.getNumberOfIrreps()]; }
62  for (int I_k=I_i; I_k<SymmInfo.getNumberOfIrreps(); I_k++){
63  int I_l = Irreps::directProd(Icenter,I_k);
64  if ((Isizes[I_k]>0)&&(Isizes[I_l]>0)){
65  if ((I_i <= I_j) && (I_j <= I_l)){
66  if (Icenter == 0){ // I_i = I_j and I_k = I_l
67  if (I_i == I_k){
68  if (allocate){ storage[Icenter][I_i][I_k] = new long long*[Isizes[I_i]*(Isizes[I_i]+1)/2]; }
69  for (int i=0; i<Isizes[I_i]; i++){
70  for (int k=i; k<Isizes[I_k]; k++){
71  if (allocate){
72  storage[Icenter][I_i][I_k][i + k*(k+1)/2] = new long long[Isizes[I_j]-i];
73  for (int j=i; j<Isizes[I_j]; j++){
74  storage[Icenter][I_i][I_k][i + k*(k+1)/2][j-i] = theTotalSize;
75  theTotalSize += Isizes[I_l] - ((i==j) ? k : j);
76  }
77  } else {
78  delete [] storage[Icenter][I_i][I_k][i + k*(k+1)/2];
79  }
80  }
81  }
82  if (!allocate){ delete [] storage[Icenter][I_i][I_k]; }
83  } else { // I_i < I_k
84  if (allocate){ storage[Icenter][I_i][I_k] = new long long*[Isizes[I_i]*Isizes[I_k]]; }
85  for (int i=0; i<Isizes[I_i]; i++){
86  for (int k=0; k<Isizes[I_k]; k++){
87  if (allocate){
88  storage[Icenter][I_i][I_k][i + k*Isizes[I_i]] = new long long[Isizes[I_j]-i];
89  for (int j=i; j<Isizes[I_j]; j++){
90  storage[Icenter][I_i][I_k][i+k*Isizes[I_i]][j-i] = theTotalSize;
91  theTotalSize += Isizes[I_l] - ((i==j) ? k : 0 );
92  }
93  } else {
94  delete [] storage[Icenter][I_i][I_k][i + k*Isizes[I_i]];
95  }
96  }
97  }
98  if (!allocate){ delete [] storage[Icenter][I_i][I_k]; }
99  }
100  } else { //Icenter !=0 ; I_i < I_j and I_k != I_l
101  if (I_i == I_k){
102  if (allocate){ storage[Icenter][I_i][I_k] = new long long*[Isizes[I_i]*(Isizes[I_i]+1)/2]; }
103  for (int i=0; i<Isizes[I_i]; i++){
104  for (int k=i; k<Isizes[I_k]; k++){
105  if (allocate){
106  storage[Icenter][I_i][I_k][i + k*(k+1)/2] = new long long[Isizes[I_j]];
107  for (int j=0; j<Isizes[I_j]; j++){
108  storage[Icenter][I_i][I_k][i + k*(k+1)/2][j] = theTotalSize;
109  theTotalSize += Isizes[I_l] - j;
110  }
111  } else {
112  delete [] storage[Icenter][I_i][I_k][i + k*(k+1)/2];
113  }
114  }
115  }
116  if (!allocate){ delete [] storage[Icenter][I_i][I_k]; }
117  } else { // I_i < I_k
118  if (allocate){ storage[Icenter][I_i][I_k] = new long long*[Isizes[I_i]*Isizes[I_k]]; }
119  for (int i=0; i<Isizes[I_i]; i++){
120  for (int k=0; k<Isizes[I_k]; k++){
121  if (allocate){
122  storage[Icenter][I_i][I_k][i + k*Isizes[I_i]] = new long long[Isizes[I_j]];
123  for (int j=0; j<Isizes[I_j]; j++){
124  storage[Icenter][I_i][I_k][i + k*Isizes[I_i]][j] = theTotalSize;
125  theTotalSize += Isizes[I_l];
126  }
127  } else {
128  delete [] storage[Icenter][I_i][I_k][i + k*Isizes[I_i]];
129  }
130  }
131  }
132  if (!allocate){ delete [] storage[Icenter][I_i][I_k]; }
133  }
134  }
135  }
136  }
137  }
138  if (!allocate){ delete [] storage[Icenter][I_i]; }
139  }
140  }
141  if (!allocate){ delete [] storage[Icenter]; }
142  }
143  if (!allocate){ delete [] storage; }
144 
145  return theTotalSize;
146 
147 }
148 
150 
151  arrayLength = calcNumberOfUniqueElements(false); //false means delete the storage!
152  delete [] theElements;
153  delete [] Isizes;
154 
155 }
156 
158 
159  for (long long count = 0; count < arrayLength; count++){ theElements[ count ] = 0.0; }
160 
161 }
162 
163 void CheMPS2::FourIndex::set(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l, const double val){
164 
165  theElements[getPointer(irrep_i, irrep_j, irrep_k, irrep_l, i, j, k, l)] = val;
166 
167 }
168 
169 void CheMPS2::FourIndex::add(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l, const double val){
170 
171  theElements[getPointer(irrep_i, irrep_j, irrep_k, irrep_l, i, j, k, l)] += val;
172 
173 }
174 
175 double CheMPS2::FourIndex::get(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l) const{
176 
177  return theElements[getPointer(irrep_i, irrep_j, irrep_k, irrep_l, i, j, k, l)];
178 
179 }
180 
181 long long CheMPS2::FourIndex::getPointer(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l) const {
182 
183  assert( Irreps::directProd( irrep_i, irrep_j ) == Irreps::directProd( irrep_k, irrep_l ) );
184 
185  if ((irrep_i <= irrep_j) && (irrep_i <= irrep_k) && (irrep_j <= irrep_l)){ // (ijkl irrep ordering)
186  return getPtrIrrepOrderOK(irrep_i,irrep_j,irrep_k,irrep_l,i,j,k,l);
187  }
188 
189  if ((irrep_j <= irrep_i) && (irrep_j <= irrep_l) && (irrep_i <= irrep_k)){ // (jilk irrep ordering)
190  return getPtrIrrepOrderOK(irrep_j,irrep_i,irrep_l,irrep_k,j,i,l,k);
191  }
192 
193  if ((irrep_k <= irrep_j) && (irrep_k <= irrep_i) && (irrep_j <= irrep_l)){ // (kjil irrep ordering)
194  return getPtrIrrepOrderOK(irrep_k,irrep_j,irrep_i,irrep_l,k,j,i,l);
195  }
196 
197  if ((irrep_j <= irrep_k) && (irrep_j <= irrep_l) && (irrep_k <= irrep_i)){ // (jkli irrep ordering)
198  return getPtrIrrepOrderOK(irrep_j,irrep_k,irrep_l,irrep_i,j,k,l,i);
199  }
200 
201  if ((irrep_i <= irrep_l) && (irrep_i <= irrep_k) && (irrep_l <= irrep_j)){ // (ilkj irrep ordering)
202  return getPtrIrrepOrderOK(irrep_i,irrep_l,irrep_k,irrep_j,i,l,k,j);
203  }
204 
205  if ((irrep_l <= irrep_i) && (irrep_l <= irrep_j) && (irrep_i <= irrep_k)){ // (lijk irrep ordering)
206  return getPtrIrrepOrderOK(irrep_l,irrep_i,irrep_j,irrep_k,l,i,j,k);
207  }
208 
209  if ((irrep_k <= irrep_l) && (irrep_k <= irrep_i) && (irrep_l <= irrep_j)){ // (klij irrep ordering)
210  return getPtrIrrepOrderOK(irrep_k,irrep_l,irrep_i,irrep_j,k,l,i,j);
211  }
212 
213  if ((irrep_l <= irrep_k) && (irrep_l <= irrep_j) && (irrep_k <= irrep_i)){ // (lkji irrep ordering)
214  return getPtrIrrepOrderOK(irrep_l,irrep_k,irrep_j,irrep_i,l,k,j,i);
215  }
216 
217  return -1;
218 
219 }
220 
221 long long CheMPS2::FourIndex::getPtrIrrepOrderOK(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l) const {
222 
223  //I_i <= I_j <= I_l and I_i <= I_k
224  int Icenter = Irreps::directProd(irrep_i,irrep_j);
225 
226  if (Icenter>0){ // I_i < I_j and I_k != I_l
227 
228  if (irrep_i == irrep_k){ //and hence I_j == I_l
229 
230  if (k>=i){
231  if (l>=j) return getPtrAllOK5(Icenter, irrep_i, irrep_k, i, j, k, l);
232  else return getPtrAllOK5(Icenter, irrep_i, irrep_k, i, l, k, j);
233  } else {
234  if (l>=j) return getPtrAllOK5(Icenter, irrep_k, irrep_i, k, j, i, l);
235  else return getPtrAllOK5(Icenter, irrep_k, irrep_i, k, l, i, j);
236  }
237 
238  } else { return storage[Icenter][irrep_i][irrep_k][i + Isizes[irrep_i]*k][j] + l; }
239 
240  } else {
241 
242  if (irrep_i == irrep_k){ // all irreps the same
243 
244  //i en j
245  if ((i < j) && (i <= k) && (j <= l)) return getPtrAllOK2(Icenter, irrep_i, irrep_k, i, j, k, l); // (ijkl ordering)
246  if ((i == j) && (i <= k) && (j <= l)){
247  if (l>=k) return getPtrAllOK1(Icenter, irrep_i, irrep_k, i, j, k, l); // (ijkl ordering)
248  else return getPtrAllOK1(Icenter, irrep_j, irrep_l, j, i, l, k); // (jilk ordering)
249  }
250  if ((j < i) && (j <= l) && (i <= k)) return getPtrAllOK2(Icenter, irrep_j, irrep_l, j, i, l, k); // (jilk ordering)
251 
252  //k en l
253  if ((k < l) && (k <= i) && (l <= j)) return getPtrAllOK2(Icenter, irrep_k, irrep_i, k, l, i, j); // (klij ordering)
254  if ((k == l) && (k <= i) && (l <= j)){
255  if (j>=i) return getPtrAllOK1(Icenter, irrep_k, irrep_i, k, l, i, j); // (klij ordering)
256  else return getPtrAllOK1(Icenter, irrep_l, irrep_j, l, k, j, i); // (lkji ordering)
257  }
258  if ((l < k) && (l <= j) && (k <= i)) return getPtrAllOK2(Icenter, irrep_l, irrep_j, l, k, j, i); // (lkji ordering)
259 
260  // k en j
261  if ((k < j) && (k <= i) && (j <= l)) return getPtrAllOK2(Icenter, irrep_k, irrep_i, k, j, i, l); // (kjil ordering)
262  if ((k == j) && (k <= i) && (j <= l)){
263  if (l>=i) return getPtrAllOK1(Icenter, irrep_k, irrep_i, k, j, i, l); // (kjil ordering)
264  else return getPtrAllOK1(Icenter, irrep_j, irrep_l, j, k, l, i); // (jkli ordering)
265  }
266  if ((j < k) && (j <= l) && (k <= i)) return getPtrAllOK2(Icenter, irrep_j, irrep_l, j, k, l, i); // (jkli ordering)
267 
268  // i en l
269  if ((i < l) && (i <= k) && (l <= j)) return getPtrAllOK2(Icenter, irrep_i, irrep_k, i, l, k, j); // (ilkj ordering)
270  if ((i == l) && (i <= k) && (l <= j)){
271  if (j>=k) return getPtrAllOK1(Icenter, irrep_i, irrep_k, i, l, k, j); // (ilkj ordering)
272  else return getPtrAllOK1(Icenter, irrep_l, irrep_j, l, i, j, k); // (lijk ordering)
273  }
274  if ((l < i) && (l <= j) && (i <= k)) return getPtrAllOK2(Icenter, irrep_l, irrep_j, l, i, j, k); // (lijk ordering)
275 
276  } else {
277 
278  if (j==i){
279  if (l>=k){ return storage[Icenter][irrep_i][irrep_k][i + Isizes[irrep_i]*k][j-i] + l-k; }
280  else { return storage[Icenter][irrep_j][irrep_l][j + Isizes[irrep_j]*l][i-j] + k-l; }
281  } else {
282  if (j>i){ return storage[Icenter][irrep_i][irrep_k][i + Isizes[irrep_i]*k][j-i] + l; }
283  else { return storage[Icenter][irrep_j][irrep_l][j + Isizes[irrep_j]*l][i-j] + k; }
284  }
285 
286  }
287 
288  }
289 
290  return -1;
291 
292 }
293 
294 int CheMPS2::FourIndex::get_irrep_size( const int irrep ) const{ return Isizes[ irrep ]; }
295 
296 long long CheMPS2::FourIndex::getPtrAllOK1(const int Icent, const int irrep_i, const int irrep_k, const int i, const int j, const int k, const int l) const{
297 
298  return storage[Icent][irrep_i][irrep_k][i + k*(k+1)/2][j-i] + l-k;
299 
300 }
301 
302 long long CheMPS2::FourIndex::getPtrAllOK2(const int Icent, const int irrep_i, const int irrep_k, const int i, const int j, const int k, const int l) const{
303 
304  return storage[Icent][irrep_i][irrep_k][i + k*(k+1)/2][j-i] + l-j;
305 
306 }
307 
308 long long CheMPS2::FourIndex::getPtrAllOK5(const int Icent, const int irrep_i, const int irrep_k, const int i, const int j, const int k, const int l) const{
309 
310  return storage[Icent][irrep_i][irrep_k][i + k*(k+1)/2][j] + l-j;
311 
312 }
313 
314 void CheMPS2::FourIndex::save(const std::string name) const{
315 
316  //The hdf5 file
317  hid_t file_id = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
318 
319  //The metadata
320  hid_t group_id = H5Gcreate(file_id, "/MetaData", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
321 
322  //The IrrepSizes
323  hsize_t dimarray = SymmInfo.getNumberOfIrreps();
324  hid_t dataspace_id = H5Screate_simple(1, &dimarray, NULL);
325  hid_t dataset_id = H5Dcreate(group_id, "IrrepSizes", H5T_STD_I32LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
326  H5Dwrite(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, Isizes);
327 
328  //Attributes
329  hid_t attribute_space_id1 = H5Screate(H5S_SCALAR);
330  hid_t attribute_id1 = H5Acreate(dataset_id, "nGroup", H5T_STD_I32LE, attribute_space_id1, H5P_DEFAULT, H5P_DEFAULT);
331  int nGroup = SymmInfo.getGroupNumber();
332  H5Awrite(attribute_id1, H5T_NATIVE_INT, &nGroup);
333 
334  hid_t attribute_space_id2 = H5Screate(H5S_SCALAR);
335  hid_t attribute_id2 = H5Acreate(dataset_id, "nIrreps", H5T_STD_I32LE, attribute_space_id2, H5P_DEFAULT, H5P_DEFAULT);
336  int nIrreps = SymmInfo.getNumberOfIrreps();
337  H5Awrite(attribute_id2, H5T_NATIVE_INT, &nIrreps);
338 
339  hid_t attribute_space_id3 = H5Screate(H5S_SCALAR);
340  hid_t attribute_id3 = H5Acreate(dataset_id, "theTotalSize", H5T_STD_I64LE, attribute_space_id3, H5P_DEFAULT, H5P_DEFAULT);
341  H5Awrite(attribute_id3, H5T_NATIVE_LLONG, &arrayLength);
342 
343  H5Aclose(attribute_id1);
344  H5Aclose(attribute_id2);
345  H5Aclose(attribute_id3);
346  H5Sclose(attribute_space_id1);
347  H5Sclose(attribute_space_id2);
348  H5Sclose(attribute_space_id3);
349 
350  H5Dclose(dataset_id);
351  H5Sclose(dataspace_id);
352 
353  H5Gclose(group_id);
354 
355  //The object itself
356  hid_t group_id7 = H5Gcreate(file_id, "/FourIndexObject", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
357 
358  hsize_t dimarray7 = arrayLength; //hsize_t is defined by default as unsigned long long, so no problem
359  hid_t dataspace_id7 = H5Screate_simple(1, &dimarray7, NULL);
360  hid_t dataset_id7 = H5Dcreate(group_id7, "Matrix elements", H5T_IEEE_F64LE, dataspace_id7, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
361  H5Dwrite(dataset_id7, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, theElements);
362 
363  H5Dclose(dataset_id7);
364  H5Sclose(dataspace_id7);
365 
366  H5Gclose(group_id7);
367 
368  H5Fclose(file_id);
369 
370 }
371 
372 void CheMPS2::FourIndex::read(const std::string name){
373 
374  //The hdf5 file
375  hid_t file_id = H5Fopen(name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
376 
377  //The metadata
378  hid_t group_id = H5Gopen(file_id, "/MetaData",H5P_DEFAULT);
379 
380  //The IrrepSizes
381  hid_t dataset_id = H5Dopen(group_id, "IrrepSizes", H5P_DEFAULT);
382 
383  //Attributes
384  hid_t attribute_id1 = H5Aopen_by_name(group_id,"IrrepSizes", "nGroup", H5P_DEFAULT, H5P_DEFAULT);
385  int nGroup;
386  H5Aread(attribute_id1, H5T_NATIVE_INT, &nGroup);
387  assert( nGroup==SymmInfo.getGroupNumber() );
388 
389  hid_t attribute_id2 = H5Aopen_by_name(group_id,"IrrepSizes", "nIrreps", H5P_DEFAULT, H5P_DEFAULT);
390  int nIrreps;
391  H5Aread(attribute_id2, H5T_NATIVE_INT, &nIrreps);
392  assert( nIrreps==SymmInfo.getNumberOfIrreps() );
393 
394  hid_t attribute_id3 = H5Aopen_by_name(group_id,"IrrepSizes", "theTotalSize", H5P_DEFAULT, H5P_DEFAULT);
395  long long theTotalSize;
396  H5Aread(attribute_id3, H5T_NATIVE_LLONG, &theTotalSize);
397  assert( theTotalSize==arrayLength );
398 
399  H5Aclose(attribute_id1);
400  H5Aclose(attribute_id2);
401  H5Aclose(attribute_id3);
402 
403  int * IsizesAgain = new int[SymmInfo.getNumberOfIrreps()];
404  H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, IsizesAgain);
405  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++){
406  assert( IsizesAgain[cnt]==Isizes[cnt] );
407  }
408  delete [] IsizesAgain;
409  H5Dclose(dataset_id);
410 
411  H5Gclose(group_id);
412 
413  std::cout << "FourIndex::read : loading " << arrayLength << " doubles." << std::endl;
414 
415  //The object itself.
416  hid_t group_id7 = H5Gopen(file_id, "/FourIndexObject", H5P_DEFAULT);
417 
418  hid_t dataset_id7 = H5Dopen(group_id7, "Matrix elements", H5P_DEFAULT);
419  H5Dread(dataset_id7, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, theElements);
420  H5Dclose(dataset_id7);
421 
422  H5Gclose(group_id7);
423 
424  H5Fclose(file_id);
425 
426  std::cout << "FourIndex::read : everything loaded!" << std::endl;
427 
428 }
429 
430 #ifdef CHEMPS2_MPI_COMPILATION
431 void CheMPS2::FourIndex::broadcast( const int ROOT ){
432 
433  assert( arrayLength <= INT_MAX ); // To be able to broadcast
434  const int size = ( int )( arrayLength );
435  if ( size > 0 ){
436  MPIchemps2::broadcast_array_double( theElements, size, ROOT );
437  }
438 
439 }
440 #endif
441 
void save(const std::string name) const
Save the FourIndex object.
Definition: FourIndex.cpp:314
void Clear()
Set all two-body matrix elements to zero.
Definition: FourIndex.cpp:157
STL namespace.
double get(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l) const
Get an element.
Definition: FourIndex.cpp:175
virtual ~FourIndex()
Destructor.
Definition: FourIndex.cpp:149
void set(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l, const double val)
Set an element.
Definition: FourIndex.cpp:163
int get_irrep_size(const int irrep) const
Get a given irrep size.
Definition: FourIndex.cpp:294
FourIndex(const int nGroup, const int *IrrepSizes)
Constructor.
Definition: FourIndex.cpp:34
void add(const int irrep_i, const int irrep_j, const int irrep_k, const int irrep_l, const int i, const int j, const int k, const int l, const double val)
Add a double to an element.
Definition: FourIndex.cpp:169
void read(const std::string name)
Load the FourIndex object.
Definition: FourIndex.cpp:372