v2DM-DOCI  1.0
TwoIndex.cpp
Go to the documentation of this file.
1 /*
2  CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3  Copyright (C) 2013, 2014 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 <cstring>
26 
27 #include "TwoIndex.h"
28 #include "MyHDF5.h"
29 
30 using namespace std;
31 
32 CheMPS2::TwoIndex::TwoIndex(const int nGroup, const int * IrrepSizes){
33 
34  SymmInfo.setGroup(nGroup);
35 
36  Isizes = new int[SymmInfo.getNumberOfIrreps()];
37  storage = new double*[SymmInfo.getNumberOfIrreps()];
38 
39  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++){
40  Isizes[cnt] = IrrepSizes[cnt];
41  if (Isizes[cnt]>0) storage[cnt] = new double[Isizes[cnt]*(Isizes[cnt]+1)/2];
42  }
43 
44 }
45 
47 {
48  SymmInfo = orig.SymmInfo;
49 
50  Isizes = new int[SymmInfo.getNumberOfIrreps()];
51  storage = new double*[SymmInfo.getNumberOfIrreps()];
52 
53  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++)
54  {
55  Isizes[cnt] = orig.Isizes[cnt];
56  if (Isizes[cnt]>0)
57  {
58  storage[cnt] = new double[Isizes[cnt]*(Isizes[cnt]+1)/2];
59  memcpy(storage[cnt], orig.storage[cnt], sizeof(double)*Isizes[cnt]*(Isizes[cnt]+1)/2);
60  }
61  }
62 }
63 
65 
66  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++) if (Isizes[cnt]>0) delete [] storage[cnt];
67  delete [] storage;
68  delete [] Isizes;
69 
70 }
71 
72 void CheMPS2::TwoIndex::set(const int irrep, const int i, const int j, const double val){
73 
74  if (i>j) storage[irrep][j + i*(i+1)/2] = val;
75  else storage[irrep][i + j*(j+1)/2] = val;
76 
77 }
78 
79 double CheMPS2::TwoIndex::get(const int irrep, const int i, const int j) const{
80 
81  if (i>j) return storage[irrep][j + i*(i+1)/2];
82  return storage[irrep][i + j*(j+1)/2];
83 
84 }
85 
86 void CheMPS2::TwoIndex::save(const std::string name) const{
87 
88  //The hdf5 file
89  hid_t file_id = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
90 
91  //The metadata
92  hid_t group_id = H5Gcreate(file_id, "/MetaData", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
93 
94  //The IrrepSizes
95  hsize_t dimarray = SymmInfo.getNumberOfIrreps();
96  hid_t dataspace_id = H5Screate_simple(1, &dimarray, NULL);
97  hid_t dataset_id = H5Dcreate(group_id, "IrrepSizes", H5T_STD_I32LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
98  H5Dwrite(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, Isizes);
99 
100  //Attributes
101  hid_t attribute_space_id1 = H5Screate(H5S_SCALAR);
102  hid_t attribute_id1 = H5Acreate(dataset_id, "nGroup", H5T_STD_I32LE, attribute_space_id1, H5P_DEFAULT, H5P_DEFAULT);
103  int nGroup = SymmInfo.getGroupNumber();
104  H5Awrite(attribute_id1, H5T_NATIVE_INT, &nGroup);
105 
106  hid_t attribute_space_id2 = H5Screate(H5S_SCALAR);
107  hid_t attribute_id2 = H5Acreate(dataset_id, "nIrreps", H5T_STD_I32LE, attribute_space_id2, H5P_DEFAULT, H5P_DEFAULT);
108  int nIrreps = SymmInfo.getNumberOfIrreps();
109  H5Awrite(attribute_id2, H5T_NATIVE_INT, &nIrreps);
110 
111  H5Aclose(attribute_id1);
112  H5Aclose(attribute_id2);
113  H5Sclose(attribute_space_id1);
114  H5Sclose(attribute_space_id2);
115 
116  H5Dclose(dataset_id);
117  H5Sclose(dataspace_id);
118 
119  H5Gclose(group_id);
120 
121  //The object itself.
122  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++){
123 
124  if (Isizes[cnt]>0) {
125 
126  std::stringstream sstream;
127  sstream << "/TwoIndex" << cnt ;
128  hid_t group_id3 = H5Gcreate(file_id, sstream.str().c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
129 
130  hsize_t dimarray3 = Isizes[cnt]*(Isizes[cnt]+1)/2;
131  hid_t dataspace_id3 = H5Screate_simple(1, &dimarray3, NULL);
132  hid_t dataset_id3 = H5Dcreate(group_id3, "Matrix elements", H5T_IEEE_F64LE, dataspace_id3, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
133  H5Dwrite(dataset_id3, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, storage[cnt]);
134 
135  H5Dclose(dataset_id3);
136  H5Sclose(dataspace_id3);
137 
138  H5Gclose(group_id3);
139 
140  }
141 
142  }
143 
144  H5Fclose(file_id);
145 
146 }
147 
148 void CheMPS2::TwoIndex::read(const std::string name){
149 
150  //The hdf5 file
151  hid_t file_id = H5Fopen(name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
152 
153  //The metadata
154  hid_t group_id = H5Gopen(file_id, "/MetaData",H5P_DEFAULT);
155 
156  //The IrrepSizes
157  hid_t dataset_id = H5Dopen(group_id, "IrrepSizes", H5P_DEFAULT);
158 
159  //Attributes
160  hid_t attribute_id1 = H5Aopen_by_name(group_id,"IrrepSizes", "nGroup", H5P_DEFAULT, H5P_DEFAULT);
161  int nGroup;
162  H5Aread(attribute_id1, H5T_NATIVE_INT, &nGroup);
163  assert( nGroup==SymmInfo.getGroupNumber() );
164 
165  hid_t attribute_id2 = H5Aopen_by_name(group_id,"IrrepSizes", "nIrreps", H5P_DEFAULT, H5P_DEFAULT);
166  int nIrreps;
167  H5Aread(attribute_id2, H5T_NATIVE_INT, &nIrreps);
168  assert( nIrreps==SymmInfo.getNumberOfIrreps() );
169 
170  H5Aclose(attribute_id1);
171  H5Aclose(attribute_id2);
172 
173  int * IsizesAgain = new int[SymmInfo.getNumberOfIrreps()];
174  H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, IsizesAgain);
175  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++){
176  assert( IsizesAgain[cnt]==Isizes[cnt] );
177  }
178  delete [] IsizesAgain;
179  H5Dclose(dataset_id);
180 
181  H5Gclose(group_id);
182 
183  //The object itself.
184  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++){
185 
186  if (Isizes[cnt]>0) {
187 
188  std::stringstream sstream;
189  sstream << "/TwoIndex" << cnt ;
190  hid_t group_id3 = H5Gopen(file_id, sstream.str().c_str(), H5P_DEFAULT);
191 
192  hid_t dataset_id3 = H5Dopen(group_id3, "Matrix elements", H5P_DEFAULT);
193  H5Dread(dataset_id3, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, storage[cnt]);
194  H5Dclose(dataset_id3);
195 
196  H5Gclose(group_id3);
197 
198  }
199 
200  }
201 
202  H5Fclose(file_id);
203 
204 }
205 
206 
207 
208 void CheMPS2::TwoIndex::save2(const std::string name) const
209 {
210  //The hdf5 file
211  hid_t file_id = H5Fopen(name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
212 
213  hid_t group_id2 = H5Gcreate(file_id, "TwoIndex", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
214 
215  //The metadata
216  hid_t group_id = H5Gcreate(group_id2, "MetaData", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
217 
218  //The IrrepSizes
219  hsize_t dimarray = SymmInfo.getNumberOfIrreps();
220  hid_t dataspace_id = H5Screate_simple(1, &dimarray, NULL);
221  hid_t dataset_id = H5Dcreate(group_id, "IrrepSizes", H5T_STD_I32LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
222  H5Dwrite(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, Isizes);
223 
224  //Attributes
225  hid_t attribute_space_id1 = H5Screate(H5S_SCALAR);
226  hid_t attribute_id1 = H5Acreate(dataset_id, "nGroup", H5T_STD_I32LE, attribute_space_id1, H5P_DEFAULT, H5P_DEFAULT);
227  int nGroup = SymmInfo.getGroupNumber();
228  H5Awrite(attribute_id1, H5T_NATIVE_INT, &nGroup);
229 
230  hid_t attribute_space_id2 = H5Screate(H5S_SCALAR);
231  hid_t attribute_id2 = H5Acreate(dataset_id, "nIrreps", H5T_STD_I32LE, attribute_space_id2, H5P_DEFAULT, H5P_DEFAULT);
232  int nIrreps = SymmInfo.getNumberOfIrreps();
233  H5Awrite(attribute_id2, H5T_NATIVE_INT, &nIrreps);
234 
235  H5Aclose(attribute_id1);
236  H5Aclose(attribute_id2);
237  H5Sclose(attribute_space_id1);
238  H5Sclose(attribute_space_id2);
239 
240  H5Dclose(dataset_id);
241  H5Sclose(dataspace_id);
242 
243  H5Gclose(group_id);
244 
245  //The object itself.
246  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++){
247 
248  if (Isizes[cnt]>0) {
249 
250  std::stringstream sstream;
251  sstream << "TwoIndex" << cnt ;
252  hid_t group_id3 = H5Gcreate(group_id2, sstream.str().c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
253 
254  hsize_t dimarray3 = Isizes[cnt]*(Isizes[cnt]+1)/2;
255  hid_t dataspace_id3 = H5Screate_simple(1, &dimarray3, NULL);
256  hid_t dataset_id3 = H5Dcreate(group_id3, "Matrix elements", H5T_IEEE_F64LE, dataspace_id3, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
257  H5Dwrite(dataset_id3, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, storage[cnt]);
258 
259  H5Dclose(dataset_id3);
260  H5Sclose(dataspace_id3);
261 
262  H5Gclose(group_id3);
263 
264  }
265 
266  }
267 
268  H5Gclose(group_id2);
269 
270  H5Fclose(file_id);
271 }
272 
273 void CheMPS2::TwoIndex::read2(const std::string name)
274 {
275 
276  //The hdf5 file
277  hid_t file_id = H5Fopen(name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
278 
279  //The metadata
280  hid_t group_id = H5Gopen(file_id, "/TwoIndex/MetaData",H5P_DEFAULT);
281 
282  //The IrrepSizes
283  hid_t dataset_id = H5Dopen(group_id, "IrrepSizes", H5P_DEFAULT);
284 
285  //Attributes
286  hid_t attribute_id1 = H5Aopen_by_name(group_id,"IrrepSizes", "nGroup", H5P_DEFAULT, H5P_DEFAULT);
287  int nGroup;
288  H5Aread(attribute_id1, H5T_NATIVE_INT, &nGroup);
289  assert( nGroup==SymmInfo.getGroupNumber() );
290 
291  hid_t attribute_id2 = H5Aopen_by_name(group_id,"IrrepSizes", "nIrreps", H5P_DEFAULT, H5P_DEFAULT);
292  int nIrreps;
293  H5Aread(attribute_id2, H5T_NATIVE_INT, &nIrreps);
294  assert( nIrreps==SymmInfo.getNumberOfIrreps() );
295 
296  H5Aclose(attribute_id1);
297  H5Aclose(attribute_id2);
298 
299  int * IsizesAgain = new int[SymmInfo.getNumberOfIrreps()];
300  H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, IsizesAgain);
301  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++){
302  assert( IsizesAgain[cnt]==Isizes[cnt] );
303  }
304  delete [] IsizesAgain;
305  H5Dclose(dataset_id);
306 
307  H5Gclose(group_id);
308 
309  //The object itself.
310  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++){
311 
312  if (Isizes[cnt]>0) {
313 
314  std::stringstream sstream;
315  sstream << "/TwoIndex/TwoIndex" << cnt ;
316  hid_t group_id3 = H5Gopen(file_id, sstream.str().c_str(), H5P_DEFAULT);
317 
318  hid_t dataset_id3 = H5Dopen(group_id3, "Matrix elements", H5P_DEFAULT);
319  H5Dread(dataset_id3, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, storage[cnt]);
320  H5Dclose(dataset_id3);
321 
322  H5Gclose(group_id3);
323 
324  }
325 
326  }
327 
328  H5Fclose(file_id);
329 }
330 
332 {
333  for (int cnt=0; cnt<SymmInfo.getNumberOfIrreps(); cnt++)
334  if (Isizes[cnt]>0)
335  memset(storage[cnt], 0, sizeof(double)*(Isizes[cnt]*(Isizes[cnt]+1)/2));
336 }
Irreps SymmInfo
Definition: TwoIndex.h:78
void reset()
set everything to zero
Definition: TwoIndex.cpp:331
virtual ~TwoIndex()
Destructor.
Definition: TwoIndex.cpp:64
STL namespace.
void read(const std::string name)
Load the TwoIndex object.
Definition: TwoIndex.cpp:148
TwoIndex(const int nGroup, const int *IrrepSizes)
Constructor.
Definition: TwoIndex.cpp:32
void read2(const std::string name)
Definition: TwoIndex.cpp:273
void set(const int irrep, const int i, const int j, const double val)
Set an element.
Definition: TwoIndex.cpp:72
double ** storage
Definition: TwoIndex.h:84
double get(const int irrep, const int i, const int j) const
Get an element.
Definition: TwoIndex.cpp:79
void save2(const std::string name) const
Definition: TwoIndex.cpp:208
void save(const std::string name) const
Save the TwoIndex object.
Definition: TwoIndex.cpp:86