diff --git a/BaseOperator.h b/BaseOperator.h index d9370dc9..4162d921 100644 --- a/BaseOperator.h +++ b/BaseOperator.h @@ -11,7 +11,12 @@ Sandeep Sharma and Garnet K.-L. Chan #include #include #include +#if __has_include("boost/bind/bind.hpp") +#include +using namespace boost::placeholders; +#else #include +#endif #include #include "SpinQuantum.h" #include "ObjectMatrix.h" diff --git a/btas/include/btas/btas.h b/btas/include/btas/btas.h index b973de81..69229666 100644 --- a/btas/include/btas/btas.h +++ b/btas/include/btas/btas.h @@ -8,7 +8,12 @@ // C++11 involves these classes, however, they're not supported by boost serialization. // They should be replaced by STL library in the future. #include +#if __has_include("boost/bind/bind.hpp") +#include +using namespace boost::placeholders; +#else #include +#endif #include #include diff --git a/genetic/GAOptimize.C b/genetic/GAOptimize.C index 074530b3..9091650c 100644 --- a/genetic/GAOptimize.C +++ b/genetic/GAOptimize.C @@ -9,7 +9,12 @@ using namespace std; #include +#if __has_include("boost/bind/bind.hpp") +#include +using namespace boost::placeholders; +#else #include +#endif #ifndef SERIAL #include diff --git a/input.C b/input.C index df25de3f..055e093e 100644 --- a/input.C +++ b/input.C @@ -187,6 +187,9 @@ void SpinAdapted::Input::initialize_defaults() m_orbformat=MOLPROFORM; m_warmup = LOCAL0; + + m_orz3rdmalt1 = -1; + m_orz3rdmalt2 = -1; } void SpinAdapted::Input::usedkey_error(string& key, string& line) { @@ -1197,6 +1200,12 @@ If 2 spins are given, the calculations of transition density matrix between wave m_reset_iterations = true; } + else if(boost::iequals(keyword, "orz3rdmalt") ) + { + m_orz3rdmalt1 = atoi(tok[1].c_str()); + m_orz3rdmalt2 = atoi(tok[2].c_str()); + } + else { pout << "Unrecognized option :: " << keyword << endl; @@ -1252,6 +1261,8 @@ If 2 spins are given, the calculations of transition density matrix between wave mpi::broadcast(world, m_load_prefix, 0); mpi::broadcast(world, m_save_prefix, 0); mpi::broadcast(world, m_calc_type, 0); + mpi::broadcast(world, m_orz3rdmalt1, 0); + mpi::broadcast(world, m_orz3rdmalt2, 0); #endif //make the scratch files diff --git a/input.h b/input.h index 5d3735d8..e5c97cb0 100644 --- a/input.h +++ b/input.h @@ -18,7 +18,11 @@ Sandeep Sharma and Garnet K.-L. Chan #include "SpinQuantum.h" #include "timer.h" #include "couplingCoeffs.h" +#if __has_include("boost/tr1/unordered_map.hpp") #include +#else +#include +#endif #include "IntegralMatrix.h" @@ -205,6 +209,9 @@ class Input { int m_kept_nevpt2_states; pair NevPrint; + int m_orz3rdmalt1; + int m_orz3rdmalt2; + friend class boost::serialization::access; template void serialize(Archive & ar, const unsigned int version) @@ -234,6 +241,7 @@ class Input { ar & n_twodot_noise & m_twodot_noise & m_twodot_gamma & m_guessState; ar & m_calc_ri_4pdm & m_store_ripdm_readable & m_nevpt2 & m_conventional_nevpt2 & m_kept_nevpt2_states & NevPrint; ar & m_act_size & m_core_size & m_virt_size & m_total_orbs & m_total_spin_orbs_symmetry & m_total_spatial_to_spin & m_total_spin_to_spatial; + ar & m_orz3rdmalt1 & m_orz3rdmalt2; } @@ -572,6 +580,8 @@ class Input { const bool &npdm_intermediate() const { return m_npdm_intermediate; } bool &npdm_multinode() { return m_npdm_multinode; } const bool &npdm_multinode() const { return m_npdm_multinode; } + int orz3rdmalt1() const { return m_orz3rdmalt1; } + int orz3rdmalt2() const { return m_orz3rdmalt2; } }; } #endif diff --git a/modules/ds0_onepdm/ds0_sweep.C b/modules/ds0_onepdm/ds0_sweep.C index 985c33a6..9aa77fd3 100644 --- a/modules/ds0_onepdm/ds0_sweep.C +++ b/modules/ds0_onepdm/ds0_sweep.C @@ -193,6 +193,8 @@ double do_one(SweepParams &sweepParams, const bool &warmUp, const bool &forward, save_onepdm_spatial_text(onepdm, state, stateB); save_onepdm_text(onepdm, state, stateB); save_onepdm_spatial_binary(onepdm, state, stateB); + + return sweepParams.get_lowest_energy()[0]; } } } diff --git a/modules/ds1_onepdm/ds1_sweep.C b/modules/ds1_onepdm/ds1_sweep.C index 16750e63..d656277b 100644 --- a/modules/ds1_onepdm/ds1_sweep.C +++ b/modules/ds1_onepdm/ds1_sweep.C @@ -192,6 +192,8 @@ double do_one(SweepParams &sweepParams, const bool &warmUp, const bool &forward, save_onepdm_spatial_text(onepdm, state, stateB); save_onepdm_text(onepdm, state, stateB); save_onepdm_spatial_binary(onepdm, state, stateB); + + return sweepParams.get_lowest_energy()[0]; } } } diff --git a/modules/npdm/fourpdm_container.C b/modules/npdm/fourpdm_container.C index 5955e9d2..c8517ac5 100644 --- a/modules/npdm/fourpdm_container.C +++ b/modules/npdm/fourpdm_container.C @@ -203,7 +203,7 @@ void Fourpdm_container::external_sort_index(const int &i, const int &j) } long outputbuffsize=sorting_buff*4; long outputbuff_position = 0; - Sortpdm::batch_index outputbuff[outputbuffsize]; + Sortpdm::batch_index *outputbuff = new Sortpdm::batch_index[outputbuffsize]; for(;;){ // select the smallest one in the current positions of different caches. int smallest = 0; @@ -230,6 +230,7 @@ void Fourpdm_container::external_sort_index(const int &i, const int &j) } fclose(outputfile); //Finish external sort of index. + delete [] outputbuff; //Clean up. for(int p=0; p< world.size();p++){ diff --git a/modules/npdm/npdm.C b/modules/npdm/npdm.C index 1f04758f..9eec765e 100644 --- a/modules/npdm/npdm.C +++ b/modules/npdm/npdm.C @@ -1,8 +1,8 @@ -/* -Developed by Sandeep Sharma and Garnet K.-L. Chan, 2012 -Copyright (c) 2012, Garnet K.-L. Chan - -This program is integrated in Molpro with the permission of +/* +Developed by Sandeep Sharma and Garnet K.-L. Chan, 2012 +Copyright (c) 2012, Garnet K.-L. Chan + +This program is integrated in Molpro with the permission of Sandeep Sharma and Garnet K.-L. Chan */ @@ -29,7 +29,7 @@ namespace Npdm { //----------------------------------------------------------------------------------------------------------------------------------------------------------- -void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, +void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, const bool &useSlater, const bool& dot_with_sys, const int state, const int stateB) { Timer timer; @@ -51,7 +51,7 @@ void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepP systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ; systemDotEnd = systemDotStart - systemDotSize; } - vector spindotsites(2); + vector spindotsites(2); spindotsites[0] = systemDotStart; spindotsites[1] = systemDotEnd; //if (useSlater) { @@ -70,7 +70,7 @@ void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepP if(dmrginp.setStateSpecific() || dmrginp.transition_diff_irrep()){ InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, state, stateB, sweepParams.get_sys_add(), dmrginp.direct(), system.get_integralIndex(), DISTRIBUTED_STORAGE, true, true); - + InitBlocks::InitNewEnvironmentBlock(environment, systemDot, newEnvironment, system, systemDot, state, stateB, sweepParams.get_sys_add(), sweepParams.get_env_add(), forward, dmrginp.direct(), @@ -79,7 +79,7 @@ void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepP else{ InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, sweepParams.current_root(), sweepParams.current_root(), sweepParams.get_sys_add(), dmrginp.direct(), system.get_integralIndex(), DISTRIBUTED_STORAGE, true, true); - + InitBlocks::InitNewEnvironmentBlock(environment, systemDot, newEnvironment, system, systemDot, sweepParams.current_root(), sweepParams.current_root(), sweepParams.get_sys_add(), sweepParams.get_env_add(), forward, dmrginp.direct(), @@ -90,21 +90,21 @@ void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepP newSystem.set_loopblock(true); system.set_loopblock(false); newEnvironment.set_loopblock(false); - InitBlocks::InitBigBlock(newSystem, newEnvironment, big); + InitBlocks::InitBigBlock(newSystem, newEnvironment, big); const int nroots = dmrginp.nroots(); std::vector solution; if(state==stateB){ solution.resize(1); DiagonalMatrix e; - GuessWave::guess_wavefunctions(solution[0], e, big, sweepParams.get_guesstype(), true, state, true, 0.0); + GuessWave::guess_wavefunctions(solution[0], e, big, sweepParams.get_guesstype(), true, state, true, 0.0); } else{ solution.resize(2); DiagonalMatrix e; - GuessWave::guess_wavefunctions(solution[0], e, big, sweepParams.get_guesstype(), true, state, true, 0.0,false); - GuessWave::guess_wavefunctions(solution[1], e, big, sweepParams.get_guesstype(), true, stateB, true, 0.0,true); + GuessWave::guess_wavefunctions(solution[0], e, big, sweepParams.get_guesstype(), true, state, true, 0.0,false); + GuessWave::guess_wavefunctions(solution[1], e, big, sweepParams.get_guesstype(), true, stateB, true, 0.0,true); } #ifndef SERIAL @@ -113,8 +113,8 @@ void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepP #endif - //GuessWave::guess_wavefunctions(solution[0], e, big, sweepParams.get_guesstype(), true, state, false, 0.0); - //GuessWave::guess_wavefunctions(solution[1], e, big, sweepParams.get_guesstype(), true, stateB, true, 0.0); + //GuessWave::guess_wavefunctions(solution[0], e, big, sweepParams.get_guesstype(), true, state, false, 0.0); + //GuessWave::guess_wavefunctions(solution[1], e, big, sweepParams.get_guesstype(), true, stateB, true, 0.0); @@ -152,7 +152,7 @@ void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepP } - + int sweepPos = sweepParams.get_block_iter(); @@ -168,13 +168,13 @@ void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepP //FIXME - //Maybe, for StateSpecific calculations, we can load rotation matrices, wavefuntions from the disk. + //Maybe, for StateSpecific calculations, we can load rotation matrices, wavefuntions from the disk. //There is no longer need to transform wavefuntions and to make rotation matrices from the density matrices. - + //FIXME - //If in the state-average pdm, different states do not share the same rotation matrices as they do in energy calculations. Making rotation matrices from - //density matrices of different states is neccessary. - + //If in the state-average pdm, different states do not share the same rotation matrices as they do in energy calculations. Making rotation matrices from + //density matrices of different states is neccessary. + //if(newSystem.get_sites().size()>1) //if (!mpigetrank()){ //LoadRotationMatrix (newSystem.get_sites(), rotateMatrix, state); @@ -187,7 +187,7 @@ void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepP #endif // Do we need to do this step for NPDM on the last sweep site? (It's not negligible cost...?) - //It crashes at the last sweep site. + //It crashes at the last sweep site. //Since it is useless, Just omit it at the last sweep site. // if( sweepParams.get_block_iter() != sweepParams.get_n_iters() - 1) { @@ -205,7 +205,7 @@ void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepP //----------------------------------------------------------------------------------------------------------------------------------------------------------- -double npdm_do_one_sweep(Npdm_driver_base &npdm_driver, SweepParams &sweepParams, const bool &warmUp, const bool &forward, +double npdm_do_one_sweep(Npdm_driver_base &npdm_driver, SweepParams &sweepParams, const bool &warmUp, const bool &forward, const bool &restart, const int &restartSize, const int state, const int stateB) { Timer sweeptimer; @@ -220,11 +220,11 @@ double npdm_do_one_sweep(Npdm_driver_base &npdm_driver, SweepParams &sweepParams // a new renormalisation sweep routine pout << ((forward) ? "\t\t\t Starting renormalisation sweep in forwards direction" : "\t\t\t Starting renormalisation sweep in backwards direction") << endl; pout << "\t\t\t ============================================================================ " << endl; - + int integralIndex = 0; if(dmrginp.setStateSpecific() || dmrginp.transition_diff_irrep()) InitBlocks::InitStartingBlock( system, forward, state, stateB, sweepParams.get_forward_starting_size(), sweepParams.get_backward_starting_size(), restartSize, restart, warmUp, integralIndex); - else + else InitBlocks::InitStartingBlock( system, forward, sweepParams.current_root(), sweepParams.current_root(), sweepParams.get_forward_starting_size(), sweepParams.get_backward_starting_size(), restartSize, restart, warmUp, integralIndex); pout << "\t\t\t Starting block is :: " << endl << system << endl; @@ -253,26 +253,26 @@ double npdm_do_one_sweep(Npdm_driver_base &npdm_driver, SweepParams &sweepParams if (dmrginp.no_transform()) sweepParams.set_guesstype() = BASIC; - else if (!warmUp && sweepParams.get_block_iter() != 0) - sweepParams.set_guesstype() = TRANSFORM; - else if (!warmUp && sweepParams.get_block_iter() == 0 && + else if (!warmUp && sweepParams.get_block_iter() != 0) + sweepParams.set_guesstype() = TRANSFORM; + else if (!warmUp && sweepParams.get_block_iter() == 0 && ((dmrginp.algorithm_method() == TWODOT_TO_ONEDOT && dmrginp.twodot_to_onedot_iter() != sweepParams.get_sweep_iter()) || dmrginp.algorithm_method() != TWODOT_TO_ONEDOT)) sweepParams.set_guesstype() = TRANSPOSE; else sweepParams.set_guesstype() = BASIC; - + //pout << "guess wave funtion type: " << sweepParams.get_guesstype()< = <\Phi_l|a^+_ja_i|\Phi_k>* // Therefore, only calculate the situations with k >= l. // for(int stateB=0; stateB<= dmrginp.nroots(); stateB++){ +const bool docalc = (dmrginp.orz3rdmalt1()==-1 && dmrginp.orz3rdmalt2()==-1); for (int state=0; state> NPDM (state,stateB)=(" << state << "," << stateB << ")" << endl; sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; Timer timerX; npdm_driver->clear(); npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0, state,stateB); ecpu = timerX.elapsedcputime();ewall=timerX.elapsedwalltime(); p3out << "\t\t\t NPDM sweep time " << ewall << " " << ecpu << endl; - } +} + } } } else { +const bool docalc = (dmrginp.orz3rdmalt1()==-1 && dmrginp.orz3rdmalt2()==-1); for (int state=0; state> NPDM (state)=(" << state << "," << ")" << endl; sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; if ( dmrginp.new_npdm_code() ) { Timer timerX; @@ -502,13 +509,14 @@ else { npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0, state,state); ecpu = timerX.elapsedcputime();ewall=timerX.elapsedwalltime(); p3out << "\t\t\t NPDM sweep time " << ewall << " " << ecpu << endl; - } + } else{ - if (npdm_order == NPDM_ONEPDM) SweepOnepdm::do_one(sweepParams, false, direction, false, 0, state); + if (npdm_order == NPDM_ONEPDM) SweepOnepdm::do_one(sweepParams, false, direction, false, 0, state); else if (npdm_order == NPDM_TWOPDM) SweepTwopdm::do_one(sweepParams, false, direction, false, 0, state); else abort(); } } +} } } @@ -518,10 +526,10 @@ else { if(transitionpdm){ for (int state=0; stateclear(); sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; @@ -570,16 +578,16 @@ else { Sweep::CanonicalizeWavefunction(sweepParams, direction, state); } // Prepare NPDM operators - dmrginp.npdm_generate() = true; + dmrginp.npdm_generate() = true; SweepGenblock::do_one(sweepParams, false, !direction, false, 0, state, state); //this will generate the cd operators - dmrginp.npdm_generate() = false; + dmrginp.npdm_generate() = false; dmrginp.set_fullrestart() = false; // Do NPDM sweep npdm_driver->clear(); npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0, state,state); } } - + } sweep_copy.savestate(direction_copy, restartsize_copy); @@ -588,4 +596,3 @@ else { } } } - diff --git a/modules/npdm/npdm_driver.C b/modules/npdm/npdm_driver.C index f4395b46..439c4cca 100644 --- a/modules/npdm/npdm_driver.C +++ b/modules/npdm/npdm_driver.C @@ -161,12 +161,12 @@ void Npdm_driver::do_parallel_intermediate_loop( const char inner, Npdm::Npdm_ex if( inner =='r') { - npdm_expectations.compute_intermediate(outerOps,dotOps,local_waves); + npdm_expectations.compute_intermediate(&outerOps,&dotOps,local_waves); } else if ( inner =='l') { - npdm_expectations.compute_intermediate(outerOps,local_waves); + npdm_expectations.compute_intermediate(&outerOps,local_waves); } else assert(false); @@ -276,7 +276,7 @@ void Npdm_driver::get_inner_Operators( const char inner, Npdm_expectations& npdm else{ boost::shared_ptr, Wavefunction> > half_waves( new std::map, Wavefunction>); - npdm_expectations.compute_intermediate(*inner_Operators[i],*dotOps,*half_waves); + npdm_expectations.compute_intermediate(&(*inner_Operators[i]),&(*dotOps),*half_waves); inner_intermediate.push_back(half_waves); } } @@ -299,7 +299,7 @@ void Npdm_driver::get_inner_Operators( const char inner, Npdm_expectations& npdm else{ boost::shared_ptr, Wavefunction> > half_waves( new std::map, Wavefunction>); - npdm_expectations.compute_intermediate(*inner_Operators[i], *half_waves); + npdm_expectations.compute_intermediate(&(*inner_Operators[i]), *half_waves); inner_intermediate.push_back(half_waves); } } diff --git a/modules/npdm/npdm_expectations.C b/modules/npdm/npdm_expectations.C index bd1fd1e4..3e53aaba 100644 --- a/modules/npdm/npdm_expectations.C +++ b/modules/npdm/npdm_expectations.C @@ -129,7 +129,7 @@ double Npdm_expectations::contract_spin_adapted_operators( int ilhs, int idot, i if ( dotOps.transpose_ ) dotOp = boost::shared_ptr( &dotOpTr, boostutils::null_deleter() ); Transposeview rhsOpTr = Transposeview(*rhsOp); if ( rhsOps.transpose_ ) rhsOp = boost::shared_ptr( &rhsOpTr, boostutils::null_deleter() ); - expectation = spinExpectation(wavefunction_0, wavefunction_1, *lhsOp, *dotOp, *rhsOp, big_); + expectation = spinExpectation(wavefunction_0, wavefunction_1, *lhsOp, *dotOp, &(*rhsOp), big_); } // X_X_0 else if ( (lhsOps.opReps_.size() > 0) && (dotOps.opReps_.size() > 0) && (rhsOps.opReps_.size() == 0) ) { @@ -137,7 +137,7 @@ double Npdm_expectations::contract_spin_adapted_operators( int ilhs, int idot, i if ( lhsOps.transpose_ ) lhsOp = boost::shared_ptr( &lhsOpTr, boostutils::null_deleter() ); Transposeview dotOpTr = Transposeview(*dotOp); if ( dotOps.transpose_ ) dotOp = boost::shared_ptr( &dotOpTr, boostutils::null_deleter() ); - expectation = spinExpectation(wavefunction_0, wavefunction_1, *lhsOp, *dotOp, *null, big_); + expectation = spinExpectation(wavefunction_0, wavefunction_1, *lhsOp, *dotOp, null, big_); } // X_0_X else if ( (lhsOps.opReps_.size() > 0) && (dotOps.opReps_.size() == 0) && (rhsOps.opReps_.size() > 0) ) { @@ -145,7 +145,7 @@ double Npdm_expectations::contract_spin_adapted_operators( int ilhs, int idot, i if ( lhsOps.transpose_ ) lhsOp = boost::shared_ptr( &lhsOpTr, boostutils::null_deleter() ); Transposeview rhsOpTr = Transposeview(*rhsOp); if ( rhsOps.transpose_ ) rhsOp = boost::shared_ptr( &rhsOpTr, boostutils::null_deleter() ); - expectation = spinExpectation(wavefunction_0, wavefunction_1, *lhsOp, *null, *rhsOp, big_); + expectation = spinExpectation(wavefunction_0, wavefunction_1, *lhsOp, *null, &(*rhsOp), big_); } // 0_X_X else if ( (lhsOps.opReps_.size() == 0) && (dotOps.opReps_.size() > 0) && (rhsOps.opReps_.size() > 0) ) { @@ -153,25 +153,25 @@ double Npdm_expectations::contract_spin_adapted_operators( int ilhs, int idot, i if ( dotOps.transpose_ ) dotOp = boost::shared_ptr( &dotOpTr, boostutils::null_deleter() ); Transposeview rhsOpTr = Transposeview(*rhsOp); if ( rhsOps.transpose_ ) rhsOp = boost::shared_ptr( &rhsOpTr, boostutils::null_deleter() ); - expectation = spinExpectation(wavefunction_0, wavefunction_1, *null, *dotOp, *rhsOp, big_); + expectation = spinExpectation(wavefunction_0, wavefunction_1, *null, *dotOp, &(*rhsOp), big_); } // X_0_0 else if ( (lhsOps.opReps_.size() > 0) && (dotOps.opReps_.size() == 0) && (rhsOps.opReps_.size() == 0) ) { Transposeview lhsOpTr = Transposeview(*lhsOp); if ( lhsOps.transpose_ ) lhsOp = boost::shared_ptr( &lhsOpTr, boostutils::null_deleter() ); - expectation = spinExpectation(wavefunction_0, wavefunction_1, *lhsOp, *null, *null, big_); + expectation = spinExpectation(wavefunction_0, wavefunction_1, *lhsOp, *null, null, big_); } // 0_X_0 else if ( (lhsOps.opReps_.size() == 0) && (dotOps.opReps_.size() > 0) && (rhsOps.opReps_.size() == 0) ) { Transposeview dotOpTr = Transposeview(*dotOp); if ( dotOps.transpose_ ) dotOp = boost::shared_ptr( &dotOpTr, boostutils::null_deleter() ); - expectation = spinExpectation(wavefunction_0, wavefunction_1, *null, *dotOp, *null, big_); + expectation = spinExpectation(wavefunction_0, wavefunction_1, *null, *dotOp, null, big_); } // 0_0_X else if ( (lhsOps.opReps_.size() == 0) && (dotOps.opReps_.size() == 0) && (rhsOps.opReps_.size() > 0) ) { Transposeview rhsOpTr = Transposeview(*rhsOp); if ( rhsOps.transpose_ ) rhsOp = boost::shared_ptr( &rhsOpTr, boostutils::null_deleter() ); - expectation = spinExpectation(wavefunction_0, wavefunction_1, *null, *null, *rhsOp, big_); + expectation = spinExpectation(wavefunction_0, wavefunction_1, *null, *null, &(*rhsOp), big_); } else abort(); @@ -640,21 +640,23 @@ void Npdm_expectations::get_op_string( NpdmSpinOps_base & rhsOps, //----------------------------------------------------------------------------------------------------------------------------------------------------------- -void Npdm_expectations::compute_intermediate( NpdmSpinOps_base & lhsOps, NpdmSpinOps_base & dotOps, std::map, Wavefunction> & waves) +void Npdm_expectations::compute_intermediate( NpdmSpinOps_base * lhsOps_ptr, NpdmSpinOps_base * dotOps_ptr, std::map, Wavefunction> & waves) { #ifndef SERIAL boost::mpi::communicator world; #endif - + NpdmSpinOps_base & lhsOps = *lhsOps_ptr; + NpdmSpinOps_base & dotOps = *dotOps_ptr; + SparseMatrix* null = 0; vector dQ = wavefunction_0.get_deltaQuantum(); assert(dQ.size()==1 ); assert(dQ[0].totalSpin.getirrep()== 0); - int lindices= &lhsOps ? lhsOps.indices_.size(): 0; - int dindices= &dotOps ? dotOps.indices_.size(): 0; + int lindices= lhsOps_ptr ? lhsOps.indices_.size(): 0; + int dindices= dotOps_ptr ? dotOps.indices_.size(): 0; int rindices= 2*npdm_order_-lindices- dindices; - int hilhs = lhsOps.opReps_.size(); - int hidot = dotOps.opReps_.size(); + int hilhs = lhsOps_ptr ? lhsOps.opReps_.size(): 0; + int hidot = dotOps_ptr ? dotOps.opReps_.size(): 0; for (int idot = 0; idot < std::max(1,hidot); ++idot) { for (int ilhs = 0; ilhs < std::max(1,hilhs); ++ilhs) { // int ls= hilhs? lhsOps.opReps_.at(ilhs)->get_deltaQuantum(0).get_s().getirrep(): 0; @@ -695,18 +697,19 @@ void Npdm_expectations::compute_intermediate( NpdmSpinOps_base & lhsOps, NpdmSpi } //----------------------------------------------------------------------------------------------------------------------------------------------------------- -void Npdm_expectations::compute_intermediate( NpdmSpinOps_base & rhsOps, std::map, Wavefunction> & waves) +void Npdm_expectations::compute_intermediate( NpdmSpinOps_base * rhsOps_ptr, std::map, Wavefunction> & waves) { #ifndef SERIAL boost::mpi::communicator world; #endif + NpdmSpinOps_base & rhsOps = *rhsOps_ptr; SparseMatrix* null = 0; vector dQ = wavefunction_1.get_deltaQuantum(); assert(dQ.size()==1 ); assert(dQ[0].totalSpin.getirrep()== 0); - int rindices= &rhsOps ? rhsOps.indices_.size(): 0; - int hirhs = rhsOps.opReps_.size(); + int rindices= rhsOps_ptr ? rhsOps.indices_.size(): 0; + int hirhs = rhsOps_ptr ? rhsOps.opReps_.size(): 0; for (int irhs = 0; irhs < hirhs; ++irhs) { int rs= !rhsOps.transpose_? rhsOps.opReps_.at(irhs)->get_deltaQuantum(0).get_s().getirrep(): (-rhsOps.opReps_.at(irhs)->get_deltaQuantum(0).get_s()).getirrep(); if(std::abs(rs)<= 2*npdm_order_-rindices ) diff --git a/modules/npdm/npdm_expectations.h b/modules/npdm/npdm_expectations.h index a5b1dfa5..249c8c68 100644 --- a/modules/npdm/npdm_expectations.h +++ b/modules/npdm/npdm_expectations.h @@ -31,9 +31,9 @@ class Npdm_expectations { get_nonspin_adapted_expectations(NpdmSpinOps_base& lhsOps, NpdmSpinOps_base& rhsOps, NpdmSpinOps_base& dotOps,std::map, Wavefunction>& leftwaves, std::map, Wavefunction>& rightwaves); void get_op_string( NpdmSpinOps_base & rhsOps, std::string& op_string); void get_op_string( NpdmSpinOps_base & lhsOps, NpdmSpinOps_base & dotOps, std::string& op_string ); - void compute_intermediate( NpdmSpinOps_base & lhsOps, NpdmSpinOps_base & dotOps, std::map, Wavefunction> & waves); + void compute_intermediate( NpdmSpinOps_base * lhsOps_ptr, NpdmSpinOps_base * dotOps_ptr, std::map, Wavefunction> & waves); - void compute_intermediate( NpdmSpinOps_base & rhsOps, std::map, Wavefunction> & waves); + void compute_intermediate( NpdmSpinOps_base * rhsOps_ptr, std::map, Wavefunction> & waves); std::vector intermediate_filenames; private: diff --git a/modules/npdm/npdm_expectations_engine.C b/modules/npdm/npdm_expectations_engine.C index 37d03665..0b3fb9ee 100644 --- a/modules/npdm/npdm_expectations_engine.C +++ b/modules/npdm/npdm_expectations_engine.C @@ -22,8 +22,9 @@ namespace Npdm{ //----------------------------------------------------------------------------------------------------------------------------------------------------------- -double spinExpectation(Wavefunction& wave1, Wavefunction& wave2, SparseMatrix& leftOp, SparseMatrix& dotOp, SparseMatrix& rightOp, const SpinBlock& big) +double spinExpectation(Wavefunction& wave1, Wavefunction& wave2, SparseMatrix& leftOp, SparseMatrix& dotOp, SparseMatrix* rightOp_ptr, const SpinBlock& big) { + SparseMatrix& rightOp = *rightOp_ptr; //calculating // do transpose specifies if we want separately. This can be avoided in some sitations if wave1 and wave2 are the same functions @@ -39,9 +40,9 @@ double spinExpectation(Wavefunction& wave1, Wavefunction& wave2, SparseMatrix& l // In spin-adapted, getirrep is the value of S // In non spin-adapted, getirrep is the value of S_z if(dmrginp.spinAdapted()) - totalspin = (&rightOp) ? rightOp.get_spin().getirrep() : 0; + totalspin = rightOp_ptr ? rightOp.get_spin().getirrep() : 0; else - totalspin = (&rightOp) ? -(rightOp.get_spin().getirrep()) : 0; + totalspin = rightOp_ptr ? -(rightOp.get_spin().getirrep()) : 0; FormLeftOp(leftBlock, leftOp, dotOp, AOp, totalspin); SpinQuantum hq(0,SpinSpace(0),IrrepSpace(0)); diff --git a/modules/npdm/npdm_expectations_engine.h b/modules/npdm/npdm_expectations_engine.h index 8968d2ad..1402a122 100644 --- a/modules/npdm/npdm_expectations_engine.h +++ b/modules/npdm/npdm_expectations_engine.h @@ -18,7 +18,7 @@ namespace Npdm{ void FormLeftOp(const SpinBlock* leftBlock, const SparseMatrix& leftOp, const SparseMatrix& dotOp, SparseMatrix& Aop, int totalspin); double DotProduct(const Wavefunction& w1, const Wavefunction& w2, const SpinBlock& big); - double spinExpectation(Wavefunction& wave1, Wavefunction& wave2, SparseMatrix &leftOp, SparseMatrix& dotOp, SparseMatrix& rightOp, const SpinBlock& big); + double spinExpectation(Wavefunction& wave1, Wavefunction& wave2, SparseMatrix &leftOp, SparseMatrix& dotOp, SparseMatrix* rightOp_ptr, const SpinBlock& big); } } diff --git a/modules/npdm/npdm_spin_transformation.h b/modules/npdm/npdm_spin_transformation.h index 8ec9038e..108bde20 100644 --- a/modules/npdm/npdm_spin_transformation.h +++ b/modules/npdm/npdm_spin_transformation.h @@ -4,8 +4,9 @@ #include #include #include +#include #include -#ifndef BOOST_1_56_0 +#if BOOST_VERSION < 105600 #include #include #include diff --git a/modules/npdm/threepdm_container.C b/modules/npdm/threepdm_container.C index d4eb0dca..acf5a74e 100644 --- a/modules/npdm/threepdm_container.C +++ b/modules/npdm/threepdm_container.C @@ -139,7 +139,7 @@ void Threepdm_container::save_spatial_npdm_text(const int &i, const int &j) } } ofs.close(); - pout << "Spatial 3PDM trace = " << trace << "\n"; + pout << "Spatial 3PDM trace = " << trace << " (state,state)=(" << i << "," << j << ")\n"; } } diff --git a/modules/twopdm/twopdm.C b/modules/twopdm/twopdm.C index 1708b8ea..3e1e7914 100644 --- a/modules/twopdm/twopdm.C +++ b/modules/twopdm/twopdm.C @@ -152,7 +152,7 @@ void compute_two_pdm_0_4_0(Wavefunction& wave1, Wavefunction& wave2, const SpinB SparseMatrix* leftOp = 0; SparseMatrix* rightOp = 0; vector expectations; - spinExpectation(wave1, wave2, *leftOp, dotop, *rightOp, big, expectations, false); + spinExpectation(wave1, wave2, leftOp, &dotop, rightOp, big, expectations, false); int ix = 2*dotindex; assign_antisymmetric(twopdm, ix, ix+1, ix+1, ix, 0.5*expectations[0]); } @@ -178,7 +178,7 @@ void compute_two_pdm_4_0_0(Wavefunction& wave1, Wavefunction& wave2, const SpinB SparseMatrix* dotOp = 0; SparseMatrix* rightOp = 0; vector expectations; - spinExpectation(wave1, wave2, leftop, *dotOp, *rightOp, big, expectations, false); + spinExpectation(wave1, wave2, &leftop, dotOp, rightOp, big, expectations, false); int ix = 2*leftindex; assign_antisymmetric(twopdm, ix, ix+1, ix+1, ix, 0.5*expectations[0]); } @@ -205,7 +205,7 @@ void compute_two_pdm_0_0_4(Wavefunction& wave1, Wavefunction& wave2, const SpinB SparseMatrix* leftOp = 0; SparseMatrix* dotop = 0; vector expectations; - spinExpectation(wave1, wave2, *leftOp, *dotop, rightop, big, expectations, false); + spinExpectation(wave1, wave2, leftOp, dotop, &rightop, big, expectations, false); int ix = 2*rightindex; assign_antisymmetric(twopdm, ix, ix+1, ix+1, ix, 0.5*expectations[0]); } @@ -245,8 +245,8 @@ void compute_two_pdm_0_2_2(Wavefunction& wave1, Wavefunction& wave2, const SpinB boost::shared_ptr rightop0 = rightBlock->get_op_rep(CRE_CRE, sq0, kx, lx); vector expectations; Transposeview rop0 = Transposeview(*rightop0), rop2 = Transposeview(*rightop2); - spinExpectation(wave1, wave2, *leftOp, *dotop0, rop0, big, expectations, false); - spinExpectation(wave1, wave2, *leftOp, *dotop2, rop2, big, expectations, false); + spinExpectation(wave1, wave2, leftOp, &(*dotop0), &rop0, big, expectations, false); + spinExpectation(wave1, wave2, leftOp, &(*dotop2), &rop2, big, expectations, false); vector indices(4,0); indices[0] = ix; indices[1] = jx; indices[2] = lx; indices[3] = kx; expectations[0]*=-1; expectations[1]*=-1; @@ -279,8 +279,8 @@ void compute_two_pdm_0_2_2(Wavefunction& wave1, Wavefunction& wave2, const SpinB int lx = rightop2->get_orbs(1); boost::shared_ptr rightop0 = rightBlock->get_op_array(CRE_DES).get_local_element(kl)[0]->getworkingrepresentation(rightBlock);//rightBlock->get_op_rep(CRE_DES_S0, kx, lx); vector expectations; - spinExpectation(wave1, wave2, *leftOp, *dotop0, *rightop0, big, expectations, false); - spinExpectation(wave1, wave2, *leftOp, *dotop2, *rightop2, big, expectations, false); + spinExpectation(wave1, wave2, leftOp, &(*dotop0), &(*rightop0), big, expectations, false); + spinExpectation(wave1, wave2, leftOp, &(*dotop2), &(*rightop2), big, expectations, false); vector indices(4,0); indices[0] = ix; indices[1] = kx; indices[2] = jx; indices[3] = lx; @@ -322,8 +322,8 @@ void compute_two_pdm_2_0_2(Wavefunction& wave1, Wavefunction& wave2, const SpinB boost::shared_ptr rightop0 = rightBlock->get_op_array(CRE_CRE).get_local_element(kl)[0]->getworkingrepresentation(rightBlock); vector expectations; Transposeview rop0 = Transposeview(*rightop0), rop2 = Transposeview(*rightop2); - spinExpectation(wave1, wave2, *leftop0, *dotop, rop0, big, expectations, false); - spinExpectation(wave1, wave2, *leftop2, *dotop, rop2, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop0), dotop, &rop0, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop2), dotop, &rop2, big, expectations, false); vector indices(4,0); indices[0] = ix; indices[1] = jx; indices[2] = lx; indices[3] = kx; expectations[0]*=-1; expectations[1]*=-1; @@ -355,8 +355,8 @@ void compute_two_pdm_2_0_2(Wavefunction& wave1, Wavefunction& wave2, const SpinB int lx = rightop2->get_orbs(1); boost::shared_ptr rightop0 = rightBlock->get_op_array(CRE_DES).get_local_element(kl)[0]->getworkingrepresentation(rightBlock);//rightBlock->get_op_rep(CRE_DES_S0, kx, lx); vector expectations; - spinExpectation(wave1, wave2, *leftop0, *dotop, *rightop0, big, expectations, false); - spinExpectation(wave1, wave2, *leftop2, *dotop, *rightop2, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop0), dotop, &(*rightop0), big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop2), dotop, &(*rightop2), big, expectations, false); vector indices(4,0); indices[0] = ix; indices[1] = kx; indices[2] = jx; indices[3] = lx; @@ -398,8 +398,8 @@ void compute_two_pdm_2_2_0(Wavefunction& wave1, Wavefunction& wave2, const SpinB boost::shared_ptr dotop0 = dotBlock->get_op_array(CRE_CRE).get_local_element(kl)[0]->getworkingrepresentation(dotBlock);//dotBlock->get_op_rep(CRE_CRE_S0, kx, lx); vector expectations; Transposeview dop0 = Transposeview(*dotop0), dop2 = Transposeview(*dotop2); - spinExpectation(wave1, wave2, *leftop0, dop0, *rightop, big, expectations, false); - spinExpectation(wave1, wave2, *leftop2, dop2, *rightop, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop0), &dop0, rightop, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop2), &dop2, rightop, big, expectations, false); vector indices(4,0); indices[0] = ix; indices[1] = jx; indices[2] = lx; indices[3] = kx; expectations[0]*=-1; expectations[1]*=-1; @@ -432,8 +432,8 @@ void compute_two_pdm_2_2_0(Wavefunction& wave1, Wavefunction& wave2, const SpinB int lx = dotop2->get_orbs(1); boost::shared_ptr dotop0 = dotBlock->get_op_array(CRE_DES).get_local_element(kl)[0]->getworkingrepresentation(dotBlock);//dotBlock->get_op_rep(CRE_DES_S0, kx, lx); vector expectations; - spinExpectation(wave1, wave2, *leftop0, *dotop0, *rightop, big, expectations, false); - spinExpectation(wave1, wave2, *leftop2, *dotop2, *rightop, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop0), &(*dotop0), rightop, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop2), &(*dotop2), rightop, big, expectations, false); vector indices(4,0); indices[0] = ix; indices[1] = kx; indices[2] = jx; indices[3] = lx; @@ -470,8 +470,8 @@ void compute_two_pdm_1_1_2(Wavefunction& wave1, Wavefunction& wave2, const SpinB boost::shared_ptr rightop0 = rightBlock->get_op_array(CRE_CRE).get_local_element(kl)[0]->getworkingrepresentation(rightBlock);//rightBlock->get_op_rep(CRE_CRE_S0, kx, lx); vector expectations; Transposeview rop0 = Transposeview(*rightop0), rop2 = Transposeview(*rightop2); - spinExpectation(wave1, wave2, *leftop, *dotop, rop0, big, expectations, false); - spinExpectation(wave1, wave2, *leftop, *dotop, rop2, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop), &(*dotop), &rop0, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop), &(*dotop), &rop2, big, expectations, false); vector indices(4,0); indices[0] = ix; indices[1] = jx; indices[2] = kx; indices[3] = lx; expectations[0]*=-1; expectations[1]*=-1; @@ -504,8 +504,8 @@ void compute_two_pdm_1_1_2(Wavefunction& wave1, Wavefunction& wave2, const SpinB boost::shared_ptr rightop0 = rightBlock->get_op_array(CRE_DES).get_local_element(kl)[0]->getworkingrepresentation(rightBlock);//rightBlock->get_op_rep(CRE_DES_S0, kx, lx); vector expectations; Transposeview tdot = Transposeview(*dotop); - spinExpectation(wave1, wave2, *leftop, tdot, *rightop0, big, expectations, true); - spinExpectation(wave1, wave2, *leftop, tdot, *rightop2, big, expectations, true); + spinExpectation(wave1, wave2, &(*leftop), &tdot, &(*rightop0), big, expectations, true); + spinExpectation(wave1, wave2, &(*leftop), &tdot, &(*rightop2), big, expectations, true); vector indices(4,0); vector expect1(2), expect2(2); expect1[0] = expectations[0]; expect1[1] = expectations[2]; @@ -556,8 +556,8 @@ void compute_two_pdm_1_2_1(Wavefunction& wave1, Wavefunction& wave2, const SpinB Transposeview tlop = Transposeview(*leftop); vector expectations; - spinExpectation(wave1, wave2, tlop, *dotop0, trop, big, expectations, false); - spinExpectation(wave1, wave2, tlop, *dotop2, trop, big, expectations, false); + spinExpectation(wave1, wave2, &tlop, &(*dotop0), &trop, big, expectations, false); + spinExpectation(wave1, wave2, &tlop, &(*dotop2), &trop, big, expectations, false); vector indices(4,0); indices[0] = jx; indices[1] = jx; indices[2] = ix; indices[3] = kx; spin_to_nonspin(indices, expectations, twopdm, D_CC_D, true); @@ -573,15 +573,15 @@ void compute_two_pdm_1_2_1(Wavefunction& wave1, Wavefunction& wave2, const SpinB Transposeview tlop = Transposeview(*leftop); vector expectations; - spinExpectation(wave1, wave2, *leftop, *dotop0, trop, big, expectations, false); - spinExpectation(wave1, wave2, *leftop, *dotop2, trop, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop), &(*dotop0), &trop, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop), &(*dotop2), &trop, big, expectations, false); vector indices(4,0); indices[0] = ix; indices[1] = jx; indices[2] = jx; indices[3] = kx; spin_to_nonspin(indices, expectations, twopdm, C_CD_D, true); expectations.resize(0); - spinExpectation(wave1, wave2, tlop, *dotop0, rightop, big, expectations, false); - spinExpectation(wave1, wave2, tlop, *dotop2, rightop, big, expectations, false); + spinExpectation(wave1, wave2, &tlop, &(*dotop0), &rightop, big, expectations, false); + spinExpectation(wave1, wave2, &tlop, &(*dotop2), &rightop, big, expectations, false); indices[0] = kx; indices[1] = jx; indices[2] = jx; indices[3] = ix; spin_to_nonspin(indices, expectations, twopdm, D_CD_C, true); @@ -622,8 +622,8 @@ void compute_two_pdm_2_1_1(Wavefunction& wave1, Wavefunction& wave2, const SpinB Transposeview tdop = Transposeview(dotop); vector expectations; - spinExpectation(wave1, wave2, *leftop0, tdop, trop, big, expectations, false); - spinExpectation(wave1, wave2, *leftop2, tdop, trop, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop0), &tdop, &trop, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop2), &tdop, &trop, big, expectations, false); vector indices(4,0); indices[0] = ix; indices[1] = ix; indices[2] = jx; indices[3] = kx; spin_to_nonspin(indices, expectations, twopdm, CC_D_D, true); @@ -648,8 +648,8 @@ void compute_two_pdm_2_1_1(Wavefunction& wave1, Wavefunction& wave2, const SpinB Transposeview tdop = Transposeview(dotop); vector expectations; - spinExpectation(wave1, wave2, *leftop0, dotop, trop, big, expectations, false); - spinExpectation(wave1, wave2, *leftop2, dotop, trop, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop0), &dotop, &trop, big, expectations, false); + spinExpectation(wave1, wave2, &(*leftop2), &dotop, &trop, big, expectations, false); vector indices(4,0); indices[0] = ix; indices[1] = jx; indices[2] = ix; indices[3] = kx; spin_to_nonspin(indices, expectations, twopdm, CD_CD, true); @@ -702,8 +702,8 @@ void compute_two_pdm_0_3_1(Wavefunction& wave1, Wavefunction& wave2, const SpinB Transposeview trop = Transposeview(rightop); vector expectations; - spinExpectation(wave1, wave2, *leftop, Dotop1, trop, big, expectations, false); - spinExpectation(wave1, wave2, *leftop, Dotop2, trop, big, expectations, false); + spinExpectation(wave1, wave2, leftop, &Dotop1, &trop, big, expectations, false); + spinExpectation(wave1, wave2, leftop, &Dotop2, &trop, big, expectations, false); vector indices(4,0); indices[0] = jx; indices[1] = jx; indices[2] = jx; indices[3] = kx; spin_to_nonspin(indices, expectations, twopdm, CC_D_D, true); @@ -789,8 +789,8 @@ void compute_two_pdm_0_3_1_notranspose(Wavefunction& wave1, Wavefunction& wave2, int kx = rightop.get_orbs(0); vector expectations; - spinExpectation(wave1, wave2, *leftop, Dotop1, rightop, big, expectations, false); - spinExpectation(wave1, wave2, *leftop, Dotop2, rightop, big, expectations, false); + spinExpectation(wave1, wave2, leftop, &Dotop1, &rightop, big, expectations, false); + spinExpectation(wave1, wave2, leftop, &Dotop2, &rightop, big, expectations, false); vector indices(4,0); indices[0] = jx; indices[1] = jx; indices[2] = jx; indices[3] = kx; spin_to_nonspin(indices, expectations, twopdm, CC_D_D, false); @@ -835,8 +835,8 @@ void compute_two_pdm_0_3_1_notranspose(Wavefunction& wave1, Wavefunction& wave2, int kx = rightop.get_orbs(0); vector expectations; - spinExpectation(wave1, wave2, *leftop, Dotop1, rightop, big, expectations, false); - spinExpectation(wave1, wave2, *leftop, Dotop2, rightop, big, expectations, false); + spinExpectation(wave1, wave2, leftop, &Dotop1, &rightop, big, expectations, false); + spinExpectation(wave1, wave2, leftop, &Dotop2, &rightop, big, expectations, false); vector indices(4,0); indices[0] = kx; indices[1] = jx; indices[2] = jx; indices[3] = jx; spin_to_nonspin(indices, expectations, twopdm, CC_D_D, false); @@ -892,8 +892,8 @@ void compute_two_pdm_3_0_1(Wavefunction& wave1, Wavefunction& wave2, const SpinB vector expectations; dotop = 0; - spinExpectation(wave1, wave2, Leftop1, *dotop, trop, big, expectations, false); - spinExpectation(wave1, wave2, Leftop2, *dotop, trop, big, expectations, false); + spinExpectation(wave1, wave2, &Leftop1, dotop, &trop, big, expectations, false); + spinExpectation(wave1, wave2, &Leftop2, dotop, &trop, big, expectations, false); vector indices(4,0); indices[0] = jx; indices[1] = jx; indices[2] = jx; indices[3] = kx; spin_to_nonspin(indices, expectations, twopdm, CC_D_D, true); @@ -938,8 +938,8 @@ void compute_two_pdm_3_1_0(Wavefunction& wave1, Wavefunction& wave2, const SpinB vector expectations; - spinExpectation(wave1, wave2, Leftop1, tdop, *rightop, big, expectations, false); - spinExpectation(wave1, wave2, Leftop2, tdop, *rightop, big, expectations, false); + spinExpectation(wave1, wave2, &Leftop1, &tdop, rightop, big, expectations, false); + spinExpectation(wave1, wave2, &Leftop2, &tdop, rightop, big, expectations, false); vector indices(4,0); indices[0] = jx; indices[1] = jx; indices[2] = jx; indices[3] = kx; @@ -1024,8 +1024,8 @@ void compute_two_pdm_1_3_0(Wavefunction& wave1, Wavefunction& wave2, const SpinB Transposeview tlop = Transposeview(leftop); vector expectations; - spinExpectation(wave1, wave2, tlop, Dotop1, *rightop, big, expectations, false); - spinExpectation(wave1, wave2, tlop, Dotop2, *rightop, big, expectations, false); + spinExpectation(wave1, wave2, &tlop, &Dotop1, rightop, big, expectations, false); + spinExpectation(wave1, wave2, &tlop, &Dotop2, rightop, big, expectations, false); vector indices(4,0); indices[0] = jx; indices[1] = jx; indices[2] = jx; indices[3] = kx; @@ -1095,8 +1095,8 @@ void compute_two_pdm_1_3(Wavefunction& wave1, Wavefunction& wave2, const SpinBlo SparseMatrix* dotop = 0; vector expectations; - spinExpectation(wave1, wave2, leftop, *dotop, Rightop1, big, expectations, false); - spinExpectation(wave1, wave2, leftop, *dotop, Rightop2, big, expectations, false); + spinExpectation(wave1, wave2, &leftop, dotop, &Rightop1, big, expectations, false); + spinExpectation(wave1, wave2, &leftop, dotop, &Rightop2, big, expectations, false); vector indices(4,0); indices[0] = kx; indices[1] = jx; indices[2] = jx; indices[3] = jx; @@ -1109,8 +1109,8 @@ void compute_two_pdm_1_3(Wavefunction& wave1, Wavefunction& wave2, const SpinBlo SparseMatrix* leftop = 0; vector expectations; - spinExpectation(wave1, wave2, *leftop, dotop, Rightop1, big, expectations, false); - spinExpectation(wave1, wave2, *leftop, dotop, Rightop2, big, expectations, false); + spinExpectation(wave1, wave2, leftop, &dotop, &Rightop1, big, expectations, false); + spinExpectation(wave1, wave2, leftop, &dotop, &Rightop2, big, expectations, false); vector indices(4,0); indices[0] = kx; indices[1] = jx; indices[2] = jx; indices[3] = jx; diff --git a/modules/twopdm/twopdm.h b/modules/twopdm/twopdm.h index 804b0c27..1e0c2430 100644 --- a/modules/twopdm/twopdm.h +++ b/modules/twopdm/twopdm.h @@ -43,8 +43,8 @@ void compute_two_pdm_3_1_0(Wavefunction& wave1, Wavefunction& wave2, const SpinB void compute_two_pdm_3_0_1(Wavefunction& wave1, Wavefunction& wave2, const SpinBlock& big, array_4d& twopdm); // done void compute_two_pdm_1_3(Wavefunction& wave1, Wavefunction& wave2, const SpinBlock& big, array_4d& twopdm);//0_1_3 and 1_0_3 // done -void spinExpectation(Wavefunction& wave1, Wavefunction& wave2, SparseMatrix &leftOp, SparseMatrix& dotOp, SparseMatrix& rightOp, const SpinBlock& big, vector& expectations, bool doTranspose); // Done -void FormLeftOp(const SpinBlock* leftBlock, const SparseMatrix& leftOp, const SparseMatrix& dotOp, SparseMatrix& Aop, int totalspin); // Done +void spinExpectation(Wavefunction& wave1, Wavefunction& wave2, SparseMatrix* leftOp_ptr, SparseMatrix* dotOp_ptr, SparseMatrix* rightOp_ptr, const SpinBlock& big, vector& expectations, bool doTranspose); // Done +void FormLeftOp(const SpinBlock* leftBlock, const SparseMatrix* leftOp_ptr, const SparseMatrix* dotOp_ptr, SparseMatrix& Aop, int totalspin); // Done void spin_to_nonspin(vector& indices, vector& coeffs, array_4d& twopdm, Oporder order, bool dotranspose); // Done void save_twopdm_text(const array_4d& twopdm, const int &i, const int &j); diff --git a/modules/twopdm/twopdm_2.C b/modules/twopdm/twopdm_2.C index 2155f2fd..6ba89d23 100644 --- a/modules/twopdm/twopdm_2.C +++ b/modules/twopdm/twopdm_2.C @@ -21,15 +21,19 @@ Sandeep Sharma and Garnet K.-L. Chan namespace SpinAdapted{ -void spinExpectation(Wavefunction& wave1, Wavefunction& wave2, SparseMatrix& leftOp, SparseMatrix& dotOp, SparseMatrix& rightOp, const SpinBlock& big, vector& expectations, bool doTranspose) + void spinExpectation(Wavefunction& wave1, Wavefunction& wave2, SparseMatrix* leftOp_ptr, SparseMatrix* dotOp_ptr, SparseMatrix* rightOp_ptr, const SpinBlock& big, vector& expectations, bool doTranspose) { + SparseMatrix& leftOp = *leftOp_ptr; + SparseMatrix& dotOp = *dotOp_ptr; + SparseMatrix& rightOp = *rightOp_ptr; + //calculating // do transpose specifies if we want separately. This can be avoided in some sitations if wave1 and wave2 are the same functions int leftindices=0, dotindices=0, rightindices=0; - leftindices = &leftOp ? leftOp.get_orbs().size() : 0; - dotindices = &dotOp ? dotOp.get_orbs().size() : 0; - rightindices = &rightOp ? rightOp.get_orbs().size() : 0; + leftindices = leftOp_ptr ? leftOp.get_orbs().size() : 0; + dotindices = dotOp_ptr ? dotOp.get_orbs().size() : 0; + rightindices = rightOp_ptr ? rightOp.get_orbs().size() : 0; int Aindices, Bindices; Aindices = leftindices+dotindices; @@ -43,10 +47,10 @@ void spinExpectation(Wavefunction& wave1, Wavefunction& wave2, SparseMatrix& lef SpinBlock* rightBlock = big.get_rightBlock(); Cre AOp; //This is just an example class - int totalspin = (&rightOp) ? rightOp.get_spin().getirrep() : 0; + int totalspin = rightOp_ptr ? rightOp.get_spin().getirrep() : 0; if (Aindices != 0) - FormLeftOp(leftBlock, leftOp, dotOp, AOp, totalspin); + FormLeftOp(leftBlock, leftOp_ptr, dotOp_ptr, AOp, totalspin); //different cases if (Aindices == 0 && Bindices == 4) @@ -74,13 +78,16 @@ void spinExpectation(Wavefunction& wave1, Wavefunction& wave2, SparseMatrix& lef } -void FormLeftOp(const SpinBlock* leftBlock, const SparseMatrix& leftOp, const SparseMatrix& dotOp, SparseMatrix& Aop, int totalspin) +void FormLeftOp(const SpinBlock* leftBlock, const SparseMatrix* leftOp_ptr, const SparseMatrix* dotOp_ptr, SparseMatrix& Aop, int totalspin) { + const SparseMatrix& leftOp = *leftOp_ptr; + const SparseMatrix& dotOp = *dotOp_ptr; + //Cre is just a class..it is not actually cre int leftindices=0, dotindices=0, rightindices=0; - leftindices = &leftOp ? leftOp.get_orbs().size() : 0; - dotindices = &dotOp ? dotOp.get_orbs().size() : 0; + leftindices = leftOp_ptr ? leftOp.get_orbs().size() : 0; + dotindices = dotOp_ptr ? dotOp.get_orbs().size() : 0; int Aindices, Bindices; Aindices = leftindices+dotindices; diff --git a/orbstring.h b/orbstring.h index 0d96f81d..88953788 100644 --- a/orbstring.h +++ b/orbstring.h @@ -32,7 +32,7 @@ class Orbstring : occ_rep (occ), sign (sn), empty (false) { // a little fix for IBMs xlC, see notes in Orbstring::distribute - for (int i = 0; i < occ_rep.size (); ++i) occ_rep [i] = abs (occ_rep [i]); + for (int i = 0; i < occ_rep.size (); ++i) occ_rep [i] = abs ((int)occ_rep [i]); } inline Orbstring (const bool* occ, const int sz, const int sn = 1) diff --git a/renormalise.C b/renormalise.C index e0d9bfd6..37e31f2a 100644 --- a/renormalise.C +++ b/renormalise.C @@ -8,7 +8,12 @@ Sandeep Sharma and Garnet K.-L. Chan #include "spinblock.h" +#if __has_include("boost/bind/bind.hpp") +#include +using namespace boost::placeholders; +#else #include +#endif #include #include #include