Contains the changed need to make the Response theory work

Author: sanshar

BaseOperator.C

@@ -213,6 +213,19 @@ void SparseMatrix::Randomise ()
 	    SpinAdapted::Randomise(operator_element(lQ, rQ));
 }
 
+double trace(const SparseMatrix& lhs)
+{
+  assert(lhs.nrows() == lhs.ncols());
+  double trace = 0.0;
+
+  for(int lQ=0;lQ<lhs.nrows();++lQ)
+    if(lhs.allowed(lQ,lQ))
+      for(int i=0;i<(lhs)(lQ,lQ).Nrows();++i)
+	trace += (lhs)(lQ,lQ)(i+1,i+1);
+
+  return trace;
+}
+
 double DotProduct(const SparseMatrix& lhs, const SparseMatrix& rhs)
 {
   double result = 0.;

BaseOperator.h

@@ -283,6 +283,7 @@ double getCommuteParity(SpinQuantum a, SpinQuantum b, SpinQuantum c);
 void Normalise(SparseMatrix& a, int* success = 0);
 void ScaleAdd(double d, const SparseMatrix& a, SparseMatrix& b);
 double DotProduct(const SparseMatrix& lhs, const SparseMatrix& rhs);
+double trace(const SparseMatrix& lhs);
 void Scale(double d, SparseMatrix& a);
 void assignloopblock(SpinBlock*& loopblock, SpinBlock*& otherblock, SpinBlock* leftSpinBlock, SpinBlock* rightSpinBlock);
 void copy(const ObjectMatrix<Matrix>& a, ObjectMatrix<Matrix>& b);

Operators.C

@@ -625,14 +625,14 @@ void SpinAdapted::DesCre::build(const SpinBlock& b)
   {
     const boost::shared_ptr<SparseMatrix> op1 = leftBlock->get_op_rep(DES, -getSpinQuantum(i), i);
     const boost::shared_ptr<SparseMatrix> op2 = rightBlock->get_op_rep(CRE, getSpinQuantum(j), j);
-    SpinAdapted::operatorfunctions::TensorProduct(rightBlock, *op2, *op1, &b, &(b.get_stateInfo()), *this, 1.0);
+    double parity = getCommuteParity(op1->get_deltaQuantum()[0], op2->get_deltaQuantum()[0], get_deltaQuantum()[0]);
+    SpinAdapted::operatorfunctions::TensorProduct(rightBlock, *op2, *op1, &b, &(b.get_stateInfo()), *this, 1.0*parity);
   }
   else if (rightBlock->get_op_array(DES).has(i))
   {
     const boost::shared_ptr<SparseMatrix> op1 = rightBlock->get_op_rep(DES, -getSpinQuantum(i), i);
     const boost::shared_ptr<SparseMatrix> op2 = leftBlock->get_op_rep(CRE, getSpinQuantum(j), j);
-    double parity = getCommuteParity(op1->get_deltaQuantum()[0], op2->get_deltaQuantum()[0], get_deltaQuantum()[0]);
-    SpinAdapted::operatorfunctions::TensorProduct(rightBlock, *op1, *op2, &b, &(b.get_stateInfo()), *this, 1.0*parity);
+    SpinAdapted::operatorfunctions::TensorProduct(rightBlock, *op1, *op2, &b, &(b.get_stateInfo()), *this, 1.0);
   }
   else
     abort();  
@@ -651,15 +651,15 @@ double SpinAdapted::DesCre::redMatrixElement(Csf c1, vector<Csf>& ladder, const
 
   TensorOp D(I, -1), C(J, 1);
 
-  TensorOp DC = D.product(C, spin, irrep);
+  TensorOp CD = C.product(D, spin, irrep);
 
 
   for (int j = 0; j < deltaQuantum.size(); ++j) {
     for (int i=0; i<ladder.size(); i++)
     {
       int index = 0; double cleb=0.0;
       if (nonZeroTensorComponent(c1, deltaQuantum[j], ladder[i], index, cleb)) {
-        std::vector<double> MatElements = calcMatrixElements(c1, DC, ladder[i]) ;
+        std::vector<double> MatElements = calcMatrixElements(c1, CD, ladder[i]) ;
         element += MatElements[index]/cleb;
         break;
       }

OverlapHelement.c

@@ -14,22 +14,30 @@ int main(int argc, char* argv []) {
 #endif
 
   ReadInputFromC(argv[1], -1);
+
   int mpsstate=0;
 
   readMPSFromDiskAndInitializeStaticVariables(mpsstate);
+
+
   test();
+  exit(0);
+  
 
   double overlap, hvalue;
 
-  long temp=1, occ;
+  long temp=1, occ=0;
 
-  //occ = temp<<63 | temp<<62 | temp<<61 | temp<<60 | temp<<59 | temp<<58 ;
-  occ = temp<<63 | temp<<62 | temp<<61 | temp<<60 | temp<<59 | temp<<58 | temp<<57 | temp<<56 | temp<<55 | temp<<54 ;
+  //11 00 11 11  11 11
+  occ = temp<<63 | temp<<62 | temp<<61 | temp<<60;// | temp<<57 | temp<<56 | temp<<55 | temp<<54 | temp<<53 | temp<<52;
+
+  //11 11 11 00  11 11
+  //occ = temp<<63 | temp<<62 | temp<<61 | temp<<60 | temp<<59 | temp<<58  | temp<<55 | temp<<54 | temp<<53 | temp<<52;
 
   evaluateOverlapAndHamiltonian(&occ, 1, &overlap, &hvalue);
 
-  printf("overlap = %10.5e  %i\n", overlap, rank);
-  printf("helement = %10.5e  %i\n", hvalue, rank);
+  printf("overlap = %15.8e  %i\n", overlap, rank);
+  printf("helement = %15.8e  %i\n", hvalue, rank);
 
 #ifndef SERIAL
   MPI_Finalize();

davidson.C

@@ -21,4 +21,16 @@ void SpinAdapted::multiply_h::operator()(Wavefunction& c, Wavefunction& v)
   block.multiplyH( c, &v, MAX_THRD);
 }
 
+SpinAdapted::multiply_h_e::multiply_h_e(const SpinBlock& b, const bool &onedot_, double E0_) : block(b), E0(E0_){}
+
+void SpinAdapted::multiply_h_e::operator()(Wavefunction& c, Wavefunction& v)
+{
+  Wavefunction oi = v; oi.Clear();
+  block.multiplyH( c, &v, MAX_THRD);
+  block.multiplyOverlap(c, &oi, MAX_THRD);
+  ScaleAdd(-E0, oi, v);
+}
+
+
+
 

davidson.h

@@ -27,6 +27,20 @@ class multiply_h : public Davidson_functor
   void operator()(Wavefunction& c, Wavefunction& v);
   const SpinBlock& get_block() {return block;}
 };
+
+//H-E
+class multiply_h_e : public Davidson_functor
+{
+ private:
+  const SpinBlock& block;
+  double E0;
+ public:
+  multiply_h_e(const SpinBlock& b, const bool &onedot_, double E0);
+  void operator()(Wavefunction& c, Wavefunction& v);
+  const SpinBlock& get_block() {return block;}
+};
+
+
 }
 
 #endif