RahulMethod for observables

src/Engine/Engine.h

@@ -30,6 +30,7 @@ Please see full open source license included in file LICENSE.
 #include "DefaultSymmetry.h"
 #include "TypeToString.h"
 #include "LabeledOperator.h"
+#include "OneOperatorSpec.h"
 
 namespace LanczosPlusPlus {
 template<typename ModelType_,
@@ -68,6 +69,9 @@ public:
 	typedef typename PsimagLite::Vector<RealType>::Type VectorRealType;
 	typedef typename PsimagLite::Vector<PsimagLite::String>::Type VectorStringType;
 	typedef LabeledOperator LabeledOperatorType;
+	typedef PsimagLite::OneOperatorSpec OneOperatorSpecType;
+	typedef typename ModelType::RahulOperatorType RahulOperatorType;
+	typedef typename ModelType::VectorRahulOperatorType VectorRahulOperatorType;
 
 	// ContF needs to support concurrency FIXME
 	static const SizeType parallelRank_ = 0;
@@ -186,6 +190,31 @@ public:
 		}
 	}
 
+	void measure(PsimagLite::String meas) const
+	{
+		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]);
+			if (site < 0)
+				err("Operator " + tokens[i] + " needs a site in brackets\n");
+			assert(i < vsites.size());
+			vsites[i] = site;
+
+			OneOperatorSpecType opspec(tokens[i]);
+
+			vops.push_back(RahulOperatorType(opspec.label, opspec.dof, opspec.transpose));
+		}
+
+		model_.rahulMethod(psiNew, vops, vsites, gsVector_);
+		const ComplexOrRealType result = gsVector_*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,

src/Engine/LanczosDriver1.h

@@ -63,6 +63,10 @@ void mainLoop3(const ModelType& model,
 	std::cout<<"Energy="<<Eg<<"\n";
 	PsimagLite::String filename = PsimagLite::basenameOf(io.filename());
 
+	const SizeType nmeas = lanczosOptions.measure.size();
+	for (SizeType i = 0; i < nmeas; ++i)
+		engine.measure(lanczosOptions.measure[i]);
+
 	for (SizeType gfi=0;gfi<lanczosOptions.gf.size();gfi++) {
 		SizeType counter = 0;
 		while (true) {

src/Engine/LanczosOptions.h

@@ -19,6 +19,7 @@ struct LanczosOptions {
 	PsimagLite::Vector<SizeType>::Type sites;
 	PsimagLite::Vector<PairSizeType>::Type spins;
 	PsimagLite::String extendedStatic;
+	PsimagLite::Vector<PsimagLite::String>::Type measure;
 
 }; // struct LanczosOptions
 }

src/Engine/ModelBase.h

@@ -23,6 +23,7 @@ Please see full open source license included in file LICENSE.
 #include "CrsMatrix.h"
 #include "BasisBase.h"
 #include "Vector.h"
+#include "RahulOperator.h"
 
 namespace LanczosPlusPlus {
 
@@ -40,6 +41,8 @@ public:
 	typedef PsimagLite::CrsMatrix<ComplexOrRealType> SparseMatrixType;
 	typedef typename PsimagLite::Vector<ComplexOrRealType>::Type VectorType;
 	typedef typename PsimagLite::Vector<SizeType>::Type VectorSizeType;
+	typedef RahulOperator<ComplexOrRealType> RahulOperatorType;
+	typedef typename PsimagLite::Vector<RahulOperatorType>::Type VectorRahulOperatorType;
 
 	virtual ~ModelBase() {}
 
@@ -82,6 +85,14 @@ public:
 
 	virtual void print(std::ostream& os) const = 0;
 
+	virtual void rahulMethod(VectorType&,
+	                         const VectorRahulOperatorType&,
+	                         const VectorSizeType&,
+	                         const VectorType&) const
+	{
+		throw PsimagLite::RuntimeError("ModelBase::rahulMethod() not impl. for this model\n");
+	}
+
 	virtual void printOperators(std::ostream&) const
 	{
 		PsimagLite::String str(__FILE__);

src/Models/HubbardOneOrbital/HubbardHelper.h

@@ -20,6 +20,9 @@ public:
 	typedef typename ModelBaseType::ComplexOrRealType ComplexOrRealType;
 	typedef typename ModelBaseType::GeometryType GeometryType;
 	typedef typename PsimagLite::Matrix<ComplexOrRealType> MatrixType;
+	typedef PsimagLite::Vector<SizeType>::Type VectorSizeType;
+	typedef typename ModelBaseType::RahulOperatorType RahulOperatorType;
+	typedef typename ModelBaseType::VectorRahulOperatorType VectorRahulOperatorType;
 
 	enum {SPIN_UP = ProgramGlobals::SPIN_UP, SPIN_DOWN = ProgramGlobals::SPIN_DOWN};
 
@@ -32,7 +35,7 @@ public:
 	    : geometry_(geometry),
 	      mp_(mp),
 	      hoppings_(geometry.numberOfSites(), geometry.numberOfSites()),
-		  hasJcoupling_(mp_.model == "SuperHubbardExtended"),
+	      hasJcoupling_(mp_.model == "SuperHubbardExtended"),
 	      hasCoulombCoupling_(mp_.model == "HubbardOneBandExtended" ||
 	                          mp_.model == "SuperHubbardExtended"),
 	      hasRashba_(mp_.model == "HubbardOneBandRashbaSOC")
@@ -123,8 +126,61 @@ public:
 		}
 	}
 
+	void rahulMethod(VectorType& psiNew,
+	                 const VectorRahulOperatorType& vops,
+	                 const VectorSizeType& vsites,
+	                 const VectorType& psi,
+	                 const BasisBaseType& basis) const
+	{
+		std::fill(psiNew.begin(), psiNew.end(), 0.0);
+		const SizeType hilbert = basis.size();
+		const SizeType nops = vops.size();
+		for (SizeType ispace = 0; ispace < hilbert; ++ispace) {
+			WordType ket1 = basis(ispace, SPIN_UP);
+			WordType ket2 = basis(ispace, SPIN_DOWN);
+			assert(ispace < psi.size());
+			ComplexOrRealType value = psi[ispace];
+			WordType ketp1 = ket1;
+			WordType ketp2 = ket2;
+			bool nonZero = false;
+			for (SizeType j = 0; j < nops; ++j) {
+				const RahulOperatorType& op = vops[j];
+				assert(j < vsites.size());
+				const SizeType site = vsites[j];
+
+				ComplexOrRealType result = 0;
+				nonZero = (op.dof() == SPIN_UP)
+				        ? applyOperator(ketp1, result, op, site)
+				        : applyOperator(ketp2, result, op, site);
+				if (!nonZero) break;
+				value *= result;
+			}
+
+			if (!nonZero) continue;
+			SizeType newI = basis.perfectIndex(ketp1, ketp2);
+			assert(newI < psiNew.size());
+			psiNew[newI] += value;
+		}
+	}
+
 private:
 
+	bool applyOperator(WordType& ketp,
+	                   ComplexOrRealType& result,
+	                   const RahulOperatorType& op,
+	                   SizeType site) const
+	{
+		const WordType ket = ketp;
+		const WordType mask = BasisType::bitmask(site);
+		WordType s = (ket & mask);
+		bool sbit = (s > 0);
+		bool sbitSaved = sbit;
+		bool nonZero = op.actOn(sbit, result);
+		if (!nonZero) return false;
+		if (sbitSaved != sbit) ketp ^= mask;
+		return true;
+	}
+
 	void calcDiagonalElements(typename PsimagLite::Vector<RealType>::Type& diag,
 	                          const BasisBaseType& basis) const
 	{
@@ -358,7 +414,7 @@ private:
 	MatrixType hoppings_;
 	MatrixType rashbaHoppings_;
 	bool hasJcoupling_;
-    bool hasCoulombCoupling_;
+	bool hasCoulombCoupling_;
 	bool hasRashba_;
 };
 

src/Models/HubbardOneOrbital/HubbardOneOrbital.h

@@ -34,6 +34,7 @@ public:
 	typedef typename BaseType::VectorType VectorType;
 	typedef typename BasisType::LabeledOperatorType LabeledOperatorType;
 	typedef HubbardHelper<BaseType, BasisType, ParametersModelType> HubbardHelperType;
+	typedef typename HubbardHelperType::VectorRahulOperatorType VectorRahulOperatorType;
 
 	static int const FERMION_SIGN = BasisType::FERMION_SIGN;
 
@@ -87,7 +88,7 @@ public:
 	                 const std::pair<SizeType,SizeType>& oldParts,
 	                 const LabeledOperator& lOperator,
 	                 SizeType spin,
-	                 SizeType orb) const
+	                 SizeType) const
 	{
 		if (lOperator.id() == LabeledOperator::Label::OPERATOR_C ||
 		        lOperator.id() == LabeledOperator::Label::OPERATOR_CDAGGER)
@@ -132,6 +133,14 @@ public:
 			printOperatorC(site,spin,os);
 	}
 
+	void rahulMethod(VectorType& psiNew,
+	                 const VectorRahulOperatorType& vops,
+	                 const VectorSizeType& vsites,
+	                 const VectorType& psi) const
+	{
+		return helper_.rahulMethod(psiNew, vops, vsites, psi, basis_);
+	}
+
 private:
 
 	void printOperatorC(SizeType site, SizeType spin, std::ostream& os) const

src/Models/HubbardOneOrbitalRashbaSOC/BasisRashbaSOC.h

@@ -53,7 +53,7 @@ public:
 
 	PairIntType parts() const
 	{
-		throw PsimagLite::RuntimeError("parts() unimplemented\n");
+		return PairIntType(0, 0); // dummy, we don't care
 	}
 
 	static const WordType& bitmask(SizeType i)

src/Models/HubbardOneOrbitalRashbaSOC/HubbardOneOrbitalRashbaSOC.h

@@ -34,6 +34,7 @@ public:
 	typedef typename BasisType::LabeledOperatorType LabeledOperatorType;
 	typedef PsimagLite::SparseRow<SparseMatrixType> SparseRowType;
 	typedef HubbardHelper<BaseType, BasisType, ParametersModelType> HubbardHelperType;
+	typedef typename HubbardHelperType::VectorRahulOperatorType VectorRahulOperatorType;
 
 	static int const FERMION_SIGN = BasisType::FERMION_SIGN;
 
@@ -84,13 +85,13 @@ public:
 		helper_.matrixVectorProduct(x, y, basis);
 	}
 
-	bool hasNewParts(std::pair<SizeType,SizeType>& newParts,
-	                 const std::pair<SizeType,SizeType>& oldParts,
-	                 const LabeledOperator& lOperator,
-	                 SizeType spin,
-	                 SizeType orb) const
+	bool hasNewParts(std::pair<SizeType,SizeType>&,
+	                 const std::pair<SizeType,SizeType>&,
+	                 const LabeledOperator&,
+	                 SizeType,
+	                 SizeType) const
 	{
-		throw std::runtime_error("hasNewParts unimplemented\n");
+		return false;
 	}
 
 	const GeometryType& geometry() const { return geometry_; }
@@ -113,6 +114,14 @@ public:
 		os<<"WARNING: printOperators unimplemented\n";
 	}
 
+	void rahulMethod(VectorType& psiNew,
+	                 const VectorRahulOperatorType& vops,
+	                 const VectorSizeType& vsites,
+	                 const VectorType& psi) const
+	{
+		return helper_.rahulMethod(psiNew, vops, vsites, psi, basis_);
+	}
+
 private:
 
 	const ParametersModelType mp_;

src/lanczos.cpp

@@ -127,7 +127,7 @@ int main(int argc,char **argv)
 	\item[-V] prints version and exits.
 	\end{itemize}
 	*/
-	while ((opt = getopt(argc, argv, "g:c:f:s:r:p:S:V")) != -1) {
+	while ((opt = getopt(argc, argv, "g:c:m:f:s:r:p:S:V")) != -1) {
 		switch (opt) {
 		case 'g':
 			lanczosOptions.gf.push_back(LabeledOperator(optarg));
@@ -138,6 +138,9 @@ int main(int argc,char **argv)
 		case 'c':
 			lanczosOptions.cicj.push_back(LabeledOperator(optarg));
 			break;
+		case 'm':
+			lanczosOptions.measure.push_back(optarg);
+			break;
 		case 's':
 			lanczosOptions.spins.clear();
 			PsimagLite::split(str, optarg, ";");