From 31ddcab69e532936dadc821ac4865f28cba05b05 Mon Sep 17 00:00:00 2001 From: Takeshi Yanai Date: Sun, 16 Oct 2016 01:05:24 +0900 Subject: [PATCH 1/8] boost_spirit: boost verson dependence is checked using BOOST_VERSION --- modules/npdm/npdm_spin_transformation.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/npdm/npdm_spin_transformation.h b/modules/npdm/npdm_spin_transformation.h index 8ec9038e..ed88e9e1 100644 --- a/modules/npdm/npdm_spin_transformation.h +++ b/modules/npdm/npdm_spin_transformation.h @@ -5,7 +5,7 @@ #include #include #include -#ifndef BOOST_1_56_0 +#if BOOST_VERSION >= 1056000 #include #include #include From c3b0f47af9c440e19d9802dc653970cec06bf5f7 Mon Sep 17 00:00:00 2001 From: Takeshi Yanai Date: Mon, 17 Oct 2016 11:04:57 +0900 Subject: [PATCH 2/8] boost/phoenix: boost verson dependence is checked using BOOST_VERSION --- modules/npdm/npdm_spin_transformation.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/modules/npdm/npdm_spin_transformation.h b/modules/npdm/npdm_spin_transformation.h index ed88e9e1..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 -#if BOOST_VERSION >= 1056000 +#if BOOST_VERSION < 105600 #include #include #include From 766606cbac1d093200fa3acda85fa093de67decd Mon Sep 17 00:00:00 2001 From: Takeshi Yanai Date: Mon, 24 Oct 2016 18:12:46 +0900 Subject: [PATCH 3/8] updates for building with clang-3.8 --- modules/npdm/fourpdm_container.C | 3 ++- orbstring.h | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) 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/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) From 8d8e3835512a647104ccf09d75c5481cf7c2fed4 Mon Sep 17 00:00:00 2001 From: Takeshi Yanai Date: Sun, 8 Jan 2017 00:38:24 +0900 Subject: [PATCH 4/8] workaround for the issue 'reference cannot be bound to dereferenced null pointer in well-defined C++ code; pointer may be assumed to always convert to true' --- modules/npdm/npdm_driver.C | 8 +-- modules/npdm/npdm_expectations.C | 35 +++++----- modules/npdm/npdm_expectations.h | 4 +- modules/npdm/npdm_expectations_engine.C | 7 +- modules/npdm/npdm_expectations_engine.h | 2 +- modules/twopdm/twopdm.C | 90 ++++++++++++------------- modules/twopdm/twopdm.h | 4 +- modules/twopdm/twopdm_2.C | 25 ++++--- 8 files changed, 93 insertions(+), 82 deletions(-) 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/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; From 3c1cb088d1979966819f0fa7441e27fcd0026c8a Mon Sep 17 00:00:00 2001 From: Takeshi Yanai Date: Thu, 26 Sep 2019 10:03:58 +0900 Subject: [PATCH 5/8] recommit --- input.C | 11 + input.C.new | 2860 +++++++++++++++++++++++++++++ input.h | 6 + input.h.new | 583 ++++++ modules/npdm/npdm.C | 129 +- modules/npdm/npdm.C.new | 598 ++++++ modules/npdm/npdm.h.new | 21 + modules/npdm/threepdm_container.C | 2 +- 8 files changed, 4148 insertions(+), 62 deletions(-) create mode 100644 input.C.new create mode 100644 input.h.new create mode 100644 modules/npdm/npdm.C.new create mode 100644 modules/npdm/npdm.h.new 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.C.new b/input.C.new new file mode 100644 index 00000000..055e093e --- /dev/null +++ b/input.C.new @@ -0,0 +1,2860 @@ +/* +Developed by Sandeep Sharma, Roberto Olivares-Amaya 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, Roberto Olivares-Amaya and Garnet K.-L. Chan +*/ + + +#include +#include +#include +#include "Symmetry.h" +#include "global.h" +#include "MatrixBLAS.h" +#include "spinblock.h" +#include "couplingCoeffs.h" +#include "genetic/GAOptimize.h" +#include "genetic/ReadIntegral.h" +#include +#include +#include +#include +#ifndef SERIAL +#include +#endif +#include +#include "fiedler.h" +#include "pario.h" + +using namespace std; + +namespace SpinAdapted { +string sym; +bool NonabelianSym; +} +void CheckFileExistence(string filename, string filetype); +void CheckFileInexistence(string filename, string filetype); + +void SpinAdapted::Input::ReadMeaningfulLine(ifstream& input, string& msg, int msgsize) +{ + bool readmore = true; + while (readmore && !input.eof()) { + char msgctr[msgsize]; + input.getline(msgctr, msgsize+1); + + msg=string(msgctr); + if(msg.size() == msgsize) { + pout << "in the process of reading line beginning with "<(world.size(),1); +#else + m_thrds_per_node = vector(1,1); +#endif + m_ham_type = QUANTUM_CHEMISTRY; + m_algorithm_type = TWODOT_TO_ONEDOT; + m_noise_type = RANDOM; + m_calc_type = DMRG; + m_solve_type = DAVIDSON; + m_stateSpecific = false; + m_implicitTranspose = true; //dont make DD just use CC^T to evaluate it + m_occupied_orbitals = -1; + m_num_Integrals = 1; + v_2.resize(1, TwoElectronArray(TwoElectronArray::restrictedNonPermSymm)); + v_1.resize(1); + coreEnergy.resize(1); + m_baseState.resize(1,0); + m_projectorState.resize(0); + m_targetState = -1; + m_guessState = 1; + m_permSymm = 2; + + m_openorbs.resize(0); + m_closedorbs.resize(0); + + m_spinAdapted = true; + m_Bogoliubov = false; + m_sys_add = 1; + m_env_add = 1; + + m_twodot_to_onedot_iter = 0; + m_integral_disk_storage_thresh = 1000; + m_max_lanczos_dimension = 5000; + + m_norbs = 0; + m_alpha = 0; + m_beta = 0; + m_hf_occ_user = ""; + + m_outputlevel = 0; + m_nquanta = 2; + m_total_symmetry_number = IrrepSpace( 0 ); + m_total_spin = 0; + m_guess_permutations = 10; + + m_direct = true; + m_nroots = 1; + m_weights.resize(m_nroots); + m_weights[0] = 1.; + + m_deflation_min_size = 2; + m_deflation_max_size = 20; + + m_transition_diff_spatial_irrep = false; + m_add_noninteracting_orbs = true; + m_bra_alpha = 0; + m_bra_beta =0; + m_non_SE = false; + m_no_transform = false; + m_do_fci = false; + m_do_npdm_ops = false; + m_do_npdm_in_core = false; + m_npdm_generate = false; + m_new_npdm_code = false; + m_do_pdm = false; + m_store_spinpdm = false; + m_spatpdm_disk_dump = false; + m_store_nonredundant_pdm =false; + m_pdm_unsorted = false; + + + m_nevpt_state_num = 0; + m_npdm_intermediate= true; + m_npdm_multinode= true; + + m_maxiter = 10; + m_oneindex_screen_tol = NUMERICAL_ZERO; + m_twoindex_screen_tol = NUMERICAL_ZERO; + + m_load_prefix = "."; + m_save_prefix = "."; + + m_maxj = 15; + m_ninej.init(m_maxj); + m_set_Sz = false; + + n_twodot_noise = 0; + m_twodot_noise = 1.0e-4; + m_twodot_gamma = 3.0e-1; + + m_sweep_tol = 1.0e-5; + m_restart = false; + m_fullrestart = false; + m_restart_warm = false; + m_backward = false; + m_reset_iterations = false; + + m_do_diis = false; + m_diis_error = 1e-2; + m_start_diis_iter = 8; + m_diis_keep_states = 6; + m_diis_error_tol = 1e-8; + m_num_spatial_orbs = 0; + m_schedule_type_default = false; + m_schedule_type_backward = false; + m_maxM = 0; + m_lastM = 500; + m_startM = 250; + + m_calc_ri_4pdm=false; + m_store_ripdm_readable=false; + m_nevpt2 = false; + m_conventional_nevpt2 = false; + m_kept_nevpt2_states = -1; + NevPrint.first=false; + NevPrint.second=0; + + //reorder options, by default it does fiedler + m_reorderType = FIEDLER; + m_reorderfile = ""; + m_gaconffile = "default"; + + m_orbformat=MOLPROFORM; + + m_warmup = LOCAL0; + + m_orz3rdmalt1 = -1; + m_orz3rdmalt2 = -1; +} + +void SpinAdapted::Input::usedkey_error(string& key, string& line) { + pout << "keyword "< usedkey(NUMKEYWORDS, -1); + int n_elec = -1; + int n_spin = -1; + int n_bra_spin =-1; //EL + sym = "c1"; + NonabelianSym = false; + std::vector orbitalfile; + + initialize_defaults(); + + if(mpigetrank() == 0) + { + pout << "Reading input file"< tok; + boost::split(tok, msg, is_any_of(", \t"), token_compress_on); + string keyword = *tok.begin(); + + if (boost::iequals(keyword, "orbs") || boost::iequals(keyword, "orbitals") || boost::iequals(keyword, " orbitals")) + { + if(usedkey[ORBS] == 0) + usedkey_error(keyword, msg); + usedkey[ORBS] = 0; + if (tok.size() < 2) { + pout << "keyword orbs should be followed by atleast a single filename and then an end line"< schd_tok; + boost::split(schd_tok, msg, is_any_of(" \t"), token_compress_on); + if (tok.size() != 1) { + if (boost::iequals(tok[1], "default")) { + m_schedule_type_default = true; + continue; + } + } + + m_sweep_iter_schedule.resize(0); + m_sweep_state_schedule.resize(0); + m_sweep_qstate_schedule.resize(0); + m_sweep_tol_schedule.resize(0); + m_sweep_noise_schedule.resize(0); + m_sweep_additional_noise_schedule.resize(0); + + int i = 0; + while(!boost::iequals(schd_tok[0], "END")) + { + + if (schd_tok.size() != 4) { + pout << "Each line of the schedule contain four entries sweep_iteration #retained states Davidson tolerance noise"<0 && m_sweep_iter_schedule[i] <= m_sweep_iter_schedule[i-1]) { + pout << "Sweep iteration at a given line should be higher than the previous sweep iteration"< and required for SOC and g-tensor calculations: + * preliminary calculations for the states should be performed without singlet embedding, otherwise the triplet excitation matrix elements will be equal to 0 +If 2 spins are given, the calculations of transition density matrix between wavefunctions with different spins (w/ or w/o different irreps) will be performed */ + { + n_bra_spin = atoi(tok[1].c_str()); + n_spin = atoi(tok[2].c_str()); + m_transition_diff_spin=true; + m_add_noninteracting_orbs = false; + m_bra_alpha = (n_elec + n_bra_spin)/2; + m_bra_beta = (n_elec - n_bra_spin)/2; + } + if ((tok.size() != 2)&&((tok.size() != 3))) { + pout << "keyword spin should be followed by a single number (or two spin numbers for SOC pdms) and then an end line"<(i)]= OneElectronArray(); +// } + for(int i= 0; i< TwoPertEnd; i++){ + //vpt_2[static_cast(i)]= TwoElectronArray(TwoElectronArray::restrictedNonPermSymm); + vpt_2[static_cast(i)]= PerturbTwoElectronArray(); + // vpt_2[i].rhf = true; + } + } + else if (boost::iequals(keyword, "restart_mps_nevpt")) { + m_calc_type = RESTART_MPS_NEVPT; + if (tok.size() !=4){ + pout << "keyword mps_nevpt should be followed by three numbers and then an end line"<(i)]= OneElectronArray(); +// } + for(int i= 0; i< TwoPertEnd; i++){ + //vpt_2[static_cast(i)]= TwoElectronArray(TwoElectronArray::restrictedNonPermSymm); + vpt_2[static_cast(i)]= PerturbTwoElectronArray(); + } + } + else if (boost::iequals(keyword, "nevpt_state_num" )) + { + if (tok.size() != 2) { + pout << "keyword nevpt_state_num should be followed by a single integer than then an end line."< pdms for SOC + else if (boost::iequals(keyword, "ds1_onepdm") || boost::iequals(keyword, "ds1_onerdm")) + m_calc_type = DS1_ONEPDM; + else if (boost::iequals(keyword, "restart_ds1_onepdm") || boost::iequals(keyword, "restart_ds1_onerdm")) + m_calc_type = RESTART_DS1_ONEPDM; +// Elvira: these are to calculate pdms to calculate S in the code for g-tensors + else if (boost::iequals(keyword, "ds0_onepdm") || boost::iequals(keyword, "ds0_onerdm")) + m_calc_type = DS0_ONEPDM; + else if (boost::iequals(keyword, "restart_ds0_onepdm") || boost::iequals(keyword, "restart_ds0_onerdm")) + m_calc_type = RESTART_DS0_ONEPDM; +// + else if (boost::iequals(keyword, "restart_tran_threepdm") || boost::iequals(keyword, "restart_tran_threerdm") || boost::iequals(keyword, "restart_tran_threerdm")) + { + m_calc_type = RESTART_T_THREEPDM; + m_restart = true; + } + else if (boost::iequals(keyword, "nevpt2") || boost::iequals(keyword, "pt2")){ + m_calc_type = NEVPT2; + m_nevpt2 = true; + m_transition_diff_spatial_irrep=false; + } + else if (boost::iequals(keyword, "ripdm")){ + m_calc_type = NEVPT2; + m_nevpt2 = false; + m_transition_diff_spatial_irrep=false; + } + else if (boost::iequals(keyword, "restart_nevpt2") || boost::iequals(keyword, "restart_pt2")){ + m_calc_type = RESTART_NEVPT2; + m_nevpt2 = true; + m_transition_diff_spatial_irrep=false; + m_restart = true; + } + else if (boost::iequals(keyword, "restart_ripdm")){ + m_calc_type = RESTART_NEVPT2; + m_nevpt2 = false; + m_transition_diff_spatial_irrep=false; + m_restart = true; + } + else if (boost::iequals(keyword, "calc_ri4pdm") || boost::iequals(keyword, "ri4pdm")) + { + m_calc_ri_4pdm = true; + } + else if (boost::iequals(keyword, "ripdm_readable")) + { + m_store_ripdm_readable = true; + } + else if (boost::iequals(keyword, "conventional_nevpt2")) + { + m_conventional_nevpt2 = true; + } + else if (boost::iequals(keyword, "M_nevpt2")){ + if (tok.size() != 2) { + pout << "keyword M_nevpt2 should be followed by a single float and then an end line."< :: iterator it = ++tok.begin(); + // the occupancies start from the second element of the tok string + if (tok.size() == 2 ) { + m_hf_occ_user = (*it).c_str(); + } + else{ + m_hf_occ_user = "manual"; + for( ; it != tok.end(); ++it ){ + int occ_i = atoi( (*it).c_str() ); + m_hf_occupancy.push_back( occ_i ); + } + } + pout << "m_hf_occ_user " << m_hf_occ_user << endl; + + } + + else if(boost::iequals(keyword, "open")) + { + vector :: iterator it = ++tok.begin(); + m_openorbs.resize(tok.size()-1,-1); + + int ii = 0; + for( ; it != tok.end(); ++it ){ + int openorb = atoi( (*it).c_str() ); + m_openorbs[ii] = openorb-1 ; + //pout << openorb<<" "; + ii++; + } + //pout << endl; + + } + else if(boost::iequals(keyword, "closed")) + { + vector :: iterator it = ++tok.begin(); + m_closedorbs.resize(tok.size()-1,-1); + //pout << "closed orbs :"; + int ii=0; + for( ; it != tok.end(); ++it ){ + int closedorb = atoi( (*it).c_str() ); + m_closedorbs[ii] = closedorb -1; + //pout << closedorb<<" "; + ii++; + } + //pout << endl; + } + else if(boost::iequals(keyword, "nroots")) + { + if(usedkey[NROOTS] == 0) + usedkey_error(keyword, msg); + usedkey[NROOTS] = 0; + std::string nroots_str; + if (tok.size() != 2) { + pout << "keyword nroots should be followed by a single integer and then an end line."< weighttoken; + boost::split(weighttoken, msg, is_any_of(" \t"), token_compress_on); + + if (!boost::iequals(weighttoken[0], "weights")) { + //manually make all the weights equal + m_weights.resize(m_nroots); + for (int i=0; i1 ) { + pout << "Currently the response code does not work with non-particle number conserving hamiltonians!!"; + abort(); + } + +//if (sym != "c1") // must be initialized even if c1 sym. + Symmetry::InitialiseTable(sym); + + if (mpigetrank() == 0) { + for (int l=0; l& oldtonew) { + string msg; int msgsize = 5000; + ReadMeaningfulLine(dumpFile, msg, msgsize); + vector tok; + std::vector diff; + boost::split(tok, msg, is_any_of(", \t"), token_compress_on); + while(msg.size() != 0) { + for (int i=0; im_norbs || oldtonew.back() < 0) { + pout << "Illegal orbital index "< tok; + int offset = m_orbformat == DMRGFORM ? 0 : 1; + std::ofstream ReorderFileOutput; + std::ifstream ReorderFileInput; + char ReorderFileName[5000]; + std::vector reorder; + + if (mpigetrank() == 0) { + ReadMeaningfulLine(dumpFile, msg, msgsize); + boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); + + if(offset != 0) + m_norbs = 2*atoi(tok[2].c_str()); + else + m_norbs = 2*atoi(tok[1].c_str()); + + m_num_spatial_orbs = 0; + m_spin_orbs_symmetry.resize(m_norbs); + m_spin_to_spatial.resize(m_norbs); + + //this is the file to which the reordering is written + sprintf(ReorderFileName, "%s%s", save_prefix().c_str(), "/RestartReorder.dat"); + } + boost::filesystem::path p(ReorderFileName); + +#ifndef SERIAL + boost::mpi::communicator world; + mpi::broadcast(world,m_restart,0); + mpi::broadcast(world,m_fullrestart,0); +#endif + + //do the reordering only if it is not a restart calculation + //if it is then just read the reorder.dat from the scratch space + if (integralIndex == 0) { + if(get_restart() || get_fullrestart() || m_calc_type == COMPRESS || m_calc_type == RESPONSE || m_calc_type == RESPONSEBW) { + if (mpigetrank() == 0) { + ReorderFileInput.open(ReorderFileName); + boost::filesystem::path ReorderFilePath(ReorderFileName); + + if(!boost::filesystem::exists(ReorderFilePath)) { + pout << "==============="<> m_reorder[i]; + ReorderFileInput.close(); + } + } + } + else { + if (mpigetrank() == 0) { + ReorderFileOutput.open(ReorderFileName); + } + + + // read the reorder file or calculate the reordering using one of the many options + if (m_reorderType == FIEDLER) { + + if (mpigetrank() == 0){ + m_reorder=get_fiedler(orbitalfile); + pout << "Fiedler-vector orbital ordering: "; + } + } + else if (m_reorderType == GAOPT) { + + ifstream gaconfFile; + + if (mpigetrank() == 0) { + if(m_gaconffile != "default") + gaconfFile.open(m_gaconffile.c_str(), ios::in); + m_reorder = get_fiedler(orbitalfile); + } + //to provide as initial guess to gaopt + m_reorder = getgaorder(gaconfFile, orbitalfile, m_reorder); + pout << "Genetic algorithm orbital ordering: "; + } + + else if (m_reorderType == MANUAL) { + if (mpigetrank() == 0) { + ifstream reorderFile(m_reorderfile.c_str()); + CheckFileExistence(m_reorderfile, "Reorder file "); + readreorderfile(reorderFile, m_reorder); + pout << "Manually provided orbital ordering: "; + } + } + else { //this is the no-reorder case + if (mpigetrank() == 0) { + m_reorder.resize(m_norbs/2); + for (int i=0; i O_unreordered(2,3) (the former is stored internally and the latter is what is given during input) + // but for this m_reorder the reorder vector below would be 3 1 2 4 and O_unreordered(1,2) -> O_reorder(3, 1) + if (mpigetrank() == 0) { + reorder.resize(m_norbs/2); + for (int i=0; i= 0 ) + ir = atoi(tok[i].c_str()) - offset; + else if (atoi(tok[i].c_str()) < -1) + ir = atoi(tok[i].c_str()) + offset; + + if (sym == "trans") ir += 1; //for translational symmetry the lowest irrep is 0 + if (sym == "lzsym") ir = atoi(tok[i].c_str()); + + + m_spin_orbs_symmetry[2*reorderOrbInd] = ir; + m_spin_orbs_symmetry[2*reorderOrbInd+1] = ir; + + if (readLine == numRead) { + m_num_spatial_orbs++; + m_spatial_to_spin.push_back(orbindex); + numRead = Symmetry::sizeofIrrep(ir); + readLine = 0; + } + m_spin_to_spatial[orbindex] = m_num_spatial_orbs-1; + m_spin_to_spatial[orbindex+1] = m_num_spatial_orbs-1; + orbindex +=2; + readLine++; + + + } + msg.resize(0); + ReadMeaningfulLine(dumpFile, msg, msgsize); + boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); + if(boost::iequals(tok[0], "IUHF")) RHF=false; + } + + + if(sym == "dinfh" ) { + m_spatial_to_spin.clear(); + for (int i=0; i= m_integral_disk_storage_thresh) // + { + for (int i=0; i orb; + for (int j=m_spatial_to_spin[i]; j& vpt2, double& coreEnergy) +{ + //TODO + //Reorder is not supported. + ifstream dumpFile; + dumpFile.open(orbitalfile.c_str(), ios::in); + + string msg; int msgsize = 5000; + vector tok; + int offset = m_orbformat == DMRGFORM ? 0 : 1; + std::ofstream ReorderFileOutput; + std::ifstream ReorderFileInput; + char ReorderFileName[5000]; + std::vector reorder(m_total_orbs); + + if (mpigetrank() == 0) { + ReadMeaningfulLine(dumpFile, msg, msgsize); + boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); + + m_norbs =m_act_size*2; + m_num_spatial_orbs = m_act_size; + m_spin_orbs_symmetry.resize(m_total_orbs*2); + m_spin_to_spatial.resize(m_total_orbs*2); + + m_total_spin_orbs_symmetry.resize(m_total_orbs*2); + m_total_spin_to_spatial.resize(m_total_orbs*2); + + //this is the file to which the reordering is written + sprintf(ReorderFileName, "%s%s", save_prefix().c_str(), "/RestartReorder.dat"); + } + + +#ifndef SERIAL + boost::mpi::communicator world; + mpi::broadcast(world,m_restart,0); + mpi::broadcast(world,m_fullrestart,0); +#endif + + //do the reordering only if it is not a restart calculation + //if it is then just read the reorder.dat from the scratch space + if(get_restart() || get_fullrestart()) { + if (mpigetrank() == 0) { + ReorderFileInput.open(ReorderFileName); + boost::filesystem::path ReorderFilePath(ReorderFileName); + + if(!boost::filesystem::exists(ReorderFilePath)) { + pout << "==============="<> m_reorder[i]; + ReorderFileInput.close(); + } + } + } + else { + if (mpigetrank() == 0) { + ReorderFileOutput.open(ReorderFileName); + } + + + // read the reorder file or calculate the reordering using one of the many options + + if (m_reorderType == FIEDLER) { + + if (mpigetrank() == 0){ + m_reorder=get_fiedler_nevpt(orbitalfile, m_act_size); + pout << "Fiedler-vector orbital ordering: "; + } + } + else if (m_reorderType == MANUAL) { + if (mpigetrank() == 0) { + ifstream reorderFile(m_reorderfile.c_str()); + CheckFileExistence(m_reorderfile, "Reorder file "); + readreorderfile(reorderFile, m_reorder); + pout << "Manually provided orbital ordering: "; + } + } + else { //this is the no-reorder case + if (mpigetrank() == 0) { + m_reorder.resize(m_act_size); + for (int i=0; i O_unreordered(2,3) (the former is stored internally and the latter is what is given during input) + // but for this m_reorder the reorder vector below would be 3 1 2 4 and O_unreordered(1,2) -> O_reorder(3, 1) + if (mpigetrank() == 0) { + reorder.resize(m_total_orbs); + for (int i=0; i= 0 ) + ir = atoi(tok[i].c_str()) - offset; + else if (atoi(tok[i].c_str()) < -1) + ir = atoi(tok[i].c_str()) + offset; + + if (sym == "trans") ir += 1; //for translational symmetry the lowest irrep is 0 + if (sym == "lzsym") ir = atoi(tok[i].c_str()); + + + m_total_spin_orbs_symmetry[2*reorderOrbInd] = ir; + m_total_spin_orbs_symmetry[2*reorderOrbInd+1] = ir; + + if (readLine == numRead) { + m_total_spatial_to_spin.push_back(orbindex); + numRead = Symmetry::sizeofIrrep(ir); + readLine = 0; + } + m_total_spin_to_spatial[orbindex] = orbindex/2; + m_total_spin_to_spatial[orbindex+1] = orbindex/2; + orbindex +=2; + readLine++; + + + } + msg.resize(0); + ReadMeaningfulLine(dumpFile, msg, msgsize); + boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); + if(boost::iequals(tok[0], "IUHF")) RHF=false; + } + + + if(sym == "dinfh" ) { + m_total_spatial_to_spin.clear(); + for (int i=0; i(m_total_spin_to_spatial.begin(),m_total_spin_to_spatial.begin()+m_act_size*2); +// m_spatial_to_spin=std::vector(m_total_spatial_to_spin.begin(),m_total_spatial_to_spin.begin()+m_act_size); +// m_spin_orbs_symmetry=std::vector(m_total_spin_orbs_symmetry.begin(),m_total_spin_orbs_symmetry.begin()+m_act_size*2); +// +// m_spatial_to_spin.push_back(m_act_size*2); +// m_spin_to_spatial.push_back(m_act_size*2); + + m_spin_to_spatial=m_total_spin_to_spatial; + m_spatial_to_spin=m_total_spatial_to_spin; + m_spin_orbs_symmetry=m_total_spin_orbs_symmetry; + + + + while (!((boost::iequals(tok[0], "&END")) || (boost::iequals(tok[0], "/")))) { + int temp; + + msg.resize(0); + ReadMeaningfulLine(dumpFile, msg, msgsize); + boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); + if(boost::iequals(tok[0], "IUHF")) RHF=false; + } + + int AOrbOffset = 0, BOrbOffset = 0; + if(!RHF) { + v1.rhf = false; + v2.rhf = false; + //for(auto i: vpt1) + // i.second.rhf = false; + vpt1.rhf = false; +// for(auto i: vpt2) +// i.second.rhf = false; +// fock.rhf = false; + } + v1.ReSize(m_total_orbs*2); + v2.ReSize(m_act_size*2); + vpt1.ReSize(m_total_orbs*2); + //for(auto i: vpt1) + // i.second.ReSize(m_total_orbs*2); + vpt2[Va].ReSize(m_total_orbs*2,m_act_size*2,m_act_size*2,m_act_size*2); + vpt2[Vi].ReSize(m_act_size*2,m_act_size*2,m_total_orbs*2,m_act_size*2); +// for(auto& i: vpt2) +// i.second.ReSize(m_total_orbs*2,m_total_orbs*2,m_total_orbs*2,m_total_orbs*2); + //fock.ReSize(m_total_orbs*2); + + + msg.resize(0); + ReadMeaningfulLine(dumpFile, msg, msgsize); //this if the first line with integrals + + if(m_total_orbs >= m_integral_disk_storage_thresh) // + { + assert(false && "Partial integral file for mps_nevpt2 is not implemented."); + } + + int i, j, k, l; + double value; + //Read zero order integral in active space; + while(msg.size() != 0) { + boost::split(tok, msg, is_any_of(" \t"), token_compress_on); + if (tok.size() != 5) { + pout << "error in reading orbital file"<(type)](2*I,2*K,2*J,2*L) = value; + } + msg.resize(0); + ReadMeaningfulLine(dumpFile, msg, msgsize); //this if the first line with integrals + } + } + + for(int type=0; type< 3 ; type++){ + msg.resize(0); + ReadMeaningfulLine(dumpFile, msg, msgsize); //this if the first line with integrals + while(msg.size() != 0) { + boost::split(tok, msg, is_any_of(" \t"), token_compress_on); + if (tok.size() != 5) { + pout << "error in reading orbital file"<(type)](2*reorder.at(i),2*reorder.at(j)) = value; +// } +// else { +// assert(false && "Only one electron integral perturbation"); +// } +// msg.resize(0); +// ReadMeaningfulLine(dumpFile, msg, msgsize); //this if the first line with integrals +// } +// +// } +// + dumpFile.close(); + } +} + +void SpinAdapted::Input::readorbitalsfile(string& orbitalfile, OneElectronArray& v1, TwoElectronArray& v2, double& coreEnergy, PairArray& vcc, CCCCArray& vcccc, CCCDArray& vcccd) { + + ifstream dumpFile; + dumpFile.open(orbitalfile.c_str(), ios::in); + + string msg; int msgsize = 5000; + vector tok; + std::ofstream ReorderFileOutput; + std::ifstream ReorderFileInput; + char ReorderFileName[5000]; + std::vector reorder; + + int offset = m_orbformat == DMRGFORM ? 0 : 1; + + if (mpigetrank() == 0) { + ReadMeaningfulLine(dumpFile, msg, msgsize); + boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); + + if(offset != 0) + m_norbs = 2*atoi(tok[2].c_str()); + else + m_norbs = 2*atoi(tok[1].c_str()); + + m_num_spatial_orbs = 0; + m_spin_orbs_symmetry.resize(m_norbs); + m_spin_to_spatial.resize(m_norbs); + + sprintf(ReorderFileName, "%s%s", save_prefix().c_str(), "/RestartReorder.dat"); + } + boost::filesystem::path p(ReorderFileName); + +#ifndef SERIAL + boost::mpi::communicator world; + mpi::broadcast(world,m_restart,0); + mpi::broadcast(world,m_fullrestart,0); +#endif + + //do the reordering only if it is not a restart calculation + //if it is then just read the reorder.dat from the scratch space + if(get_restart() || get_fullrestart()) { + if (mpigetrank() == 0) { + ReorderFileInput.open(ReorderFileName); + if(!boost::filesystem::exists(p)) { + pout << "==============="<> m_reorder[i]; + ReorderFileInput.close(); + } + } + } + else { + if (mpigetrank() == 0) { + ReorderFileOutput.open(ReorderFileName); + } + + // read the reorder file or calculate the reordering using one of the many options + if (m_reorderType == FIEDLER) { + if (mpigetrank() == 0) { + m_reorder=get_fiedler_bcs(orbitalfile); + pout << "Fiedler-vector orbital ordering: "; + } + } else if (m_reorderType == GAOPT) { + if (mpigetrank() == 0) { + pout << "GAOPT orbital ordering for BCS calculation not implemented"; + } + abort(); + } else if (m_reorderType == MANUAL) { + if (mpigetrank() == 0) { + ifstream reorderFile(m_reorderfile.c_str()); + CheckFileExistence(m_reorderfile, "Reorder file "); + readreorderfile(reorderFile, m_reorder); + pout << "Manually provided orbital ordering: "; + } + } else { //this is the no-reorder case + if (mpigetrank() == 0) { + m_reorder.resize(m_norbs/2); + for (int i=0; i O_unreordered(2,3) (the former is stored internally and the latter is what is given during input) + // but for this m_reorder the reorder vector below would be 3 1 2 4 and O_unreordered(1,2) -> O_reorder(3, 1) + if (mpigetrank() == 0) { + reorder.resize(m_norbs/2); + for (int i=0; i= 0 ) + ir = atoi(tok[i].c_str()) - offset; + else if (atoi(tok[i].c_str()) < -1) + ir = atoi(tok[i].c_str()) + offset; + + if (sym == "trans") ir += 1; //for translational symmetry the lowest irrep is 0 + if (sym == "lzsym") ir = atoi(tok[i].c_str()); + + m_spin_orbs_symmetry[2*reorderOrbInd] = ir; + m_spin_orbs_symmetry[2*reorderOrbInd+1] = ir; + + if (readLine == numRead) { + m_num_spatial_orbs++; + m_spatial_to_spin.push_back(orbindex); + numRead = Symmetry::sizeofIrrep(ir); + readLine = 0; + } + m_spin_to_spatial[orbindex] = m_num_spatial_orbs-1; + m_spin_to_spatial[orbindex+1] = m_num_spatial_orbs-1; + orbindex +=2; + readLine++; + } + msg.resize(0); + ReadMeaningfulLine(dumpFile, msg, msgsize); + boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); + if(boost::iequals(tok[0], "IUHF")) + RHF=false; + } + + if(sym == "dinfh" ) { + m_spatial_to_spin.clear(); + for (int i=0; i= m_integral_disk_storage_thresh) // + { + for (int i=0; i orb; + for (int j=m_spatial_to_spin[i]; j SpinAdapted::Input::get_fiedler(string& dumpname){ + Matrix fiedler; + ifstream dumpFile; + dumpFile.open(dumpname.c_str(), ios::in); + genetic::ReadIntegral(dumpFile, fiedler); + dumpFile.close(); + SymmetricMatrix fiedler_sym; + fiedler_sym << fiedler; + std::vector findices = fiedler_reorder(fiedler_sym); + return findices; +} + +std::vector SpinAdapted::Input::get_fiedler_nevpt(string& dumpname, int nact){ + Matrix fiedler; + ifstream dumpFile; + dumpFile.open(dumpname.c_str(), ios::in); + genetic::ReadIntegral_nevpt(dumpFile, fiedler, nact); + dumpFile.close(); + SymmetricMatrix fiedler_sym; + fiedler_sym << fiedler; + std::vector findices = fiedler_reorder(fiedler_sym); + return findices; +} + +std::vector SpinAdapted::Input::get_fiedler_bcs(string& dumpname) { + Matrix fiedler; + ifstream dumpFile; + dumpFile.open(dumpname.c_str(), ios::in); + genetic::ReadIntegral_BCS(dumpFile, fiedler); + dumpFile.close(); + SymmetricMatrix fiedler_sym; + fiedler_sym << fiedler; + std::vector findices = fiedler_reorder(fiedler_sym); + return findices; +} + +std::vector SpinAdapted::Input::getgaorder(ifstream& gaconfFile, string& orbitalfile, std::vector& fiedlerorder) +{ + ifstream dumpFile; dumpFile.open(orbitalfile.c_str()); + return genetic::gaordering(gaconfFile, dumpFile, fiedlerorder).Gen().Sequence(); +} + +#ifdef MOLPRO +void SpinAdapted::Input::writeSummaryForMolpro() +{ + pout << setw(50) << "Total number of orbitals : " ; + pout << m_norbs/2 << endl; + pout << setw(50) << "Symmetry of targeted wavefunctions : " ; + pout << m_alpha + m_beta << ":" << m_alpha-m_beta << ":" << m_total_symmetry_number.getirrep()+1 << endl; + pout << setw(50) << "Number of wavefunctions targeted : " ; + pout << m_nroots << endl; + if (m_nroots >1) { + pout << setw(50) << "The weights of the wavefunctions : "; + for (int i=0; i1) { +// printf("%-50s : ", "The weights of the wavefunctions"); +// for (int i=0; i 1) + { + pout << setw(50) << left << "The weights of the wavefunctions" << " : "; + for (int i = 0; i < m_nroots; ++i) pout << left << setprecision(2) << setw(10) << scientific << m_weights[i]; + pout << endl; + } + pout << setw(50) << left << "Symmetry of the molecule" << " : " + << sym.c_str() << endl; + if (sym != "c1") + { + pout << setw(50) << left << "Irreducible representations of the orbitals" << " : "; + for (int i = 0; i < m_spin_orbs_symmetry.size(); i+=2) + { + pout << Symmetry::stringOfIrrep(m_spin_orbs_symmetry[i]) << " "; + } + pout << endl; + } + + pout << endl << "Schedule" << endl; + pout << "--------" << endl; + pout << setw(10) << left << "Iter" << " : " + << setw(20) << left << "# States" << " " + << setw(20) << left << "Davidson_tol" << " " + << setw(20) << left << "Random_noise" << endl; + + for (int i = 0; i < m_sweep_iter_schedule.size(); ++i) + { + pout << setw(10) << left << m_sweep_iter_schedule[i] << " : " + << setw(20) << left << m_sweep_state_schedule[i] << " " + << setw(20) << left << scientific << m_sweep_tol_schedule[i] << " " + << setw(20) << left << scientific << m_sweep_noise_schedule[i] << endl; + } + if (m_algorithm_type == TWODOT_TO_ONEDOT) + { + pout << setw(50) << left << "Switching from twodot to onedot algorithm" << " : " + << m_twodot_to_onedot_iter << endl; + } + pout << setw(50) << left << "Maximum sweep iterations" << " : " << m_maxiter << endl; + +#ifndef SERIAL + } +#endif +} + +void SpinAdapted::Input::performSanityTest() +{ +#ifndef SERIAL + if (mpigetrank() == 0) { +#endif + if (m_num_spatial_orbs <= 0) { + pout << "total number of orbitals has to be a positive number"< 200) { + pout << "Number of orbitals cannot be greater than 130"< 10000) { + pout << "default schedule only works for maxM less than 10000"< m_startM) { + pout << "lastM is larger than startM" << endl; + pout << "Make sure you specify a lastM larger than " << lastM << endl; + pout << "or specify a larger startM " << endl; + abort(); + } + } + else { + // This needs to go at the end of sanity check, or move up schedule + if (m_maxM != 0) { + pout << "With detailed schedule a non-zero maxM should not be specified"<= m_maxiter) { + pout << "Switch from twodot to onedot algorithm cannot happen after maxiter"< openCopy, closedCopy, reorder; + openCopy.resize(m_openorbs.size(), -1); + closedCopy.resize(m_closedorbs.size(), -1); + + for (int i=0; i hf_occupancy_tmp(m_num_spatial_orbs*2,0); + + if (m_Bogoliubov) { + // overwrite hf_occ_user option, since initial guess is always vacuum in BCS case + m_hf_occupancy.assign(m_norbs, 0); + } + else if (m_hf_occ_user == "manual") { + //check if n_orbs is correct and if n_elec is correct + if (m_hf_occupancy.size() != m_num_spatial_orbs ) { + pout << "ERROR: The length of user-defined HF occupancies does not match the number of orbitals " << endl; + pout << "Length of occupancies is: " << m_hf_occupancy.size() << ", and the number of orbitals is: " << m_num_spatial_orbs<< endl; + abort(); + } + int UserElectrons = 0; + for (int i=0; i ele_map; +// for( int i = 0; i < m_norbs/2; ++i ){ +// ele_map.insert( pair( v_1[0](2*i, 2*i), i ) ); +// } + +// multimap :: iterator it_alpha = ele_map.begin(); +// for( int i = 0; i < m_alpha; ++i ){ +// int ia = it_alpha->second; +// hf_occupancy_tmp.at( 2*ia ) = 1; +// if (i < m_beta) +// hf_occupancy_tmp.at( 2*ia+1 ) = 1; +// ++it_alpha; +// } + hf_occupancy_tmp = this->hfOccGenerator_(); + } + else { + pout << "currently other options besides manual, integral and canonical are not implemented."< reorder(m_num_spatial_orbs); + for (int i=0; i= dmrginp.sweep_iter_schedule()[iter]) + current = iter; + + } + + int nroots = m_nroots; + if (m_noise_type == EXCITEDSTATE && m_sweep_additional_noise_schedule[current] != 0.0 && nroots == 1) + nroots++; + + if (setStateSpecific()) + return 1; + + return nroots; +} + +std::vector SpinAdapted::Input::weights(int sweep_iter) const +{ + int iter; + int current = 0; + for (iter = 0; iter < dmrginp.sweep_iter_schedule().size(); ++iter) + { + if (sweep_iter >= dmrginp.sweep_iter_schedule()[iter]) + current = iter; + + } + + if (setStateSpecific()) + return std::vector(1,1.0); + + int nroots = this->nroots(sweep_iter); + std::vector weights(nroots); + for (int i=0; i< m_nroots; i++) + weights[i] = m_weights[i]; + if (nroots == m_nroots+1) + weights[nroots-1] = m_sweep_additional_noise_schedule[current]; + if (nroots != m_nroots+1 && nroots != m_nroots) + { + pout << "Something wrong with the number of weights of wavefunctions!"<& sites = b.get_sites(); + int n=0, s=0; + for (int i=0; i NevPrint; + int m_orz3rdmalt1; + int m_orz3rdmalt2; + friend class boost::serialization::access; template void serialize(Archive & ar, const unsigned int version) @@ -234,6 +237,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 +576,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/input.h.new b/input.h.new new file mode 100644 index 00000000..50e4a35b --- /dev/null +++ b/input.h.new @@ -0,0 +1,583 @@ +/* +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 +*/ + + +#ifndef SPIN_INPUT_HEADER_H +#define SPIN_INPUT_HEADER_H +#include +#include +#include +#include +#include +#include "IrrepSpace.h" +#include "SpinQuantum.h" +#include "timer.h" +#include "couplingCoeffs.h" +#include +#include "IntegralMatrix.h" + + +namespace SpinAdapted{ +class SpinBlock; +class OneElectronArray; +class TwoElectronArray; +//typedef OneElectronArray PerturbOneElectronArray; +//typedef TwoElectronArray PerturbTwoElectronArray; +class PerturbTwoElectronArray; +class PairArray; +class CCCCArray; +class CCCDArray; +enum OnePerturbType{ Va_1=0, Vi_1, Vai_1, OnePertEnd}; +enum TwoPerturbType{Va=0,Vi,Vab,Vij,Vai,Vabi,Vaij,Vabij,TwoPertEnd}; +//enum OnePerturbType : int; +//enum TwoPerturbType : int; + +enum WarmUpTypes {WILSON, LOCAL0, LOCAL2, LOCAL3, LOCAL4}; +enum hamTypes {QUANTUM_CHEMISTRY, HUBBARD, BCS, HEISENBERG}; +enum solveTypes {LANCZOS, DAVIDSON, CONJUGATE_GRADIENT}; +enum algorithmTypes {ONEDOT, TWODOT, TWODOT_TO_ONEDOT}; +enum noiseTypes {RANDOM, EXCITEDSTATE}; +enum calcType {DMRG, ONEPDM, TWOPDM, THREEPDM, FOURPDM, NEVPT2PDM, RESTART_TWOPDM, + RESTART_ONEPDM, RESTART_THREEPDM, RESTART_FOURPDM, RESTART_NEVPT2PDM, TINYCALC, FCI, + EXCITEDDMRG, CALCOVERLAP, CALCHAMILTONIAN, COMPRESS, RESPONSE, RESPONSEBW, + TRANSITION_ONEPDM, TRANSITION_TWOPDM, TRANSITION_THREEPDM, RESTART_T_ONEPDM, RESTART_T_TWOPDM, RESTART_T_THREEPDM, + NEVPT2,RESTART_NEVPT2, MPS_NEVPT, RESTART_MPS_NEVPT, + DS1_ONEPDM, RESTART_DS1_ONEPDM, DS0_ONEPDM, RESTART_DS0_ONEPDM}; +enum orbitalFormat{MOLPROFORM, DMRGFORM}; +enum reorderType{FIEDLER, GAOPT, MANUAL, NOREORDER}; +enum keywords{ORBS, LASTM, STARTM, MAXM, REORDER, HF_OCC, SCHEDULE, SYM, NELECS, SPIN, IRREP, + MAXJ, PREFIX, NROOTS, DOCD, DEFLATION_MAX_SIZE, MAXITER, BASENERGY, + SCREEN_TOL, ODOT, SWEEP_TOL, OUTPUTLEVEL, NONSPINADAPTED, BOGOLIUBOV, TWODOT_NOISE, WARMUP, NPDM_INTERMEDIATE, NPDM_NO_INTERMEDIATE, NPDM_MULTINODE, NPDM_NO_MULTINODE, SPECIFICPDM, NUMKEYWORDS}; + +class Input { + + private: + std::vector m_thrds_per_node; + int m_norbs; + int m_alpha; + int m_beta; +// + int m_bra_alpha; + int m_bra_beta; + int n_bra_spin; +// + int m_Sz; + bool m_spinAdapted; + bool m_Bogoliubov; + int m_permSymm; + + IrrepSpace m_total_symmetry_number; + IrrepSpace m_bra_symmetry_number;// This is used when bra and ket have different spatial symmetry irrep; + // It is only used for transition density matrix calculations. + bool m_transition_diff_spatial_irrep; + bool m_transition_diff_spin =false; //Elvira: This is used when bra and ket have different spins + SpinQuantum m_bra_molecule_quantum; // It is only to calculate pdms for SOC + + SpinQuantum m_molecule_quantum; + int m_total_spin; + int m_guess_permutations; + bool m_stateSpecific; //when targetting excited states we switch from state + //average to statespecific + int m_occupied_orbitals; + + vector m_openorbs; + vector m_closedorbs; + vector m_baseState; + vector m_projectorState; + int m_targetState; + int m_guessState; + + std::vector m_hf_occupancy; + std::string m_hf_occ_user; + std::vector m_weights; + + std::vector m_sweep_iter_schedule; + std::vector m_sweep_state_schedule; + std::vector m_sweep_qstate_schedule; + std::vector m_sweep_tol_schedule; + std::vector m_sweep_noise_schedule; + std::vector m_sweep_additional_noise_schedule; + bool m_schedule_type_default; + bool m_schedule_type_backward; + int m_lastM; + int m_startM; + int m_maxM; + int m_integral_disk_storage_thresh; + int m_num_Integrals; + + bool m_do_diis; + double m_diis_error; + int m_start_diis_iter; + int m_diis_keep_states; + double m_diis_error_tol; + + calcType m_calc_type; + noiseTypes m_noise_type; + hamTypes m_ham_type; + WarmUpTypes m_warmup; + int m_nroots; + solveTypes m_solve_type; + bool m_do_deriv; + bool m_do_fci; + bool m_do_pdm; + bool m_do_npdm_ops; + bool m_do_npdm_in_core; + bool m_npdm_generate; + bool m_new_npdm_code; + bool m_store_spinpdm; + bool m_spatpdm_disk_dump; + bool m_pdm_unsorted; + bool m_store_nonredundant_pdm; + bool m_npdm_intermediate; + std::vector m_specificpdm; + bool m_npdm_multinode; + bool m_set_Sz; + int m_maxiter; + double m_oneindex_screen_tol; + double m_twoindex_screen_tol; + bool m_no_transform; + bool m_add_noninteracting_orbs; + bool m_non_SE; + int m_nquanta; + int m_sys_add; + int m_env_add; + + int m_deflation_min_size; + int m_deflation_max_size; + + algorithmTypes m_algorithm_type; + int m_twodot_to_onedot_iter; + std::vector< std::map > m_quantaToKeep; + std::string m_save_prefix; + std::string m_load_prefix; + bool m_direct; + std::vector m_orbenergies; + + int m_maxj; + ninejCoeffs m_ninej; + int m_max_lanczos_dimension; + + int n_twodot_noise; + double m_twodot_noise; + double m_twodot_gamma; + + double m_sweep_tol; + bool m_restart; + bool m_backward; + bool m_fullrestart; + bool m_restart_warm; + bool m_reset_iterations; + bool m_implicitTranspose; + + + std::vector m_spin_vector; + std::vector m_spin_orbs_symmetry; + int m_num_spatial_orbs; + std::vector m_spatial_to_spin; + std::vector m_spin_to_spatial; + + int m_act_size; + int m_core_size; + int m_virt_size; + int m_total_orbs; + int m_nevpt_state_num; + std::vector m_total_spin_orbs_symmetry; + std::vector m_total_spatial_to_spin; + std::vector m_total_spin_to_spatial; + + int m_outputlevel; + orbitalFormat m_orbformat; + + int m_reorderType; + string m_reorderfile; + std::vector m_reorder;//this can be manual, fiedler, gaopt or noreorder + string m_gaconffile; + + bool m_calc_ri_4pdm; + bool m_store_ripdm_readable; + bool m_nevpt2; + bool m_conventional_nevpt2; + 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) + { + ar & m_thrds_per_node & m_spinAdapted & m_Bogoliubov & m_stateSpecific & m_implicitTranspose & m_num_Integrals; + ar & m_norbs & m_alpha & m_beta & m_solve_type & m_Sz & m_set_Sz & m_baseState& m_projectorState& m_targetState; + ar & m_spin_vector & m_spin_orbs_symmetry & m_guess_permutations & m_nroots & m_weights & m_hf_occ_user & m_hf_occupancy; + ar & m_sweep_iter_schedule & m_sweep_state_schedule & m_sweep_qstate_schedule & m_sweep_tol_schedule & m_sweep_noise_schedule &m_sweep_additional_noise_schedule & m_reorder; + ar & m_molecule_quantum & m_total_symmetry_number & m_total_spin & m_orbenergies & m_add_noninteracting_orbs; + ar & m_bra_symmetry_number & m_permSymm & m_openorbs & m_closedorbs; + ar & m_bra_molecule_quantum & m_bra_symmetry_number & m_add_noninteracting_orbs; //diff spins case + ar & m_bra_alpha & m_bra_beta & n_bra_spin; //diff spins case + ar & m_transition_diff_spin; + ar & m_non_SE; +// ar & m_save_prefix & m_load_prefix & m_direct & m_max_lanczos_dimension; + ar & m_direct & m_max_lanczos_dimension; + ar & m_deflation_min_size & m_deflation_max_size & m_outputlevel & m_reorderfile; + ar & m_algorithm_type & m_twodot_to_onedot_iter & m_orbformat ; + ar & m_nquanta & m_sys_add & m_env_add & m_do_fci & m_no_transform ; + ar & m_do_pdm & m_do_npdm_ops & m_do_npdm_in_core & m_npdm_generate & m_new_npdm_code & m_specificpdm; + ar & m_transition_diff_spatial_irrep & m_occupied_orbitals; + ar & m_store_spinpdm &m_spatpdm_disk_dump & m_pdm_unsorted & m_npdm_intermediate & m_npdm_multinode; + ar & m_maxj & m_ninej & m_maxiter & m_do_deriv & m_oneindex_screen_tol & m_twoindex_screen_tol & m_quantaToKeep & m_noise_type; + ar & m_sweep_tol & m_restart & m_backward & m_fullrestart & m_restart_warm & m_reset_iterations & m_calc_type & m_ham_type & m_warmup; + ar & m_do_diis & m_diis_error & m_start_diis_iter & m_diis_keep_states & m_diis_error_tol & m_num_spatial_orbs; + ar & m_spatial_to_spin & m_spin_to_spatial & m_maxM & m_schedule_type_backward & m_schedule_type_default &m_integral_disk_storage_thresh; + 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; + } + + + void initialize_defaults(); + + std::vector hfOccGenerator_ (); + + public: + //Input() : m_ninej(ninejCoeffs::getinstance()){} + Input() {} + Input (const std::string& config_name); + // ROA + void initCumulTimer() + { + ddscreen = boost::shared_ptr (new cumulTimer()); + cdscreen = boost::shared_ptr (new cumulTimer()); + dscreen = boost::shared_ptr (new cumulTimer()); + buildcsfops = boost::shared_ptr (new cumulTimer()); + sysdotmake = boost::shared_ptr (new cumulTimer()); + initnewsystem = boost::shared_ptr (new cumulTimer()); + guessgenT = boost::shared_ptr (new cumulTimer()); + multiplierT = boost::shared_ptr (new cumulTimer()); + operrotT = boost::shared_ptr (new cumulTimer()); + davidsonT = boost::shared_ptr (new cumulTimer()); + rotmatrixT = boost::shared_ptr (new cumulTimer()); + blockdavid = boost::shared_ptr (new cumulTimer()); + datatransfer = boost::shared_ptr (new cumulTimer()); + hmultiply = boost::shared_ptr (new cumulTimer()); + oneelecT = boost::shared_ptr (new cumulTimer()); + twoelecT = boost::shared_ptr (new cumulTimer()); + makeopsT = boost::shared_ptr (new cumulTimer()); + collectqT = boost::shared_ptr (new cumulTimer()); + opallocate = boost::shared_ptr (new cumulTimer()); + opcatenate = boost::shared_ptr (new cumulTimer()); + oprelease = boost::shared_ptr (new cumulTimer()); + opequateT = boost::shared_ptr (new cumulTimer()); + justmultiply = boost::shared_ptr (new cumulTimer()); + spinrotation = boost::shared_ptr (new cumulTimer()); + otherrotation = boost::shared_ptr (new cumulTimer()); + solvewf = boost::shared_ptr (new cumulTimer()); + postwfrearrange = boost::shared_ptr (new cumulTimer()); + couplingcoeff = boost::shared_ptr (new cumulTimer()); + buildsumblock = boost::shared_ptr (new cumulTimer()); + buildblockops = boost::shared_ptr (new cumulTimer()); + addnoise = boost::shared_ptr (new cumulTimer()); + s0time = boost::shared_ptr (new cumulTimer()); + s1time = boost::shared_ptr (new cumulTimer()); + s2time = boost::shared_ptr (new cumulTimer()); + blockintegrals = boost::shared_ptr (new cumulTimer()); + blocksites = boost::shared_ptr (new cumulTimer()); + statetensorproduct = boost::shared_ptr (new cumulTimer()); + statecollectquanta = boost::shared_ptr (new cumulTimer()); + builditeratorsT = boost::shared_ptr (new cumulTimer()); + diskio = boost::shared_ptr (new cumulTimer()); + } + void writeSummary(); +#ifdef MOLPRO + void writeSummaryForMolpro(); +#endif + void performSanityTest(); + void generateDefaultSchedule(); + void readorbitalsfile(string& dumpFile, OneElectronArray& v1, TwoElectronArray& v2, double& coreEnergy, int integralIndex); + void readorbitalsfile(string& dumpFile, OneElectronArray& v1, TwoElectronArray& v2, OneElectronArray& vpt1, std::map& vpt2, double& coreEnergy); + void readorbitalsfile(string& dumpFile, OneElectronArray& v1, TwoElectronArray& v2, double& coreEnergy, PairArray& vcc, CCCCArray& vcccc, CCCDArray& vcccd); + int getNumIntegrals() { return m_num_Integrals;} + void readreorderfile(ifstream& dumpFile, std::vector& reorder); + std::vector getgaorder(ifstream& gaconfFile, string& orbitalfile, std::vector& fiedlerorder); + std::vector get_fiedler(string& dumpname); + std::vector get_fiedler_nevpt(string& dumpname, int nact); + std::vector get_fiedler_bcs(string& dumpname); + void usedkey_error(string& key, string& line); + void makeInitialHFGuess(); + static void ReadMeaningfulLine(ifstream&, string&, int); + + boost::shared_ptr dscreen; + boost::shared_ptr cdscreen; + boost::shared_ptr ddscreen; + boost::shared_ptr buildcsfops; + boost::shared_ptr sysdotmake; + boost::shared_ptr initnewsystem; + boost::shared_ptr diskio; + boost::shared_ptr builditeratorsT; + boost::shared_ptr blockintegrals; + boost::shared_ptr blocksites; + boost::shared_ptr statetensorproduct; + boost::shared_ptr statecollectquanta; + boost::shared_ptr guessgenT; + boost::shared_ptr multiplierT; + boost::shared_ptr operrotT; + boost::shared_ptr davidsonT; + boost::shared_ptr rotmatrixT; + boost::shared_ptr blockdavid; + boost::shared_ptr datatransfer; + boost::shared_ptr hmultiply; + boost::shared_ptr oneelecT; + boost::shared_ptr twoelecT; + boost::shared_ptr makeopsT; + boost::shared_ptr collectqT; + boost::shared_ptr opallocate; + boost::shared_ptr opcatenate; + boost::shared_ptr oprelease; + boost::shared_ptr opequateT; + boost::shared_ptr justmultiply; + boost::shared_ptr spinrotation; + boost::shared_ptr otherrotation; + boost::shared_ptr solvewf; + boost::shared_ptr postwfrearrange; + boost::shared_ptr couplingcoeff; + boost::shared_ptr buildsumblock; + boost::shared_ptr buildblockops; + boost::shared_ptr addnoise; + boost::shared_ptr s0time; + boost::shared_ptr s1time; + boost::shared_ptr s2time; + + std::vector& get_openorbs() { return m_openorbs;} + std::vector& get_closedorbs() { return m_closedorbs;} + const std::vector& baseStates() const {return m_baseState;} + const int& targetState() const {return m_targetState;} + const int& guessState() const {return m_guessState;} + int& setGuessState() {return m_guessState;} + const std::vector& projectorStates() const {return m_projectorState;} + + std::vector& baseStates() {return m_baseState;} + int& targetState() {return m_targetState;} + std::vector& projectorStates() {return m_projectorState;} + + void reorderOpenAndClosed(); + const int& num_occupied_orbitals() const {return m_occupied_orbitals;} + const bool& doimplicitTranspose() const {return m_implicitTranspose;} + bool& setimplicitTranspose() {return m_implicitTranspose;} + const bool& setStateSpecific() const {return m_stateSpecific;} + bool& setStateSpecific() {return m_stateSpecific;} + const orbitalFormat& orbformat() const {return m_orbformat;} + const int& outputlevel() const {return m_outputlevel;} + int& setOutputlevel() {return m_outputlevel;} + const vector& spatial_to_spin() const {return m_spatial_to_spin;} + int spatial_to_spin(int i) const {return m_spatial_to_spin.at(i);} + const vector& spin_to_spatial() const {return m_spin_to_spatial;} + const double& diis_error_tol() const {return m_diis_error_tol;} + const bool& do_diis() const {return m_do_diis;} + const double& diis_error() const {return m_diis_error;} + const int& start_diis_iter() const {return m_start_diis_iter;} + const int& diis_keep_states() const {return m_diis_keep_states;} + + bool use_partial_two_integrals() const {return (m_norbs/2 >= m_integral_disk_storage_thresh);} + bool& set_fullrestart() {return m_fullrestart;} + const bool& get_fullrestart() const {return m_fullrestart;} + const bool& get_backward() const {return m_backward;} + const double& get_sweep_tol() const {return m_sweep_tol;} + const int& get_twodot_method() const {return n_twodot_noise;} + const double& get_twodot_noise() const {return m_twodot_noise;} + double& set_twodot_noise() {return m_twodot_noise;} + const double& get_twodot_gamma() const {return m_twodot_gamma;} + double& set_twodot_gamma() {return m_twodot_gamma;} + const bool& get_restart() const {return m_restart;} + const bool& get_restart_warm() const {return m_restart_warm;} + const bool& get_reset_iterations() const {return m_reset_iterations;} + const ninejCoeffs& get_ninej() const {return m_ninej;} + const hamTypes &hamiltonian() const {return m_ham_type;} + const WarmUpTypes &warmup() const {return m_warmup;} + const int &guess_permutations() const { return m_guess_permutations; } + const int &max_lanczos_dimension() const {return m_max_lanczos_dimension;} + std::vector thrds_per_node() const { return m_thrds_per_node; } + const calcType &calc_type() const { return m_calc_type; } + calcType &calc_type() { return m_calc_type; } + const solveTypes &solve_method() const { return m_solve_type; } + const noiseTypes &noise_type() const {return m_noise_type;} + const bool &set_Sz() const {return m_set_Sz;} + const algorithmTypes &algorithm_method() const { return m_algorithm_type; } + algorithmTypes &set_algorithm_method() { return m_algorithm_type; } + int twodot_to_onedot_iter() const { return m_twodot_to_onedot_iter; } + std::vector< std::map >& get_quantaToKeep() { return m_quantaToKeep;} + const std::vector &hf_occupancy() const { return m_hf_occupancy; } + const std::vector &spin_orbs_symmetry() const { return m_spin_orbs_symmetry; } + std::vector weights(int sweep_iter) const;// { return m_weights; } + std::vector weights() const { return m_weights; } + const std::vector &sweep_iter_schedule() const { return m_sweep_iter_schedule; } + const std::vector &sweep_state_schedule() const { return m_sweep_state_schedule; } + const std::vector &sweep_qstate_schedule() const { return m_sweep_qstate_schedule; } + const std::vector &sweep_tol_schedule() const { return m_sweep_tol_schedule; } + const std::vector &sweep_noise_schedule() const { return m_sweep_noise_schedule; } + const std::vector &sweep_additional_noise_schedule() const { return m_sweep_additional_noise_schedule; } + std::vector &set_sweep_noise_schedule() { return m_sweep_noise_schedule; } + std::vector &set_sweep_additional_noise_schedule() { return m_sweep_additional_noise_schedule; } + int& Sz() {return m_Sz;} + int nroots(int sweep_iter) const; + int nroots() const {return m_nroots;} + int real_particle_number() const { return (m_alpha + m_beta);} + int total_particle_number() const { if(!m_add_noninteracting_orbs) return (m_alpha + m_beta); else return (2*m_alpha); } + int bra_particle_number() const { if(!m_add_noninteracting_orbs) return (m_bra_alpha + m_bra_beta); else return (2*m_bra_alpha); } // diff spins case + bool calc_ri_4pdm() const {return m_calc_ri_4pdm;} + bool store_ripdm_readable() const {return m_store_ripdm_readable;} + bool nevpt2() const {return m_nevpt2;} + bool read_higherpdm() const {return m_conventional_nevpt2;} + int kept_nevpt2_states() const {return m_kept_nevpt2_states;} + bool Print() const {return NevPrint.first;} + int PrintIndex() const{return NevPrint.second;} + void SetPrint(bool p, int i=0){NevPrint.first = p;NevPrint.second=i;} + const SpinSpace total_spin_number() const { if (!m_add_noninteracting_orbs) return SpinSpace(m_alpha - m_beta); else return SpinSpace(0); } + const SpinSpace bra_spin_number() const { if (!m_add_noninteracting_orbs) return SpinSpace(m_bra_alpha - m_bra_beta); else return SpinSpace(0); } // diff spins case + int last_site() const + { + if(m_spinAdapted) + { + if(m_calc_type == MPS_NEVPT) + return m_act_size; + else + return m_num_spatial_orbs; + } + else + { + if(m_calc_type == MPS_NEVPT) + return 2*m_act_size; + else + return 2*m_num_spatial_orbs; + } + } + const bool &no_transform() const { return m_no_transform; } + const int &deflation_min_size() const { return m_deflation_min_size; } + const bool &direct() const { return m_direct; } + const int &deflation_max_size() const { return m_deflation_max_size; } + const IrrepSpace &total_symmetry_number() const { return m_total_symmetry_number; } + const IrrepSpace &bra_symmetry_number() const { return m_bra_symmetry_number; } + const SpinQuantum &molecule_quantum() const { return m_molecule_quantum; } + const SpinQuantum &bra_molecule_quantum() const { return m_bra_molecule_quantum; } // diff spins case + const int &sys_add() const { return m_sys_add; } + const bool &add_noninteracting_orbs() const {return m_add_noninteracting_orbs;} + bool &add_noninteracting_orbs() {return m_add_noninteracting_orbs;} + const bool &non_SE() {return m_non_SE;} //for calculations without singlet embedding + const int &nquanta() const { return m_nquanta; } + const int &env_add() const { return m_env_add; } + const bool &do_fci() const { return m_do_fci; } + const int &max_iter() const { return m_maxiter; } + const double &oneindex_screen_tol() const { return m_oneindex_screen_tol; } + double &oneindex_screen_tol() { return m_oneindex_screen_tol; } + const double &twoindex_screen_tol() const { return m_twoindex_screen_tol; } + double &twoindex_screen_tol() { return m_twoindex_screen_tol; } + const int &total_spin() const {return m_total_spin;} + const std::vector &spin_vector() const { return m_spin_vector; } + const std::string &save_prefix() const { return m_save_prefix; } + const std::string &load_prefix() const { return m_load_prefix; } + SpinQuantum& set_molecule_quantum() {return m_molecule_quantum;} + + bool transition_diff_spin() { + return m_transition_diff_spin; + } + + SpinQuantum effective_molecule_quantum() { + if (!m_add_noninteracting_orbs) + return m_molecule_quantum; + else + return SpinQuantum(m_molecule_quantum.particleNumber + m_molecule_quantum.totalSpin.getirrep(), SpinSpace(0), m_molecule_quantum.orbitalSymmetry); + //return SpinQuantum(total_particle_number() + total_spin_number().getirrep(), SpinSpace(0), total_symmetry_number()); + } + +// + SpinQuantum bra_quantum() { + if (!m_add_noninteracting_orbs) + { + if (m_transition_diff_spin) { + return SpinQuantum(bra_particle_number(), SpinSpace(m_bra_alpha - m_bra_beta), bra_symmetry_number()); + } + else + return SpinQuantum(total_particle_number(), SpinSpace(m_alpha - m_beta), bra_symmetry_number()); + } + else { + if (m_transition_diff_spin) { + return SpinQuantum(bra_particle_number()+bra_spin_number().getirrep(), SpinSpace(0), bra_symmetry_number()); + } + else + return SpinQuantum(total_particle_number() + total_spin_number().getirrep(), SpinSpace(0), bra_symmetry_number()); } +} + + + + vector effective_molecule_quantum_vec() { + vector q; + if (!m_Bogoliubov) q.push_back(effective_molecule_quantum()); + else { + SpinQuantum q_max = effective_molecule_quantum(); + for (int n = 0; n <= q_max.get_n(); n+=2) { + q.push_back(SpinQuantum(n, q_max.get_s(), q_max.get_symm())); + } + } + return q; + } + vector bra_quantum_vec() { + vector q; + if (!m_Bogoliubov) q.push_back(bra_quantum()); + else { + SpinQuantum q_max = bra_quantum(); + for (int n = 0; n <= q_max.get_n(); n+=2) { + q.push_back(SpinQuantum(n, q_max.get_s(), q_max.get_symm())); + } + } + return q; + } + bool transition_diff_irrep(){ + return m_transition_diff_spatial_irrep; + } + std::vector& get_orbenergies() {return m_orbenergies;} + int getHFQuanta(const SpinBlock& b) const; + const bool &do_pdm() const {return m_do_pdm;} + bool &do_pdm() {return m_do_pdm;} + const bool &do_npdm_ops() const {return m_do_npdm_ops;} + bool &do_npdm_ops() {return m_do_npdm_ops;} + const bool &do_npdm_in_core() const {return m_do_npdm_in_core;} + bool &do_npdm_in_core() {return m_do_npdm_in_core;} + bool new_npdm_code() const{ + if( m_do_pdm) return m_new_npdm_code; + else return false; + } + void set_new_npdm_code(){ m_new_npdm_code= true;} + const bool &npdm_generate() const { return m_npdm_generate;} + bool &npdm_generate() { return m_npdm_generate;} + const bool &store_spinpdm() const {return m_store_spinpdm;} + bool &store_spinpdm() {return m_store_spinpdm;} + const bool &spatpdm_disk_dump() const {return m_spatpdm_disk_dump;} + bool &spatpdm_disk_dump() {return m_spatpdm_disk_dump;} + const bool &pdm_unsorted() const {return m_pdm_unsorted;} + bool &pdm_unsorted(){return m_pdm_unsorted;} + const std::vector &specificpdm() const {return m_specificpdm;} + std::vector &specificpdm() {return m_specificpdm;} + const bool &store_nonredundant_pdm() const { return m_store_nonredundant_pdm;} + bool &store_nonredundant_pdm() { return m_store_nonredundant_pdm;} + int slater_size() const {return m_norbs;} + const std::vector &reorder_vector() {return m_reorder;} + const int &act_size() const { return m_act_size;} + const int &core_size() const { return m_core_size;} + const int &virt_size() const { return m_virt_size;} + const int &total_size() const { return m_total_orbs;} + bool spinAdapted() {return m_spinAdapted;} + const int &nevpt_state_num() const {return m_nevpt_state_num;} + bool &npdm_intermediate() { return m_npdm_intermediate; } + 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/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.C.new b/modules/npdm/npdm.C.new new file mode 100644 index 00000000..9eec765e --- /dev/null +++ b/modules/npdm/npdm.C.new @@ -0,0 +1,598 @@ +/* +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 +*/ + +#include "npdm.h" +#include "sweep.h" +#include "sweepgenblock.h" +#include "density.h" +#include "sweeponepdm.h" // For legacy version of 1pdm +#include "sweeptwopdm.h" // For legacy version of 2pdm +#include "npdm_driver.h" +#include "nevpt2_npdm_driver.h" +#include "pario.h" +#include "ds1_sweeponepdm.h" //EL +#include "ds0_sweeponepdm.h" //EL + + +void dmrg(double sweep_tol); +void restart(double sweep_tol, bool reset_iter); +void ReadInput(char* conf); +void fullrestartGenblock(); + +namespace SpinAdapted { +namespace Npdm { + +//----------------------------------------------------------------------------------------------------------------------------------------------------------- + +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; + //mcheck("at the start of block and decimate"); + // figure out if we are going forward or backwards + dmrginp.guessgenT -> start(); + bool forward = (system.get_sites() [0] == 0); + SpinBlock systemDot; + SpinBlock envDot; + int systemDotStart, systemDotEnd; + int systemDotSize = sweepParams.get_sys_add() - 1; + if (forward) + { + systemDotStart = dmrginp.spinAdapted() ? *system.get_sites().rbegin () + 1 : (*system.get_sites().rbegin ())/2 + 1 ; + systemDotEnd = systemDotStart + systemDotSize; + } + else + { + systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ; + systemDotEnd = systemDotStart - systemDotSize; + } + vector spindotsites(2); + spindotsites[0] = systemDotStart; + spindotsites[1] = systemDotEnd; + //if (useSlater) { + systemDot = SpinBlock(systemDotStart, systemDotEnd, system.get_integralIndex(), true); + //SpinBlock::store(true, systemDot.get_sites(), systemDot); + //} + //else + //SpinBlock::restore(true, spindotsites, systemDot); + SpinBlock environment, environmentDot, newEnvironment; + + int environmentDotStart, environmentDotEnd, environmentStart, environmentEnd; + + const int nexact = forward ? sweepParams.get_forward_starting_size() : sweepParams.get_backward_starting_size(); + + system.addAdditionalCompOps(); + 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(), + sweepParams.get_onedot(), nexact, useSlater, environment.get_integralIndex(), true, true, true); + } + 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(), + sweepParams.get_onedot(), nexact, useSlater, environment.get_integralIndex(), true, true, true); + + } + SpinBlock big; + newSystem.set_loopblock(true); + system.set_loopblock(false); + newEnvironment.set_loopblock(false); + 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); + + } + 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); + } + +#ifndef SERIAL + mpi::communicator world; + mpi::broadcast(world, solution, 0); +#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); + + + + //bra and ket rotation matrices are calculated from different density matrices. + + + std::vector rotateMatrix; + std::vector rotateMatrixB; + + if(state!=stateB){ + + DensityMatrix tracedMatrix(newSystem.get_braStateInfo()); + tracedMatrix.allocate(newSystem.get_braStateInfo()); + tracedMatrix.makedensitymatrix(std::vector(1,solution[0]), big, std::vector(1,1.0), 0.0, 0.0, false); + rotateMatrix.clear(); + + DensityMatrix tracedMatrixB(newSystem.get_ketStateInfo()); + tracedMatrixB.allocate(newSystem.get_ketStateInfo()); + tracedMatrixB.makedensitymatrix(std::vector(1,solution[1]), big, std::vector(1,1.0), 0.0, 0.0, false); + rotateMatrixB.clear(); + if (!mpigetrank()){ + double error = makeRotateMatrix(tracedMatrix, rotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates()); + error = makeRotateMatrix(tracedMatrixB, rotateMatrixB, sweepParams.get_keep_states(), sweepParams.get_keep_qstates()); + } + } + else{ + DensityMatrix tracedMatrix(newSystem.get_stateInfo()); + tracedMatrix.allocate(newSystem.get_stateInfo()); + tracedMatrix.makedensitymatrix(std::vector(1,solution[0]), big, std::vector(1,1.0), 0.0, 0.0, false); + rotateMatrix.clear(); + if (!mpigetrank()){ + double error = makeRotateMatrix(tracedMatrix, rotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates()); + } + + } + + + + + + int sweepPos = sweepParams.get_block_iter(); + int endPos = sweepParams.get_n_iters()-1; + npdm_driver.compute_npdm_elements(solution, big, sweepPos, endPos); + SaveRotationMatrix (newSystem.get_sites(), rotateMatrix, state); + solution[0].SaveWavefunctionInfo (big.get_braStateInfo(), big.get_leftBlock()->get_sites(), state); + + if(state!=stateB){ + SaveRotationMatrix (newSystem.get_sites(), rotateMatrixB, stateB); + solution[1].SaveWavefunctionInfo (big.get_ketStateInfo(), big.get_leftBlock()->get_sites(), stateB); + } + + + //FIXME + //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(newSystem.get_sites().size()>1) + //if (!mpigetrank()){ + //LoadRotationMatrix (newSystem.get_sites(), rotateMatrix, state); + //LoadRotationMatrix (newSystem.get_sites(), rotateMatrixB, stateB); + //} + #ifndef SERIAL + mpi::broadcast(world,rotateMatrix,0); + if(state!=stateB) + mpi::broadcast(world,rotateMatrixB,0); + #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. + //Since it is useless, Just omit it at the last sweep site. + // if( sweepParams.get_block_iter() != sweepParams.get_n_iters() - 1) + { + if(state!=stateB) + newSystem.transform_operators(rotateMatrix,rotateMatrixB); + else + newSystem.transform_operators(rotateMatrix); + } + + //newSystem.transform_operators(rotateMatrix,rotateMatrixB); + ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime(); + p3out << "NPDM block and decimate and compute elements " << ewall << " " << ecpu << endl; + +} + +//----------------------------------------------------------------------------------------------------------------------------------------------------------- + +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; + pout.precision(12); + SpinBlock system; + const int nroots = dmrginp.nroots(); + std::vector finalEnergy(nroots,0.); + std::vector finalEnergy_spins(nroots,0.); + double finalError = 0.; + + sweepParams.set_sweep_parameters(); + // 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 + 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; + + if (!restart) sweepParams.set_block_iter() = 0; + if(dmrginp.setStateSpecific() || dmrginp.transition_diff_irrep()){ + if (!restart) SpinBlock::store (forward, system.get_sites(), system, state, stateB ); // if restart, just restoring an existing block -- + } + else{ + if (!restart) SpinBlock::store (forward, system.get_sites(), system, sweepParams.current_root(), sweepParams.current_root() ); // if restart, just restoring an existing block -- + } + + sweepParams.savestate(forward, system.get_sites().size()); + bool dot_with_sys = true; + + // Loop over all block sites + for (; sweepParams.get_block_iter() < sweepParams.get_n_iters(); ) { + Timer timer; + + pout << "\n\t\t\t Block Iteration :: " << sweepParams.get_block_iter() << endl; + pout << "\t\t\t ----------------------------" << endl; + if (forward) { p1out << "\t\t\t Current direction is :: Forwards " << endl; } + else { p1out << "\t\t\t Current direction is :: Backwards " << endl; } + + //if (SHOW_MORE) pout << "system block" << endl << system << endl; + + 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 && + ((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()< npdm_driver; + //if ( (dmrginp.hamiltonian() == QUANTUM_CHEMISTRY) && dmrginp.spinAdapted() ) { + if ( (dmrginp.hamiltonian() == QUANTUM_CHEMISTRY) || (dmrginp.hamiltonian() == BCS)) { + //By default, new_npdm_code is false. + //For npdm_order 1 or 2. new_npdm_code is determined by default or manual setting. + //For the other situation, only old or new code is suitable. + if(npdm_order == NPDM_PAIRMATRIX || npdm_order == NPDM_THREEPDM || npdm_order == NPDM_FOURPDM || npdm_order == NPDM_NEVPT2 || transitionpdm == true || dmrginp.spinAdapted() == false || dmrginp.setStateSpecific()) + dmrginp.set_new_npdm_code(); + + if(dmrginp.new_npdm_code()){ + if (npdm_order == NPDM_ONEPDM) npdm_driver = boost::shared_ptr( new Onepdm_driver( dmrginp.last_site() ) ); + else if ((npdm_order == NPDM_DS0)||(npdm_order == NPDM_DS1)) npdm_driver = boost::shared_ptr( new Onepdm_driver( dmrginp.last_site() ) ); + else if (npdm_order == NPDM_TWOPDM) npdm_driver = boost::shared_ptr( new Twopdm_driver( dmrginp.last_site() ) ); + else if (npdm_order == NPDM_THREEPDM) npdm_driver = boost::shared_ptr( new Threepdm_driver( dmrginp.last_site() ) ); + else if (npdm_order == NPDM_FOURPDM) npdm_driver = boost::shared_ptr( new Fourpdm_driver( dmrginp.last_site() ) ); + else if (npdm_order == NPDM_NEVPT2) npdm_driver = boost::shared_ptr( new Nevpt2_npdm_driver( dmrginp.last_site() ) ); + else if (npdm_order == NPDM_PAIRMATRIX) npdm_driver = boost::shared_ptr( new Pairpdm_driver( dmrginp.last_site() ) ); + else abort(); + } + } + +if (dS) { + for (int state=0; stateclear(); + sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; + if (npdm_order == NPDM_DS1) ds1_onepdm::do_one(sweepParams, false, !direction, false, 0, state, stateB); + else if (npdm_order == NPDM_DS0) ds0_onepdm::do_one(sweepParams, false, !direction, false, 0, state, stateB); + else abort(); + } + } +} +else { + if (dmrginp.specificpdm().size()!=0) + { + Timer timer; + dmrginp.set_fullrestart() = true; + sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; + dmrginp.npdm_generate() = true; + if ( !dmrginp.setStateSpecific()) + SweepGenblock::do_one(sweepParams, false, !direction, false, 0, -1, -1); //this will generate the cd operators + else if (dmrginp.specificpdm().size()==1) + SweepGenblock::do_one(sweepParams, false, !direction, false, 0, dmrginp.specificpdm()[0], dmrginp.specificpdm()[0]); //this will generate the cd operators + else if (dmrginp.specificpdm().size()==2) + SweepGenblock::do_one(sweepParams, false, !direction, false, 0, dmrginp.specificpdm()[0], dmrginp.specificpdm()[1]); //this will generate the cd operators + else abort(); + dmrginp.npdm_generate() = false; + ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime(); + p3out << "\t\t\t NPDM SweepGenblock time " << ewall << " " << ecpu << endl; + dmrginp.set_fullrestart() = false; + + sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; + Timer timerX; + npdm_driver->clear(); + if (dmrginp.specificpdm().size()==1) + npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0,dmrginp.specificpdm()[0],dmrginp.specificpdm()[0]); + else if (dmrginp.specificpdm().size()==2) + npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0,dmrginp.specificpdm()[0],dmrginp.specificpdm()[1]); + else abort(); + ecpu = timerX.elapsedcputime();ewall=timerX.elapsedwalltime(); + p3out << "\t\t\t NPDM sweep time " << ewall << " " << ecpu << endl; + return; + } + + if(dmrginp.transition_diff_irrep()){ + // It is used when bra and ket has different spatial irrep + // For now, only the transtion pdm between two wavefuntions( 1 as bra and 0 as ket) are calculation + // If the spatial irrep information is stored in wavefuntions, transition pdm among several wavefunctions( i for bra and j for ket, there are n(n-1)/2 kinds of situations.) are possible. + for (int state=0; stateclear(); + 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 if( !dmrginp.setStateSpecific()){ + Timer timer; + dmrginp.set_fullrestart() = true; + dmrginp.npdm_generate() = true; + sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; + SweepGenblock::do_one(sweepParams, false, !direction, false, 0, -1, -1); //this will generate the cd operators + dmrginp.npdm_generate() = false; + ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime(); + p3out << "\t\t\t NPDM SweepGenblock time " << ewall << " " << ecpu << endl; + dmrginp.set_fullrestart() = false; + + +//pout << "(orz3rdmalt1,orz3rdmalt2) = (" << dmrginp.orz3rdmalt1() << "," << dmrginp.orz3rdmalt2() << ")" << endl; + if(transitionpdm){ + // <\Phi_k|a^+_ia_j|\Phi_l> = <\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; + npdm_driver->clear(); + 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); + else if (npdm_order == NPDM_TWOPDM) SweepTwopdm::do_one(sweepParams, false, direction, false, 0, state); + else abort(); + } + } +} + } + + } + + else { + // state-specific + if(transitionpdm){ + for (int state=0; stateclear(); + sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; + 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{ + for (int state=0; stateclear(); + npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0, state,state); + } + } + +} + + sweep_copy.savestate(direction_copy, restartsize_copy); + +} +} +} +} diff --git a/modules/npdm/npdm.h.new b/modules/npdm/npdm.h.new new file mode 100644 index 00000000..c7e43ec0 --- /dev/null +++ b/modules/npdm/npdm.h.new @@ -0,0 +1,21 @@ +/* +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 +*/ + +#ifndef NPDM_HEADER +#define NPDM_HEADER + +namespace SpinAdapted{ +enum NpdmOrder{NPDM_NEVPT2, NPDM_ONEPDM, NPDM_TWOPDM, NPDM_THREEPDM, NPDM_FOURPDM, NPDM_PAIRMATRIX, NPDM_OVERLAP, NPDM_EMPTY, NPDM_DS0, NPDM_DS1}; +namespace Npdm{ + + void npdm(NpdmOrder npdm_order, bool restartpdm=false, bool transitionpdm=false, bool dS = false); + +} +} + +#endif 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"; } } From 168b8b05882edda04a5fdad717f13c643176f309 Mon Sep 17 00:00:00 2001 From: Takeshi Yanai Date: Thu, 26 Sep 2019 10:14:52 +0900 Subject: [PATCH 6/8] recommit --- input.C.new | 2860 --------------------------------------- input.h.new | 583 -------- modules/npdm/npdm.C.new | 598 -------- modules/npdm/npdm.h.new | 21 - 4 files changed, 4062 deletions(-) delete mode 100644 input.C.new delete mode 100644 input.h.new delete mode 100644 modules/npdm/npdm.C.new delete mode 100644 modules/npdm/npdm.h.new diff --git a/input.C.new b/input.C.new deleted file mode 100644 index 055e093e..00000000 --- a/input.C.new +++ /dev/null @@ -1,2860 +0,0 @@ -/* -Developed by Sandeep Sharma, Roberto Olivares-Amaya 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, Roberto Olivares-Amaya and Garnet K.-L. Chan -*/ - - -#include -#include -#include -#include "Symmetry.h" -#include "global.h" -#include "MatrixBLAS.h" -#include "spinblock.h" -#include "couplingCoeffs.h" -#include "genetic/GAOptimize.h" -#include "genetic/ReadIntegral.h" -#include -#include -#include -#include -#ifndef SERIAL -#include -#endif -#include -#include "fiedler.h" -#include "pario.h" - -using namespace std; - -namespace SpinAdapted { -string sym; -bool NonabelianSym; -} -void CheckFileExistence(string filename, string filetype); -void CheckFileInexistence(string filename, string filetype); - -void SpinAdapted::Input::ReadMeaningfulLine(ifstream& input, string& msg, int msgsize) -{ - bool readmore = true; - while (readmore && !input.eof()) { - char msgctr[msgsize]; - input.getline(msgctr, msgsize+1); - - msg=string(msgctr); - if(msg.size() == msgsize) { - pout << "in the process of reading line beginning with "<(world.size(),1); -#else - m_thrds_per_node = vector(1,1); -#endif - m_ham_type = QUANTUM_CHEMISTRY; - m_algorithm_type = TWODOT_TO_ONEDOT; - m_noise_type = RANDOM; - m_calc_type = DMRG; - m_solve_type = DAVIDSON; - m_stateSpecific = false; - m_implicitTranspose = true; //dont make DD just use CC^T to evaluate it - m_occupied_orbitals = -1; - m_num_Integrals = 1; - v_2.resize(1, TwoElectronArray(TwoElectronArray::restrictedNonPermSymm)); - v_1.resize(1); - coreEnergy.resize(1); - m_baseState.resize(1,0); - m_projectorState.resize(0); - m_targetState = -1; - m_guessState = 1; - m_permSymm = 2; - - m_openorbs.resize(0); - m_closedorbs.resize(0); - - m_spinAdapted = true; - m_Bogoliubov = false; - m_sys_add = 1; - m_env_add = 1; - - m_twodot_to_onedot_iter = 0; - m_integral_disk_storage_thresh = 1000; - m_max_lanczos_dimension = 5000; - - m_norbs = 0; - m_alpha = 0; - m_beta = 0; - m_hf_occ_user = ""; - - m_outputlevel = 0; - m_nquanta = 2; - m_total_symmetry_number = IrrepSpace( 0 ); - m_total_spin = 0; - m_guess_permutations = 10; - - m_direct = true; - m_nroots = 1; - m_weights.resize(m_nroots); - m_weights[0] = 1.; - - m_deflation_min_size = 2; - m_deflation_max_size = 20; - - m_transition_diff_spatial_irrep = false; - m_add_noninteracting_orbs = true; - m_bra_alpha = 0; - m_bra_beta =0; - m_non_SE = false; - m_no_transform = false; - m_do_fci = false; - m_do_npdm_ops = false; - m_do_npdm_in_core = false; - m_npdm_generate = false; - m_new_npdm_code = false; - m_do_pdm = false; - m_store_spinpdm = false; - m_spatpdm_disk_dump = false; - m_store_nonredundant_pdm =false; - m_pdm_unsorted = false; - - - m_nevpt_state_num = 0; - m_npdm_intermediate= true; - m_npdm_multinode= true; - - m_maxiter = 10; - m_oneindex_screen_tol = NUMERICAL_ZERO; - m_twoindex_screen_tol = NUMERICAL_ZERO; - - m_load_prefix = "."; - m_save_prefix = "."; - - m_maxj = 15; - m_ninej.init(m_maxj); - m_set_Sz = false; - - n_twodot_noise = 0; - m_twodot_noise = 1.0e-4; - m_twodot_gamma = 3.0e-1; - - m_sweep_tol = 1.0e-5; - m_restart = false; - m_fullrestart = false; - m_restart_warm = false; - m_backward = false; - m_reset_iterations = false; - - m_do_diis = false; - m_diis_error = 1e-2; - m_start_diis_iter = 8; - m_diis_keep_states = 6; - m_diis_error_tol = 1e-8; - m_num_spatial_orbs = 0; - m_schedule_type_default = false; - m_schedule_type_backward = false; - m_maxM = 0; - m_lastM = 500; - m_startM = 250; - - m_calc_ri_4pdm=false; - m_store_ripdm_readable=false; - m_nevpt2 = false; - m_conventional_nevpt2 = false; - m_kept_nevpt2_states = -1; - NevPrint.first=false; - NevPrint.second=0; - - //reorder options, by default it does fiedler - m_reorderType = FIEDLER; - m_reorderfile = ""; - m_gaconffile = "default"; - - m_orbformat=MOLPROFORM; - - m_warmup = LOCAL0; - - m_orz3rdmalt1 = -1; - m_orz3rdmalt2 = -1; -} - -void SpinAdapted::Input::usedkey_error(string& key, string& line) { - pout << "keyword "< usedkey(NUMKEYWORDS, -1); - int n_elec = -1; - int n_spin = -1; - int n_bra_spin =-1; //EL - sym = "c1"; - NonabelianSym = false; - std::vector orbitalfile; - - initialize_defaults(); - - if(mpigetrank() == 0) - { - pout << "Reading input file"< tok; - boost::split(tok, msg, is_any_of(", \t"), token_compress_on); - string keyword = *tok.begin(); - - if (boost::iequals(keyword, "orbs") || boost::iequals(keyword, "orbitals") || boost::iequals(keyword, " orbitals")) - { - if(usedkey[ORBS] == 0) - usedkey_error(keyword, msg); - usedkey[ORBS] = 0; - if (tok.size() < 2) { - pout << "keyword orbs should be followed by atleast a single filename and then an end line"< schd_tok; - boost::split(schd_tok, msg, is_any_of(" \t"), token_compress_on); - if (tok.size() != 1) { - if (boost::iequals(tok[1], "default")) { - m_schedule_type_default = true; - continue; - } - } - - m_sweep_iter_schedule.resize(0); - m_sweep_state_schedule.resize(0); - m_sweep_qstate_schedule.resize(0); - m_sweep_tol_schedule.resize(0); - m_sweep_noise_schedule.resize(0); - m_sweep_additional_noise_schedule.resize(0); - - int i = 0; - while(!boost::iequals(schd_tok[0], "END")) - { - - if (schd_tok.size() != 4) { - pout << "Each line of the schedule contain four entries sweep_iteration #retained states Davidson tolerance noise"<0 && m_sweep_iter_schedule[i] <= m_sweep_iter_schedule[i-1]) { - pout << "Sweep iteration at a given line should be higher than the previous sweep iteration"< and required for SOC and g-tensor calculations: - * preliminary calculations for the states should be performed without singlet embedding, otherwise the triplet excitation matrix elements will be equal to 0 -If 2 spins are given, the calculations of transition density matrix between wavefunctions with different spins (w/ or w/o different irreps) will be performed */ - { - n_bra_spin = atoi(tok[1].c_str()); - n_spin = atoi(tok[2].c_str()); - m_transition_diff_spin=true; - m_add_noninteracting_orbs = false; - m_bra_alpha = (n_elec + n_bra_spin)/2; - m_bra_beta = (n_elec - n_bra_spin)/2; - } - if ((tok.size() != 2)&&((tok.size() != 3))) { - pout << "keyword spin should be followed by a single number (or two spin numbers for SOC pdms) and then an end line"<(i)]= OneElectronArray(); -// } - for(int i= 0; i< TwoPertEnd; i++){ - //vpt_2[static_cast(i)]= TwoElectronArray(TwoElectronArray::restrictedNonPermSymm); - vpt_2[static_cast(i)]= PerturbTwoElectronArray(); - // vpt_2[i].rhf = true; - } - } - else if (boost::iequals(keyword, "restart_mps_nevpt")) { - m_calc_type = RESTART_MPS_NEVPT; - if (tok.size() !=4){ - pout << "keyword mps_nevpt should be followed by three numbers and then an end line"<(i)]= OneElectronArray(); -// } - for(int i= 0; i< TwoPertEnd; i++){ - //vpt_2[static_cast(i)]= TwoElectronArray(TwoElectronArray::restrictedNonPermSymm); - vpt_2[static_cast(i)]= PerturbTwoElectronArray(); - } - } - else if (boost::iequals(keyword, "nevpt_state_num" )) - { - if (tok.size() != 2) { - pout << "keyword nevpt_state_num should be followed by a single integer than then an end line."< pdms for SOC - else if (boost::iequals(keyword, "ds1_onepdm") || boost::iequals(keyword, "ds1_onerdm")) - m_calc_type = DS1_ONEPDM; - else if (boost::iequals(keyword, "restart_ds1_onepdm") || boost::iequals(keyword, "restart_ds1_onerdm")) - m_calc_type = RESTART_DS1_ONEPDM; -// Elvira: these are to calculate pdms to calculate S in the code for g-tensors - else if (boost::iequals(keyword, "ds0_onepdm") || boost::iequals(keyword, "ds0_onerdm")) - m_calc_type = DS0_ONEPDM; - else if (boost::iequals(keyword, "restart_ds0_onepdm") || boost::iequals(keyword, "restart_ds0_onerdm")) - m_calc_type = RESTART_DS0_ONEPDM; -// - else if (boost::iequals(keyword, "restart_tran_threepdm") || boost::iequals(keyword, "restart_tran_threerdm") || boost::iequals(keyword, "restart_tran_threerdm")) - { - m_calc_type = RESTART_T_THREEPDM; - m_restart = true; - } - else if (boost::iequals(keyword, "nevpt2") || boost::iequals(keyword, "pt2")){ - m_calc_type = NEVPT2; - m_nevpt2 = true; - m_transition_diff_spatial_irrep=false; - } - else if (boost::iequals(keyword, "ripdm")){ - m_calc_type = NEVPT2; - m_nevpt2 = false; - m_transition_diff_spatial_irrep=false; - } - else if (boost::iequals(keyword, "restart_nevpt2") || boost::iequals(keyword, "restart_pt2")){ - m_calc_type = RESTART_NEVPT2; - m_nevpt2 = true; - m_transition_diff_spatial_irrep=false; - m_restart = true; - } - else if (boost::iequals(keyword, "restart_ripdm")){ - m_calc_type = RESTART_NEVPT2; - m_nevpt2 = false; - m_transition_diff_spatial_irrep=false; - m_restart = true; - } - else if (boost::iequals(keyword, "calc_ri4pdm") || boost::iequals(keyword, "ri4pdm")) - { - m_calc_ri_4pdm = true; - } - else if (boost::iequals(keyword, "ripdm_readable")) - { - m_store_ripdm_readable = true; - } - else if (boost::iequals(keyword, "conventional_nevpt2")) - { - m_conventional_nevpt2 = true; - } - else if (boost::iequals(keyword, "M_nevpt2")){ - if (tok.size() != 2) { - pout << "keyword M_nevpt2 should be followed by a single float and then an end line."< :: iterator it = ++tok.begin(); - // the occupancies start from the second element of the tok string - if (tok.size() == 2 ) { - m_hf_occ_user = (*it).c_str(); - } - else{ - m_hf_occ_user = "manual"; - for( ; it != tok.end(); ++it ){ - int occ_i = atoi( (*it).c_str() ); - m_hf_occupancy.push_back( occ_i ); - } - } - pout << "m_hf_occ_user " << m_hf_occ_user << endl; - - } - - else if(boost::iequals(keyword, "open")) - { - vector :: iterator it = ++tok.begin(); - m_openorbs.resize(tok.size()-1,-1); - - int ii = 0; - for( ; it != tok.end(); ++it ){ - int openorb = atoi( (*it).c_str() ); - m_openorbs[ii] = openorb-1 ; - //pout << openorb<<" "; - ii++; - } - //pout << endl; - - } - else if(boost::iequals(keyword, "closed")) - { - vector :: iterator it = ++tok.begin(); - m_closedorbs.resize(tok.size()-1,-1); - //pout << "closed orbs :"; - int ii=0; - for( ; it != tok.end(); ++it ){ - int closedorb = atoi( (*it).c_str() ); - m_closedorbs[ii] = closedorb -1; - //pout << closedorb<<" "; - ii++; - } - //pout << endl; - } - else if(boost::iequals(keyword, "nroots")) - { - if(usedkey[NROOTS] == 0) - usedkey_error(keyword, msg); - usedkey[NROOTS] = 0; - std::string nroots_str; - if (tok.size() != 2) { - pout << "keyword nroots should be followed by a single integer and then an end line."< weighttoken; - boost::split(weighttoken, msg, is_any_of(" \t"), token_compress_on); - - if (!boost::iequals(weighttoken[0], "weights")) { - //manually make all the weights equal - m_weights.resize(m_nroots); - for (int i=0; i1 ) { - pout << "Currently the response code does not work with non-particle number conserving hamiltonians!!"; - abort(); - } - -//if (sym != "c1") // must be initialized even if c1 sym. - Symmetry::InitialiseTable(sym); - - if (mpigetrank() == 0) { - for (int l=0; l& oldtonew) { - string msg; int msgsize = 5000; - ReadMeaningfulLine(dumpFile, msg, msgsize); - vector tok; - std::vector diff; - boost::split(tok, msg, is_any_of(", \t"), token_compress_on); - while(msg.size() != 0) { - for (int i=0; im_norbs || oldtonew.back() < 0) { - pout << "Illegal orbital index "< tok; - int offset = m_orbformat == DMRGFORM ? 0 : 1; - std::ofstream ReorderFileOutput; - std::ifstream ReorderFileInput; - char ReorderFileName[5000]; - std::vector reorder; - - if (mpigetrank() == 0) { - ReadMeaningfulLine(dumpFile, msg, msgsize); - boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); - - if(offset != 0) - m_norbs = 2*atoi(tok[2].c_str()); - else - m_norbs = 2*atoi(tok[1].c_str()); - - m_num_spatial_orbs = 0; - m_spin_orbs_symmetry.resize(m_norbs); - m_spin_to_spatial.resize(m_norbs); - - //this is the file to which the reordering is written - sprintf(ReorderFileName, "%s%s", save_prefix().c_str(), "/RestartReorder.dat"); - } - boost::filesystem::path p(ReorderFileName); - -#ifndef SERIAL - boost::mpi::communicator world; - mpi::broadcast(world,m_restart,0); - mpi::broadcast(world,m_fullrestart,0); -#endif - - //do the reordering only if it is not a restart calculation - //if it is then just read the reorder.dat from the scratch space - if (integralIndex == 0) { - if(get_restart() || get_fullrestart() || m_calc_type == COMPRESS || m_calc_type == RESPONSE || m_calc_type == RESPONSEBW) { - if (mpigetrank() == 0) { - ReorderFileInput.open(ReorderFileName); - boost::filesystem::path ReorderFilePath(ReorderFileName); - - if(!boost::filesystem::exists(ReorderFilePath)) { - pout << "==============="<> m_reorder[i]; - ReorderFileInput.close(); - } - } - } - else { - if (mpigetrank() == 0) { - ReorderFileOutput.open(ReorderFileName); - } - - - // read the reorder file or calculate the reordering using one of the many options - if (m_reorderType == FIEDLER) { - - if (mpigetrank() == 0){ - m_reorder=get_fiedler(orbitalfile); - pout << "Fiedler-vector orbital ordering: "; - } - } - else if (m_reorderType == GAOPT) { - - ifstream gaconfFile; - - if (mpigetrank() == 0) { - if(m_gaconffile != "default") - gaconfFile.open(m_gaconffile.c_str(), ios::in); - m_reorder = get_fiedler(orbitalfile); - } - //to provide as initial guess to gaopt - m_reorder = getgaorder(gaconfFile, orbitalfile, m_reorder); - pout << "Genetic algorithm orbital ordering: "; - } - - else if (m_reorderType == MANUAL) { - if (mpigetrank() == 0) { - ifstream reorderFile(m_reorderfile.c_str()); - CheckFileExistence(m_reorderfile, "Reorder file "); - readreorderfile(reorderFile, m_reorder); - pout << "Manually provided orbital ordering: "; - } - } - else { //this is the no-reorder case - if (mpigetrank() == 0) { - m_reorder.resize(m_norbs/2); - for (int i=0; i O_unreordered(2,3) (the former is stored internally and the latter is what is given during input) - // but for this m_reorder the reorder vector below would be 3 1 2 4 and O_unreordered(1,2) -> O_reorder(3, 1) - if (mpigetrank() == 0) { - reorder.resize(m_norbs/2); - for (int i=0; i= 0 ) - ir = atoi(tok[i].c_str()) - offset; - else if (atoi(tok[i].c_str()) < -1) - ir = atoi(tok[i].c_str()) + offset; - - if (sym == "trans") ir += 1; //for translational symmetry the lowest irrep is 0 - if (sym == "lzsym") ir = atoi(tok[i].c_str()); - - - m_spin_orbs_symmetry[2*reorderOrbInd] = ir; - m_spin_orbs_symmetry[2*reorderOrbInd+1] = ir; - - if (readLine == numRead) { - m_num_spatial_orbs++; - m_spatial_to_spin.push_back(orbindex); - numRead = Symmetry::sizeofIrrep(ir); - readLine = 0; - } - m_spin_to_spatial[orbindex] = m_num_spatial_orbs-1; - m_spin_to_spatial[orbindex+1] = m_num_spatial_orbs-1; - orbindex +=2; - readLine++; - - - } - msg.resize(0); - ReadMeaningfulLine(dumpFile, msg, msgsize); - boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); - if(boost::iequals(tok[0], "IUHF")) RHF=false; - } - - - if(sym == "dinfh" ) { - m_spatial_to_spin.clear(); - for (int i=0; i= m_integral_disk_storage_thresh) // - { - for (int i=0; i orb; - for (int j=m_spatial_to_spin[i]; j& vpt2, double& coreEnergy) -{ - //TODO - //Reorder is not supported. - ifstream dumpFile; - dumpFile.open(orbitalfile.c_str(), ios::in); - - string msg; int msgsize = 5000; - vector tok; - int offset = m_orbformat == DMRGFORM ? 0 : 1; - std::ofstream ReorderFileOutput; - std::ifstream ReorderFileInput; - char ReorderFileName[5000]; - std::vector reorder(m_total_orbs); - - if (mpigetrank() == 0) { - ReadMeaningfulLine(dumpFile, msg, msgsize); - boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); - - m_norbs =m_act_size*2; - m_num_spatial_orbs = m_act_size; - m_spin_orbs_symmetry.resize(m_total_orbs*2); - m_spin_to_spatial.resize(m_total_orbs*2); - - m_total_spin_orbs_symmetry.resize(m_total_orbs*2); - m_total_spin_to_spatial.resize(m_total_orbs*2); - - //this is the file to which the reordering is written - sprintf(ReorderFileName, "%s%s", save_prefix().c_str(), "/RestartReorder.dat"); - } - - -#ifndef SERIAL - boost::mpi::communicator world; - mpi::broadcast(world,m_restart,0); - mpi::broadcast(world,m_fullrestart,0); -#endif - - //do the reordering only if it is not a restart calculation - //if it is then just read the reorder.dat from the scratch space - if(get_restart() || get_fullrestart()) { - if (mpigetrank() == 0) { - ReorderFileInput.open(ReorderFileName); - boost::filesystem::path ReorderFilePath(ReorderFileName); - - if(!boost::filesystem::exists(ReorderFilePath)) { - pout << "==============="<> m_reorder[i]; - ReorderFileInput.close(); - } - } - } - else { - if (mpigetrank() == 0) { - ReorderFileOutput.open(ReorderFileName); - } - - - // read the reorder file or calculate the reordering using one of the many options - - if (m_reorderType == FIEDLER) { - - if (mpigetrank() == 0){ - m_reorder=get_fiedler_nevpt(orbitalfile, m_act_size); - pout << "Fiedler-vector orbital ordering: "; - } - } - else if (m_reorderType == MANUAL) { - if (mpigetrank() == 0) { - ifstream reorderFile(m_reorderfile.c_str()); - CheckFileExistence(m_reorderfile, "Reorder file "); - readreorderfile(reorderFile, m_reorder); - pout << "Manually provided orbital ordering: "; - } - } - else { //this is the no-reorder case - if (mpigetrank() == 0) { - m_reorder.resize(m_act_size); - for (int i=0; i O_unreordered(2,3) (the former is stored internally and the latter is what is given during input) - // but for this m_reorder the reorder vector below would be 3 1 2 4 and O_unreordered(1,2) -> O_reorder(3, 1) - if (mpigetrank() == 0) { - reorder.resize(m_total_orbs); - for (int i=0; i= 0 ) - ir = atoi(tok[i].c_str()) - offset; - else if (atoi(tok[i].c_str()) < -1) - ir = atoi(tok[i].c_str()) + offset; - - if (sym == "trans") ir += 1; //for translational symmetry the lowest irrep is 0 - if (sym == "lzsym") ir = atoi(tok[i].c_str()); - - - m_total_spin_orbs_symmetry[2*reorderOrbInd] = ir; - m_total_spin_orbs_symmetry[2*reorderOrbInd+1] = ir; - - if (readLine == numRead) { - m_total_spatial_to_spin.push_back(orbindex); - numRead = Symmetry::sizeofIrrep(ir); - readLine = 0; - } - m_total_spin_to_spatial[orbindex] = orbindex/2; - m_total_spin_to_spatial[orbindex+1] = orbindex/2; - orbindex +=2; - readLine++; - - - } - msg.resize(0); - ReadMeaningfulLine(dumpFile, msg, msgsize); - boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); - if(boost::iequals(tok[0], "IUHF")) RHF=false; - } - - - if(sym == "dinfh" ) { - m_total_spatial_to_spin.clear(); - for (int i=0; i(m_total_spin_to_spatial.begin(),m_total_spin_to_spatial.begin()+m_act_size*2); -// m_spatial_to_spin=std::vector(m_total_spatial_to_spin.begin(),m_total_spatial_to_spin.begin()+m_act_size); -// m_spin_orbs_symmetry=std::vector(m_total_spin_orbs_symmetry.begin(),m_total_spin_orbs_symmetry.begin()+m_act_size*2); -// -// m_spatial_to_spin.push_back(m_act_size*2); -// m_spin_to_spatial.push_back(m_act_size*2); - - m_spin_to_spatial=m_total_spin_to_spatial; - m_spatial_to_spin=m_total_spatial_to_spin; - m_spin_orbs_symmetry=m_total_spin_orbs_symmetry; - - - - while (!((boost::iequals(tok[0], "&END")) || (boost::iequals(tok[0], "/")))) { - int temp; - - msg.resize(0); - ReadMeaningfulLine(dumpFile, msg, msgsize); - boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); - if(boost::iequals(tok[0], "IUHF")) RHF=false; - } - - int AOrbOffset = 0, BOrbOffset = 0; - if(!RHF) { - v1.rhf = false; - v2.rhf = false; - //for(auto i: vpt1) - // i.second.rhf = false; - vpt1.rhf = false; -// for(auto i: vpt2) -// i.second.rhf = false; -// fock.rhf = false; - } - v1.ReSize(m_total_orbs*2); - v2.ReSize(m_act_size*2); - vpt1.ReSize(m_total_orbs*2); - //for(auto i: vpt1) - // i.second.ReSize(m_total_orbs*2); - vpt2[Va].ReSize(m_total_orbs*2,m_act_size*2,m_act_size*2,m_act_size*2); - vpt2[Vi].ReSize(m_act_size*2,m_act_size*2,m_total_orbs*2,m_act_size*2); -// for(auto& i: vpt2) -// i.second.ReSize(m_total_orbs*2,m_total_orbs*2,m_total_orbs*2,m_total_orbs*2); - //fock.ReSize(m_total_orbs*2); - - - msg.resize(0); - ReadMeaningfulLine(dumpFile, msg, msgsize); //this if the first line with integrals - - if(m_total_orbs >= m_integral_disk_storage_thresh) // - { - assert(false && "Partial integral file for mps_nevpt2 is not implemented."); - } - - int i, j, k, l; - double value; - //Read zero order integral in active space; - while(msg.size() != 0) { - boost::split(tok, msg, is_any_of(" \t"), token_compress_on); - if (tok.size() != 5) { - pout << "error in reading orbital file"<(type)](2*I,2*K,2*J,2*L) = value; - } - msg.resize(0); - ReadMeaningfulLine(dumpFile, msg, msgsize); //this if the first line with integrals - } - } - - for(int type=0; type< 3 ; type++){ - msg.resize(0); - ReadMeaningfulLine(dumpFile, msg, msgsize); //this if the first line with integrals - while(msg.size() != 0) { - boost::split(tok, msg, is_any_of(" \t"), token_compress_on); - if (tok.size() != 5) { - pout << "error in reading orbital file"<(type)](2*reorder.at(i),2*reorder.at(j)) = value; -// } -// else { -// assert(false && "Only one electron integral perturbation"); -// } -// msg.resize(0); -// ReadMeaningfulLine(dumpFile, msg, msgsize); //this if the first line with integrals -// } -// -// } -// - dumpFile.close(); - } -} - -void SpinAdapted::Input::readorbitalsfile(string& orbitalfile, OneElectronArray& v1, TwoElectronArray& v2, double& coreEnergy, PairArray& vcc, CCCCArray& vcccc, CCCDArray& vcccd) { - - ifstream dumpFile; - dumpFile.open(orbitalfile.c_str(), ios::in); - - string msg; int msgsize = 5000; - vector tok; - std::ofstream ReorderFileOutput; - std::ifstream ReorderFileInput; - char ReorderFileName[5000]; - std::vector reorder; - - int offset = m_orbformat == DMRGFORM ? 0 : 1; - - if (mpigetrank() == 0) { - ReadMeaningfulLine(dumpFile, msg, msgsize); - boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); - - if(offset != 0) - m_norbs = 2*atoi(tok[2].c_str()); - else - m_norbs = 2*atoi(tok[1].c_str()); - - m_num_spatial_orbs = 0; - m_spin_orbs_symmetry.resize(m_norbs); - m_spin_to_spatial.resize(m_norbs); - - sprintf(ReorderFileName, "%s%s", save_prefix().c_str(), "/RestartReorder.dat"); - } - boost::filesystem::path p(ReorderFileName); - -#ifndef SERIAL - boost::mpi::communicator world; - mpi::broadcast(world,m_restart,0); - mpi::broadcast(world,m_fullrestart,0); -#endif - - //do the reordering only if it is not a restart calculation - //if it is then just read the reorder.dat from the scratch space - if(get_restart() || get_fullrestart()) { - if (mpigetrank() == 0) { - ReorderFileInput.open(ReorderFileName); - if(!boost::filesystem::exists(p)) { - pout << "==============="<> m_reorder[i]; - ReorderFileInput.close(); - } - } - } - else { - if (mpigetrank() == 0) { - ReorderFileOutput.open(ReorderFileName); - } - - // read the reorder file or calculate the reordering using one of the many options - if (m_reorderType == FIEDLER) { - if (mpigetrank() == 0) { - m_reorder=get_fiedler_bcs(orbitalfile); - pout << "Fiedler-vector orbital ordering: "; - } - } else if (m_reorderType == GAOPT) { - if (mpigetrank() == 0) { - pout << "GAOPT orbital ordering for BCS calculation not implemented"; - } - abort(); - } else if (m_reorderType == MANUAL) { - if (mpigetrank() == 0) { - ifstream reorderFile(m_reorderfile.c_str()); - CheckFileExistence(m_reorderfile, "Reorder file "); - readreorderfile(reorderFile, m_reorder); - pout << "Manually provided orbital ordering: "; - } - } else { //this is the no-reorder case - if (mpigetrank() == 0) { - m_reorder.resize(m_norbs/2); - for (int i=0; i O_unreordered(2,3) (the former is stored internally and the latter is what is given during input) - // but for this m_reorder the reorder vector below would be 3 1 2 4 and O_unreordered(1,2) -> O_reorder(3, 1) - if (mpigetrank() == 0) { - reorder.resize(m_norbs/2); - for (int i=0; i= 0 ) - ir = atoi(tok[i].c_str()) - offset; - else if (atoi(tok[i].c_str()) < -1) - ir = atoi(tok[i].c_str()) + offset; - - if (sym == "trans") ir += 1; //for translational symmetry the lowest irrep is 0 - if (sym == "lzsym") ir = atoi(tok[i].c_str()); - - m_spin_orbs_symmetry[2*reorderOrbInd] = ir; - m_spin_orbs_symmetry[2*reorderOrbInd+1] = ir; - - if (readLine == numRead) { - m_num_spatial_orbs++; - m_spatial_to_spin.push_back(orbindex); - numRead = Symmetry::sizeofIrrep(ir); - readLine = 0; - } - m_spin_to_spatial[orbindex] = m_num_spatial_orbs-1; - m_spin_to_spatial[orbindex+1] = m_num_spatial_orbs-1; - orbindex +=2; - readLine++; - } - msg.resize(0); - ReadMeaningfulLine(dumpFile, msg, msgsize); - boost::split(tok, msg, is_any_of("=, \t"), token_compress_on); - if(boost::iequals(tok[0], "IUHF")) - RHF=false; - } - - if(sym == "dinfh" ) { - m_spatial_to_spin.clear(); - for (int i=0; i= m_integral_disk_storage_thresh) // - { - for (int i=0; i orb; - for (int j=m_spatial_to_spin[i]; j SpinAdapted::Input::get_fiedler(string& dumpname){ - Matrix fiedler; - ifstream dumpFile; - dumpFile.open(dumpname.c_str(), ios::in); - genetic::ReadIntegral(dumpFile, fiedler); - dumpFile.close(); - SymmetricMatrix fiedler_sym; - fiedler_sym << fiedler; - std::vector findices = fiedler_reorder(fiedler_sym); - return findices; -} - -std::vector SpinAdapted::Input::get_fiedler_nevpt(string& dumpname, int nact){ - Matrix fiedler; - ifstream dumpFile; - dumpFile.open(dumpname.c_str(), ios::in); - genetic::ReadIntegral_nevpt(dumpFile, fiedler, nact); - dumpFile.close(); - SymmetricMatrix fiedler_sym; - fiedler_sym << fiedler; - std::vector findices = fiedler_reorder(fiedler_sym); - return findices; -} - -std::vector SpinAdapted::Input::get_fiedler_bcs(string& dumpname) { - Matrix fiedler; - ifstream dumpFile; - dumpFile.open(dumpname.c_str(), ios::in); - genetic::ReadIntegral_BCS(dumpFile, fiedler); - dumpFile.close(); - SymmetricMatrix fiedler_sym; - fiedler_sym << fiedler; - std::vector findices = fiedler_reorder(fiedler_sym); - return findices; -} - -std::vector SpinAdapted::Input::getgaorder(ifstream& gaconfFile, string& orbitalfile, std::vector& fiedlerorder) -{ - ifstream dumpFile; dumpFile.open(orbitalfile.c_str()); - return genetic::gaordering(gaconfFile, dumpFile, fiedlerorder).Gen().Sequence(); -} - -#ifdef MOLPRO -void SpinAdapted::Input::writeSummaryForMolpro() -{ - pout << setw(50) << "Total number of orbitals : " ; - pout << m_norbs/2 << endl; - pout << setw(50) << "Symmetry of targeted wavefunctions : " ; - pout << m_alpha + m_beta << ":" << m_alpha-m_beta << ":" << m_total_symmetry_number.getirrep()+1 << endl; - pout << setw(50) << "Number of wavefunctions targeted : " ; - pout << m_nroots << endl; - if (m_nroots >1) { - pout << setw(50) << "The weights of the wavefunctions : "; - for (int i=0; i1) { -// printf("%-50s : ", "The weights of the wavefunctions"); -// for (int i=0; i 1) - { - pout << setw(50) << left << "The weights of the wavefunctions" << " : "; - for (int i = 0; i < m_nroots; ++i) pout << left << setprecision(2) << setw(10) << scientific << m_weights[i]; - pout << endl; - } - pout << setw(50) << left << "Symmetry of the molecule" << " : " - << sym.c_str() << endl; - if (sym != "c1") - { - pout << setw(50) << left << "Irreducible representations of the orbitals" << " : "; - for (int i = 0; i < m_spin_orbs_symmetry.size(); i+=2) - { - pout << Symmetry::stringOfIrrep(m_spin_orbs_symmetry[i]) << " "; - } - pout << endl; - } - - pout << endl << "Schedule" << endl; - pout << "--------" << endl; - pout << setw(10) << left << "Iter" << " : " - << setw(20) << left << "# States" << " " - << setw(20) << left << "Davidson_tol" << " " - << setw(20) << left << "Random_noise" << endl; - - for (int i = 0; i < m_sweep_iter_schedule.size(); ++i) - { - pout << setw(10) << left << m_sweep_iter_schedule[i] << " : " - << setw(20) << left << m_sweep_state_schedule[i] << " " - << setw(20) << left << scientific << m_sweep_tol_schedule[i] << " " - << setw(20) << left << scientific << m_sweep_noise_schedule[i] << endl; - } - if (m_algorithm_type == TWODOT_TO_ONEDOT) - { - pout << setw(50) << left << "Switching from twodot to onedot algorithm" << " : " - << m_twodot_to_onedot_iter << endl; - } - pout << setw(50) << left << "Maximum sweep iterations" << " : " << m_maxiter << endl; - -#ifndef SERIAL - } -#endif -} - -void SpinAdapted::Input::performSanityTest() -{ -#ifndef SERIAL - if (mpigetrank() == 0) { -#endif - if (m_num_spatial_orbs <= 0) { - pout << "total number of orbitals has to be a positive number"< 200) { - pout << "Number of orbitals cannot be greater than 130"< 10000) { - pout << "default schedule only works for maxM less than 10000"< m_startM) { - pout << "lastM is larger than startM" << endl; - pout << "Make sure you specify a lastM larger than " << lastM << endl; - pout << "or specify a larger startM " << endl; - abort(); - } - } - else { - // This needs to go at the end of sanity check, or move up schedule - if (m_maxM != 0) { - pout << "With detailed schedule a non-zero maxM should not be specified"<= m_maxiter) { - pout << "Switch from twodot to onedot algorithm cannot happen after maxiter"< openCopy, closedCopy, reorder; - openCopy.resize(m_openorbs.size(), -1); - closedCopy.resize(m_closedorbs.size(), -1); - - for (int i=0; i hf_occupancy_tmp(m_num_spatial_orbs*2,0); - - if (m_Bogoliubov) { - // overwrite hf_occ_user option, since initial guess is always vacuum in BCS case - m_hf_occupancy.assign(m_norbs, 0); - } - else if (m_hf_occ_user == "manual") { - //check if n_orbs is correct and if n_elec is correct - if (m_hf_occupancy.size() != m_num_spatial_orbs ) { - pout << "ERROR: The length of user-defined HF occupancies does not match the number of orbitals " << endl; - pout << "Length of occupancies is: " << m_hf_occupancy.size() << ", and the number of orbitals is: " << m_num_spatial_orbs<< endl; - abort(); - } - int UserElectrons = 0; - for (int i=0; i ele_map; -// for( int i = 0; i < m_norbs/2; ++i ){ -// ele_map.insert( pair( v_1[0](2*i, 2*i), i ) ); -// } - -// multimap :: iterator it_alpha = ele_map.begin(); -// for( int i = 0; i < m_alpha; ++i ){ -// int ia = it_alpha->second; -// hf_occupancy_tmp.at( 2*ia ) = 1; -// if (i < m_beta) -// hf_occupancy_tmp.at( 2*ia+1 ) = 1; -// ++it_alpha; -// } - hf_occupancy_tmp = this->hfOccGenerator_(); - } - else { - pout << "currently other options besides manual, integral and canonical are not implemented."< reorder(m_num_spatial_orbs); - for (int i=0; i= dmrginp.sweep_iter_schedule()[iter]) - current = iter; - - } - - int nroots = m_nroots; - if (m_noise_type == EXCITEDSTATE && m_sweep_additional_noise_schedule[current] != 0.0 && nroots == 1) - nroots++; - - if (setStateSpecific()) - return 1; - - return nroots; -} - -std::vector SpinAdapted::Input::weights(int sweep_iter) const -{ - int iter; - int current = 0; - for (iter = 0; iter < dmrginp.sweep_iter_schedule().size(); ++iter) - { - if (sweep_iter >= dmrginp.sweep_iter_schedule()[iter]) - current = iter; - - } - - if (setStateSpecific()) - return std::vector(1,1.0); - - int nroots = this->nroots(sweep_iter); - std::vector weights(nroots); - for (int i=0; i< m_nroots; i++) - weights[i] = m_weights[i]; - if (nroots == m_nroots+1) - weights[nroots-1] = m_sweep_additional_noise_schedule[current]; - if (nroots != m_nroots+1 && nroots != m_nroots) - { - pout << "Something wrong with the number of weights of wavefunctions!"<& sites = b.get_sites(); - int n=0, s=0; - for (int i=0; i -#include -#include -#include -#include -#include "IrrepSpace.h" -#include "SpinQuantum.h" -#include "timer.h" -#include "couplingCoeffs.h" -#include -#include "IntegralMatrix.h" - - -namespace SpinAdapted{ -class SpinBlock; -class OneElectronArray; -class TwoElectronArray; -//typedef OneElectronArray PerturbOneElectronArray; -//typedef TwoElectronArray PerturbTwoElectronArray; -class PerturbTwoElectronArray; -class PairArray; -class CCCCArray; -class CCCDArray; -enum OnePerturbType{ Va_1=0, Vi_1, Vai_1, OnePertEnd}; -enum TwoPerturbType{Va=0,Vi,Vab,Vij,Vai,Vabi,Vaij,Vabij,TwoPertEnd}; -//enum OnePerturbType : int; -//enum TwoPerturbType : int; - -enum WarmUpTypes {WILSON, LOCAL0, LOCAL2, LOCAL3, LOCAL4}; -enum hamTypes {QUANTUM_CHEMISTRY, HUBBARD, BCS, HEISENBERG}; -enum solveTypes {LANCZOS, DAVIDSON, CONJUGATE_GRADIENT}; -enum algorithmTypes {ONEDOT, TWODOT, TWODOT_TO_ONEDOT}; -enum noiseTypes {RANDOM, EXCITEDSTATE}; -enum calcType {DMRG, ONEPDM, TWOPDM, THREEPDM, FOURPDM, NEVPT2PDM, RESTART_TWOPDM, - RESTART_ONEPDM, RESTART_THREEPDM, RESTART_FOURPDM, RESTART_NEVPT2PDM, TINYCALC, FCI, - EXCITEDDMRG, CALCOVERLAP, CALCHAMILTONIAN, COMPRESS, RESPONSE, RESPONSEBW, - TRANSITION_ONEPDM, TRANSITION_TWOPDM, TRANSITION_THREEPDM, RESTART_T_ONEPDM, RESTART_T_TWOPDM, RESTART_T_THREEPDM, - NEVPT2,RESTART_NEVPT2, MPS_NEVPT, RESTART_MPS_NEVPT, - DS1_ONEPDM, RESTART_DS1_ONEPDM, DS0_ONEPDM, RESTART_DS0_ONEPDM}; -enum orbitalFormat{MOLPROFORM, DMRGFORM}; -enum reorderType{FIEDLER, GAOPT, MANUAL, NOREORDER}; -enum keywords{ORBS, LASTM, STARTM, MAXM, REORDER, HF_OCC, SCHEDULE, SYM, NELECS, SPIN, IRREP, - MAXJ, PREFIX, NROOTS, DOCD, DEFLATION_MAX_SIZE, MAXITER, BASENERGY, - SCREEN_TOL, ODOT, SWEEP_TOL, OUTPUTLEVEL, NONSPINADAPTED, BOGOLIUBOV, TWODOT_NOISE, WARMUP, NPDM_INTERMEDIATE, NPDM_NO_INTERMEDIATE, NPDM_MULTINODE, NPDM_NO_MULTINODE, SPECIFICPDM, NUMKEYWORDS}; - -class Input { - - private: - std::vector m_thrds_per_node; - int m_norbs; - int m_alpha; - int m_beta; -// - int m_bra_alpha; - int m_bra_beta; - int n_bra_spin; -// - int m_Sz; - bool m_spinAdapted; - bool m_Bogoliubov; - int m_permSymm; - - IrrepSpace m_total_symmetry_number; - IrrepSpace m_bra_symmetry_number;// This is used when bra and ket have different spatial symmetry irrep; - // It is only used for transition density matrix calculations. - bool m_transition_diff_spatial_irrep; - bool m_transition_diff_spin =false; //Elvira: This is used when bra and ket have different spins - SpinQuantum m_bra_molecule_quantum; // It is only to calculate pdms for SOC - - SpinQuantum m_molecule_quantum; - int m_total_spin; - int m_guess_permutations; - bool m_stateSpecific; //when targetting excited states we switch from state - //average to statespecific - int m_occupied_orbitals; - - vector m_openorbs; - vector m_closedorbs; - vector m_baseState; - vector m_projectorState; - int m_targetState; - int m_guessState; - - std::vector m_hf_occupancy; - std::string m_hf_occ_user; - std::vector m_weights; - - std::vector m_sweep_iter_schedule; - std::vector m_sweep_state_schedule; - std::vector m_sweep_qstate_schedule; - std::vector m_sweep_tol_schedule; - std::vector m_sweep_noise_schedule; - std::vector m_sweep_additional_noise_schedule; - bool m_schedule_type_default; - bool m_schedule_type_backward; - int m_lastM; - int m_startM; - int m_maxM; - int m_integral_disk_storage_thresh; - int m_num_Integrals; - - bool m_do_diis; - double m_diis_error; - int m_start_diis_iter; - int m_diis_keep_states; - double m_diis_error_tol; - - calcType m_calc_type; - noiseTypes m_noise_type; - hamTypes m_ham_type; - WarmUpTypes m_warmup; - int m_nroots; - solveTypes m_solve_type; - bool m_do_deriv; - bool m_do_fci; - bool m_do_pdm; - bool m_do_npdm_ops; - bool m_do_npdm_in_core; - bool m_npdm_generate; - bool m_new_npdm_code; - bool m_store_spinpdm; - bool m_spatpdm_disk_dump; - bool m_pdm_unsorted; - bool m_store_nonredundant_pdm; - bool m_npdm_intermediate; - std::vector m_specificpdm; - bool m_npdm_multinode; - bool m_set_Sz; - int m_maxiter; - double m_oneindex_screen_tol; - double m_twoindex_screen_tol; - bool m_no_transform; - bool m_add_noninteracting_orbs; - bool m_non_SE; - int m_nquanta; - int m_sys_add; - int m_env_add; - - int m_deflation_min_size; - int m_deflation_max_size; - - algorithmTypes m_algorithm_type; - int m_twodot_to_onedot_iter; - std::vector< std::map > m_quantaToKeep; - std::string m_save_prefix; - std::string m_load_prefix; - bool m_direct; - std::vector m_orbenergies; - - int m_maxj; - ninejCoeffs m_ninej; - int m_max_lanczos_dimension; - - int n_twodot_noise; - double m_twodot_noise; - double m_twodot_gamma; - - double m_sweep_tol; - bool m_restart; - bool m_backward; - bool m_fullrestart; - bool m_restart_warm; - bool m_reset_iterations; - bool m_implicitTranspose; - - - std::vector m_spin_vector; - std::vector m_spin_orbs_symmetry; - int m_num_spatial_orbs; - std::vector m_spatial_to_spin; - std::vector m_spin_to_spatial; - - int m_act_size; - int m_core_size; - int m_virt_size; - int m_total_orbs; - int m_nevpt_state_num; - std::vector m_total_spin_orbs_symmetry; - std::vector m_total_spatial_to_spin; - std::vector m_total_spin_to_spatial; - - int m_outputlevel; - orbitalFormat m_orbformat; - - int m_reorderType; - string m_reorderfile; - std::vector m_reorder;//this can be manual, fiedler, gaopt or noreorder - string m_gaconffile; - - bool m_calc_ri_4pdm; - bool m_store_ripdm_readable; - bool m_nevpt2; - bool m_conventional_nevpt2; - 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) - { - ar & m_thrds_per_node & m_spinAdapted & m_Bogoliubov & m_stateSpecific & m_implicitTranspose & m_num_Integrals; - ar & m_norbs & m_alpha & m_beta & m_solve_type & m_Sz & m_set_Sz & m_baseState& m_projectorState& m_targetState; - ar & m_spin_vector & m_spin_orbs_symmetry & m_guess_permutations & m_nroots & m_weights & m_hf_occ_user & m_hf_occupancy; - ar & m_sweep_iter_schedule & m_sweep_state_schedule & m_sweep_qstate_schedule & m_sweep_tol_schedule & m_sweep_noise_schedule &m_sweep_additional_noise_schedule & m_reorder; - ar & m_molecule_quantum & m_total_symmetry_number & m_total_spin & m_orbenergies & m_add_noninteracting_orbs; - ar & m_bra_symmetry_number & m_permSymm & m_openorbs & m_closedorbs; - ar & m_bra_molecule_quantum & m_bra_symmetry_number & m_add_noninteracting_orbs; //diff spins case - ar & m_bra_alpha & m_bra_beta & n_bra_spin; //diff spins case - ar & m_transition_diff_spin; - ar & m_non_SE; -// ar & m_save_prefix & m_load_prefix & m_direct & m_max_lanczos_dimension; - ar & m_direct & m_max_lanczos_dimension; - ar & m_deflation_min_size & m_deflation_max_size & m_outputlevel & m_reorderfile; - ar & m_algorithm_type & m_twodot_to_onedot_iter & m_orbformat ; - ar & m_nquanta & m_sys_add & m_env_add & m_do_fci & m_no_transform ; - ar & m_do_pdm & m_do_npdm_ops & m_do_npdm_in_core & m_npdm_generate & m_new_npdm_code & m_specificpdm; - ar & m_transition_diff_spatial_irrep & m_occupied_orbitals; - ar & m_store_spinpdm &m_spatpdm_disk_dump & m_pdm_unsorted & m_npdm_intermediate & m_npdm_multinode; - ar & m_maxj & m_ninej & m_maxiter & m_do_deriv & m_oneindex_screen_tol & m_twoindex_screen_tol & m_quantaToKeep & m_noise_type; - ar & m_sweep_tol & m_restart & m_backward & m_fullrestart & m_restart_warm & m_reset_iterations & m_calc_type & m_ham_type & m_warmup; - ar & m_do_diis & m_diis_error & m_start_diis_iter & m_diis_keep_states & m_diis_error_tol & m_num_spatial_orbs; - ar & m_spatial_to_spin & m_spin_to_spatial & m_maxM & m_schedule_type_backward & m_schedule_type_default &m_integral_disk_storage_thresh; - 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; - } - - - void initialize_defaults(); - - std::vector hfOccGenerator_ (); - - public: - //Input() : m_ninej(ninejCoeffs::getinstance()){} - Input() {} - Input (const std::string& config_name); - // ROA - void initCumulTimer() - { - ddscreen = boost::shared_ptr (new cumulTimer()); - cdscreen = boost::shared_ptr (new cumulTimer()); - dscreen = boost::shared_ptr (new cumulTimer()); - buildcsfops = boost::shared_ptr (new cumulTimer()); - sysdotmake = boost::shared_ptr (new cumulTimer()); - initnewsystem = boost::shared_ptr (new cumulTimer()); - guessgenT = boost::shared_ptr (new cumulTimer()); - multiplierT = boost::shared_ptr (new cumulTimer()); - operrotT = boost::shared_ptr (new cumulTimer()); - davidsonT = boost::shared_ptr (new cumulTimer()); - rotmatrixT = boost::shared_ptr (new cumulTimer()); - blockdavid = boost::shared_ptr (new cumulTimer()); - datatransfer = boost::shared_ptr (new cumulTimer()); - hmultiply = boost::shared_ptr (new cumulTimer()); - oneelecT = boost::shared_ptr (new cumulTimer()); - twoelecT = boost::shared_ptr (new cumulTimer()); - makeopsT = boost::shared_ptr (new cumulTimer()); - collectqT = boost::shared_ptr (new cumulTimer()); - opallocate = boost::shared_ptr (new cumulTimer()); - opcatenate = boost::shared_ptr (new cumulTimer()); - oprelease = boost::shared_ptr (new cumulTimer()); - opequateT = boost::shared_ptr (new cumulTimer()); - justmultiply = boost::shared_ptr (new cumulTimer()); - spinrotation = boost::shared_ptr (new cumulTimer()); - otherrotation = boost::shared_ptr (new cumulTimer()); - solvewf = boost::shared_ptr (new cumulTimer()); - postwfrearrange = boost::shared_ptr (new cumulTimer()); - couplingcoeff = boost::shared_ptr (new cumulTimer()); - buildsumblock = boost::shared_ptr (new cumulTimer()); - buildblockops = boost::shared_ptr (new cumulTimer()); - addnoise = boost::shared_ptr (new cumulTimer()); - s0time = boost::shared_ptr (new cumulTimer()); - s1time = boost::shared_ptr (new cumulTimer()); - s2time = boost::shared_ptr (new cumulTimer()); - blockintegrals = boost::shared_ptr (new cumulTimer()); - blocksites = boost::shared_ptr (new cumulTimer()); - statetensorproduct = boost::shared_ptr (new cumulTimer()); - statecollectquanta = boost::shared_ptr (new cumulTimer()); - builditeratorsT = boost::shared_ptr (new cumulTimer()); - diskio = boost::shared_ptr (new cumulTimer()); - } - void writeSummary(); -#ifdef MOLPRO - void writeSummaryForMolpro(); -#endif - void performSanityTest(); - void generateDefaultSchedule(); - void readorbitalsfile(string& dumpFile, OneElectronArray& v1, TwoElectronArray& v2, double& coreEnergy, int integralIndex); - void readorbitalsfile(string& dumpFile, OneElectronArray& v1, TwoElectronArray& v2, OneElectronArray& vpt1, std::map& vpt2, double& coreEnergy); - void readorbitalsfile(string& dumpFile, OneElectronArray& v1, TwoElectronArray& v2, double& coreEnergy, PairArray& vcc, CCCCArray& vcccc, CCCDArray& vcccd); - int getNumIntegrals() { return m_num_Integrals;} - void readreorderfile(ifstream& dumpFile, std::vector& reorder); - std::vector getgaorder(ifstream& gaconfFile, string& orbitalfile, std::vector& fiedlerorder); - std::vector get_fiedler(string& dumpname); - std::vector get_fiedler_nevpt(string& dumpname, int nact); - std::vector get_fiedler_bcs(string& dumpname); - void usedkey_error(string& key, string& line); - void makeInitialHFGuess(); - static void ReadMeaningfulLine(ifstream&, string&, int); - - boost::shared_ptr dscreen; - boost::shared_ptr cdscreen; - boost::shared_ptr ddscreen; - boost::shared_ptr buildcsfops; - boost::shared_ptr sysdotmake; - boost::shared_ptr initnewsystem; - boost::shared_ptr diskio; - boost::shared_ptr builditeratorsT; - boost::shared_ptr blockintegrals; - boost::shared_ptr blocksites; - boost::shared_ptr statetensorproduct; - boost::shared_ptr statecollectquanta; - boost::shared_ptr guessgenT; - boost::shared_ptr multiplierT; - boost::shared_ptr operrotT; - boost::shared_ptr davidsonT; - boost::shared_ptr rotmatrixT; - boost::shared_ptr blockdavid; - boost::shared_ptr datatransfer; - boost::shared_ptr hmultiply; - boost::shared_ptr oneelecT; - boost::shared_ptr twoelecT; - boost::shared_ptr makeopsT; - boost::shared_ptr collectqT; - boost::shared_ptr opallocate; - boost::shared_ptr opcatenate; - boost::shared_ptr oprelease; - boost::shared_ptr opequateT; - boost::shared_ptr justmultiply; - boost::shared_ptr spinrotation; - boost::shared_ptr otherrotation; - boost::shared_ptr solvewf; - boost::shared_ptr postwfrearrange; - boost::shared_ptr couplingcoeff; - boost::shared_ptr buildsumblock; - boost::shared_ptr buildblockops; - boost::shared_ptr addnoise; - boost::shared_ptr s0time; - boost::shared_ptr s1time; - boost::shared_ptr s2time; - - std::vector& get_openorbs() { return m_openorbs;} - std::vector& get_closedorbs() { return m_closedorbs;} - const std::vector& baseStates() const {return m_baseState;} - const int& targetState() const {return m_targetState;} - const int& guessState() const {return m_guessState;} - int& setGuessState() {return m_guessState;} - const std::vector& projectorStates() const {return m_projectorState;} - - std::vector& baseStates() {return m_baseState;} - int& targetState() {return m_targetState;} - std::vector& projectorStates() {return m_projectorState;} - - void reorderOpenAndClosed(); - const int& num_occupied_orbitals() const {return m_occupied_orbitals;} - const bool& doimplicitTranspose() const {return m_implicitTranspose;} - bool& setimplicitTranspose() {return m_implicitTranspose;} - const bool& setStateSpecific() const {return m_stateSpecific;} - bool& setStateSpecific() {return m_stateSpecific;} - const orbitalFormat& orbformat() const {return m_orbformat;} - const int& outputlevel() const {return m_outputlevel;} - int& setOutputlevel() {return m_outputlevel;} - const vector& spatial_to_spin() const {return m_spatial_to_spin;} - int spatial_to_spin(int i) const {return m_spatial_to_spin.at(i);} - const vector& spin_to_spatial() const {return m_spin_to_spatial;} - const double& diis_error_tol() const {return m_diis_error_tol;} - const bool& do_diis() const {return m_do_diis;} - const double& diis_error() const {return m_diis_error;} - const int& start_diis_iter() const {return m_start_diis_iter;} - const int& diis_keep_states() const {return m_diis_keep_states;} - - bool use_partial_two_integrals() const {return (m_norbs/2 >= m_integral_disk_storage_thresh);} - bool& set_fullrestart() {return m_fullrestart;} - const bool& get_fullrestart() const {return m_fullrestart;} - const bool& get_backward() const {return m_backward;} - const double& get_sweep_tol() const {return m_sweep_tol;} - const int& get_twodot_method() const {return n_twodot_noise;} - const double& get_twodot_noise() const {return m_twodot_noise;} - double& set_twodot_noise() {return m_twodot_noise;} - const double& get_twodot_gamma() const {return m_twodot_gamma;} - double& set_twodot_gamma() {return m_twodot_gamma;} - const bool& get_restart() const {return m_restart;} - const bool& get_restart_warm() const {return m_restart_warm;} - const bool& get_reset_iterations() const {return m_reset_iterations;} - const ninejCoeffs& get_ninej() const {return m_ninej;} - const hamTypes &hamiltonian() const {return m_ham_type;} - const WarmUpTypes &warmup() const {return m_warmup;} - const int &guess_permutations() const { return m_guess_permutations; } - const int &max_lanczos_dimension() const {return m_max_lanczos_dimension;} - std::vector thrds_per_node() const { return m_thrds_per_node; } - const calcType &calc_type() const { return m_calc_type; } - calcType &calc_type() { return m_calc_type; } - const solveTypes &solve_method() const { return m_solve_type; } - const noiseTypes &noise_type() const {return m_noise_type;} - const bool &set_Sz() const {return m_set_Sz;} - const algorithmTypes &algorithm_method() const { return m_algorithm_type; } - algorithmTypes &set_algorithm_method() { return m_algorithm_type; } - int twodot_to_onedot_iter() const { return m_twodot_to_onedot_iter; } - std::vector< std::map >& get_quantaToKeep() { return m_quantaToKeep;} - const std::vector &hf_occupancy() const { return m_hf_occupancy; } - const std::vector &spin_orbs_symmetry() const { return m_spin_orbs_symmetry; } - std::vector weights(int sweep_iter) const;// { return m_weights; } - std::vector weights() const { return m_weights; } - const std::vector &sweep_iter_schedule() const { return m_sweep_iter_schedule; } - const std::vector &sweep_state_schedule() const { return m_sweep_state_schedule; } - const std::vector &sweep_qstate_schedule() const { return m_sweep_qstate_schedule; } - const std::vector &sweep_tol_schedule() const { return m_sweep_tol_schedule; } - const std::vector &sweep_noise_schedule() const { return m_sweep_noise_schedule; } - const std::vector &sweep_additional_noise_schedule() const { return m_sweep_additional_noise_schedule; } - std::vector &set_sweep_noise_schedule() { return m_sweep_noise_schedule; } - std::vector &set_sweep_additional_noise_schedule() { return m_sweep_additional_noise_schedule; } - int& Sz() {return m_Sz;} - int nroots(int sweep_iter) const; - int nroots() const {return m_nroots;} - int real_particle_number() const { return (m_alpha + m_beta);} - int total_particle_number() const { if(!m_add_noninteracting_orbs) return (m_alpha + m_beta); else return (2*m_alpha); } - int bra_particle_number() const { if(!m_add_noninteracting_orbs) return (m_bra_alpha + m_bra_beta); else return (2*m_bra_alpha); } // diff spins case - bool calc_ri_4pdm() const {return m_calc_ri_4pdm;} - bool store_ripdm_readable() const {return m_store_ripdm_readable;} - bool nevpt2() const {return m_nevpt2;} - bool read_higherpdm() const {return m_conventional_nevpt2;} - int kept_nevpt2_states() const {return m_kept_nevpt2_states;} - bool Print() const {return NevPrint.first;} - int PrintIndex() const{return NevPrint.second;} - void SetPrint(bool p, int i=0){NevPrint.first = p;NevPrint.second=i;} - const SpinSpace total_spin_number() const { if (!m_add_noninteracting_orbs) return SpinSpace(m_alpha - m_beta); else return SpinSpace(0); } - const SpinSpace bra_spin_number() const { if (!m_add_noninteracting_orbs) return SpinSpace(m_bra_alpha - m_bra_beta); else return SpinSpace(0); } // diff spins case - int last_site() const - { - if(m_spinAdapted) - { - if(m_calc_type == MPS_NEVPT) - return m_act_size; - else - return m_num_spatial_orbs; - } - else - { - if(m_calc_type == MPS_NEVPT) - return 2*m_act_size; - else - return 2*m_num_spatial_orbs; - } - } - const bool &no_transform() const { return m_no_transform; } - const int &deflation_min_size() const { return m_deflation_min_size; } - const bool &direct() const { return m_direct; } - const int &deflation_max_size() const { return m_deflation_max_size; } - const IrrepSpace &total_symmetry_number() const { return m_total_symmetry_number; } - const IrrepSpace &bra_symmetry_number() const { return m_bra_symmetry_number; } - const SpinQuantum &molecule_quantum() const { return m_molecule_quantum; } - const SpinQuantum &bra_molecule_quantum() const { return m_bra_molecule_quantum; } // diff spins case - const int &sys_add() const { return m_sys_add; } - const bool &add_noninteracting_orbs() const {return m_add_noninteracting_orbs;} - bool &add_noninteracting_orbs() {return m_add_noninteracting_orbs;} - const bool &non_SE() {return m_non_SE;} //for calculations without singlet embedding - const int &nquanta() const { return m_nquanta; } - const int &env_add() const { return m_env_add; } - const bool &do_fci() const { return m_do_fci; } - const int &max_iter() const { return m_maxiter; } - const double &oneindex_screen_tol() const { return m_oneindex_screen_tol; } - double &oneindex_screen_tol() { return m_oneindex_screen_tol; } - const double &twoindex_screen_tol() const { return m_twoindex_screen_tol; } - double &twoindex_screen_tol() { return m_twoindex_screen_tol; } - const int &total_spin() const {return m_total_spin;} - const std::vector &spin_vector() const { return m_spin_vector; } - const std::string &save_prefix() const { return m_save_prefix; } - const std::string &load_prefix() const { return m_load_prefix; } - SpinQuantum& set_molecule_quantum() {return m_molecule_quantum;} - - bool transition_diff_spin() { - return m_transition_diff_spin; - } - - SpinQuantum effective_molecule_quantum() { - if (!m_add_noninteracting_orbs) - return m_molecule_quantum; - else - return SpinQuantum(m_molecule_quantum.particleNumber + m_molecule_quantum.totalSpin.getirrep(), SpinSpace(0), m_molecule_quantum.orbitalSymmetry); - //return SpinQuantum(total_particle_number() + total_spin_number().getirrep(), SpinSpace(0), total_symmetry_number()); - } - -// - SpinQuantum bra_quantum() { - if (!m_add_noninteracting_orbs) - { - if (m_transition_diff_spin) { - return SpinQuantum(bra_particle_number(), SpinSpace(m_bra_alpha - m_bra_beta), bra_symmetry_number()); - } - else - return SpinQuantum(total_particle_number(), SpinSpace(m_alpha - m_beta), bra_symmetry_number()); - } - else { - if (m_transition_diff_spin) { - return SpinQuantum(bra_particle_number()+bra_spin_number().getirrep(), SpinSpace(0), bra_symmetry_number()); - } - else - return SpinQuantum(total_particle_number() + total_spin_number().getirrep(), SpinSpace(0), bra_symmetry_number()); } -} - - - - vector effective_molecule_quantum_vec() { - vector q; - if (!m_Bogoliubov) q.push_back(effective_molecule_quantum()); - else { - SpinQuantum q_max = effective_molecule_quantum(); - for (int n = 0; n <= q_max.get_n(); n+=2) { - q.push_back(SpinQuantum(n, q_max.get_s(), q_max.get_symm())); - } - } - return q; - } - vector bra_quantum_vec() { - vector q; - if (!m_Bogoliubov) q.push_back(bra_quantum()); - else { - SpinQuantum q_max = bra_quantum(); - for (int n = 0; n <= q_max.get_n(); n+=2) { - q.push_back(SpinQuantum(n, q_max.get_s(), q_max.get_symm())); - } - } - return q; - } - bool transition_diff_irrep(){ - return m_transition_diff_spatial_irrep; - } - std::vector& get_orbenergies() {return m_orbenergies;} - int getHFQuanta(const SpinBlock& b) const; - const bool &do_pdm() const {return m_do_pdm;} - bool &do_pdm() {return m_do_pdm;} - const bool &do_npdm_ops() const {return m_do_npdm_ops;} - bool &do_npdm_ops() {return m_do_npdm_ops;} - const bool &do_npdm_in_core() const {return m_do_npdm_in_core;} - bool &do_npdm_in_core() {return m_do_npdm_in_core;} - bool new_npdm_code() const{ - if( m_do_pdm) return m_new_npdm_code; - else return false; - } - void set_new_npdm_code(){ m_new_npdm_code= true;} - const bool &npdm_generate() const { return m_npdm_generate;} - bool &npdm_generate() { return m_npdm_generate;} - const bool &store_spinpdm() const {return m_store_spinpdm;} - bool &store_spinpdm() {return m_store_spinpdm;} - const bool &spatpdm_disk_dump() const {return m_spatpdm_disk_dump;} - bool &spatpdm_disk_dump() {return m_spatpdm_disk_dump;} - const bool &pdm_unsorted() const {return m_pdm_unsorted;} - bool &pdm_unsorted(){return m_pdm_unsorted;} - const std::vector &specificpdm() const {return m_specificpdm;} - std::vector &specificpdm() {return m_specificpdm;} - const bool &store_nonredundant_pdm() const { return m_store_nonredundant_pdm;} - bool &store_nonredundant_pdm() { return m_store_nonredundant_pdm;} - int slater_size() const {return m_norbs;} - const std::vector &reorder_vector() {return m_reorder;} - const int &act_size() const { return m_act_size;} - const int &core_size() const { return m_core_size;} - const int &virt_size() const { return m_virt_size;} - const int &total_size() const { return m_total_orbs;} - bool spinAdapted() {return m_spinAdapted;} - const int &nevpt_state_num() const {return m_nevpt_state_num;} - bool &npdm_intermediate() { return m_npdm_intermediate; } - 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/npdm/npdm.C.new b/modules/npdm/npdm.C.new deleted file mode 100644 index 9eec765e..00000000 --- a/modules/npdm/npdm.C.new +++ /dev/null @@ -1,598 +0,0 @@ -/* -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 -*/ - -#include "npdm.h" -#include "sweep.h" -#include "sweepgenblock.h" -#include "density.h" -#include "sweeponepdm.h" // For legacy version of 1pdm -#include "sweeptwopdm.h" // For legacy version of 2pdm -#include "npdm_driver.h" -#include "nevpt2_npdm_driver.h" -#include "pario.h" -#include "ds1_sweeponepdm.h" //EL -#include "ds0_sweeponepdm.h" //EL - - -void dmrg(double sweep_tol); -void restart(double sweep_tol, bool reset_iter); -void ReadInput(char* conf); -void fullrestartGenblock(); - -namespace SpinAdapted { -namespace Npdm { - -//----------------------------------------------------------------------------------------------------------------------------------------------------------- - -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; - //mcheck("at the start of block and decimate"); - // figure out if we are going forward or backwards - dmrginp.guessgenT -> start(); - bool forward = (system.get_sites() [0] == 0); - SpinBlock systemDot; - SpinBlock envDot; - int systemDotStart, systemDotEnd; - int systemDotSize = sweepParams.get_sys_add() - 1; - if (forward) - { - systemDotStart = dmrginp.spinAdapted() ? *system.get_sites().rbegin () + 1 : (*system.get_sites().rbegin ())/2 + 1 ; - systemDotEnd = systemDotStart + systemDotSize; - } - else - { - systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ; - systemDotEnd = systemDotStart - systemDotSize; - } - vector spindotsites(2); - spindotsites[0] = systemDotStart; - spindotsites[1] = systemDotEnd; - //if (useSlater) { - systemDot = SpinBlock(systemDotStart, systemDotEnd, system.get_integralIndex(), true); - //SpinBlock::store(true, systemDot.get_sites(), systemDot); - //} - //else - //SpinBlock::restore(true, spindotsites, systemDot); - SpinBlock environment, environmentDot, newEnvironment; - - int environmentDotStart, environmentDotEnd, environmentStart, environmentEnd; - - const int nexact = forward ? sweepParams.get_forward_starting_size() : sweepParams.get_backward_starting_size(); - - system.addAdditionalCompOps(); - 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(), - sweepParams.get_onedot(), nexact, useSlater, environment.get_integralIndex(), true, true, true); - } - 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(), - sweepParams.get_onedot(), nexact, useSlater, environment.get_integralIndex(), true, true, true); - - } - SpinBlock big; - newSystem.set_loopblock(true); - system.set_loopblock(false); - newEnvironment.set_loopblock(false); - 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); - - } - 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); - } - -#ifndef SERIAL - mpi::communicator world; - mpi::broadcast(world, solution, 0); -#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); - - - - //bra and ket rotation matrices are calculated from different density matrices. - - - std::vector rotateMatrix; - std::vector rotateMatrixB; - - if(state!=stateB){ - - DensityMatrix tracedMatrix(newSystem.get_braStateInfo()); - tracedMatrix.allocate(newSystem.get_braStateInfo()); - tracedMatrix.makedensitymatrix(std::vector(1,solution[0]), big, std::vector(1,1.0), 0.0, 0.0, false); - rotateMatrix.clear(); - - DensityMatrix tracedMatrixB(newSystem.get_ketStateInfo()); - tracedMatrixB.allocate(newSystem.get_ketStateInfo()); - tracedMatrixB.makedensitymatrix(std::vector(1,solution[1]), big, std::vector(1,1.0), 0.0, 0.0, false); - rotateMatrixB.clear(); - if (!mpigetrank()){ - double error = makeRotateMatrix(tracedMatrix, rotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates()); - error = makeRotateMatrix(tracedMatrixB, rotateMatrixB, sweepParams.get_keep_states(), sweepParams.get_keep_qstates()); - } - } - else{ - DensityMatrix tracedMatrix(newSystem.get_stateInfo()); - tracedMatrix.allocate(newSystem.get_stateInfo()); - tracedMatrix.makedensitymatrix(std::vector(1,solution[0]), big, std::vector(1,1.0), 0.0, 0.0, false); - rotateMatrix.clear(); - if (!mpigetrank()){ - double error = makeRotateMatrix(tracedMatrix, rotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates()); - } - - } - - - - - - int sweepPos = sweepParams.get_block_iter(); - int endPos = sweepParams.get_n_iters()-1; - npdm_driver.compute_npdm_elements(solution, big, sweepPos, endPos); - SaveRotationMatrix (newSystem.get_sites(), rotateMatrix, state); - solution[0].SaveWavefunctionInfo (big.get_braStateInfo(), big.get_leftBlock()->get_sites(), state); - - if(state!=stateB){ - SaveRotationMatrix (newSystem.get_sites(), rotateMatrixB, stateB); - solution[1].SaveWavefunctionInfo (big.get_ketStateInfo(), big.get_leftBlock()->get_sites(), stateB); - } - - - //FIXME - //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(newSystem.get_sites().size()>1) - //if (!mpigetrank()){ - //LoadRotationMatrix (newSystem.get_sites(), rotateMatrix, state); - //LoadRotationMatrix (newSystem.get_sites(), rotateMatrixB, stateB); - //} - #ifndef SERIAL - mpi::broadcast(world,rotateMatrix,0); - if(state!=stateB) - mpi::broadcast(world,rotateMatrixB,0); - #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. - //Since it is useless, Just omit it at the last sweep site. - // if( sweepParams.get_block_iter() != sweepParams.get_n_iters() - 1) - { - if(state!=stateB) - newSystem.transform_operators(rotateMatrix,rotateMatrixB); - else - newSystem.transform_operators(rotateMatrix); - } - - //newSystem.transform_operators(rotateMatrix,rotateMatrixB); - ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime(); - p3out << "NPDM block and decimate and compute elements " << ewall << " " << ecpu << endl; - -} - -//----------------------------------------------------------------------------------------------------------------------------------------------------------- - -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; - pout.precision(12); - SpinBlock system; - const int nroots = dmrginp.nroots(); - std::vector finalEnergy(nroots,0.); - std::vector finalEnergy_spins(nroots,0.); - double finalError = 0.; - - sweepParams.set_sweep_parameters(); - // 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 - 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; - - if (!restart) sweepParams.set_block_iter() = 0; - if(dmrginp.setStateSpecific() || dmrginp.transition_diff_irrep()){ - if (!restart) SpinBlock::store (forward, system.get_sites(), system, state, stateB ); // if restart, just restoring an existing block -- - } - else{ - if (!restart) SpinBlock::store (forward, system.get_sites(), system, sweepParams.current_root(), sweepParams.current_root() ); // if restart, just restoring an existing block -- - } - - sweepParams.savestate(forward, system.get_sites().size()); - bool dot_with_sys = true; - - // Loop over all block sites - for (; sweepParams.get_block_iter() < sweepParams.get_n_iters(); ) { - Timer timer; - - pout << "\n\t\t\t Block Iteration :: " << sweepParams.get_block_iter() << endl; - pout << "\t\t\t ----------------------------" << endl; - if (forward) { p1out << "\t\t\t Current direction is :: Forwards " << endl; } - else { p1out << "\t\t\t Current direction is :: Backwards " << endl; } - - //if (SHOW_MORE) pout << "system block" << endl << system << endl; - - 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 && - ((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()< npdm_driver; - //if ( (dmrginp.hamiltonian() == QUANTUM_CHEMISTRY) && dmrginp.spinAdapted() ) { - if ( (dmrginp.hamiltonian() == QUANTUM_CHEMISTRY) || (dmrginp.hamiltonian() == BCS)) { - //By default, new_npdm_code is false. - //For npdm_order 1 or 2. new_npdm_code is determined by default or manual setting. - //For the other situation, only old or new code is suitable. - if(npdm_order == NPDM_PAIRMATRIX || npdm_order == NPDM_THREEPDM || npdm_order == NPDM_FOURPDM || npdm_order == NPDM_NEVPT2 || transitionpdm == true || dmrginp.spinAdapted() == false || dmrginp.setStateSpecific()) - dmrginp.set_new_npdm_code(); - - if(dmrginp.new_npdm_code()){ - if (npdm_order == NPDM_ONEPDM) npdm_driver = boost::shared_ptr( new Onepdm_driver( dmrginp.last_site() ) ); - else if ((npdm_order == NPDM_DS0)||(npdm_order == NPDM_DS1)) npdm_driver = boost::shared_ptr( new Onepdm_driver( dmrginp.last_site() ) ); - else if (npdm_order == NPDM_TWOPDM) npdm_driver = boost::shared_ptr( new Twopdm_driver( dmrginp.last_site() ) ); - else if (npdm_order == NPDM_THREEPDM) npdm_driver = boost::shared_ptr( new Threepdm_driver( dmrginp.last_site() ) ); - else if (npdm_order == NPDM_FOURPDM) npdm_driver = boost::shared_ptr( new Fourpdm_driver( dmrginp.last_site() ) ); - else if (npdm_order == NPDM_NEVPT2) npdm_driver = boost::shared_ptr( new Nevpt2_npdm_driver( dmrginp.last_site() ) ); - else if (npdm_order == NPDM_PAIRMATRIX) npdm_driver = boost::shared_ptr( new Pairpdm_driver( dmrginp.last_site() ) ); - else abort(); - } - } - -if (dS) { - for (int state=0; stateclear(); - sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; - if (npdm_order == NPDM_DS1) ds1_onepdm::do_one(sweepParams, false, !direction, false, 0, state, stateB); - else if (npdm_order == NPDM_DS0) ds0_onepdm::do_one(sweepParams, false, !direction, false, 0, state, stateB); - else abort(); - } - } -} -else { - if (dmrginp.specificpdm().size()!=0) - { - Timer timer; - dmrginp.set_fullrestart() = true; - sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; - dmrginp.npdm_generate() = true; - if ( !dmrginp.setStateSpecific()) - SweepGenblock::do_one(sweepParams, false, !direction, false, 0, -1, -1); //this will generate the cd operators - else if (dmrginp.specificpdm().size()==1) - SweepGenblock::do_one(sweepParams, false, !direction, false, 0, dmrginp.specificpdm()[0], dmrginp.specificpdm()[0]); //this will generate the cd operators - else if (dmrginp.specificpdm().size()==2) - SweepGenblock::do_one(sweepParams, false, !direction, false, 0, dmrginp.specificpdm()[0], dmrginp.specificpdm()[1]); //this will generate the cd operators - else abort(); - dmrginp.npdm_generate() = false; - ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime(); - p3out << "\t\t\t NPDM SweepGenblock time " << ewall << " " << ecpu << endl; - dmrginp.set_fullrestart() = false; - - sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; - Timer timerX; - npdm_driver->clear(); - if (dmrginp.specificpdm().size()==1) - npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0,dmrginp.specificpdm()[0],dmrginp.specificpdm()[0]); - else if (dmrginp.specificpdm().size()==2) - npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0,dmrginp.specificpdm()[0],dmrginp.specificpdm()[1]); - else abort(); - ecpu = timerX.elapsedcputime();ewall=timerX.elapsedwalltime(); - p3out << "\t\t\t NPDM sweep time " << ewall << " " << ecpu << endl; - return; - } - - if(dmrginp.transition_diff_irrep()){ - // It is used when bra and ket has different spatial irrep - // For now, only the transtion pdm between two wavefuntions( 1 as bra and 0 as ket) are calculation - // If the spatial irrep information is stored in wavefuntions, transition pdm among several wavefunctions( i for bra and j for ket, there are n(n-1)/2 kinds of situations.) are possible. - for (int state=0; stateclear(); - 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 if( !dmrginp.setStateSpecific()){ - Timer timer; - dmrginp.set_fullrestart() = true; - dmrginp.npdm_generate() = true; - sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; - SweepGenblock::do_one(sweepParams, false, !direction, false, 0, -1, -1); //this will generate the cd operators - dmrginp.npdm_generate() = false; - ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime(); - p3out << "\t\t\t NPDM SweepGenblock time " << ewall << " " << ecpu << endl; - dmrginp.set_fullrestart() = false; - - -//pout << "(orz3rdmalt1,orz3rdmalt2) = (" << dmrginp.orz3rdmalt1() << "," << dmrginp.orz3rdmalt2() << ")" << endl; - if(transitionpdm){ - // <\Phi_k|a^+_ia_j|\Phi_l> = <\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; - npdm_driver->clear(); - 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); - else if (npdm_order == NPDM_TWOPDM) SweepTwopdm::do_one(sweepParams, false, direction, false, 0, state); - else abort(); - } - } -} - } - - } - - else { - // state-specific - if(transitionpdm){ - for (int state=0; stateclear(); - sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy; - 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{ - for (int state=0; stateclear(); - npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0, state,state); - } - } - -} - - sweep_copy.savestate(direction_copy, restartsize_copy); - -} -} -} -} diff --git a/modules/npdm/npdm.h.new b/modules/npdm/npdm.h.new deleted file mode 100644 index c7e43ec0..00000000 --- a/modules/npdm/npdm.h.new +++ /dev/null @@ -1,21 +0,0 @@ -/* -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 -*/ - -#ifndef NPDM_HEADER -#define NPDM_HEADER - -namespace SpinAdapted{ -enum NpdmOrder{NPDM_NEVPT2, NPDM_ONEPDM, NPDM_TWOPDM, NPDM_THREEPDM, NPDM_FOURPDM, NPDM_PAIRMATRIX, NPDM_OVERLAP, NPDM_EMPTY, NPDM_DS0, NPDM_DS1}; -namespace Npdm{ - - void npdm(NpdmOrder npdm_order, bool restartpdm=false, bool transitionpdm=false, bool dS = false); - -} -} - -#endif From 60abaf9ab08f2160b5b1bd9aa53cf7ebb0fd4c39 Mon Sep 17 00:00:00 2001 From: Takeshi Yanai Date: Fri, 5 Mar 2021 20:36:48 +0900 Subject: [PATCH 7/8] related to compatibility to icpc --- BaseOperator.h | 5 +++++ btas/include/btas/btas.h | 5 +++++ genetic/GAOptimize.C | 5 +++++ input.h | 4 ++++ renormalise.C | 5 +++++ 5 files changed, 24 insertions(+) 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.h b/input.h index 50e4a35b..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" 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 From 65d4f79eab1fb6058dbeeb9178c4f2228981ad2c Mon Sep 17 00:00:00 2001 From: Jakub Chalupsky Date: Mon, 28 Aug 2023 14:16:14 +0200 Subject: [PATCH 8/8] return statements added --- modules/ds0_onepdm/ds0_sweep.C | 2 ++ modules/ds1_onepdm/ds1_sweep.C | 2 ++ 2 files changed, 4 insertions(+) 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]; } } }