CheMPS2
DMRGSCFrotations.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 <math.h>
22 #include <algorithm>
23 #include <assert.h>
24 #include <string>
25 
26 #include "DMRGSCFrotations.h"
27 #include "Lapack.h"
28 
29 using std::min;
30 using std::max;
31 using std::string;
32 
33 void CheMPS2::DMRGSCFrotations::fetch( double * eri, const FourIndex * ORIG_VMAT, const int irrep1, const int irrep2, const int irrep3, const int irrep4, DMRGSCFindices * idx, const int start, const int stop, const bool pack ){
34 
35  if ( pack ){
36 
37  assert( irrep1 == irrep2 );
38  assert( irrep3 == irrep4 );
39 
40  const int NORB12 = idx->getNORB( irrep1 );
41  const int NORB34 = idx->getNORB( irrep3 );
42 
43  int counter = 0; // counter = cnt3 + ( cnt4 * ( cnt4 + 1 )) / 2
44  for ( int cnt4 = 0; cnt4 < NORB34; cnt4++ ){
45  for ( int cnt3 = 0; cnt3 <= cnt4; cnt3++ ){
46  if (( start <= counter ) && ( counter < stop )){
47  for ( int cnt2 = 0; cnt2 < NORB12; cnt2++ ){
48  for ( int cnt1 = 0; cnt1 < NORB12; cnt1++ ){
49  eri[ cnt1 + NORB12 * ( cnt2 + NORB12 * ( counter - start ) ) ]
50  = ORIG_VMAT->get( irrep1, irrep3, irrep2, irrep4, cnt1, cnt3, cnt2, cnt4 );
51  // Indices (12) and indices (34) are Coulomb pairs
52  }
53  }
54  }
55  counter++;
56  }
57  }
58 
59  } else {
60 
61  assert( Irreps::directProd( irrep1, irrep2 ) == Irreps::directProd( irrep3, irrep4 ) );
62 
63  const int NORB1 = idx->getNORB( irrep1 );
64  const int NORB2 = idx->getNORB( irrep2 );
65  const int NORB3 = idx->getNORB( irrep3 );
66  const int NORB4 = idx->getNORB( irrep4 );
67 
68  int counter = 0; // counter = cnt3 + NORB3 * cnt4
69  for ( int cnt4 = 0; cnt4 < NORB4; cnt4++ ){
70  for ( int cnt3 = 0; cnt3 < NORB3; cnt3++ ){
71  if (( start <= counter ) && ( counter < stop )){
72  for ( int cnt2 = 0; cnt2 < NORB2; cnt2++ ){
73  for ( int cnt1 = 0; cnt1 < NORB1; cnt1++ ){
74  eri[ cnt1 + NORB1 * ( cnt2 + NORB2 * ( counter - start ) ) ]
75  = ORIG_VMAT->get( irrep1, irrep3, irrep2, irrep4, cnt1, cnt3, cnt2, cnt4 );
76  // Indices (12) and indices (34) are Coulomb pairs
77  }
78  }
79  }
80  counter++;
81  }
82  }
83 
84  }
85 
86 }
87 
88 void CheMPS2::DMRGSCFrotations::write( double * eri, FourIndex * NEW_VMAT, DMRGSCFintegrals * ROT_TEI, const char space1, const char space2, const char space3, const char space4, const int irrep1, const int irrep2, const int irrep3, const int irrep4, DMRGSCFindices * idx, const int start, const int stop, const bool pack ){
89 
90  bool written = false;
91  if (( space1 == space2 ) && ( space1 == space3 ) && ( space1 == space4 )){ // All four spaces equal
92 
93  if ( pack ){
94 
95  assert( irrep1 == irrep2 );
96  assert( irrep3 == irrep4 );
97 
98  const int NEW12 = dimension( idx, irrep1, space1 );
99  const int NEW34 = dimension( idx, irrep3, space3 );
100  const int SIZE = stop - start;
101 
102  int counter = 0; // counter = cnt1 + ( cnt2 * ( cnt2 + 1 )) / 2
103  for ( int cnt2 = 0; cnt2 < NEW12; cnt2++ ){
104  for ( int cnt1 = 0; cnt1 <= cnt2; cnt1++ ){
105  if (( start <= counter ) && ( counter < stop )){
106  for ( int cnt4 = 0; cnt4 < NEW34; cnt4++ ){
107  for ( int cnt3 = 0; cnt3 <= cnt4; cnt3++ ){
108  NEW_VMAT->set( irrep1, irrep3, irrep2, irrep4, cnt1, cnt3, cnt2, cnt4,
109  eri[ ( counter - start ) + SIZE * ( cnt3 + NEW34 * cnt4 ) ] );
110  // Indices (12) and indices (34) are Coulomb pairs
111  }
112  }
113  }
114  counter++;
115  }
116  }
117  written = true;
118 
119  } else {
120 
121  assert( Irreps::directProd( irrep1, irrep2 ) == Irreps::directProd( irrep3, irrep4 ) );
122 
123  const int NEW1 = dimension( idx, irrep1, space1 );
124  const int NEW2 = dimension( idx, irrep2, space2 );
125  const int NEW3 = dimension( idx, irrep3, space3 );
126  const int NEW4 = dimension( idx, irrep4, space4 );
127  const int SIZE = stop - start;
128 
129  int counter = 0; // counter = cnt1 + NEW1 * cnt2
130  for ( int cnt2 = 0; cnt2 < NEW2; cnt2++ ){
131  for ( int cnt1 = 0; cnt1 < NEW1; cnt1++ ){
132  if (( start <= counter ) && ( counter < stop )){
133  for ( int cnt4 = 0; cnt4 < NEW4; cnt4++ ){
134  for ( int cnt3 = 0; cnt3 < NEW3; cnt3++ ){
135  NEW_VMAT->set( irrep1, irrep3, irrep2, irrep4, cnt1, cnt3, cnt2, cnt4,
136  eri[ ( counter - start ) + SIZE * ( cnt3 + NEW3 * cnt4 ) ] );
137  // Indices (12) and indices (34) are Coulomb pairs
138  }
139  }
140  }
141  counter++;
142  }
143  }
144  written = true;
145 
146  }
147  }
148 
149  if (( space1 == 'C' ) && ( space2 == 'C' ) && ( space3 == 'F' ) && ( space4 == 'F' )){
150 
151  if ( pack ){
152 
153  assert( irrep1 == irrep2 );
154  assert( irrep3 == irrep4 );
155 
156  const int NEW12 = dimension( idx, irrep1, space1 );
157  const int NEW34 = dimension( idx, irrep3, space3 );
158  const int SIZE = stop - start;
159 
160  int counter = 0; // counter = cnt1 + ( cnt2 * ( cnt2 + 1 )) / 2
161  for ( int cnt2 = 0; cnt2 < NEW12; cnt2++ ){
162  for ( int cnt1 = 0; cnt1 <= cnt2; cnt1++ ){
163  if (( start <= counter ) && ( counter < stop )){
164  for ( int cnt4 = 0; cnt4 < NEW34; cnt4++ ){
165  for ( int cnt3 = 0; cnt3 <= cnt4; cnt3++ ){
166  ROT_TEI->set_coulomb( irrep1, irrep2, irrep3, irrep4, cnt1, cnt2, cnt3, cnt4,
167  eri[ ( counter - start ) + SIZE * ( cnt3 + NEW34 * cnt4 ) ] );
168  // Indices (12) and indices (34) are Coulomb pairs
169  }
170  }
171  }
172  counter++;
173  }
174  }
175  written = true;
176 
177  } else {
178 
179  assert( Irreps::directProd( irrep1, irrep2 ) == Irreps::directProd( irrep3, irrep4 ) );
180 
181  const int NEW1 = dimension( idx, irrep1, space1 );
182  const int NEW2 = dimension( idx, irrep2, space2 );
183  const int NEW3 = dimension( idx, irrep3, space3 );
184  const int NEW4 = dimension( idx, irrep4, space4 );
185  const int SIZE = stop - start;
186 
187  int counter = 0; // counter = cnt1 + NEW1 * cnt2
188  for ( int cnt2 = 0; cnt2 < NEW2; cnt2++ ){
189  for ( int cnt1 = 0; cnt1 < NEW1; cnt1++ ){
190  if (( start <= counter ) && ( counter < stop )){
191  for ( int cnt4 = 0; cnt4 < NEW4; cnt4++ ){
192  for ( int cnt3 = 0; cnt3 < NEW3; cnt3++ ){
193  ROT_TEI->set_coulomb( irrep1, irrep2, irrep3, irrep4, cnt1, cnt2, cnt3, cnt4,
194  eri[ ( counter - start ) + SIZE * ( cnt3 + NEW3 * cnt4 ) ] );
195  // Indices (12) and indices (34) are Coulomb pairs
196  }
197  }
198  }
199  counter++;
200  }
201  }
202  written = true;
203 
204  }
205  }
206 
207  if (( space1 == 'C' ) && ( space2 == 'V' ) && ( space3 == 'C' ) && ( space4 == 'V' )){ // ( C V | C V )
208 
209  assert( pack == false );
210  assert( Irreps::directProd( irrep1, irrep2 ) == Irreps::directProd( irrep3, irrep4 ) );
211 
212  const int NEW1 = dimension( idx, irrep1, space1 );
213  const int NEW2 = dimension( idx, irrep2, space2 );
214  const int NEW3 = dimension( idx, irrep3, space3 );
215  const int NEW4 = dimension( idx, irrep4, space4 );
216  const int JUMP2 = jump( idx, irrep2, space2 );
217  const int JUMP4 = jump( idx, irrep4, space4 );
218  const int SIZE = stop - start;
219 
220  int counter = 0; // counter = cnt1 + NEW1 * cnt2
221  for ( int cnt2 = 0; cnt2 < NEW2; cnt2++ ){
222  for ( int cnt1 = 0; cnt1 < NEW1; cnt1++ ){
223  if (( start <= counter ) && ( counter < stop )){
224  for ( int cnt4 = 0; cnt4 < NEW4; cnt4++ ){
225  for ( int cnt3 = 0; cnt3 < NEW3; cnt3++ ){
226  ROT_TEI->set_exchange( irrep1, irrep3, irrep2, irrep4, cnt1, cnt3, JUMP2 + cnt2, JUMP4 + cnt4,
227  eri[ ( counter - start ) + SIZE * ( cnt3 + NEW3 * cnt4 ) ] );
228  // Indices (12) and indices (34) are Coulomb pairs
229  }
230  }
231  }
232  counter++;
233  }
234  }
235  written = true;
236 
237  }
238 
239  assert( written == true );
240 
241 }
242 
243 void CheMPS2::DMRGSCFrotations::blockwise_first( double * origin, double * target, int orig1, int dim2, const int dim34, double * umat1, int new1, int lda1 ){
244 
245  char notrans = 'N';
246  double one = 1.0;
247  double set = 0.0;
248  int right_dim = dim2 * dim34;
249  dgemm_( &notrans, &notrans, &new1, &right_dim, &orig1, &one, umat1, &lda1, origin, &orig1, &set, target, &new1 );
250 
251 }
252 
253 void CheMPS2::DMRGSCFrotations::blockwise_second( double * origin, double * target, int dim1, int orig2, const int dim34, double * umat2, int new2, int lda2 ){
254 
255  char trans = 'T';
256  char notrans = 'N';
257  double one = 1.0;
258  double set = 0.0;
259  const int jump_old = dim1 * orig2;
260  const int jump_new = dim1 * new2;
261  #pragma omp parallel for schedule(static)
262  for ( int index = 0; index < dim34; index++ ){
263  dgemm_( &notrans, &trans, &dim1, &new2, &orig2, &one, origin + jump_old * index, &dim1, umat2, &lda2, &set, target + jump_new * index, &dim1 );
264  }
265 
266 }
267 
268 void CheMPS2::DMRGSCFrotations::blockwise_third( double * origin, double * target, const int dim12, int orig3, const int dim4, double * umat3, int new3, int lda3 ){
269 
270  char trans = 'T';
271  char notrans = 'N';
272  double one = 1.0;
273  double set = 0.0;
274  int left_dim = dim12;
275  const int jump_old = dim12 * orig3;
276  const int jump_new = dim12 * new3;
277  #pragma omp parallel for schedule(static)
278  for ( int index = 0; index < dim4; index++ ){
279  dgemm_( &notrans, &trans, &left_dim, &new3, &orig3, &one, origin + jump_old * index, &left_dim, umat3, &lda3, &set, target + jump_new * index, &left_dim );
280  }
281 
282 }
283 
284 void CheMPS2::DMRGSCFrotations::blockwise_fourth( double * origin, double * target, const int dim12, int dim3, int orig4, double * umat4, int new4, int lda4 ){
285 
286  char trans = 'T';
287  char notrans = 'N';
288  double one = 1.0;
289  double set = 0.0;
290  int left_dim = dim12 * dim3;
291  dgemm_( &notrans, &trans, &left_dim, &new4, &orig4, &one, origin, &left_dim, umat4, &lda4, &set, target, &left_dim );
292 
293 }
294 
295 void CheMPS2::DMRGSCFrotations::open_file( hid_t * file_id, hid_t * dspc_id, hid_t * dset_id, const int first, const int second, const string filename ){
296 
297  file_id[0] = H5Fcreate( filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT );
298  hsize_t fdim_h5[] = { (hsize_t) second, (hsize_t) first }; // C is row major: [ col + ncol * row ] is assumed
299  dspc_id[0] = H5Screate_simple( 2, fdim_h5, NULL );
300  dset_id[0] = H5Dcreate( file_id[0], "storage", H5T_NATIVE_DOUBLE, dspc_id[0], H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT );
301 
302 }
303 
304 void CheMPS2::DMRGSCFrotations::close_file( hid_t file_id, hid_t dspc_id, hid_t dset_id ){
305 
306  H5Dclose( dset_id );
307  H5Sclose( dspc_id );
308  H5Fclose( file_id );
309 
310 }
311 
312 void CheMPS2::DMRGSCFrotations::write_file( hid_t dspc_id, hid_t dset_id, double * eri, const int start, const int size, const int first_write ){
313 
314  hsize_t stride_h5[] = { 1, 1 };
315  hsize_t count_h5[] = { 1, 1 };
316  hsize_t start_h5[] = { (hsize_t) start, 0 };
317  hsize_t block_h5[] = { (hsize_t) size, (hsize_t) first_write };
318  H5Sselect_hyperslab( dspc_id, H5S_SELECT_SET, start_h5, stride_h5, count_h5, block_h5 );
319  hsize_t mem_h5 = size * first_write; // Should be OK to multiply as integers as it is smaller than mem_size
320  hid_t mem_id = H5Screate_simple( 1, &mem_h5, NULL );
321  H5Dwrite( dset_id, H5T_NATIVE_DOUBLE, mem_id, dspc_id, H5P_DEFAULT, eri );
322  H5Sclose( mem_id );
323 
324 }
325 
326 void CheMPS2::DMRGSCFrotations::read_file( hid_t dspc_id, hid_t dset_id, double * eri, const int start, const int size, const int second_read ){
327 
328  hsize_t stride_h5[] = { 1, 1 };
329  hsize_t count_h5[] = { 1, 1 };
330  hsize_t start_h5[] = { 0, (hsize_t) start };
331  hsize_t block_h5[] = { (hsize_t) second_read, (hsize_t) size };
332  H5Sselect_hyperslab( dspc_id, H5S_SELECT_SET, start_h5, stride_h5, count_h5, block_h5 );
333  hsize_t mem_h5 = second_read * size; // Should be OK to multiply as integers as it is smaller than mem_size
334  hid_t mem_id = H5Screate_simple( 1, &mem_h5, NULL );
335  H5Dread( dset_id, H5T_NATIVE_DOUBLE, mem_id, dspc_id, H5P_DEFAULT, eri );
336  H5Sclose( mem_id );
337 
338 }
339 
340 int CheMPS2::DMRGSCFrotations::dimension( DMRGSCFindices * idx, const int irrep, const char space ){
341 
342  if ( space == 'O' ){ return idx->getNOCC( irrep ); }
343  if ( space == 'A' ){ return idx->getNDMRG( irrep ); }
344  if ( space == 'V' ){ return idx->getNVIRT( irrep ); }
345  if ( space == 'C' ){ return ( idx->getNOCC( irrep ) + idx->getNDMRG( irrep ) ); }
346  if ( space == 'F' ){ return idx->getNORB( irrep ); }
347  return -1;
348 
349 }
350 
351 int CheMPS2::DMRGSCFrotations::jump( DMRGSCFindices * idx, const int irrep, const char space ){
352 
353  if ( space == 'A' ){ return idx->getNOCC( irrep ); }
354  if ( space == 'V' ){ return idx->getNOCC( irrep ) + idx->getNDMRG( irrep ); }
355  return 0; // O, F, C
356 
357 }
358 
359 void CheMPS2::DMRGSCFrotations::unpackage_second( double * mem1, double * mem2, const int SIZE, const int ORIG ){
360 
361  #pragma omp parallel for schedule(static)
362  for ( int cnt4 = 0; cnt4 < ORIG; cnt4++ ){
363  for ( int cnt3 = 0; cnt3 < ORIG; cnt3++ ){
364  const int combined = (( cnt3 < cnt4 ) ? ( cnt3 + ( cnt4 * ( cnt4 + 1 )) / 2 )
365  : ( cnt4 + ( cnt3 * ( cnt3 + 1 )) / 2 ));
366  #pragma omp simd
367  for ( int cnt12 = 0; cnt12 < SIZE; cnt12++ ){
368  mem2[ cnt12 + SIZE * ( cnt3 + ORIG * cnt4 ) ] = mem1[ cnt12 + SIZE * combined ];
369  }
370  }
371  }
372 
373 }
374 
375 void CheMPS2::DMRGSCFrotations::package_first( double * mem1, double * mem2, const int NEW, const int PACKED, const int SIZE ){
376 
377  #pragma omp parallel for schedule(static)
378  for ( int cnt34 = 0; cnt34 < SIZE; cnt34++ ){
379  for ( int cnt2 = 0; cnt2 < NEW; cnt2++ ){
380  #pragma omp simd
381  for ( int cnt1 = 0; cnt1 <= cnt2; cnt1++ ){
382  mem2[ cnt1 + ( cnt2 * ( cnt2 + 1 ))/2 + PACKED * cnt34 ] = mem1[ cnt1 + NEW * ( cnt2 + NEW * cnt34 ) ];
383  }
384  }
385  }
386 
387 }
388 
389 void CheMPS2::DMRGSCFrotations::rotate( const FourIndex * ORIG_VMAT, FourIndex * NEW_VMAT, DMRGSCFintegrals * ROT_TEI, const char space1, const char space2, const char space3, const char space4, DMRGSCFindices * idx, DMRGSCFunitary * umat, double * mem1, double * mem2, const int mem_size, const string filename ){
390 
391  /* Matrix elements ( 1 2 | 3 4 ) */
392 
393  assert(( space1 == 'O' ) || ( space1 == 'A' ) || ( space1 == 'V' ) || ( space1 == 'C' ) || ( space1 == 'F' ));
394  assert(( space2 == 'O' ) || ( space2 == 'A' ) || ( space2 == 'V' ) || ( space2 == 'C' ) || ( space2 == 'F' ));
395  assert(( space3 == 'O' ) || ( space3 == 'A' ) || ( space3 == 'V' ) || ( space3 == 'C' ) || ( space3 == 'F' ));
396  assert(( space4 == 'O' ) || ( space4 == 'A' ) || ( space4 == 'V' ) || ( space4 == 'C' ) || ( space4 == 'F' ));
397 
398  const int num_irreps = idx->getNirreps();
399  const bool equal12 = ( space1 == space2 );
400  const bool equal34 = ( space3 == space4 );
401  const bool eightfold = (( space1 == space3 ) && ( space2 == space4 ));
402 
403  for ( int irrep1 = 0; irrep1 < num_irreps; irrep1++ ){
404  for ( int irrep2 = (( equal12 ) ? irrep1 : 0 ); irrep2 < num_irreps; irrep2++ ){ // irrep2 >= irrep1 if space1 == space2
405  const int product_symm = Irreps::directProd( irrep1, irrep2 );
406  for ( int irrep3 = (( eightfold ) ? irrep1 : 0 ); irrep3 < num_irreps; irrep3++ ){
407  const int irrep4 = Irreps::directProd( product_symm, irrep3 );
408  if ( irrep4 >= (( equal34 ) ? irrep3 : 0 ) ){ // irrep4 >= irrep3 if space3 == space4
409 
410  const int NEW1 = dimension( idx, irrep1, space1 );
411  const int NEW2 = dimension( idx, irrep2, space2 );
412  const int NEW3 = dimension( idx, irrep3, space3 );
413  const int NEW4 = dimension( idx, irrep4, space4 );
414 
415  if (( NEW1 > 0 ) && ( NEW2 > 0 ) && ( NEW3 > 0 ) && ( NEW4 > 0 )){
416 
417  const int ORIG1 = idx->getNORB( irrep1 );
418  const int ORIG2 = idx->getNORB( irrep2 );
419  const int ORIG3 = idx->getNORB( irrep3 );
420  const int ORIG4 = idx->getNORB( irrep4 );
421 
422  double * umat1 = umat->getBlock( irrep1 ) + jump( idx, irrep1, space1 );
423  double * umat2 = umat->getBlock( irrep2 ) + jump( idx, irrep2, space2 );
424  double * umat3 = umat->getBlock( irrep3 ) + jump( idx, irrep3, space3 );
425  double * umat4 = umat->getBlock( irrep4 ) + jump( idx, irrep4, space4 );
426 
427  const int block_size1 = mem_size / ( ORIG1 * ORIG2 ); // Floor of amount of times orig( first ) fits in mem_size
428  const int block_size2 = mem_size / ( ORIG3 * ORIG4 ); // Floor of amount of times orig( second ) fits in mem_size
429  assert( block_size1 > 0 );
430  assert( block_size2 > 0 );
431 
432  const bool pack_first = (( equal12 ) && ( irrep1 == irrep2 ));
433  const bool pack_second = (( equal34 ) && ( irrep3 == irrep4 ));
434  const int first_size = (( pack_first ) ? ( NEW1 * ( NEW1 + 1 )) / 2 : NEW1 * NEW2 );
435  const int second_size = (( pack_second ) ? ( ORIG3 * ( ORIG3 + 1 )) / 2 : ORIG3 * ORIG4 );
436 
437  const bool io_free = (( block_size1 >= second_size ) && ( block_size2 >= first_size ));
438  hid_t file_id, dspc_id, dset_id;
439 
440  if ( io_free == false ){
441  assert( filename.compare( "edmistonruedenberg" ) != 0 );
442  open_file( &file_id, &dspc_id, &dset_id, first_size, second_size, filename );
443  }
444 
445  // First half transformation
446  int start = 0;
447  while ( start < second_size ){
448  const int stop = min( start + block_size1, second_size );
449  const int size = stop - start;
450  fetch( mem1, ORIG_VMAT, irrep1, irrep2, irrep3, irrep4, idx, start, stop, pack_second );
451  blockwise_first( mem1, mem2, ORIG1, ORIG2, size, umat1, NEW1, ORIG1 );
452  blockwise_second( mem2, mem1, NEW1, ORIG2, size, umat2, NEW2, ORIG2 );
453  if ( pack_first ){
454  package_first( mem1, mem2, NEW1, first_size, size );
455  double * temp = mem1;
456  mem1 = mem2;
457  mem2 = temp;
458  }
459  if ( io_free == false ){ write_file( dspc_id, dset_id, mem1, start, size, first_size ); }
460  start += size;
461  }
462  assert( start == second_size );
463 
464  // Do the second half transformation
465  start = 0;
466  while ( start < first_size ){
467  const int stop = min( start + block_size2, first_size );
468  const int size = stop - start;
469  if ( io_free == false ){ read_file( dspc_id, dset_id, mem1, start, size, second_size ); }
470  if ( pack_second ){
471  unpackage_second( mem1, mem2, size, ORIG3 );
472  double * temp = mem1;
473  mem1 = mem2;
474  mem2 = temp;
475  }
476  blockwise_fourth( mem1, mem2, size, ORIG3, ORIG4, umat4, NEW4, ORIG4 );
477  blockwise_third( mem2, mem1, size, ORIG3, NEW4, umat3, NEW3, ORIG3 );
478  write( mem1, NEW_VMAT, ROT_TEI, space1, space2, space3, space4, irrep1, irrep2, irrep3, irrep4, idx, start, stop, pack_first );
479  start += size;
480  }
481  assert( start == first_size );
482  if ( io_free == false ){ close_file( file_id, dspc_id, dset_id ); }
483  }
484  }
485  }
486  }
487  }
488 
489 }
490 
491 
CheMPS2::DMRGSCFindices::getNORB
int getNORB(const int irrep) const
Get the number of orbitals for an irrep.
Definition: DMRGSCFindices.cpp:98
CheMPS2::DMRGSCFmatrix::getBlock
double * getBlock(const int irrep)
Get a matrix block.
Definition: DMRGSCFmatrix.cpp:87
CheMPS2::DMRGSCFrotations::rotate
static void rotate(const FourIndex *ORIG_VMAT, FourIndex *NEW_VMAT, DMRGSCFintegrals *ROT_TEI, const char space1, const char space2, const char space3, const char space4, DMRGSCFindices *idx, DMRGSCFunitary *umat, double *mem1, double *mem2, const int mem_size, const string filename)
Fill the rotated two-body matrix elements for the space. If the blocks become too large...
Definition: DMRGSCFrotations.cpp:389
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&#39;s convent...
Definition: Irreps.h:123
CheMPS2::DMRGSCFindices::getNirreps
int getNirreps() const
Get the number of irreps.
Definition: DMRGSCFindices.cpp:96