1 /*
2 CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3 Copyright (C) 2013-2018 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 <iostream>
22 #include <fstream>
23 #include <string>
24 #include <math.h>
25 #include <algorithm>
26 #include <assert.h>
27
28 #include "CASSCF.h"
29 #include "Lapack.h"
30 #include "Special.h"
31 #include "MPIchemps2.h"
32
33 using std::string;
34 using std::ifstream;
35 using std::cout;
36 using std::endl;
37 using std::max;
38
CASSCF(Hamiltonian * ham_in,int * docc,int * socc,int * nocc,int * ndmrg,int * nvirt,const string new_tmp_folder)39 CheMPS2::CASSCF::CASSCF( Hamiltonian * ham_in, int * docc, int * socc, int * nocc, int * ndmrg, int * nvirt, const string new_tmp_folder ){
40
41 #ifdef CHEMPS2_MPI_COMPILATION
42 const bool am_i_master = ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER );
43 #else
44 const bool am_i_master = true;
45 #endif
46
47 NUCL_ORIG = ham_in->getEconst();
48 TMAT_ORIG = ham_in->getTmat();
49 VMAT_ORIG = ham_in->getVmat();
50
51 L = ham_in->getL();
52 SymmInfo.setGroup( ham_in->getNGroup() );
53 num_irreps = SymmInfo.getNumberOfIrreps();
54 successful_solve = false;
55
56 if (( am_i_master ) && ( docc != NULL ) && ( socc != NULL )){
57 cout << "DOCC = [ ";
58 for ( int irrep = 0; irrep < num_irreps-1; irrep++ ){ cout << docc[ irrep ] << " , "; }
59 cout << docc[ num_irreps - 1 ] << " ]" << endl;
60 cout << "SOCC = [ ";
61 for ( int irrep = 0; irrep < num_irreps-1; irrep++ ){ cout << socc[ irrep ] << " , "; }
62 cout << socc[ num_irreps - 1 ] << " ]" << endl;
63 }
64
65 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
66 const int norb_in = nocc[ irrep ] + ndmrg[ irrep ] + nvirt[ irrep ];
67 const int norb_ham = VMAT_ORIG->get_irrep_size( irrep );
68 if (( norb_ham != norb_in ) && ( am_i_master )){
69 cout << "CASSCF::CASSCF : nocc[" << irrep << "] + ndmrg[" << irrep << "] + nvirt[" << irrep << "] = " << norb_in
70 << " and in the Hamiltonian norb[" << irrep << "] = " << norb_ham << "." << endl;
71 }
72 assert( norb_ham == norb_in );
73 }
74
75 iHandler = new DMRGSCFindices( L, SymmInfo.getGroupNumber(), nocc, ndmrg, nvirt );
76 unitary = new DMRGSCFunitary( iHandler );
77
78 // Allocate space for the DMRG 1DM and 2DM
79 nOrbDMRG = iHandler->getDMRGcumulative( num_irreps );
80 DMRG1DM = new double[nOrbDMRG * nOrbDMRG];
81 DMRG2DM = new double[nOrbDMRG * nOrbDMRG * nOrbDMRG * nOrbDMRG];
82
83 // To store the F-matrix and Q-matrix(occ,act)
84 theFmatrix = new DMRGSCFmatrix( iHandler ); theFmatrix->clear();
85 theQmatOCC = new DMRGSCFmatrix( iHandler ); theQmatOCC->clear();
86 theQmatACT = new DMRGSCFmatrix( iHandler ); theQmatACT->clear();
87 theQmatWORK= new DMRGSCFmatrix( iHandler ); theQmatWORK->clear();
88 theTmatrix = new DMRGSCFmatrix( iHandler ); theTmatrix->clear();
89
90 if ( am_i_master ){
91 if (( docc != NULL ) && ( socc != NULL )){ checkHF( docc, socc ); } // Print the MO info. This requires the iHandler to be created...
92 iHandler->Print();
93 cout << "DMRGSCF::setupStart : Number of variables in the x-matrix = " << unitary->getNumVariablesX() << endl;
94 }
95
96 this->tmp_folder = new_tmp_folder;
97
98 }
99
~CASSCF()100 CheMPS2::CASSCF::~CASSCF(){
101
102 delete [] DMRG1DM;
103 delete [] DMRG2DM;
104
105 //The following objects depend on iHandler: delete them first
106 delete theFmatrix;
107 delete theQmatOCC;
108 delete theQmatACT;
109 delete theQmatWORK;
110 delete theTmatrix;
111 delete unitary;
112
113 delete iHandler;
114
115 }
116
get_num_irreps()117 int CheMPS2::CASSCF::get_num_irreps(){ return num_irreps; }
118
copy2DMover(TwoDM * theDMRG2DM,const int LAS,double * two_dm)119 void CheMPS2::CASSCF::copy2DMover( TwoDM * theDMRG2DM, const int LAS, double * two_dm ){
120
121 for ( int i1 = 0; i1 < LAS; i1++ ){
122 for ( int i2 = 0; i2 < LAS; i2++ ){
123 for ( int i3 = 0; i3 < LAS; i3++ ){
124 for ( int i4 = 0; i4 < LAS; i4++ ){
125 // The assignment has been changed to an addition for state-averaged calculations!
126 two_dm[ i1 + LAS * ( i2 + LAS * ( i3 + LAS * i4 ) ) ] += theDMRG2DM->getTwoDMA_HAM( i1, i2, i3, i4 );
127 }
128 }
129 }
130 }
131
132 }
133
setDMRG1DM(const int num_elec,const int LAS,double * one_dm,double * two_dm)134 void CheMPS2::CASSCF::setDMRG1DM( const int num_elec, const int LAS, double * one_dm, double * two_dm ){
135
136 #ifdef CHEMPS2_MPI_COMPILATION
137 const bool am_i_master = ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER );
138 #else
139 const bool am_i_master = true;
140 #endif
141
142 if ( am_i_master ){
143 const double prefactor = 1.0 / ( num_elec - 1 );
144 for ( int cnt1 = 0; cnt1 < LAS; cnt1++ ){
145 for ( int cnt2 = cnt1; cnt2 < LAS; cnt2++ ){
146 double value = 0.0;
147 for ( int sum = 0; sum < LAS; sum++ ){ value += two_dm[ cnt1 + LAS * ( sum + LAS * ( cnt2 + LAS * sum ) ) ]; }
148 one_dm[ cnt1 + LAS * cnt2 ] = prefactor * value;
149 one_dm[ cnt2 + LAS * cnt1 ] = one_dm[ cnt1 + LAS * cnt2 ];
150 }
151 }
152 }
153
154 #ifdef CHEMPS2_MPI_COMPILATION
155 MPIchemps2::broadcast_array_double( one_dm, LAS * LAS, MPI_CHEMPS2_MASTER );
156 #endif
157
158 }
159
fillLocalizedOrbitalRotations(DMRGSCFunitary * umat,DMRGSCFindices * idx,double * eigenvecs)160 void CheMPS2::CASSCF::fillLocalizedOrbitalRotations( DMRGSCFunitary * umat, DMRGSCFindices * idx, double * eigenvecs ){
161
162 const int n_irreps = idx->getNirreps();
163 const int tot_dmrg = idx->getDMRGcumulative( n_irreps );
164 const int size = tot_dmrg * tot_dmrg;
165
166 for ( int cnt = 0; cnt < size; cnt++ ){ eigenvecs[ cnt ] = 0.0; }
167 int passed = 0;
168 for ( int irrep = 0; irrep < n_irreps; irrep++ ){
169
170 const int NDMRG = idx->getNDMRG( irrep );
171 if ( NDMRG > 0 ){
172
173 double * Ublock = umat->getBlock( irrep );
174 double * Eblock = eigenvecs + passed * ( 1 + tot_dmrg );
175
176 for ( int row = 0; row < NDMRG; row++ ){
177 for ( int col = 0; col < NDMRG; col++ ){
178 Eblock[ row + tot_dmrg * col ] = Ublock[ col + NDMRG * row ]; //Eigs = Unit^T
179 }
180 }
181
182 }
183
184 passed += NDMRG;
185
186 }
187
188 }
189
rotateOldToNew(DMRGSCFmatrix * myMatrix)190 void CheMPS2::CASSCF::rotateOldToNew( DMRGSCFmatrix * myMatrix ){
191
192 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
193
194 int NORB = iHandler->getNORB( irrep );
195 if ( NORB > 0 ){
196 double * Umat = unitary->getBlock( irrep );
197 double * work = theQmatWORK->getBlock( irrep );
198 double * block = myMatrix->getBlock( irrep );
199 double one = 1.0;
200 double set = 0.0;
201 char trans = 'T';
202 char notrans = 'N';
203 dgemm_( ¬rans, ¬rans, &NORB, &NORB, &NORB, &one, Umat, &NORB, block, &NORB, &set, work, &NORB );
204 dgemm_( ¬rans, &trans, &NORB, &NORB, &NORB, &one, work, &NORB, Umat, &NORB, &set, block, &NORB );
205 }
206 }
207
208 }
209
constructCoulombAndExchangeMatrixInOrigIndices(DMRGSCFmatrix * density,DMRGSCFmatrix * result)210 void CheMPS2::CASSCF::constructCoulombAndExchangeMatrixInOrigIndices( DMRGSCFmatrix * density, DMRGSCFmatrix * result ){
211
212 for ( int irrepQ = 0; irrepQ < num_irreps; irrepQ++ ){
213
214 const int linearsizeQ = iHandler->getNORB( irrepQ );
215 const int triangsizeQ = ( linearsizeQ * ( linearsizeQ + 1 ) ) / 2;
216 int myindices[ 2 ];
217
218 #pragma omp parallel for schedule(static)
219 for ( int combinedindex = 0; combinedindex < triangsizeQ; combinedindex++ ){
220
221 Special::invert_triangle_two( combinedindex, myindices );
222 const int rowQ = myindices[ 0 ];
223 const int colQ = myindices[ 1 ];
224
225 double value = 0.0;
226
227 for ( int irrepN = 0; irrepN < num_irreps; irrepN++ ){
228 const int linearsizeN = iHandler->getNORB( irrepN );
229 for ( int rowN = 0; rowN < linearsizeN; rowN++ ){
230
231 value += density->get( irrepN, rowN, rowN ) * ( VMAT_ORIG->get( irrepQ, irrepN, irrepQ, irrepN, rowQ, rowN, colQ, rowN )
232 - 0.5 * VMAT_ORIG->get( irrepQ, irrepQ, irrepN, irrepN, rowQ, colQ, rowN, rowN ) );
233
234 for ( int colN = rowN + 1; colN < linearsizeN; colN++ ){
235
236 value += density->get( irrepN, rowN, colN ) * ( 2 * VMAT_ORIG->get( irrepQ, irrepN, irrepQ, irrepN, rowQ, rowN, colQ, colN )
237 - 0.5 * VMAT_ORIG->get( irrepQ, irrepQ, irrepN, irrepN, rowQ, colQ, rowN, colN )
238 - 0.5 * VMAT_ORIG->get( irrepQ, irrepQ, irrepN, irrepN, rowQ, colQ, colN, rowN ) );
239
240 }
241 }
242 }
243
244 result->set( irrepQ, rowQ, colQ, value );
245 result->set( irrepQ, colQ, rowQ, value );
246
247 }
248 }
249
250 }
251
buildQmatOCC()252 void CheMPS2::CASSCF::buildQmatOCC(){
253
254 #ifdef CHEMPS2_MPI_COMPILATION
255 const bool am_i_master = ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER );
256 #else
257 const bool am_i_master = true;
258 #endif
259
260 if ( am_i_master ){
261 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
262
263 int NORB = iHandler->getNORB( irrep );
264 if ( NORB > 0 ){
265 int NOCC = iHandler->getNOCC( irrep );
266 double two = 2.0;
267 double set = 0.0;
268 char trans = 'T';
269 char notrans = 'N';
270 double * Umat = unitary->getBlock( irrep );
271 double * work = theQmatWORK->getBlock( irrep );
272 dgemm_( &trans, ¬rans, &NORB, &NORB, &NOCC, &two, Umat, &NORB, Umat, &NORB, &set, work, &NORB );
273 }
274 }
275 constructCoulombAndExchangeMatrixInOrigIndices( theQmatWORK, theQmatOCC );
276 rotateOldToNew( theQmatOCC );
277 }
278
279 #ifdef CHEMPS2_MPI_COMPILATION
280 theQmatOCC->broadcast( MPI_CHEMPS2_MASTER );
281 #endif
282
283 }
284
buildQmatACT()285 void CheMPS2::CASSCF::buildQmatACT(){
286
287 #ifdef CHEMPS2_MPI_COMPILATION
288 const bool am_i_master = ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER );
289 #else
290 const bool am_i_master = true;
291 #endif
292
293 if ( am_i_master ){
294 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
295
296 int NORB = iHandler->getNORB( irrep );
297 if ( NORB > 0 ){
298 int NACT = iHandler->getNDMRG( irrep );
299 double one = 1.0;
300 double set = 0.0;
301 char trans = 'T';
302 char notrans = 'N';
303 double * Umat = unitary->getBlock( irrep ) + iHandler->getNOCC( irrep );
304 double * work = theQmatWORK->getBlock( irrep );
305 double * work2 = theQmatACT->getBlock( irrep );
306 double * RDM = DMRG1DM + iHandler->getDMRGcumulative( irrep ) * ( 1 + nOrbDMRG );
307 dgemm_( &trans, ¬rans, &NORB, &NACT, &NACT, &one, Umat, &NORB, RDM, &nOrbDMRG, &set, work2, &NORB );
308 dgemm_( ¬rans, ¬rans, &NORB, &NORB, &NACT, &one, work2, &NORB, Umat, &NORB, &set, work, &NORB );
309 }
310 }
311 constructCoulombAndExchangeMatrixInOrigIndices( theQmatWORK, theQmatACT );
312 rotateOldToNew( theQmatACT );
313 }
314
315 #ifdef CHEMPS2_MPI_COMPILATION
316 theQmatACT->broadcast( MPI_CHEMPS2_MASTER );
317 #endif
318
319 }
320
buildTmatrix()321 void CheMPS2::CASSCF::buildTmatrix(){
322
323 #ifdef CHEMPS2_MPI_COMPILATION
324 const bool am_i_master = ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER );
325 #else
326 const bool am_i_master = true;
327 #endif
328
329 if ( am_i_master ){
330 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
331 const int NumORB = iHandler->getNORB( irrep );
332 for ( int row = 0; row < NumORB; row++ ){
333 for ( int col = 0; col < NumORB; col++ ){
334 theTmatrix->set( irrep, row, col, TMAT_ORIG->get( irrep, row, col ) );
335 }
336 }
337 }
338 rotateOldToNew( theTmatrix );
339 }
340
341 #ifdef CHEMPS2_MPI_COMPILATION
342 theTmatrix->broadcast( MPI_CHEMPS2_MASTER );
343 #endif
344
345 }
346
fillConstAndTmatDMRG(Hamiltonian * HamDMRG) const347 void CheMPS2::CASSCF::fillConstAndTmatDMRG( Hamiltonian * HamDMRG ) const{
348
349 //Constant part of the energy
350 double constant = NUCL_ORIG;
351 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
352 const int NOCC = iHandler->getNOCC( irrep );
353 for ( int occ = 0; occ < NOCC; occ++ ){
354 constant += ( 2 * theTmatrix->get( irrep, occ, occ )
355 + theQmatOCC->get( irrep, occ, occ ) );
356 }
357 }
358 HamDMRG->setEconst( constant );
359
360 //One-body terms: diagonal in the irreps
361 for ( int irrep = 0; irrep < num_irreps; irrep++ ){
362 const int JUMP = iHandler->getDMRGcumulative( irrep );
363 const int NACT = iHandler->getNDMRG( irrep );
364 const int NOCC = iHandler->getNOCC( irrep );
365 for ( int cnt1 = 0; cnt1 < NACT; cnt1++ ){
366 for ( int cnt2 = cnt1; cnt2 < NACT; cnt2++ ){
367 HamDMRG->setTmat( JUMP + cnt1, JUMP + cnt2, ( theTmatrix->get( irrep, NOCC + cnt1, NOCC + cnt2 )
368 + theQmatOCC->get( irrep, NOCC + cnt1, NOCC + cnt2 ) ) );
369 }
370 }
371 }
372
373 }
374
copy_active(double * origin,DMRGSCFmatrix * result,const DMRGSCFindices * idx,const bool one_rdm)375 void CheMPS2::CASSCF::copy_active( double * origin, DMRGSCFmatrix * result, const DMRGSCFindices * idx, const bool one_rdm ){
376
377 result->clear();
378
379 const int n_irreps = idx->getNirreps();
380 const int tot_dmrg = idx->getDMRGcumulative( n_irreps );
381
382 int passed = 0;
383 for ( int irrep = 0; irrep < n_irreps; irrep++ ){
384
385 const int NOCC = idx->getNOCC( irrep );
386 const int NACT = idx->getNDMRG( irrep );
387
388 if ( one_rdm ){
389 for ( int orb = 0; orb < NOCC; orb++ ){
390 result->set( irrep, orb, orb, 2.0 );
391 }
392 }
393
394 for ( int row = 0; row < NACT; row++ ){
395 for ( int col = 0; col < NACT; col++ ){
396 result->set( irrep, NOCC + row, NOCC + col, origin[ passed + row + tot_dmrg * ( passed + col ) ] );
397 }
398 }
399
400 passed += NACT;
401 }
402
403 assert( passed == tot_dmrg );
404
405 }
406
copy_active(const DMRGSCFmatrix * origin,double * result,const DMRGSCFindices * idx)407 void CheMPS2::CASSCF::copy_active( const DMRGSCFmatrix * origin, double * result, const DMRGSCFindices * idx ){
408
409 const int n_irreps = idx->getNirreps();
410 const int tot_dmrg = idx->getDMRGcumulative( n_irreps );
411
412 for ( int cnt = 0; cnt < tot_dmrg * tot_dmrg; cnt++ ){ result[ cnt ] = 0.0; }
413
414 int passed = 0;
415 for ( int irrep = 0; irrep < n_irreps; irrep++ ){
416
417 const int NOCC = idx->getNOCC( irrep );
418 const int NACT = idx->getNDMRG( irrep );
419
420 for ( int row = 0; row < NACT; row++ ){
421 for ( int col = 0; col < NACT; col++ ){
422 result[ passed + row + tot_dmrg * ( passed + col ) ] = origin->get( irrep, NOCC + row, NOCC + col );
423 }
424 }
425
426 passed += NACT;
427 }
428
429 assert( passed == tot_dmrg );
430
431 }
432
block_diagonalize(const char space,const DMRGSCFmatrix * Mat,DMRGSCFunitary * Umat,double * work1,double * work2,const DMRGSCFindices * idx,const bool invert,double * two_dm,double * three_dm,double * contract)433 void CheMPS2::CASSCF::block_diagonalize( const char space, const DMRGSCFmatrix * Mat, DMRGSCFunitary * Umat, double * work1, double * work2, const DMRGSCFindices * idx, const bool invert, double * two_dm, double * three_dm, double * contract ){
434
435 #ifdef CHEMPS2_MPI_COMPILATION
436 const bool am_i_master = ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER );
437 #else
438 const bool am_i_master = true;
439 #endif
440
441 const int n_irreps = idx->getNirreps();
442 const int tot_dmrg = idx->getDMRGcumulative( n_irreps );
443
444 if ( am_i_master ){
445
446 for ( int irrep = 0; irrep < n_irreps; irrep++ ){
447
448 int NTOTAL = idx->getNORB( irrep );
449 int NROTATE = (( space == 'O' ) ? idx->getNOCC( irrep ) : (( space == 'A' ) ? idx->getNDMRG( irrep ) : idx->getNVIRT( irrep )));
450 int NSHIFT = (( space == 'O' ) ? 0 : (( space == 'A' ) ? idx->getNOCC( irrep ) : NTOTAL - NROTATE ));
451 int NJUMP = idx->getDMRGcumulative( irrep );
452 if ( NROTATE > 1 ){
453
454 // Diagonalize the relevant block of Mat
455 {
456 for ( int row = 0; row < NROTATE; row++ ){
457 for ( int col = 0; col < NROTATE; col++ ){
458 work1[ row + NROTATE * col ] = Mat->get( irrep, NSHIFT + row, NSHIFT + col );
459 }
460 }
461 char jobz = 'V';
462 char uplo = 'U';
463 int info;
464 int size = max( 3 * NROTATE - 1, NROTATE * NROTATE );
465 dsyev_( &jobz, &uplo, &NROTATE, work1, &NROTATE, work2 + size, work2, &size, &info );
466 }
467
468 // Invert the order (large to small)
469 if ( invert ){
470 for ( int col = 0; col < NROTATE / 2; col++ ){
471 for ( int row = 0; row < NROTATE; row++ ){
472 const double temp = work1[ row + NROTATE * ( NROTATE - 1 - col ) ];
473 work1[ row + NROTATE * ( NROTATE - 1 - col ) ] = work1[ row + NROTATE * col ];
474 work1[ row + NROTATE * col ] = temp;
475 }
476 }
477 }
478
479 // Adjust the u-matrix accordingly
480 double * umatrix = Umat->getBlock( irrep ) + NSHIFT;
481 for ( int row = 0; row < NROTATE; row++ ){
482 for ( int col = 0; col < NTOTAL; col++ ){
483 work2[ row + NROTATE * col ] = umatrix[ row + NTOTAL * col ];
484 }
485 }
486 char trans = 'T';
487 char notrans = 'N';
488 double one = 1.0;
489 double set = 0.0;
490 dgemm_( &trans, ¬rans, &NROTATE, &NTOTAL, &NROTATE, &one, work1, &NROTATE, work2, &NROTATE, &set, umatrix, &NTOTAL );
491
492 // Adjust the two_dm, three_dm, and contract objects accordingly
493 if ( space == 'A' ){
494 if ( two_dm != NULL ){ rotate_active_space_object( 4, two_dm, work2, work1, tot_dmrg, NJUMP, NROTATE ); }
495 if ( three_dm != NULL ){ rotate_active_space_object( 6, three_dm, work2, work1, tot_dmrg, NJUMP, NROTATE ); }
496 if ( contract != NULL ){ rotate_active_space_object( 6, contract, work2, work1, tot_dmrg, NJUMP, NROTATE ); }
497 }
498 }
499 }
500 }
501
502 #ifdef CHEMPS2_MPI_COMPILATION
503 Umat->broadcast( MPI_CHEMPS2_MASTER );
504 if ( space == 'A' ){
505 if ( two_dm != NULL ){ MPIchemps2::broadcast_array_double( two_dm, tot_dmrg * tot_dmrg * tot_dmrg * tot_dmrg, MPI_CHEMPS2_MASTER ); }
506 if ( three_dm != NULL ){ MPIchemps2::broadcast_array_double( three_dm, tot_dmrg * tot_dmrg * tot_dmrg * tot_dmrg * tot_dmrg * tot_dmrg, MPI_CHEMPS2_MASTER ); }
507 if ( contract != NULL ){ MPIchemps2::broadcast_array_double( contract, tot_dmrg * tot_dmrg * tot_dmrg * tot_dmrg * tot_dmrg * tot_dmrg, MPI_CHEMPS2_MASTER ); }
508 }
509 #endif
510
511 }
512
rotate_active_space_object(const int num_indices,double * object,double * work,double * rotation,const int LAS,const int NJUMP,const int NROTATE)513 void CheMPS2::CASSCF::rotate_active_space_object( const int num_indices, double * object, double * work, double * rotation, const int LAS, const int NJUMP, const int NROTATE ){
514
515 assert( num_indices >= 2 );
516 assert( num_indices <= 6 );
517
518 int power[] = { 1,
519 LAS,
520 LAS * LAS,
521 LAS * LAS * LAS,
522 LAS * LAS * LAS * LAS,
523 LAS * LAS * LAS * LAS * LAS,
524 LAS * LAS * LAS * LAS * LAS * LAS };
525
526 for ( int rot_index = num_indices - 1; rot_index >= 0; rot_index-- ){
527 for ( int block = 0; block < power[ num_indices - 1 - rot_index ]; block++ ){
528 double * mat = object + power[ rot_index ] * NJUMP + power[ rot_index + 1 ] * block;
529 int ROTDIM = NROTATE;
530 char notrans = 'N';
531 double one = 1.0;
532 double set = 0.0;
533 dgemm_( ¬rans, ¬rans, power + rot_index, &ROTDIM, &ROTDIM, &one, mat, power + rot_index, rotation, &ROTDIM, &set, work, power + rot_index );
534 int size = power[ rot_index ] * NROTATE;
535 int inc1 = 1;
536 dcopy_( &size, work, &inc1, mat, &inc1 );
537 }
538 }
539
540 }
541
542
543