C++ Problem::SetupReorderD2h方法代码示例

double CheMPS2::CASSCF::caspt2( const int Nelectrons, const int TwoS, const int Irrep, ConvergenceScheme * OptScheme, const int rootNum, DMRGSCFoptions * scf_options, const double IPEA, const double IMAG, const bool PSEUDOCANONICAL, const bool CHECKPOINT, const bool CUMULANT ){

   #ifdef CHEMPS2_MPI_COMPILATION
      const bool am_i_master = ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER );
   #else
      const bool am_i_master = true;
   #endif

   const int num_elec = Nelectrons - 2 * iHandler->getNOCCsum();
   assert( num_elec >= 0 );

   if ( CASPT2::vector_length( iHandler ) == 0 ){
      if ( am_i_master ){
         cout << "CheMPS2::CASSCF::caspt2 : There are no CASPT2 excitations between the CORE, ACTIVE, and VIRTUAL orbital spaces." << endl;
      }
      return 0.0;
   }

   //Determine the maximum NORB(irrep) and the max_block_size for the ERI orbital rotation
   const int maxlinsize      = iHandler->getNORBmax();
   const long long fullsize  = ((long long) maxlinsize ) * ((long long) maxlinsize ) * ((long long) maxlinsize ) * ((long long) maxlinsize );
   const string tmp_filename = tmp_folder + "/" + CheMPS2::DMRGSCF_eri_storage_name;
   const int dmrgsize_power4 = nOrbDMRG * nOrbDMRG * nOrbDMRG * nOrbDMRG;
   //For (ERI rotation, update unitary, block diagonalize, orbital localization)
   DMRGSCFintegrals * theRotatedTEI = new DMRGSCFintegrals( iHandler );
   const int temp_work_size = (( fullsize > CheMPS2::DMRGSCF_max_mem_eri_tfo ) ? CheMPS2::DMRGSCF_max_mem_eri_tfo : fullsize );
   const int work_mem_size  = max( max( temp_work_size , maxlinsize * maxlinsize * 4 ) , dmrgsize_power4 );
   const int tot_dmrg_power6 = dmrgsize_power4 * nOrbDMRG * nOrbDMRG;
   double * mem1 = new double[ work_mem_size ];
   double * mem2 = new double[ ( PSEUDOCANONICAL ) ? work_mem_size : max( work_mem_size, tot_dmrg_power6 ) ];

   // Rotate to pseudocanonical orbitals
   if ( PSEUDOCANONICAL ){
      assert( successful_solve ); // DMRG1RDM needs to be set by CASSCF::solve for buildQmatACT()
      buildTmatrix();
      buildQmatOCC();
      buildQmatACT();
      construct_fock( theFmatrix, theTmatrix, theQmatOCC, theQmatACT, iHandler );
      block_diagonalize( 'O', theFmatrix, unitary, mem1, mem2, iHandler, false, NULL, NULL, NULL );
      block_diagonalize( 'A', theFmatrix, unitary, mem1, mem2, iHandler, false, NULL, NULL, NULL );
      block_diagonalize( 'V', theFmatrix, unitary, mem1, mem2, iHandler, false, NULL, NULL, NULL );
   } else {
      if ( successful_solve == false ){
         assert( scf_options->getStoreUnitary() );
         if ( am_i_master ){
            struct stat file_info;
            const int file_stat = stat( (scf_options->getUnitaryStorageName()).c_str(), &file_info );
            assert( file_stat == 0 );
            unitary->loadU( scf_options->getUnitaryStorageName() );
         }
         #ifdef CHEMPS2_MPI_COMPILATION
         unitary->broadcast( MPI_CHEMPS2_MASTER );
         #endif
      }
   }

   // Fill active space Hamiltonian
   Hamiltonian * HamAS = new Hamiltonian( nOrbDMRG, SymmInfo.getGroupNumber(), iHandler->getIrrepOfEachDMRGorbital() );
   Problem * Prob = new Problem( HamAS, TwoS, num_elec, Irrep );
   Prob->SetupReorderD2h(); // Doesn't matter if the group isn't D2h, Prob checks it.
   buildTmatrix();
   buildQmatOCC();
   fillConstAndTmatDMRG( HamAS );
   if ( am_i_master ){
      DMRGSCFrotations::rotate( VMAT_ORIG, HamAS->getVmat(), NULL, 'A', 'A', 'A', 'A', iHandler, unitary, mem1, mem2, work_mem_size, tmp_filename );
   }
   #ifdef CHEMPS2_MPI_COMPILATION
   HamAS->getVmat()->broadcast( MPI_CHEMPS2_MASTER );
   #endif

   double E_CASSCF = 0.0;
   double * three_dm = new double[ tot_dmrg_power6 ];
   double * contract = new double[ tot_dmrg_power6 ];
   for ( int cnt = 0; cnt < tot_dmrg_power6; cnt++ ){ contract[ cnt ] = 0.0; }

   int next_hamorb1 = 0;
   int next_hamorb2 = 0;
   const bool make_checkpt = (( CUMULANT == false ) && ( CHECKPOINT ));
   bool checkpt_loaded = false;
   if ( make_checkpt ){
      assert(( OptScheme != NULL ) || ( rootNum > 1 ));
      checkpt_loaded = read_f4rdm_checkpoint( CheMPS2::DMRGSCF_f4rdm_name, &next_hamorb1, &next_hamorb2, tot_dmrg_power6, contract );
   }

   // Solve the active space problem
   if (( OptScheme == NULL ) && ( rootNum == 1 )){ // Do FCI

      if ( am_i_master ){
         const int nalpha = ( num_elec + TwoS ) / 2;
         const int nbeta  = ( num_elec - TwoS ) / 2;
         const double workmem = 1000.0; // 1GB
         const int verbose = 2;
         CheMPS2::FCI * theFCI = new CheMPS2::FCI( HamAS, nalpha, nbeta, Irrep, workmem, verbose );
         double * inoutput = new double[ theFCI->getVecLength(0) ];
         theFCI->ClearVector( theFCI->getVecLength(0), inoutput );
         inoutput[ theFCI->LowestEnergyDeterminant() ] = 1.0;
         E_CASSCF = theFCI->GSDavidson( inoutput );
         theFCI->Fill2RDM( inoutput, DMRG2DM );                     // 2-RDM
         theFCI->Fill3RDM( inoutput, three_dm );                    // 3-RDM
         setDMRG1DM( num_elec, nOrbDMRG, DMRG1DM, DMRG2DM );        // 1-RDM
//.........这里部分代码省略.........

double CheMPS2::CASSCF::solve( const int Nelectrons, const int TwoS, const int Irrep, ConvergenceScheme * OptScheme, const int rootNum, DMRGSCFoptions * scf_options ){

   const int num_elec = Nelectrons - 2 * iHandler->getNOCCsum();
   assert( num_elec >= 0 );
   assert(( OptScheme != NULL ) || (( OptScheme == NULL ) && ( rootNum == 1 )));

   // Convergence variables
   double gradNorm = 1.0;
   double updateNorm = 1.0;
   double * gradient = new double[ unitary->getNumVariablesX() ];
   for ( int cnt = 0; cnt < unitary->getNumVariablesX(); cnt++ ){ gradient[ cnt ] = 0.0; }
   double * diis_vec = NULL;
   double Energy = 1e8;

   // The CheMPS2::Problem for the inner DMRG calculation
   Hamiltonian * HamDMRG = new Hamiltonian(nOrbDMRG, SymmInfo.getGroupNumber(), iHandler->getIrrepOfEachDMRGorbital());
   Problem * Prob = new Problem(HamDMRG, TwoS, num_elec, Irrep);
   Prob->SetupReorderD2h(); //Doesn't matter if the group isn't D2h, Prob checks it.

   // Determine the maximum NORB(irrep) and the max_block_size for the ERI orbital rotation
   const int maxlinsize      = iHandler->getNORBmax();
   const long long fullsize  = ((long long) maxlinsize ) * ((long long) maxlinsize ) * ((long long) maxlinsize ) * ((long long) maxlinsize );
   const string tmp_filename = CheMPS2::defaultTMPpath + "/" + CheMPS2::DMRGSCF_eri_storage_name;
   const int dmrgsize_power4 = nOrbDMRG * nOrbDMRG * nOrbDMRG * nOrbDMRG;
   //For (ERI rotation, update unitary, block diagonalize, orbital localization)
   const int temp_work_size = (( fullsize > CheMPS2::DMRGSCF_max_mem_eri_tfo ) ? CheMPS2::DMRGSCF_max_mem_eri_tfo : fullsize );
   const int work_mem_size = max( max( temp_work_size , maxlinsize * maxlinsize * 4 ) , dmrgsize_power4 );
   double * mem1 = new double[ work_mem_size ];
   double * mem2 = new double[ work_mem_size ];

   //The two-body rotator and Edmiston-Ruedenberg active space localizer
   EdmistonRuedenberg * theLocalizer = NULL;
   if ( scf_options->getWhichActiveSpace() == 2 ){ theLocalizer = new EdmistonRuedenberg( HamDMRG->getVmat(), iHandler->getGroupNumber() ); }

   //Load unitary from disk
   if ( scf_options->getStoreUnitary() ){
      struct stat file_info;
      int master_stat = stat( (scf_options->getUnitaryStorageName()).c_str(), &file_info );
      if ( master_stat == 0 ){ unitary->loadU( scf_options->getUnitaryStorageName() ); }
   }

   //Load DIIS from disk
   DIIS * diis = NULL;
   if (( scf_options->getDoDIIS() ) && ( scf_options->getStoreDIIS() )){
      struct stat file_info;
      int master_stat = stat( (scf_options->getDIISStorageName()).c_str(), &file_info );
      if ( master_stat == 0 ){
         const int diis_vec_size = iHandler->getROTparamsize();
         diis = new DIIS( diis_vec_size, unitary->getNumVariablesX(), scf_options->getNumDIISVecs() );
         diis->loadDIIS( scf_options->getDIISStorageName() );
         diis_vec = new double[ diis_vec_size ];
      }
   }

   int nIterations = 0;

   /*******************************
   ***   Actual DMRGSCF loops   ***
   *******************************/
   while (( gradNorm > scf_options->getGradientThreshold() ) && ( nIterations < scf_options->getMaxIterations() )){

      nIterations++;

      //Update the unitary transformation
      if ( unitary->getNumVariablesX() > 0 ){

         unitary->updateUnitary( mem1, mem2, gradient, true, true ); //multiply = compact = true

         if (( scf_options->getDoDIIS() ) && ( updateNorm <= scf_options->getDIISGradientBranch() )){
            if ( scf_options->getWhichActiveSpace() == 1 ){
               cout << "DMRGSCF::solve : DIIS has started. Active space not rotated to NOs anymore!" << endl;
            }
            if ( scf_options->getWhichActiveSpace() == 2 ){
               cout << "DMRGSCF::solve : DIIS has started. Active space not rotated to localized orbitals anymore!" << endl;
            }
            if ( diis == NULL ){
               const int diis_vec_size = iHandler->getROTparamsize();
               diis = new DIIS( diis_vec_size, unitary->getNumVariablesX(), scf_options->getNumDIISVecs() );
               diis_vec = new double[ diis_vec_size ];
               unitary->makeSureAllBlocksDetOne( mem1, mem2 );
            }
            unitary->getLog( diis_vec, mem1, mem2 );
            diis->appendNew( gradient, diis_vec );
            diis->calculateParam( diis_vec );
            unitary->updateUnitary( mem1, mem2, diis_vec, false, false ); //multiply = compact = false
         }
      }
      if (( scf_options->getStoreUnitary() ) && ( gradNorm != 1.0 )){ unitary->saveU( scf_options->getUnitaryStorageName() ); }
      if (( scf_options->getStoreDIIS() ) && ( updateNorm != 1.0 ) && ( diis != NULL )){ diis->saveDIIS( scf_options->getDIISStorageName() ); }
      int master_diis = (( diis != NULL ) ? 1 : 0 );

      //Fill HamDMRG
      buildQmatOCC();
      buildTmatrix();
      fillConstAndTmatDMRG( HamDMRG );
      DMRGSCFVmatRotations::rotate( VMAT_ORIG, HamDMRG->getVmat(), NULL, 'A', 'A', 'A', 'A', iHandler, unitary, mem1, mem2, work_mem_size, tmp_filename );

      //Localize the active space and reorder the orbitals within each irrep based on the exchange matrix
      if (( scf_options->getWhichActiveSpace() == 2 ) && ( master_diis == 0 )){ //When the DIIS has started: stop
         theLocalizer->Optimize(mem1, mem2, scf_options->getStartLocRandom()); //Default EDMISTONRUED_gradThreshold and EDMISTONRUED_maxIter used
//.........这里部分代码省略.........