Auxiliary fields

lib/qcd/utils/CovariantLaplacian.h

@@ -144,6 +144,25 @@ class LaplacianAdjointField: public Metric<typename Impl::Field> {
     } 
   }
 
+  // separating this temporarily
+  void MDeriv(const GaugeField& left, const GaugeField& right,
+              GaugeField& der) {
+    // in is anti-hermitian
+    RealD factor = -kappa / (double(4 * Nd));
+
+    for (int mu = 0; mu < Nd; mu++) {
+      GaugeLinkField der_mu(der._grid);
+      der_mu = zero;
+      for (int nu = 0; nu < Nd; nu++) {
+        GaugeLinkField left_nu = PeekIndex<LorentzIndex>(left, nu);
+        GaugeLinkField right_nu = PeekIndex<LorentzIndex>(right, nu);
+        der_mu += U[mu] * Cshift(left_nu, mu, 1) * adj(U[mu]) * right_nu;
+        der_mu += U[mu] * Cshift(right_nu, mu, 1) * adj(U[mu]) * left_nu;
+      }
+      PokeIndex<LorentzIndex>(der, -factor * der_mu, mu);
+    }
+  }
+
   void Minv(const GaugeField& in, GaugeField& inverted){
     HermitianLinearOperator<LaplacianAdjointField<Impl>,GaugeField> HermOp(*this);
     Solver(HermOp, in, inverted);
@@ -157,6 +176,14 @@ class LaplacianAdjointField: public Metric<typename Impl::Field> {
     P = Gp; 
   }
 
+  void MInvSquareRoot(GaugeField& P){
+    GaugeField Gp(P._grid);
+    HermitianLinearOperator<LaplacianAdjointField<Impl>,GaugeField> HermOp(*this);
+    ConjugateGradientMultiShift<GaugeField> msCG(param.MaxIter,PowerInvHalf);
+    msCG(HermOp,P,Gp);
+    P = Gp; 
+  }
+
 
 
  private:

lib/qcd/utils/Metric.h

@@ -40,7 +40,9 @@ public:
   virtual void M(const Field&, Field&)     = 0;
   virtual void Minv(const Field&, Field&)  = 0;
   virtual void MSquareRoot(Field&) = 0;
+  virtual void MInvSquareRoot(Field&) = 0;
   virtual void MDeriv(const Field&, Field&) = 0;
+  virtual void MDeriv(const Field&, const Field&, Field&) = 0;
 };
 
 
@@ -58,9 +60,15 @@ public:
   virtual void MSquareRoot(Field& P){
     // do nothing
   }
+  virtual void MInvSquareRoot(Field& P){
+    // do nothing
+  }
   virtual void MDeriv(const Field& in, Field& out){
     out = zero;
   }
+  virtual void MDeriv(const Field& left, const Field& right, Field& out){
+    out = zero;
+  }
 
 };
 
@@ -76,18 +84,37 @@ public:
   Metric<MomentaField>& M;
   MomentaField Mom;
 
-  GeneralisedMomenta(GridBase* grid, Metric<MomentaField>& M): M(M), Mom(grid){}
+  // Auxiliary fields
+  // not hard coded but inherit the type from the metric
+  // created Nd new fields
+  // hide these in the metric?
+  //typedef Lattice<iVector<iScalar<iMatrix<vComplex, Nc> >, Nd/2 > > AuxiliaryMomentaType;
+  MomentaField AuxMom;
+  MomentaField AuxField;
+
+  GeneralisedMomenta(GridBase* grid, Metric<MomentaField>& M): M(M), Mom(grid), AuxMom(grid), AuxField(grid){}
 
   // Correct
   void MomentaDistribution(GridParallelRNG& pRNG){
     // Generate a distribution for
-    // 1/2 P^dag G P
+    // P^dag G P
     // where G = M^-1
 
     // Generate gaussian momenta
     Implementation::generate_momenta(Mom, pRNG);
     // Modify the distribution with the metric
     M.MSquareRoot(Mom);
+
+    if (0) {
+      // Auxiliary momenta
+      // do nothing if trivial, so hide in the metric
+      MomentaField AuxMomTemp(Mom._grid);
+      Implementation::generate_momenta(AuxMom, pRNG);
+      Implementation::generate_momenta(AuxField, pRNG);
+      // Modify the distribution with the metric
+      // Aux^dag M Aux
+      M.MInvSquareRoot(AuxMom);  // AuxMom = M^{-1/2} AuxMomTemp
+    }
   }
 
   // Correct
@@ -99,17 +126,34 @@ public:
     Hloc = zero;
     for (int mu = 0; mu < Nd; mu++) {
       // This is not very general
-      // hide in the operators
+      // hide in the metric
       auto Mom_mu = PeekIndex<LorentzIndex>(Mom, mu);
       auto inv_mu = PeekIndex<LorentzIndex>(inv, mu);
       Hloc += trace(Mom_mu * inv_mu);
     }
+
+    if (0) {
+      // Auxiliary Fields
+      // hide in the metric
+      M.M(AuxMom, inv);
+      for (int mu = 0; mu < Nd; mu++) {
+        // This is not very general
+        // hide in the operators
+        auto inv_mu = PeekIndex<LorentzIndex>(inv, mu);
+        auto am_mu = PeekIndex<LorentzIndex>(AuxMom, mu);
+        auto af_mu = PeekIndex<LorentzIndex>(AuxField, mu);
+        Hloc += trace(am_mu * inv_mu);// p M p
+        Hloc += trace(af_mu * af_mu);
+      }
+    }
+
     Complex Hsum = sum(Hloc);
     return Hsum.real();
   }
 
   // Correct
   void DerivativeU(MomentaField& in, MomentaField& der){
+
     // Compute the derivative of the kinetic term
     // with respect to the gauge field
     MomentaField MDer(in._grid);
@@ -118,16 +162,54 @@ public:
     M.Minv(in, X);  // X = G in
     M.MDeriv(X, MDer);  // MDer = U * dS/dU
     der = Implementation::projectForce(MDer);  // Ta if gauge fields
+    
   }
 
+  void AuxiliaryFieldsDerivative(MomentaField& der){
+    der = zero;
+    if (0){
+    // Auxiliary fields
+    MomentaField der_temp(der._grid);
+    MomentaField X(der._grid);
+    X=zero;
+    //M.M(AuxMom, X); // X = M Aux
+    // Two derivative terms
+    // the Mderiv need separation of left and right terms
+    M.MDeriv(AuxMom, der); 
+
+
+    // this one should not be necessary (identical to the previous one)
+    //M.MDeriv(X, AuxMom, der_temp); der += der_temp;
+
+    der = -1.0*Implementation::projectForce(der);
+    }
+  }
 
-  //   
   void DerivativeP(MomentaField& der){
     der = zero;
     M.Minv(Mom, der);
+    // is the projection necessary here?
+    // no for fields in the algebra
     der = Implementation::projectForce(der); 
   }
 
+  void update_auxiliary_momenta(RealD ep){
+    if(0){
+      AuxMom -= ep * AuxField;
+    }
+  }
+
+  void update_auxiliary_fields(RealD ep){
+    if (0) {
+      MomentaField tmp(AuxMom._grid);
+      MomentaField tmp2(AuxMom._grid);
+      M.M(AuxMom, tmp);
+      // M.M(tmp, tmp2);
+      AuxField += ep * tmp;  // M^2 AuxMom
+      // factor of 2?
+    }
+  }
+
 };