27 #include "FourIndex.h" 30 #include "MPIchemps2.h" 36 SymmInfo.setGroup(nGroup);
38 Isizes =
new int[SymmInfo.getNumberOfIrreps()];
39 for (
int Icenter=0; Icenter<SymmInfo.getNumberOfIrreps(); Icenter++){
40 Isizes[Icenter] = IrrepSizes[Icenter];
43 arrayLength = calcNumberOfUniqueElements(
true);
44 theElements =
new double[arrayLength];
50 long long CheMPS2::FourIndex::calcNumberOfUniqueElements(
const bool allocate){
53 long long theTotalSize = 0;
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)){
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++){
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);
78 delete [] storage[Icenter][I_i][I_k][i + k*(k+1)/2];
82 if (!allocate){
delete [] storage[Icenter][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++){
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 );
94 delete [] storage[Icenter][I_i][I_k][i + k*Isizes[I_i]];
98 if (!allocate){
delete [] storage[Icenter][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++){
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;
112 delete [] storage[Icenter][I_i][I_k][i + k*(k+1)/2];
116 if (!allocate){
delete [] storage[Icenter][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++){
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];
128 delete [] storage[Icenter][I_i][I_k][i + k*Isizes[I_i]];
132 if (!allocate){
delete [] storage[Icenter][I_i][I_k]; }
138 if (!allocate){
delete [] storage[Icenter][I_i]; }
141 if (!allocate){
delete [] storage[Icenter]; }
143 if (!allocate){
delete [] storage; }
151 arrayLength = calcNumberOfUniqueElements(
false);
152 delete [] theElements;
159 for (
long long count = 0; count < arrayLength; count++){ theElements[ count ] = 0.0; }
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){
165 theElements[getPointer(irrep_i, irrep_j, irrep_k, irrep_l, i, j, k, l)] = val;
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){
171 theElements[getPointer(irrep_i, irrep_j, irrep_k, irrep_l, i, j, k, l)] += val;
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{
177 return theElements[getPointer(irrep_i, irrep_j, irrep_k, irrep_l, i, j, k, l)];
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 {
183 assert( Irreps::directProd( irrep_i, irrep_j ) == Irreps::directProd( irrep_k, irrep_l ) );
185 if ((irrep_i <= irrep_j) && (irrep_i <= irrep_k) && (irrep_j <= irrep_l)){
186 return getPtrIrrepOrderOK(irrep_i,irrep_j,irrep_k,irrep_l,i,j,k,l);
189 if ((irrep_j <= irrep_i) && (irrep_j <= irrep_l) && (irrep_i <= irrep_k)){
190 return getPtrIrrepOrderOK(irrep_j,irrep_i,irrep_l,irrep_k,j,i,l,k);
193 if ((irrep_k <= irrep_j) && (irrep_k <= irrep_i) && (irrep_j <= irrep_l)){
194 return getPtrIrrepOrderOK(irrep_k,irrep_j,irrep_i,irrep_l,k,j,i,l);
197 if ((irrep_j <= irrep_k) && (irrep_j <= irrep_l) && (irrep_k <= irrep_i)){
198 return getPtrIrrepOrderOK(irrep_j,irrep_k,irrep_l,irrep_i,j,k,l,i);
201 if ((irrep_i <= irrep_l) && (irrep_i <= irrep_k) && (irrep_l <= irrep_j)){
202 return getPtrIrrepOrderOK(irrep_i,irrep_l,irrep_k,irrep_j,i,l,k,j);
205 if ((irrep_l <= irrep_i) && (irrep_l <= irrep_j) && (irrep_i <= irrep_k)){
206 return getPtrIrrepOrderOK(irrep_l,irrep_i,irrep_j,irrep_k,l,i,j,k);
209 if ((irrep_k <= irrep_l) && (irrep_k <= irrep_i) && (irrep_l <= irrep_j)){
210 return getPtrIrrepOrderOK(irrep_k,irrep_l,irrep_i,irrep_j,k,l,i,j);
213 if ((irrep_l <= irrep_k) && (irrep_l <= irrep_j) && (irrep_k <= irrep_i)){
214 return getPtrIrrepOrderOK(irrep_l,irrep_k,irrep_j,irrep_i,l,k,j,i);
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 {
224 int Icenter = Irreps::directProd(irrep_i,irrep_j);
228 if (irrep_i == irrep_k){
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);
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);
238 }
else {
return storage[Icenter][irrep_i][irrep_k][i + Isizes[irrep_i]*k][j] + l; }
242 if (irrep_i == irrep_k){
245 if ((i < j) && (i <= k) && (j <= l))
return getPtrAllOK2(Icenter, irrep_i, irrep_k, i, j, k, l);
246 if ((i == j) && (i <= k) && (j <= l)){
247 if (l>=k)
return getPtrAllOK1(Icenter, irrep_i, irrep_k, i, j, k, l);
248 else return getPtrAllOK1(Icenter, irrep_j, irrep_l, j, i, l, k);
250 if ((j < i) && (j <= l) && (i <= k))
return getPtrAllOK2(Icenter, irrep_j, irrep_l, j, i, l, k);
253 if ((k < l) && (k <= i) && (l <= j))
return getPtrAllOK2(Icenter, irrep_k, irrep_i, k, l, i, j);
254 if ((k == l) && (k <= i) && (l <= j)){
255 if (j>=i)
return getPtrAllOK1(Icenter, irrep_k, irrep_i, k, l, i, j);
256 else return getPtrAllOK1(Icenter, irrep_l, irrep_j, l, k, j, i);
258 if ((l < k) && (l <= j) && (k <= i))
return getPtrAllOK2(Icenter, irrep_l, irrep_j, l, k, j, i);
261 if ((k < j) && (k <= i) && (j <= l))
return getPtrAllOK2(Icenter, irrep_k, irrep_i, k, j, i, l);
262 if ((k == j) && (k <= i) && (j <= l)){
263 if (l>=i)
return getPtrAllOK1(Icenter, irrep_k, irrep_i, k, j, i, l);
264 else return getPtrAllOK1(Icenter, irrep_j, irrep_l, j, k, l, i);
266 if ((j < k) && (j <= l) && (k <= i))
return getPtrAllOK2(Icenter, irrep_j, irrep_l, j, k, l, i);
269 if ((i < l) && (i <= k) && (l <= j))
return getPtrAllOK2(Icenter, irrep_i, irrep_k, i, l, k, j);
270 if ((i == l) && (i <= k) && (l <= j)){
271 if (j>=k)
return getPtrAllOK1(Icenter, irrep_i, irrep_k, i, l, k, j);
272 else return getPtrAllOK1(Icenter, irrep_l, irrep_j, l, i, j, k);
274 if ((l < i) && (l <= j) && (i <= k))
return getPtrAllOK2(Icenter, irrep_l, irrep_j, l, i, j, k);
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; }
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; }
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{
298 return storage[Icent][irrep_i][irrep_k][i + k*(k+1)/2][j-i] + l-k;
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{
304 return storage[Icent][irrep_i][irrep_k][i + k*(k+1)/2][j-i] + l-j;
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{
310 return storage[Icent][irrep_i][irrep_k][i + k*(k+1)/2][j] + l-j;
317 hid_t file_id = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
320 hid_t group_id = H5Gcreate(file_id,
"/MetaData", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
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);
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);
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);
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);
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);
350 H5Dclose(dataset_id);
351 H5Sclose(dataspace_id);
356 hid_t group_id7 = H5Gcreate(file_id,
"/FourIndexObject", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
358 hsize_t dimarray7 = arrayLength;
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);
363 H5Dclose(dataset_id7);
364 H5Sclose(dataspace_id7);
375 hid_t file_id = H5Fopen(name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
378 hid_t group_id = H5Gopen(file_id,
"/MetaData",H5P_DEFAULT);
381 hid_t dataset_id = H5Dopen(group_id,
"IrrepSizes", H5P_DEFAULT);
384 hid_t attribute_id1 = H5Aopen_by_name(group_id,
"IrrepSizes",
"nGroup", H5P_DEFAULT, H5P_DEFAULT);
386 H5Aread(attribute_id1, H5T_NATIVE_INT, &nGroup);
387 assert( nGroup==SymmInfo.getGroupNumber() );
389 hid_t attribute_id2 = H5Aopen_by_name(group_id,
"IrrepSizes",
"nIrreps", H5P_DEFAULT, H5P_DEFAULT);
391 H5Aread(attribute_id2, H5T_NATIVE_INT, &nIrreps);
392 assert( nIrreps==SymmInfo.getNumberOfIrreps() );
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 );
399 H5Aclose(attribute_id1);
400 H5Aclose(attribute_id2);
401 H5Aclose(attribute_id3);
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] );
408 delete [] IsizesAgain;
409 H5Dclose(dataset_id);
413 std::cout <<
"FourIndex::read : loading " << arrayLength <<
" doubles." << std::endl;
416 hid_t group_id7 = H5Gopen(file_id,
"/FourIndexObject", H5P_DEFAULT);
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);
426 std::cout <<
"FourIndex::read : everything loaded!" << std::endl;
430 #ifdef CHEMPS2_MPI_COMPILATION 431 void CheMPS2::FourIndex::broadcast(
const int ROOT ){
433 assert( arrayLength <= INT_MAX );
434 const int size = ( int )( arrayLength );
436 MPIchemps2::broadcast_array_double( theElements, size, ROOT );
void save(const std::string name) const
Save the FourIndex object.
void Clear()
Set all two-body matrix elements to zero.
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.
virtual ~FourIndex()
Destructor.
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.
int get_irrep_size(const int irrep) const
Get a given irrep size.
FourIndex(const int nGroup, const int *IrrepSizes)
Constructor.
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.
void read(const std::string name)
Load the FourIndex object.