From f3f24b3017eea98336c1c88783d51c7fd9c767c3 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Thu, 8 Nov 2018 12:55:25 +0000
Subject: [PATCH] Optional Twisted BC's added, in "DoubleStore" for WilsonImpl.
 Untested but doesn't affect answers when twists are all zero. The zero is the
 default behaviour for ImplParams.

---
 Grid/qcd/action/ActionParams.h                |  7 +++++--
 Grid/qcd/action/fermion/FermionOperatorImpl.h | 20 ++++++++++++++++---
 2 files changed, 22 insertions(+), 5 deletions(-)

diff --git a/Grid/qcd/action/ActionParams.h b/Grid/qcd/action/ActionParams.h
index d25b60a9..76868abc 100644
--- a/Grid/qcd/action/ActionParams.h
+++ b/Grid/qcd/action/ActionParams.h
@@ -44,12 +44,15 @@ namespace QCD {
   
   struct WilsonImplParams {
     bool overlapCommsCompute;
+    std::vector<Real> twist_n_2pi_L;
     std::vector<Complex> boundary_phases;
     WilsonImplParams() : overlapCommsCompute(false) {
       boundary_phases.resize(Nd, 1.0);
+      twist_n_2pi_L.resize(Nd, 0.0);
     };
-    WilsonImplParams(const std::vector<Complex> phi)
-      : boundary_phases(phi), overlapCommsCompute(false) {}
+    WilsonImplParams(const std::vector<Complex> phi) : boundary_phases(phi), overlapCommsCompute(false) {
+      twist_n_2pi_L.resize(Nd, 0.0);
+    }
   };
 
   struct StaggeredImplParams {
diff --git a/Grid/qcd/action/fermion/FermionOperatorImpl.h b/Grid/qcd/action/fermion/FermionOperatorImpl.h
index b0ffa90e..721004e1 100644
--- a/Grid/qcd/action/fermion/FermionOperatorImpl.h
+++ b/Grid/qcd/action/fermion/FermionOperatorImpl.h
@@ -240,16 +240,30 @@ namespace QCD {
       GaugeLinkField tmp(GaugeGrid);
 
       Lattice<iScalar<vInteger> > coor(GaugeGrid);
+      ////////////////////////////////////////////////////
+      // apply any boundary phase or twists
+      ////////////////////////////////////////////////////
       for (int mu = 0; mu < Nd; mu++) {
 
-	      auto pha = Params.boundary_phases[mu];
-	      scalar_type phase( real(pha),imag(pha) );
+	////////// boundary phase /////////////
+	auto pha = Params.boundary_phases[mu];
+	scalar_type phase( real(pha),imag(pha) );
 
-        int Lmu = GaugeGrid->GlobalDimensions()[mu] - 1;
+	int L   = GaugeGrid->GlobalDimensions()[mu];
+        int Lmu = L - 1;
 
         LatticeCoordinate(coor, mu);
 
         U = PeekIndex<LorentzIndex>(Umu, mu);
+
+	// apply any twists
+	RealD theta = Params.twist_n_2pi_L[mu] * 2*M_PI / L;
+	if ( theta != 0.0) { 
+	  scalar_type twphase(::cos(theta),::sin(theta));
+	  U = twphase*U;
+	  std::cout << GridLogMessage << " Twist ["<<mu<<"] "<< Params.twist_n_2pi_L[mu]<< " phase"<<phase <<std::endl;
+	}
+
         tmp = where(coor == Lmu, phase * U, U);
         PokeIndex<LorentzIndex>(Uds, tmp, mu);