v2DM-DOCI  1.0
Hamiltonian.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 <string>
24 #include <cstring>
25 #include <fstream>
26 
27 #include "Irreps.h"
28 #include "TwoIndex.h"
29 #include "FourIndex.h"
30 #include "Hamiltonian.h"
31 #include "MyHDF5.h"
32 
33 using std::cout;
34 using std::endl;
35 using std::string;
36 using std::ifstream;
37 
38 CheMPS2::Hamiltonian::Hamiltonian(const int Norbitals, const int nGroup, const int * OrbIrreps){
39 
40  L = Norbitals;
41  assert( nGroup>=0 );
42  assert( nGroup<=7 );
43  SymmInfo.setGroup(nGroup);
44 
45  orb2irrep.reset(new int[L]);
46  orb2indexSy.reset(new int[L]);
47  int nIrreps = SymmInfo.getNumberOfIrreps();
48  irrep2num_orb.reset(new int[nIrreps]);
49  for (int cnt=0; cnt<nIrreps; cnt++) irrep2num_orb[cnt] = 0;
50  for (int cnt=0; cnt<L; cnt++){
51  assert( OrbIrreps[cnt]>=0 );
52  assert( OrbIrreps[cnt]<nIrreps );
53  orb2irrep[cnt] = OrbIrreps[cnt];
54  orb2indexSy[cnt] = irrep2num_orb[orb2irrep[cnt]];
55  irrep2num_orb[orb2irrep[cnt]]++;
56  }
57 
58  Tmat.reset(new TwoIndex(SymmInfo.getGroupNumber(),irrep2num_orb.get()));
59  Vmat.reset(new FourIndex(SymmInfo.getGroupNumber(),irrep2num_orb.get()));
60 
61  Ne = 0;
62 }
63 
64 CheMPS2::Hamiltonian::Hamiltonian(const Hamiltonian &orig)
65 {
66  L = orig.L;
67  Ne = orig.Ne;
68  SymmInfo = orig.SymmInfo;
69  Econst = orig.Econst;
70 
71  orb2irrep.reset(new int[L]);
72  orb2indexSy.reset(new int[L]);
73  int nIrreps = SymmInfo.getNumberOfIrreps();
74  irrep2num_orb.reset(new int[nIrreps]);
75 
76  memcpy(orb2irrep.get(), orig.orb2irrep.get(), sizeof(int)*L);
77  memcpy(orb2indexSy.get(), orig.orb2indexSy.get(), sizeof(int)*L);
78  memcpy(irrep2num_orb.get(), orig.irrep2num_orb.get(), sizeof(int)*nIrreps);
79 
80  Tmat.reset(new TwoIndex(*orig.Tmat));
81  Vmat.reset(new FourIndex(*orig.Vmat));
82 }
83 
84 CheMPS2::Hamiltonian& CheMPS2::Hamiltonian::operator=(const CheMPS2::Hamiltonian &orig)
85 {
86  L = orig.L;
87  Ne = orig.Ne;
88  SymmInfo = orig.SymmInfo;
89  Econst = orig.Econst;
90 
91  orb2irrep.reset(new int[L]);
92  orb2indexSy.reset(new int[L]);
93  int nIrreps = SymmInfo.getNumberOfIrreps();
94  irrep2num_orb.reset(new int[nIrreps]);
95 
96  memcpy(orb2irrep.get(), orig.orb2irrep.get(), sizeof(int)*L);
97  memcpy(orb2indexSy.get(), orig.orb2indexSy.get(), sizeof(int)*L);
98  memcpy(irrep2num_orb.get(), orig.irrep2num_orb.get(), sizeof(int)*nIrreps);
99 
100  Tmat.reset(new TwoIndex(*orig.Tmat));
101  Vmat.reset(new FourIndex(*orig.Vmat));
102 
103  return *this;
104 }
105 
106 CheMPS2::Hamiltonian::Hamiltonian(const string file_psi4text){
107 
108  CreateAndFillFromPsi4dump( file_psi4text );
109 
110 }
111 
112 CheMPS2::Hamiltonian::Hamiltonian(const bool fileh5, const string main_file, const string file_tmat, const string file_vmat){
113 
114  Ne = 0;
115 
116  if (fileh5){
117  CreateAndFillFromH5( main_file, file_tmat, file_vmat );
118  } else {
119  CreateAndFillFromPsi4dump( main_file );
120  }
121 
122 }
123 
124 /*
125 CheMPS2::Hamiltonian::~Hamiltonian(){
126 
127  delete [] orb2irrep;
128  delete [] orb2indexSy;
129  delete [] irrep2num_orb;
130  delete Tmat;
131  delete Vmat;
132 
133 }
134 */
135 
136 int CheMPS2::Hamiltonian::getL() const{ return L; }
137 
138 int CheMPS2::Hamiltonian::getNGroup() const{ return SymmInfo.getGroupNumber(); }
139 
140 int CheMPS2::Hamiltonian::getOrbitalIrrep(const int nOrb) const{ return orb2irrep[nOrb]; }
141 
142 void CheMPS2::Hamiltonian::setEconst(const double val){ Econst = val; }
143 
144 double CheMPS2::Hamiltonian::getEconst() const{ return Econst; }
145 
146 void CheMPS2::Hamiltonian::setTmat(const int index1, const int index2, const double val){
147 
148  assert( orb2irrep[index1]==orb2irrep[index2] );
149  Tmat->set(orb2irrep[index1], orb2indexSy[index1], orb2indexSy[index2], val);
150 
151 }
152 
153 double CheMPS2::Hamiltonian::getTmat(const int index1, const int index2) const{
154 
155  if (orb2irrep[index1]==orb2irrep[index2]){
156  return Tmat->get(orb2irrep[index1], orb2indexSy[index1], orb2indexSy[index2]);
157  }
158 
159  return 0.0;
160 
161 }
162 
163 void CheMPS2::Hamiltonian::setVmat(const int index1, const int index2, const int index3, const int index4, const double val){
164 
165  assert( Irreps::directProd(orb2irrep[index1],orb2irrep[index2]) == Irreps::directProd(orb2irrep[index3],orb2irrep[index4]) );
166  Vmat->set(orb2irrep[index1], orb2irrep[index2], orb2irrep[index3], orb2irrep[index4], orb2indexSy[index1], orb2indexSy[index2], orb2indexSy[index3], orb2indexSy[index4], val);
167 
168 }
169 
170 void CheMPS2::Hamiltonian::addToVmat(const int index1, const int index2, const int index3, const int index4, const double val){
171 
172  assert( Irreps::directProd(orb2irrep[index1],orb2irrep[index2]) == Irreps::directProd(orb2irrep[index3],orb2irrep[index4]) );
173  Vmat->add(orb2irrep[index1], orb2irrep[index2], orb2irrep[index3], orb2irrep[index4], orb2indexSy[index1], orb2indexSy[index2], orb2indexSy[index3], orb2indexSy[index4], val);
174 
175 }
176 
177 double CheMPS2::Hamiltonian::getVmat(const int index1, const int index2, const int index3, const int index4) const{
178 
179  if ( Irreps::directProd(orb2irrep[index1],orb2irrep[index2]) == Irreps::directProd(orb2irrep[index3],orb2irrep[index4]) ){
180  return Vmat->get(orb2irrep[index1], orb2irrep[index2], orb2irrep[index3], orb2irrep[index4], orb2indexSy[index1], orb2indexSy[index2], orb2indexSy[index3], orb2indexSy[index4]);
181  }
182 
183  return 0.0;
184 
185 }
186 
187 void CheMPS2::Hamiltonian::save(const string file_parent, const string file_tmat, const string file_vmat) const{
188 
189  Tmat->save(file_tmat);
190  Vmat->save(file_vmat);
191 
192  //The hdf5 file
193  hid_t file_id = H5Fcreate(file_parent.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
194 
195  //The data
196  hid_t group_id = H5Gcreate(file_id, "/Data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
197 
198  //The chain length
199  hsize_t dimarray = 1;
200  hid_t dataspace_id = H5Screate_simple(1, &dimarray, NULL);
201  hid_t dataset_id = H5Dcreate(group_id, "L", H5T_STD_I32LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
202  H5Dwrite(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &L);
203 
204  //The group number
205  hsize_t dimarray2 = 1;
206  hid_t dataspace_id2 = H5Screate_simple(1, &dimarray2, NULL);
207  hid_t dataset_id2 = H5Dcreate(group_id, "nGroup", H5T_STD_I32LE, dataspace_id2, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
208  int nGroup = SymmInfo.getGroupNumber();
209  H5Dwrite(dataset_id2, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &nGroup);
210 
211  //orb2irrep
212  hsize_t dimarray3 = L;
213  hid_t dataspace_id3 = H5Screate_simple(1, &dimarray3, NULL);
214  hid_t dataset_id3 = H5Dcreate(group_id, "orb2irrep", H5T_STD_I32LE, dataspace_id3, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
215  H5Dwrite(dataset_id3, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, orb2irrep.get());
216 
217  //Econst
218  hsize_t dimarray4 = 1;
219  hid_t dataspace_id4 = H5Screate_simple(1, &dimarray4, NULL);
220  hid_t dataset_id4 = H5Dcreate(group_id, "Econst", H5T_IEEE_F64LE, dataspace_id4, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
221  H5Dwrite(dataset_id4, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &Econst);
222 
223  H5Dclose(dataset_id);
224  H5Sclose(dataspace_id);
225  H5Dclose(dataset_id2);
226  H5Sclose(dataspace_id2);
227  H5Dclose(dataset_id3);
228  H5Sclose(dataspace_id3);
229  H5Dclose(dataset_id4);
230  H5Sclose(dataspace_id4);
231 
232  H5Gclose(group_id);
233 
234  H5Fclose(file_id);
235 
236 }
237 
238 void CheMPS2::Hamiltonian::read(const string file_parent, const string file_tmat, const string file_vmat){
239 
240  Tmat->read(file_tmat);
241  Vmat->read(file_vmat);
242 
243  //The hdf5 file
244  hid_t file_id = H5Fopen(file_parent.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
245 
246  //The data
247  hid_t group_id = H5Gopen(file_id, "/Data",H5P_DEFAULT);
248 
249  //The chain length
250  hid_t dataset_id1 = H5Dopen(group_id, "L", H5P_DEFAULT);
251  int Lagain;
252  H5Dread(dataset_id1, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &Lagain);
253  assert( L==Lagain );
254 
255  //The group number
256  hid_t dataset_id2 = H5Dopen(group_id, "nGroup", H5P_DEFAULT);
257  int nGroup;
258  H5Dread(dataset_id2, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &nGroup);
259  assert( nGroup==SymmInfo.getGroupNumber() );
260 
261  //orb2irrep
262  hid_t dataset_id3 = H5Dopen(group_id, "orb2irrep", H5P_DEFAULT);
263  int * orb2irrepAgain = new int[L];
264  H5Dread(dataset_id3, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, orb2irrepAgain);
265  for (int cnt=0; cnt<L; cnt++){
266  assert( orb2irrep[cnt]==orb2irrepAgain[cnt] );
267  }
268  delete [] orb2irrepAgain;
269 
270  //Econst
271  hid_t dataset_id4 = H5Dopen(group_id, "Econst", H5P_DEFAULT);
272  H5Dread(dataset_id4, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &Econst);
273 
274  H5Dclose(dataset_id1);
275  H5Dclose(dataset_id2);
276  H5Dclose(dataset_id3);
277  H5Dclose(dataset_id4);
278 
279  H5Gclose(group_id);
280 
281  H5Fclose(file_id);
282 
283  if (CheMPS2::HAMILTONIAN_debugPrint) debugcheck();
284 
285 }
286 
287 void CheMPS2::Hamiltonian::CreateAndFillFromH5(const string file_parent, const string file_tmat, const string file_vmat){
288 
289  //The hdf5 file
290  hid_t file_id = H5Fopen(file_parent.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
291 
292  //The data
293  hid_t group_id = H5Gopen(file_id, "/Data",H5P_DEFAULT);
294 
295  //The number of orbitals: Set L directly!
296  hid_t dataset_id1 = H5Dopen(group_id, "L", H5P_DEFAULT);
297  H5Dread(dataset_id1, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &L);
298 
299  //The group number: SymmInfo.setGroup()
300  hid_t dataset_id2 = H5Dopen(group_id, "nGroup", H5P_DEFAULT);
301  int nGroup_LOADH5;
302  H5Dread(dataset_id2, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &nGroup_LOADH5);
303  SymmInfo.setGroup( nGroup_LOADH5 );
304 
305 
306  hid_t dataset_id5 = H5Dopen(group_id, "nelectrons", H5P_DEFAULT);
307  H5Dread(dataset_id5, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &Ne);
308 
309  //The irrep of each orbital: Create and fill orb2irrep directly!
310  hid_t dataset_id3 = H5Dopen(group_id, "orb2irrep", H5P_DEFAULT);
311  orb2irrep.reset(new int[L]);
312  H5Dread(dataset_id3, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, orb2irrep.get());
313 
314  H5Dclose(dataset_id1);
315  H5Dclose(dataset_id2);
316  H5Dclose(dataset_id3);
317  H5Dclose(dataset_id5);
318 
319  H5Gclose(group_id);
320 
321  H5Fclose(file_id);
322 
323  orb2indexSy.reset(new int[L]);
324  int nIrreps = SymmInfo.getNumberOfIrreps();
325  irrep2num_orb.reset(new int[nIrreps]);
326 
327  for (int irrep = 0; irrep < nIrreps; irrep++){ irrep2num_orb[ irrep ] = 0; }
328  for (int orb = 0; orb < L; orb++){
329  orb2indexSy[ orb ] = irrep2num_orb[ orb2irrep[ orb ] ];
330  irrep2num_orb[ orb2irrep[ orb ] ]++;
331  }
332 
333  Tmat.reset(new TwoIndex(SymmInfo.getGroupNumber(),irrep2num_orb.get()));
334  Vmat.reset(new FourIndex(SymmInfo.getGroupNumber(),irrep2num_orb.get()));
335 
336  read(file_parent, file_tmat, file_vmat);
337 
338 }
339 
340 //Works for the file mointegrals/mointegrals.cc_PRINT which can be used as a plugin in psi4 beta5
341 void CheMPS2::Hamiltonian::CreateAndFillFromPsi4dump(const string filename){
342 
343  string line, part;
344  int pos;
345 
346  ifstream inputfile(filename.c_str());
347 
348  //First go to the start of the integral dump.
349  bool stop = false;
350  string start = "**** Molecular Integrals For CheMPS Start Here";
351  do{
352  getline(inputfile,line);
353  pos = line.find(start);
354  if (pos==0) stop = true;
355  } while (!stop);
356 
357  //Get the group name and convert it to the group number
358  getline(inputfile,line);
359  pos = line.find("=");
360  part = line.substr(pos+2,line.size()-pos-3);
361  int nGroup = 0;
362  stop = false;
363  do {
364  if (part.compare(SymmInfo.getGroupName(nGroup))==0) stop = true;
365  else nGroup += 1;
366  } while (!stop);
367  SymmInfo.setGroup(nGroup);
368  //cout << "The group was found to be " << SymmInfo.getGroupName() << " ." << endl;
369 
370  //This line says how many irreps there are: skip.
371  getline(inputfile,line);
372 
373  //This line contains the nuclear energy part.
374  getline(inputfile,line);
375  pos = line.find("=");
376  part = line.substr(pos+2,line.size()-pos-3);
377  Econst = atof(part.c_str());
378 
379  //This line contains the number of MO's.
380  getline(inputfile,line);
381  pos = line.find("=");
382  part = line.substr(pos+2,line.size()-pos-3);
383  L = atoi(part.c_str());
384 
385  //This line contains only text
386  getline(inputfile,line);
387 
388  //This line contains the irrep numbers --> allocate, read in & set
389  getline(inputfile,line);
390 
391  orb2irrep.reset(new int[L]);
392  orb2indexSy.reset(new int[L]);
393  int nIrreps = SymmInfo.getNumberOfIrreps();
394  irrep2num_orb.reset(new int[nIrreps]);
395 
396  pos = 0;
397  do {
398  orb2irrep[pos] = atoi(line.substr(2*pos,1).c_str());
399  pos++;
400  } while (2*pos < (int)line.size()-1);
401 
402  for (int cnt=0; cnt<nIrreps; cnt++) irrep2num_orb[cnt] = 0;
403  for (int cnt=0; cnt<L; cnt++){
404  orb2indexSy[cnt] = irrep2num_orb[orb2irrep[cnt]];
405  irrep2num_orb[orb2irrep[cnt]]++;
406  }
407  Tmat.reset(new TwoIndex(SymmInfo.getGroupNumber(),irrep2num_orb.get()));
408  Vmat.reset(new FourIndex(SymmInfo.getGroupNumber(),irrep2num_orb.get()));
409 
410  //Skip three lines --> number of double occupations, single occupations and test line
411  getline(inputfile,line);
412  getline(inputfile,line);
413  getline(inputfile,line);
414 
415  //Read in one-electron integrals
416  getline(inputfile,line);
417  int pos2, index1, index2;
418  double value;
419  while( (line.substr(0,1)).compare("*")!=0 ){
420 
421  pos = 0;
422  pos2 = line.find(" ",pos);
423  index1 = atoi(line.substr(pos,pos2-pos).c_str());
424 
425  pos = pos2+1;
426  pos2 = line.find(" ",pos);
427  index2 = atoi(line.substr(pos,pos2-pos).c_str());
428 
429  value = atof(line.substr(pos2+1,line.size()-pos2-2).c_str());
430 
431  setTmat(index1,index2,value);
432 
433  getline(inputfile,line);
434 
435  }
436 
437  //Read in two-electron integrals --> in file: chemical notation; in Vmat: physics notation
438  getline(inputfile,line);
439  int index3, index4;
440  while( (line.substr(0,1)).compare("*")!=0 ){
441 
442  pos = 0;
443  pos2 = line.find(" ",pos);
444  index1 = atoi(line.substr(pos,pos2-pos).c_str());
445 
446  pos = pos2+1;
447  pos2 = line.find(" ",pos);
448  index2 = atoi(line.substr(pos,pos2-pos).c_str());
449 
450  pos = pos2+1;
451  pos2 = line.find(" ",pos);
452  index3 = atoi(line.substr(pos,pos2-pos).c_str());
453 
454  pos = pos2+1;
455  pos2 = line.find(" ",pos);
456  index4 = atoi(line.substr(pos,pos2-pos).c_str());
457 
458  value = atof(line.substr(pos2+1,line.size()-pos2-2).c_str());
459 
460  setVmat(index1, index3, index2, index4, value);
461 
462  getline(inputfile,line);
463 
464  }
465 
466  if (CheMPS2::HAMILTONIAN_debugPrint) debugcheck();
467 
468  inputfile.close();
469 
470 }
471 
472 void CheMPS2::Hamiltonian::debugcheck() const{
473 
474  cout << "Econst = " << Econst << endl;
475 
476  double test = 0.0;
477  double test2 = 0.0;
478  double test3 = 0.0;
479  for (int i=0; i<L; i++){
480  test3 += getTmat(i,i);
481  for (int j=0; j<L; j++){
482  test += getTmat(i,j);
483  if (i<=j) test2 += getTmat(i,j);
484  }
485  }
486  cout << "1-electron integrals: Trace : " << test3 << endl;
487  cout << "1-electron integrals: Sum over all elements : " << test << endl;
488  cout << "1-electron integrals: Sum over Tij with i<=j : " << test2 << endl;
489 
490  test = 0.0;
491  test2 = 0.0;
492  test3 = 0.0;
493  for (int i=0; i<L; i++){
494  test3 += getVmat(i,i,i,i);
495  for (int j=0; j<L; j++){
496  for (int k=0; k<L; k++){
497  for (int l=0; l<L; l++){
498  test += getVmat(i,j,k,l);
499  if ((i<=j) && (j<=k) && (k<=l)) test2 += getVmat(i,j,k,l);
500  }
501  }
502  }
503  }
504  cout << "2-electron integrals: Trace : " << test3 << endl;
505  cout << "2-electron integrals: Sum over all elements : " << test << endl;
506  cout << "2-electron integrals: Sum over Vijkl with i<=j<=k<=l : " << test2 << endl;
507 
508 }
509 
510 
511 int CheMPS2::Hamiltonian::getNe() const
512 {
513  return Ne;
514 }
515 
516 void CheMPS2::Hamiltonian::save2(const string filename) const
517 {
518  //The hdf5 file
519  hid_t file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
520 
521  //The data
522  hid_t group_id = H5Gcreate(file_id, "/Data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
523 
524  //The chain length
525  hsize_t dimarray = 1;
526  hid_t dataspace_id = H5Screate_simple(1, &dimarray, NULL);
527  hid_t dataset_id = H5Dcreate(group_id, "L", H5T_STD_I32LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
528  H5Dwrite(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &L);
529 
530  //The group number
531  hsize_t dimarray2 = 1;
532  hid_t dataspace_id2 = H5Screate_simple(1, &dimarray2, NULL);
533  hid_t dataset_id2 = H5Dcreate(group_id, "nGroup", H5T_STD_I32LE, dataspace_id2, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
534  int nGroup = SymmInfo.getGroupNumber();
535  H5Dwrite(dataset_id2, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &nGroup);
536 
537  //orb2irrep
538  hsize_t dimarray3 = L;
539  hid_t dataspace_id3 = H5Screate_simple(1, &dimarray3, NULL);
540  hid_t dataset_id3 = H5Dcreate(group_id, "orb2irrep", H5T_STD_I32LE, dataspace_id3, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
541  H5Dwrite(dataset_id3, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, orb2irrep.get());
542 
543  //Econst
544  hsize_t dimarray4 = 1;
545  hid_t dataspace_id4 = H5Screate_simple(1, &dimarray4, NULL);
546  hid_t dataset_id4 = H5Dcreate(group_id, "Econst", H5T_IEEE_F64LE, dataspace_id4, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
547  H5Dwrite(dataset_id4, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &Econst);
548 
549  hid_t dataspace_id5 = H5Screate(H5S_SCALAR);
550  hid_t dataset_id5 = H5Dcreate(group_id, "nelectrons", H5T_STD_I32LE, dataspace_id5, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
551  H5Dwrite(dataset_id5, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &Ne);
552 
553  H5Dclose(dataset_id);
554  H5Sclose(dataspace_id);
555  H5Dclose(dataset_id2);
556  H5Sclose(dataspace_id2);
557  H5Dclose(dataset_id3);
558  H5Sclose(dataspace_id3);
559  H5Dclose(dataset_id4);
560  H5Sclose(dataspace_id4);
561  H5Dclose(dataset_id5);
562  H5Sclose(dataspace_id5);
563 
564  H5Gclose(group_id);
565 
566  H5Fclose(file_id);
567 
568  Tmat->save2(filename);
569  Vmat->save2(filename);
570 }
571 
572 void CheMPS2::Hamiltonian::read2(const string filename)
573 {
574  Tmat->read2(filename);
575  Vmat->read2(filename);
576 
577  //The hdf5 file
578  hid_t file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
579 
580  //The data
581  hid_t group_id = H5Gopen(file_id, "/Data",H5P_DEFAULT);
582 
583  //The chain length
584  hid_t dataset_id1 = H5Dopen(group_id, "L", H5P_DEFAULT);
585  int Lagain;
586  H5Dread(dataset_id1, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &Lagain);
587  assert( L==Lagain );
588 
589  //The group number
590  hid_t dataset_id2 = H5Dopen(group_id, "nGroup", H5P_DEFAULT);
591  int nGroup;
592  H5Dread(dataset_id2, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &nGroup);
593  assert( nGroup==SymmInfo.getGroupNumber() );
594 
595  hid_t dataset_id5 = H5Dopen(group_id, "nelectrons", H5P_DEFAULT);
596  H5Dread(dataset_id5, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &Ne);
597 
598  //orb2irrep
599  hid_t dataset_id3 = H5Dopen(group_id, "orb2irrep", H5P_DEFAULT);
600  int * orb2irrepAgain = new int[L];
601  H5Dread(dataset_id3, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, orb2irrepAgain);
602  for (int cnt=0; cnt<L; cnt++){
603  assert( orb2irrep[cnt]==orb2irrepAgain[cnt] );
604  }
605  delete [] orb2irrepAgain;
606 
607  //Econst
608  hid_t dataset_id4 = H5Dopen(group_id, "Econst", H5P_DEFAULT);
609  H5Dread(dataset_id4, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &Econst);
610 
611  H5Dclose(dataset_id1);
612  H5Dclose(dataset_id2);
613  H5Dclose(dataset_id3);
614  H5Dclose(dataset_id4);
615  H5Dclose(dataset_id5);
616 
617  H5Gclose(group_id);
618 
619  H5Fclose(file_id);
620 
621  if (CheMPS2::HAMILTONIAN_debugPrint) debugcheck();
622 
623 }
624 
625 void CheMPS2::Hamiltonian::setNe(int ne)
626 {
627  this->Ne = ne;
628 }
629 
630 
631 CheMPS2::Hamiltonian CheMPS2::Hamiltonian::CreateFromH5(const string filename)
632 {
633  //The hdf5 file
634  hid_t file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
635 
636  //The data
637  hid_t group_id = H5Gopen(file_id, "/Data",H5P_DEFAULT);
638 
639  //The number of orbitals: Set L directly!
640  int L;
641  hid_t dataset_id1 = H5Dopen(group_id, "L", H5P_DEFAULT);
642  H5Dread(dataset_id1, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &L);
643 
644  //The group number: SymmInfo.setGroup()
645  hid_t dataset_id2 = H5Dopen(group_id, "nGroup", H5P_DEFAULT);
646  int nGroup_LOADH5;
647  H5Dread(dataset_id2, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &nGroup_LOADH5);
648 
649  //The irrep of each orbital: Create and fill orb2irrep directly!
650  hid_t dataset_id3 = H5Dopen(group_id, "orb2irrep", H5P_DEFAULT);
651  std::unique_ptr<int []> orb2irrep(new int[L]);
652  H5Dread(dataset_id3, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, orb2irrep.get());
653 
654  H5Dclose(dataset_id1);
655  H5Dclose(dataset_id2);
656  H5Dclose(dataset_id3);
657 
658  H5Gclose(group_id);
659 
660  H5Fclose(file_id);
661 
662  Hamiltonian my_obj(L, nGroup_LOADH5, orb2irrep.get());
663 
664  my_obj.read2(filename);
665 
666  return my_obj;
667 }
668 
669 void CheMPS2::Hamiltonian::reset()
670 {
671  Tmat->reset();
672  Vmat->reset();
673 }
CheMPS2::Hamiltonian::CreateFromH5
static Hamiltonian CreateFromH5(const string filename)
Definition: Hamiltonian.cpp:631
CheMPS2::Hamiltonian::reset
void reset()
set everything to zero
Definition: Hamiltonian.cpp:669
CheMPS2::Irreps::setGroup
bool setGroup(const int nGroup)
Set the group.
Definition: Irreps.cpp:50
CheMPS2::Hamiltonian::orb2indexSy
std::unique_ptr< int[]> orb2indexSy
Definition: Hamiltonian.h:193
CheMPS2::Irreps::getNumberOfIrreps
int getNumberOfIrreps() const
Get the number of irreps for the currently activated group.
Definition: Irreps.cpp:94
CheMPS2::Hamiltonian::setTmat
void setTmat(const int index1, const int index2, const double val)
Set a Tmat element.
Definition: Hamiltonian.cpp:146
CheMPS2::HAMILTONIAN_debugPrint
const bool HAMILTONIAN_debugPrint
Definition: Options.h:55
CheMPS2::Hamiltonian::getOrbitalIrrep
int getOrbitalIrrep(const int nOrb) const
Get an orbital irrep number.
Definition: Hamiltonian.cpp:140
CheMPS2::Hamiltonian::getVmat
double getVmat(const int index1, const int index2, const int index3, const int index4) const
Get a Vmat element.
Definition: Hamiltonian.cpp:177
CheMPS2::Hamiltonian::debugcheck
void debugcheck() const
Debug check certain elements and sums.
Definition: Hamiltonian.cpp:472
CheMPS2::Hamiltonian::setVmat
void setVmat(const int index1, const int index2, const int index3, const int index4, const double val)
Set a Vmat element.
Definition: Hamiltonian.cpp:163
CheMPS2::Hamiltonian::save
void save(const string file_parent=HAMILTONIAN_ParentStorageName, const string file_tmat=HAMILTONIAN_TmatStorageName, const string file_vmat=HAMILTONIAN_VmatStorageName) const
Save the Hamiltonian.
Definition: Hamiltonian.cpp:187
CheMPS2::Hamiltonian::setEconst
void setEconst(const double val)
Set the constant energy.
Definition: Hamiltonian.cpp:142
CheMPS2::Irreps::getGroupNumber
int getGroupNumber() const
Get the group number.
Definition: Irreps.cpp:66
CheMPS2::Hamiltonian::Vmat
std::unique_ptr< FourIndex > Vmat
Definition: Hamiltonian.h:199
CheMPS2::Hamiltonian::read
void read(const string file_parent=HAMILTONIAN_ParentStorageName, const string file_tmat=HAMILTONIAN_TmatStorageName, const string file_vmat=HAMILTONIAN_VmatStorageName)
Load the Hamiltonian.
Definition: Hamiltonian.cpp:238
CheMPS2::Hamiltonian::save2
void save2(const string filename) const
Definition: Hamiltonian.cpp:516
CheMPS2::Hamiltonian::getTmat
double getTmat(const int index1, const int index2) const
Get a Tmat element.
Definition: Hamiltonian.cpp:153
CheMPS2::Hamiltonian::Tmat
std::unique_ptr< TwoIndex > Tmat
Definition: Hamiltonian.h:196
CheMPS2::Hamiltonian::CreateAndFillFromH5
void CreateAndFillFromH5(const string file_parent, const string file_tmat, const string file_vmat)
Definition: Hamiltonian.cpp:287
CheMPS2::Hamiltonian::addToVmat
void addToVmat(const int index1, const int index2, const int index3, const int index4, const double val)
Add to Vmat element.
Definition: Hamiltonian.cpp:170
CheMPS2::Hamiltonian::CreateAndFillFromPsi4dump
void CreateAndFillFromPsi4dump(const string filename)
Definition: Hamiltonian.cpp:341
CheMPS2::Irreps::directProd
static int directProd(const int Irrep1, const int Irrep2)
Get the direct product of the irreps with numbers Irrep1 and Irrep2: a bitwise XOR for Psi4's convent...
Definition: Irreps.h:117
CheMPS2::Hamiltonian::operator=
Hamiltonian & operator=(const Hamiltonian &)
Definition: Hamiltonian.cpp:84
CheMPS2::Hamiltonian::orb2irrep
std::unique_ptr< int[]> orb2irrep
Definition: Hamiltonian.h:187
CheMPS2::Hamiltonian::read2
void read2(const string filename)
Definition: Hamiltonian.cpp:572
CheMPS2::Hamiltonian::getNGroup
int getNGroup() const
Get the group number.
Definition: Hamiltonian.cpp:138
CheMPS2::Hamiltonian::getEconst
double getEconst() const
Get the constant energy.
Definition: Hamiltonian.cpp:144
CheMPS2::Hamiltonian::Hamiltonian
Hamiltonian(const int Norbitals, const int nGroup, const int *OrbIrreps)
Constructor.
Definition: Hamiltonian.cpp:38
CheMPS2::Hamiltonian::getL
int getL() const
Get the number of orbitals.
Definition: Hamiltonian.cpp:136
CheMPS2::Hamiltonian::irrep2num_orb
std::unique_ptr< int[]> irrep2num_orb
Definition: Hamiltonian.h:190