Excited states

Say
Excited=5
to get the gs and the next 4 excited states.
This label is optional and defaults to 0 (which means g.s.)
Only energies are printed for now.

src/Engine/DefaultSymmetry.h

@@ -40,6 +40,7 @@ public:
 	typedef GeometryType_ GeometryType;
 	typedef PsimagLite::CrsMatrix<ComplexOrRealType> SparseMatrixType;
 	typedef typename PsimagLite::Vector<ComplexOrRealType>::Type VectorType;
+	typedef typename PsimagLite::Vector<VectorType>::Type VectorVectorType;
 
 	DefaultSymmetry(const BasisType&,
 	                const GeometryType&,
@@ -97,9 +98,7 @@ public:
 		throw std::runtime_error("DefaultSymmetry: cannot call transformMatrix\n");
 	}
 
-	void transformGs(VectorType&,SizeType)
-	{
-	}
+	void transform(VectorVectorType&, SizeType) {}
 
 	SizeType sectors() const { return 1; }
 

src/Engine/Engine.h

@@ -55,6 +55,7 @@ public:
 	typedef PsimagLite::Random48<RealType> RandomType;
 	typedef PsimagLite::ParametersForSolver<RealType> ParametersForSolverType;
 	typedef typename PsimagLite::Vector<ComplexOrRealType>::Type VectorType;
+	typedef typename PsimagLite::Vector<VectorType>::Type VectorVectorType;
 	typedef PsimagLite::LanczosSolver<ParametersForSolverType,
 	InternalProductType,
 	VectorType> LanczosSolverType;
@@ -79,26 +80,32 @@ public:
 
 	enum {PLUS,MINUS};
 
-	Engine(const ModelType& model,
-	       SizeType,
-	       InputType& io)
+	Engine(const ModelType& model, InputType& io)
 	    : model_(model),
 	      progress_("Engine"),
 	      io_(io),
 	      options_("")
 	{
 		io_.readline(options_,"SolverOptions=");
-		computeGroundState();
+		SizeType excited = 0;
+
+		try {
+			io_.readline(excited, "Excited=");
+		} catch (std::exception&) {}
+
+		computeAllStatesBelow(excited);
 	}
 
-	RealType gsEnergy() const
+	RealType energies(SizeType ind) const
 	{
-		return gsEnergy_;
+		assert(ind < energies_.size());
+		return energies_[ind];
 	}
 
-	const VectorType& eigenvector() const
+	const VectorType& eigenvector(SizeType ind) const
 	{
-		return gsVector_;
+		assert(ind < vectors_.size());
+		return vectors_[ind];
 	}
 
 	//! Calc Green function G(isite,jsite)  (still diagonal in spin)
@@ -138,6 +145,9 @@ public:
 			throw std::runtime_error(str.c_str());
 		}
 
+		assert(0 < vectors_.size());
+		const VectorType& gsVector = vectors_[0];
+
 		typedef typename ContinuedFractionCollectionType::ContinuedFractionType
 		        ContinuedFractionType;
 
@@ -164,7 +174,7 @@ public:
 			VectorType modifVector;
 			getModifiedState(modifVector,
 			                 lOperator,
-			                 gsVector_,
+			                 gsVector,
 			                 *basisNew,
 			                 type,
 			                 isite,
@@ -190,13 +200,16 @@ public:
 		}
 	}
 
-	void measure(PsimagLite::String meas) const
+	void measure(const VectorStringType& braOpKet) const
 	{
+		if (braOpKet.size() != 3)
+			err("LanczosDriver1: Only dressed brakets allowed (FATAL ERROR)\n");
+
+		const PsimagLite::String meas = braOpKet[1];
 		VectorStringType tokens;
 		PsimagLite::split(tokens, meas, ";");
 		const SizeType n = tokens.size();
 		VectorRahulOperatorType vops;
-		VectorType psiNew(gsVector_.size());
 		VectorSizeType vsites(n);
 		for (SizeType i = 0; i < n; ++i) {
 			int site = OneOperatorSpecType::extractSiteIfAny(tokens[i]);
@@ -210,21 +223,31 @@ public:
 			vops.push_back(RahulOperatorType(opspec.label, opspec.dof, opspec.transpose));
 		}
 
-		model_.rahulMethod(psiNew, vops, vsites, gsVector_);
-		const ComplexOrRealType result = gsVector_*psiNew;
+		SizeType ketIndex = 0; // use braOpKet[2] to derive this one FIXME TODO
+		checkBraOrKet(braOpKet[2], ketIndex);
+		const VectorType& ketVector = vectors_[ketIndex];
+		VectorType psiNew(ketVector.size());
+
+		model_.rahulMethod(psiNew, vops, vsites, ketVector);
+
+		SizeType braIndex = 0; // use braOpKet[0] to derive this one FIXME TODO
+		checkBraOrKet(braOpKet[0], braIndex);
+		const VectorType& braVector = vectors_[braIndex];
+
+		const ComplexOrRealType result = braVector*psiNew;
 		std::cout<<"<gs|"<<meas<<"|gs> = "<<result<<"\n";
 	}
 
 	void twoPoint(PsimagLite::Matrix<typename VectorType::value_type>& result,
 	              const LabeledOperatorType& lOperator,
 	              const PsimagLite::Vector<PairType>::Type& spins,
-	              const PairType& orbs) const
+	              const PairType& orbs,
+	              const PairType& braAndKet) const
 	{
 		for (SizeType i=0;i<spins.size();i++) {
 			std::cout<<"spins="<<spins[i].first<<" "<<spins[i].second<<"\n";
-			twoPoint(result,lOperator,spins[i],orbs);
+			twoPoint(result, lOperator, spins[i], orbs, braAndKet);
 		}
-
 	}
 
 	/* PSIDOC TwoPointCorrelations
@@ -233,7 +256,8 @@ public:
 	void twoPoint(PsimagLite::Matrix<typename VectorType::value_type>& result,
 	              const LabeledOperatorType& lOperator,
 	              const PairType& spins,
-	              const PairType& orbs) const
+	              const PairType& orbs,
+	              const PairType& braAndKet) const
 	{
 		const BasisType* basisNew = 0;
 		PairType oldParts = model_.basis().parts();
@@ -260,6 +284,10 @@ public:
 			basisNew = &model_.basis();
 		}
 
+		checkBraOrKet("bra", braAndKet.first);
+		checkBraOrKet("ket", braAndKet.second);
+		const VectorType& braVector = vectors_[braAndKet.first];
+		const VectorType& ketVector = vectors_[braAndKet.second];
 		SizeType total =result.n_row();
 
 		for (SizeType isite=0;isite<total;isite++)
@@ -276,7 +304,7 @@ public:
 			accModifiedState(modifVector1,
 			                 lOperator,
 			                 *basisNew,
-			                 gsVector_,
+			                 ketVector,
 			                 isite,
 			                 spins.first,
 			                 orbs.first,
@@ -287,7 +315,7 @@ public:
 				accModifiedState(modifVector2,
 				                 lOperator,
 				                 *basisNew,
-				                 gsVector_,
+				                 braVector,
 				                 jsite,
 				                 spins.second,
 				                 orbs.second,
@@ -303,9 +331,11 @@ public:
 	ComplexOrRealType manyPoint(const VectorSizeType& sites,
 	                            const PsimagLite::Vector<LabeledOperatorType>::Type& what,
 	                            const VectorSizeType& spins,
-	                            const VectorSizeType& orbs) const
+	                            const VectorSizeType& orbs,
+	                            const PairType& braAndKet) const
 	{
-		VectorType tmpVector = gsVector_;
+		checkBraOrKet("ket", braAndKet.second);
+		VectorType tmpVector = vectors_[braAndKet.second];
 		const BasisType* basisOld = &(model_.basis());
 		RealType isign = 1.0;
 		PairType oldParts = model_.basis().parts();
@@ -342,7 +372,10 @@ public:
 		oldParts = model_.basis().parts();
 		if (oldParts != newParts) return 0.0;
 
-		return gsVector_*tmpVector;
+		checkBraOrKet("bra", braAndKet.first);
+		const VectorType& braVector = vectors_[braAndKet.first];
+
+		return braVector*tmpVector;
 	}
 
 private:
@@ -523,50 +556,57 @@ private:
 		                  isign);
 	}
 
-	void computeGroundState()
+	void computeAllStatesBelow(SizeType excited)
 	{
+		const SizeType excitedPlusOne = excited + 1;
 		SpecialSymmetryType rs(model_.basis(),model_.geometry(),options_);
 		InternalProductType hamiltonian(model_,rs);
 		ParametersForSolverType params(io_,"Lanczos");
 		LanczosSolverType lanczosSolver(hamiltonian,params);
 
-		gsEnergy_ = 1e10;
 		SizeType offset = model_.size();
 		SizeType currentOffset = 0;
 
-		for (SizeType i=0;i<rs.sectors();i++) {
+		for (SizeType i = 0; i < rs.sectors(); ++i) {
 			hamiltonian.specialSymmetrySector(i);
 			SizeType n = hamiltonian.rows();
 			if (n == 0) continue;
 			VectorType initial(n);
 			PsimagLite::fillRandom(initial);
-			VectorType gsVector1(n);
-			RealType gsEnergy1 = 0;
+			VectorVectorType zs(excitedPlusOne, VectorType(n));
+			VectorRealType eigs(excitedPlusOne);
 
 			try {
-				lanczosSolver.computeOneState(gsEnergy1, gsVector1, initial, 0);
+				lanczosSolver.computeAllStatesBelow(eigs, zs, initial, excitedPlusOne);
 			} catch (std::exception&) {
 
 				std::cerr<<"Engine: Lanczos Solver failed ";
 				std::cerr<<" trying exact diagonalization...\n";
-				VectorRealType eigs(hamiltonian.rows());
+				VectorRealType eigs2(n);
 				MatrixType fm;
-				hamiltonian.fullDiag(eigs,fm);
-				for (SizeType j = 0; j < eigs.size(); ++j)
-					gsVector1[j] = fm(j,0);
-				gsEnergy1 = eigs[0];
-				std::cout<<"Found lowest eigenvalue= "<<gsEnergy1<<"\n";
+				hamiltonian.fullDiag(eigs2, fm);
+				for (SizeType k = 0; k < excitedPlusOne; ++k) {
+					for (SizeType j = 0; j < n; ++j)
+						zs[k][j] = fm(j, k);
+					eigs[k] = eigs2[k];
+				}
 			}
 
-			if (gsEnergy1<gsEnergy_) {
-				gsVector_=gsVector1;
-				gsEnergy_=gsEnergy1;
+			if (eigs[0] < energies_[0]) {
+				for (SizeType j = 0; j < excitedPlusOne; ++j) {
+					vectors_[j] = zs[j];
+					energies_[j] = eigs[j];
+				}
+
 				offset = currentOffset;
 			}
-			currentOffset +=  gsVector1.size();
+
+			currentOffset +=  zs[0].size();
 		}
-		rs.transformGs(gsVector_,offset);
-		std::cout<<"#GSNorm="<<PsimagLite::real(gsVector_*gsVector_)<<"\n";
+
+		rs.transform(vectors_, offset);
+
+		printEnergiesAndNorms();
 	}
 
 	template<typename ContinuedFractionType>
@@ -596,15 +636,34 @@ private:
 		RealType diagonalFactor = (isDiagonal) ? 1 : 0.5;
 		s2 *= diagonalFactor;
 
-		cf.set(ab, gsEnergy_, PsimagLite::real(weight*s2), s);
+		assert(0 < energies_.size());
+		const RealType gsEnergy = energies_[0];
+		cf.set(ab, gsEnergy, PsimagLite::real(weight*s2), s);
+	}
+
+	void checkBraOrKet(PsimagLite::String name, SizeType ind) const
+	{
+		if (ind < vectors_.size()) return;
+
+		err("Wrong " + name + " FATAL ERROR\n");
+	}
+
+	void printEnergiesAndNorms() const
+	{
+		const SizeType excited = energies_.size();
+		assert(excited == vectors_.size());
+		for (SizeType i = 0; i < excited; ++i) {
+			const RealType val = PsimagLite::real(vectors_[i]*vectors_[i]);
+			std::cout<<"E["<<i<<"]="<<energies_[i]<<" norm="<<val<<"\n";
+		}
 	}
 
 	const ModelType& model_;
 	PsimagLite::ProgressIndicator progress_;
 	InputType& io_;
 	PsimagLite::String options_;
-	RealType gsEnergy_;
-	VectorType gsVector_;
+	VectorRealType energies_;
+	VectorVectorType vectors_;
 }; // class ContinuedFraction
 } // namespace Dmrg
 

src/Engine/LanczosDriver.h

@@ -39,7 +39,9 @@ template<typename ModelType>
 SizeType maxOrbitals(const ModelType& model);
 
 template<typename EngineType>
-void extendedStatic(PsimagLite::String manypoint, const EngineType& engine);
+void extendedStatic(PsimagLite::String manypoint,
+                    const EngineType& engine,
+                    const typename EngineType::PairType& braAndKet);
 
 template<typename ModelType,
          typename SpecialSymmetryType,

src/Engine/LanczosDriver1.h

@@ -15,7 +15,9 @@ SizeType maxOrbitals(const ModelType& model)
 }
 
 template<typename EngineType>
-void extendedStatic(PsimagLite::String manypoint, const EngineType& engine)
+void extendedStatic(PsimagLite::String manypoint,
+                    const EngineType& engine,
+                    const typename EngineType::PairType& braAndKet)
 {
 	typedef typename EngineType::VectorSizeType VectorSizeType;
 	PsimagLite::Vector<PsimagLite::String>::Type str;
@@ -38,7 +40,7 @@ void extendedStatic(PsimagLite::String manypoint, const EngineType& engine)
 	}
 
 	std::cout<<"<gs|"<<manypoint<<"|gs>=";
-	std::cout<<engine.manyPoint(sites,whats,spins,orbs);
+	std::cout<<engine.manyPoint(sites, whats, spins, orbs, braAndKet);
 	std::cout<<"\n";
 }
 
@@ -56,10 +58,10 @@ void mainLoop3(const ModelType& model,
 	typedef PsimagLite::Vector<PsimagLite::String>::Type VectorStringType;
 
 	const GeometryType& geometry = model.geometry();
-	EngineType engine(model,geometry.numberOfSites(),io);
+	EngineType engine(model, io);
 
 	//! get the g.s.:
-	RealType Eg = engine.gsEnergy();
+	RealType Eg = engine.energies(0);
 	std::cout.precision(8);
 	std::cout<<"Energy="<<Eg<<"\n";
 	PsimagLite::String filename = PsimagLite::basenameOf(io.filename());
@@ -72,10 +74,7 @@ void mainLoop3(const ModelType& model,
 		for (SizeType j = 0; j < ntokens; ++j) {
 			VectorStringType braOpKet;
 			PsimagLite::split(braOpKet, tokens[j], "|");
-			const SizeType ind = (braOpKet.size() == 3) ? 1 : 0;
-			if (braOpKet.size() != 1 && ind == 0)
-				err("Wrong braket\n");
-			engine.measure(braOpKet[ind]);
+			engine.measure(braOpKet);
 		}
 	}
 
@@ -133,7 +132,8 @@ void mainLoop3(const ModelType& model,
 				engine.twoPoint(cicjMatrix,
 				                lanczosOptions.cicj[cicji],
 				                lanczosOptions.spins,
-				                std::pair<SizeType,SizeType>(orb1,orb2));
+				                std::pair<SizeType,SizeType>(orb1,orb2),
+				                std::pair<SizeType,SizeType>(0, 0));
 				std::cout<<cicjMatrix;
 			}
 		}
@@ -141,7 +141,7 @@ void mainLoop3(const ModelType& model,
 
 	if (lanczosOptions.split >= 0) {
 		LanczosPlusPlus::ReducedDensityMatrix<ModelType> reducedDensityMatrix(model,
-		                                                                      engine.eigenvector(),
+		                                                                      engine.eigenvector(0),
 		                                                                      lanczosOptions.split);
 		reducedDensityMatrix.printAll(std::cout);
 	}
@@ -150,7 +150,7 @@ void mainLoop3(const ModelType& model,
 		PsimagLite::Vector<PsimagLite::String>::Type str;
 		PsimagLite::split(str, lanczosOptions.extendedStatic, ",");
 		for (SizeType i = 0; i < str.size(); ++i)
-			extendedStatic(str[i],engine);
+			extendedStatic(str[i], engine, typename EngineType::PairType(0, 0));
 	}
 }
 

src/Engine/ProgramGlobals.h

@@ -84,6 +84,7 @@ DISCLOSED WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS.
 #include <climits>
 #include "Vector.h"
 #include "BitManip.h"
+#include "CrsMatrix.h"
 
 namespace LanczosPlusPlus {
 
@@ -185,6 +186,26 @@ struct ProgramGlobals {
 		os<<"\n--------------\n";
 	}
 
+	template<typename SomeVectorType>
+	static typename PsimagLite::EnableIf<PsimagLite::IsVectorLike<SomeVectorType>::True,
+	void>::Type
+	transform(SomeVectorType& gs,
+	          SizeType offset,
+	          SomeVectorType& gstmp,
+	          const PsimagLite::CrsMatrix<typename SomeVectorType::value_type>& tr)
+	{
+		for (SizeType i=0;i<gs.size();i++) {
+			assert(i+offset<gstmp.size());
+			gstmp[i+offset]=gs[i];
+		}
+
+		PsimagLite::CrsMatrix<typename SomeVectorType::value_type> rT;
+		transposeConjugate(rT, tr);
+		gs.clear();
+		gs.resize(tr.rows());
+		multiply(gs,rT,gstmp);
+	}
+
 }; // ProgramGlobals
 
 } // namespace LanczosPlusPlus

src/Engine/ReflectionSymmetry.h

@@ -22,6 +22,7 @@ Please see full open source license included in file LICENSE.
 #include "ProgressIndicator.h"
 #include "CrsMatrix.h"
 #include "Vector.h"
+#include "ProgramGlobals.h"
 
 namespace LanczosPlusPlus {
 
@@ -60,6 +61,7 @@ public:
 	typedef PsimagLite::CrsMatrix<ComplexOrRealType> SparseMatrixType;
 	typedef typename PsimagLite::Vector<ComplexOrRealType>::Type VectorType;
 	typedef typename PsimagLite::Vector<RealType>::Type VectorRealType;
+	typedef typename PsimagLite::Vector<VectorType>::Type VectorVectorType;
 
 	ReflectionSymmetry(const BasisType& basis,
 	                   const GeometryType& geometry,
@@ -150,19 +152,14 @@ public:
 		split(matrix1[0],matrix1[1],matrix2);
 	}
 
-	void transformGs(VectorType& gs,SizeType offset)
+	void transform(VectorVectorType& zs,SizeType offset)
 	{
-		VectorType gstmp(transform_.rows(),0);
+		VectorType gstmp(transform_.rows());
 
-		for (SizeType i=0;i<gs.size();i++) {
-			assert(i+offset<gstmp.size());
-			gstmp[i+offset]=gs[i];
-		}
-		SparseMatrixType rT;
-		transposeConjugate(rT,transform_);
-		gs.clear();
-		gs.resize(transform_.rows());
-		multiply(gs,rT,gstmp);
+		const SizeType excitedPlusOne = zs.size();
+
+		for (SizeType i = 0; i < excitedPlusOne; ++i)
+			ProgramGlobals::transform(zs[i], offset, gstmp, transform_);
 	}
 
 	SizeType sectors() const { return 2; }

src/Engine/TranslationSymmetry.h

@@ -23,6 +23,7 @@ Please see full open source license included in file LICENSE.
 #include "CrsMatrix.h"
 #include "Vector.h"
 #include "SparseVector.h"
+#include "ProgramGlobals.h"
 
 namespace LanczosPlusPlus {
 
@@ -184,6 +185,7 @@ public:
 	typedef PsimagLite::CrsMatrix<ComplexOrRealType> SparseMatrixType;
 	typedef typename PsimagLite::Vector<ComplexOrRealType>::Type VectorType;
 	typedef typename PsimagLite::Vector<RealType>::Type VectorRealType;
+	typedef typename PsimagLite::Vector<VectorType>::Type VectorVectorType;
 
 	TranslationSymmetry(const BasisType& basis,
 	                    const GeometryType& geometry,
@@ -265,19 +267,15 @@ public:
 		split(matrix1,matrix2);
 	}
 
-	void transformGs(VectorType& gs,SizeType offset)
+	void transform(VectorVectorType& zs, SizeType offset)
 	{
 		VectorType gstmp(transform_.rows(),0);
 
-		for (SizeType i=0;i<gs.size();i++) {
-			assert(i+offset<gstmp.size());
-			gstmp[i+offset]=gs[i];
+		const SizeType excitedPlusOne = zs.size();
+
+		for (SizeType i = 0; i < excitedPlusOne; ++i) {
+			ProgramGlobals::transform(zs[i], offset, gstmp, transform_);
 		}
-		SparseMatrixType rT;
-		transposeConjugate(rT,transform_);
-		gs.clear();
-		gs.resize(transform_.rows());
-		multiply(gs,rT,gstmp);
 	}
 
 	SizeType sectors() const { return kspace_.size(); }
@@ -394,6 +392,20 @@ private:
 		}
 	}
 
+	void transform(VectorType& gs,SizeType offset, VectorType& gstmp)
+	{
+		for (SizeType i=0;i<gs.size();i++) {
+			assert(i+offset<gstmp.size());
+			gstmp[i+offset]=gs[i];
+		}
+
+		SparseMatrixType rT;
+		transposeConjugate(rT,transform_);
+		gs.clear();
+		gs.resize(transform_.rows());
+		multiply(gs,rT,gstmp);
+	}
+
 	PsimagLite::ProgressIndicator progress_;
 	SparseMatrixType transform_;
 	KspaceType kspace_;

src/Version.h

 #ifndef LANCZOSPP_VERSION
-#define LANCZOSPP_VERSION "1.68"
+#define LANCZOSPP_VERSION "1.69"
 #endif