From 886fb44128ff923c9c2a27bed6e4b53ff82f95ee Mon Sep 17 00:00:00 2001
From: Luchang Jin <ljin.luchang@gmail.com>
Date: Sat, 18 May 2024 02:23:42 -0400
Subject: [PATCH] add topo-measure-gpt by Christoph

---
 qlat/bin/topo-measure-gpt             |  3 +
 qlat/qlat/__main__.py                 |  3 +
 qlat/qlat/scripts/meson.build         |  1 +
 qlat/qlat/scripts/topo_measure.py     |  4 ++
 qlat/qlat/scripts/topo_measure_gpt.py | 85 +++++++++++++++++++++++++++
 5 files changed, 96 insertions(+)
 create mode 100755 qlat/bin/topo-measure-gpt
 create mode 100644 qlat/qlat/scripts/topo_measure_gpt.py

diff --git a/qlat/bin/topo-measure-gpt b/qlat/bin/topo-measure-gpt
new file mode 100755
index 000000000..3b09145b0
--- /dev/null
+++ b/qlat/bin/topo-measure-gpt
@@ -0,0 +1,3 @@
+#!/usr/bin/env python3
+
+import qlat.scripts.topo_measure_gpt
diff --git a/qlat/qlat/__main__.py b/qlat/qlat/__main__.py
index 1fa53d1ac..db73fc57a 100644
--- a/qlat/qlat/__main__.py
+++ b/qlat/qlat/__main__.py
@@ -12,6 +12,7 @@
     python3 -m qlat fields-properly-truncate ...
     python3 -m qlat gauge-fix-coulomb ...
     python3 -m qlat topo-measure ...
+    python3 -m qlat topo-measure-gpt ...
 """.strip()
 
 if len(sys.argv) < 2:
@@ -41,6 +42,8 @@
     from .scripts import gauge_fix_coulomb
 elif action == "topo-measure":
     from .scripts import topo_measure
+elif action == "topo-measure-gpt":
+    from .scripts import topo_measure_gpt
 else:
     print(usage_message)
     exit()
diff --git a/qlat/qlat/scripts/meson.build b/qlat/qlat/scripts/meson.build
index 289b14ea4..900b0cfc6 100644
--- a/qlat/qlat/scripts/meson.build
+++ b/qlat/qlat/scripts/meson.build
@@ -9,6 +9,7 @@ py_fns = [
   'gauge_fix_coulomb.py',
   'qlat_config.py',
   'topo_measure.py',
+  'topo_measure_gpt.py',
   ]
 py_files = files(py_fns)
 
diff --git a/qlat/qlat/scripts/topo_measure.py b/qlat/qlat/scripts/topo_measure.py
index c69aa3929..de0fc3347 100644
--- a/qlat/qlat/scripts/topo_measure.py
+++ b/qlat/qlat/scripts/topo_measure.py
@@ -3,6 +3,10 @@
 import numpy as np
 from pprint import pformat
 
+q.displayln_info("Topological charge measurement with Qlattice")
+q.displayln_info("by Luchang Jin")
+q.displayln_info("2024/01/25")
+
 if len(sys.argv) == 1:
     q.displayln_info("Usage: topo-measure [ --source source_config ] [ --output output.pickle ] [ --show-topo-terms ] [ --density-field-path path_for_density_field ]")
 
diff --git a/qlat/qlat/scripts/topo_measure_gpt.py b/qlat/qlat/scripts/topo_measure_gpt.py
new file mode 100644
index 000000000..9bde4bb45
--- /dev/null
+++ b/qlat/qlat/scripts/topo_measure_gpt.py
@@ -0,0 +1,85 @@
+import gpt as g
+import numpy as np
+import sys, os
+
+g.message("Topological charge measurement with GPT")
+g.message("by Christoph Lehner")
+g.message("2024/05/17")
+
+config = g.default.get_single("--config",None)
+
+def smear(Uin, c1, epsilon):
+    Usm = g.copy(Uin)
+    if c1 == 0.0:
+        action = g.qcd.gauge.action.wilson(2.0 * Uin[0].otype.shape[0])
+    else:
+        action = g.qcd.gauge.action.improved_with_rectangle(2.0 * Uin[0].otype.shape[0], c1)
+    g.algorithms.integrator.euler(
+        Usm, lambda: [g(-u) for u in action.gradient(Usm, Usm)], 1.0
+    )(epsilon)
+    return Usm
+
+smear_info_list = [
+    [ 20, 0.05, 0.0 ],
+    [ 20, 0.05, 0.0 ],
+    [ 20, 0.05, 0.0 ],
+    [ 50, 0.01, -1.4008 ],
+    [ 50, 0.01, -1.4008 ],
+    [ 50, 0.01, -1.4008 ],
+    [ 50, 0.01, -1.4008 ]
+]
+
+U = g.load(config)
+
+vol3d = float(np.prod(U[0].grid.gdimensions[0:3]))
+
+w = g.corr_io.writer(f"{config}.gluonic_luchang")
+output = g.gpt_io.writer(f"{config}.gluonic_fields_luchang", mpi=[2,2,2,4])
+
+# first plaquette
+g.message("Plaquette")
+P = g.slice(g.qcd.gauge.rectangle(U, 1, 1, field=True) / vol3d, 3)
+w.write("P", P)
+
+fQ = open(f"{config}.Q_luchang","wt")
+fE = open(f"{config}.E_luchang","wt")
+
+tau = 0.0
+U_wf = U
+c = {}
+for nsmear, epsilon, c1 in smear_info_list:
+
+    for i in range(nsmear):
+        tau += epsilon
+        U_wf = smear(U_wf, c1, epsilon)
+        g.message("%g" % tau, g.qcd.gauge.plaquette(U_wf))
+
+    g.message("Field Strength")
+    E_field = g.qcd.gauge.energy_density(U_wf, field=True)
+    E = g.slice(E_field / vol3d, 3)
+    w.write("E(%g)" % tau, E)
+
+    E = sum(E).real / len(E)
+    t2E = tau**2 * E
+    g.message("t2E = ", t2E)
+
+    g.message("Topology")
+    Q_field = g.qcd.gauge.topological_charge_5LI(U_wf, cache=c, field=True)
+    Q = g.slice(Q_field / vol3d, 3)
+    w.write("Q(%g)" % tau, Q)
+
+    Q = sum(Q).real / len(Q)
+    fQ.write("%g %.15g\n" % (tau, Q))
+    fQ.flush()
+
+    output.write({
+        "t" : tau,
+        "Q" : Q_field,
+        "E" : E_field
+    })
+    output.flush()
+
+    fE.write("%g %.15g\n" % (tau, E))
+    fE.flush()
+
+exit()