Merge pull request #70 from issp-center-dev/temperatures

Output temperatures in Kelvin (revert to the old behavior)
issp-center-dev · Oct 23, 2024 · 9ba4cc9 · 9ba4cc9
2 parents 1ae0bf9 + 08d7c1b
commit 9ba4cc9
Show file tree

Hide file tree

Showing 7 changed files with 65 additions and 37 deletions.
diff --git a/abics/sampling/mc_mpi.py b/abics/sampling/mc_mpi.py
@@ -210,6 +210,10 @@ class ParallelMC(object):
     "MPI size"
     kTs: list[float]
     "Temperatures"
+    T2E: float
+    "Temperature to energy conversion factor"
+    E2T: float
+    "Energy to temperature conversion factor"
     nreplicas: int
     "The number of replicas"
     model: Model
@@ -230,6 +234,7 @@ def __init__(
         configs,
         kTs,
         write_node=True,
+        T2E=1.0,
     ):
         """
 
@@ -245,14 +250,18 @@ def __init__(
             Configurations
         kTs: list
             Temperature list
+        T2E: float
+            Temperature to energy conversion factor
         """
         self.comm = comm
         self.rank = self.comm.Get_rank()
         self.procs = self.comm.Get_size()
-        self.kTs = kTs
         self.model = model
         self.nreplicas = len(configs)
         self.write_node = write_node
+        self.T2E = T2E
+        self.E2T = 1.0 / T2E
+        self.kTs = [T2E * T for T in kTs]
 
         ## mycalc.kT and mycalc.config should be set later
         myconfig = configs[0]

diff --git a/abics/sampling/pamc.py b/abics/sampling/pamc.py
@@ -99,6 +99,7 @@ def from_dict(cls, d):
             kTend = d["kTend"]
             kTnum = d["kTnum"]
             params.kTs = np.linspace(kTstart, kTend, kTnum)
+
         if "nsteps_between_anneal" in d:
             params.nsteps = d["nsteps_between_anneal"] * kTnum
             if "nsteps" in d:
@@ -157,7 +158,8 @@ def __init__(
         model: Model,
         configs,
         kTs: np.ndarray,
-        write_node=True,
+        write_node: bool=True,
+        T2E: float=1.0,
     ):
         """
 
@@ -179,8 +181,8 @@ def __init__(
         if kTs[0] < kTs[-1]:
             logger.warning("Warning: kTs[0] < kTs[-1]. abICS reverses kTs.")
             kTs = kTs[::-1]
-        super().__init__(comm, MCalgo, model, configs, kTs, write_node=write_node)
-        self.mycalc.kT = kTs[0]
+        super().__init__(comm, MCalgo, model, configs, kTs, write_node=write_node, T2E=T2E)
+        self.mycalc.kT = self.kTs[0]
         self.mycalc.config = configs[self.rank]
         self.betas = 1.0 / np.array(kTs)
         self.obs_save = []
@@ -360,7 +362,8 @@ def run(
                 self.anneal(self.mycalc.energy)
                 logger.info(
                     "--Anneal from T={} to {}".format(
-                        self.kTs[self.Tindex - 1], self.kTs[self.Tindex]
+                        self.E2T*self.kTs[self.Tindex - 1],
+                        self.E2T*self.kTs[self.Tindex]
                     )
                 )
                 if self.Tindex % resample_frequency == 0:
@@ -389,7 +392,7 @@ def run(
             self.acceptance_ratios.append(naccepted / ntrials)
             with open("acceptance_ratio.dat", "w") as f:
                 for T, ar in zip(self.kTs, self.acceptance_ratios):
-                    f.write(f"{T} {ar}\n")
+                    f.write(f"{self.E2T*T} {ar}\n")
             self.Tindex += 1
             self.isamples_offsets[self.Tindex] = nsample
         output.close()
@@ -421,7 +424,7 @@ def postproc(self):
         isamples_offsets = self.isamples_offsets
         comm = self.comm
         postproc(
-            kTs, obs_save, dlogz, logweight_history, obsnames, isamples_offsets, comm
+            kTs, obs_save, dlogz, logweight_history, obsnames, isamples_offsets, comm, E2T=self.E2T
         )
 
 
@@ -433,6 +436,7 @@ def postproc(
     obsnames,
     isamples_offsets,
     comm,
+    E2T:float=1.0,
 ):
     procs = comm.Get_size()
     rank = comm.Get_rank()
@@ -505,7 +509,7 @@ def postproc(
                 f.write(f"# $7: ERROR of <{oname}^2> - <{oname}>^2\n")
 
                 for iT in range(nT):
-                    f.write(f"{kTs[iT]}")
+                    f.write(f"{E2T*kTs[iT]}")
                     for j in range(3):
                         f.write(f" {o_mean[iT, 3*iobs+j]} {o_err[iT, 3*iobs+j]}")
                     f.write("\n")
@@ -523,5 +527,5 @@ def postproc(
             for i, (dlz, dlz_e) in enumerate(zip(dlogZ, dlogZ_err)):
                 F += dlz
                 dF += dlz_e
-                f.write(f"{kTs[i]} {F} {dF} {dlz} {dlz_e}\n")
+                f.write(f"{E2T*kTs[i]} {F} {dF} {dlz} {dlz_e}\n")
     comm.Barrier()
diff --git a/abics/sampling/rxmc.py b/abics/sampling/rxmc.py
@@ -135,7 +135,8 @@ def __init__(
         model: Model,
         configs,
         kTs: list[float],
-        write_node=True,
+        write_node: bool=True,
+        T2E: float=1.0,
     ):
         """
 
@@ -154,10 +155,10 @@ def __init__(
         subdirs: boolean
             If true, working directory for this rank is made
         """
-        super().__init__(comm, MCalgo, model, configs, kTs, write_node=write_node)
-        self.mycalc.kT = kTs[self.rank]
+        super().__init__(comm, MCalgo, model, configs, kTs, write_node=write_node, T2E=T2E)
+        self.mycalc.kT = self.kTs[self.rank]
         self.mycalc.config = configs[self.rank]
-        self.betas = 1.0 / np.array(kTs)
+        self.betas = 1.0 / np.array(self.kTs)
         self.rank_to_T = np.arange(0, self.procs, 1, dtype=np.int64)
         self.float_buffer = np.array(0.0, dtype=np.float64)
         self.int_buffer = np.array(0, dtype=np.int64)
@@ -360,9 +361,9 @@ def run(
         with open("acceptance_ratio.dat", "w") as f:
             for T, acc, trial in zip(self.kTs, self.naccepted, self.ntrials):
                 if trial > 0:
-                    f.write(f"{T} {acc/trial}\n")
+                    f.write(f"{self.E2T*T} {acc/trial}\n")
                 else:
-                    f.write(f"{T} {acc}/{trial}\n")
+                    f.write(f"{self.E2T*T} {acc}/{trial}\n")
 
         if nsample != 0:
             obs_buffer = np.empty(obs.shape)
@@ -406,7 +407,7 @@ def postproc(self, throw_out: int | float):
         kTs = self.kTs
         comm = self.comm
         obsnames = self.obsnames
-        postproc(obs_save, Trank_hist, kT_hist, kTs, comm, obsnames, throw_out)
+        postproc(obs_save, Trank_hist, kT_hist, kTs, comm, obsnames, throw_out, E2T=self.E2T)
 
     def __merge_obs(self):
         if hasattr(self, "obs_save0"):
@@ -427,7 +428,9 @@ def jackknife(X: np.ndarray) -> np.ndarray:
 
 
 def postproc(obs_save, Trank_hist, kT_hist, kTs, comm,
-             obsnames, throw_out: int | float):
+             obsnames, throw_out: int | float,
+             E2T: float = 1.0,
+             ):
     assert throw_out >= 0
     rank = comm.Get_rank()
     nT = comm.Get_size()
@@ -513,7 +516,7 @@ def postproc(obs_save, Trank_hist, kT_hist, kTs, comm,
                 f.write(f"# $6: <{oname}^2> - <{oname}>^2\n")
                 f.write(f"# $7: ERROR of <{oname}^2> - <{oname}>^2\n")
                 for iT in range(nT):
-                    f.write(f"{kTs[iT]}")
+                    f.write(f"{E2T*kTs[iT]}")
                     for itype in range(ntype):
                         f.write(f" {obs_all[iT, itype, iobs]}")
                     f.write("\n")
@@ -527,17 +530,17 @@ def postproc(obs_save, Trank_hist, kT_hist, kTs, comm,
             F = 0.0
             dF = 0.0
             if kTs[0] <= kTs[-1]:
-                f.write(f"{kTs[-1]} {F} {dF} {0.0} {0.0}\n")
+                f.write(f"{E2T*kTs[-1]} {F} {dF} {0.0} {0.0}\n")
                 for iT in np.arange(1, nT)[::-1]:
                     dlz, dlz_e = dlogz[iT]
                     F += dlz
                     dF += dlz_e
-                    f.write(f"{kTs[iT-1]} {F} {dF} {dlz} {dlz_e}\n")
+                    f.write(f"{E2T*kTs[iT-1]} {F} {dF} {dlz} {dlz_e}\n")
             else:
-                f.write(f"{kTs[0]} {F} {dF} {0.0} {0.0}\n")
+                f.write(f"{E2T*kTs[0]} {F} {dF} {0.0} {0.0}\n")
                 for iT in np.arange(0, nT - 1):
                     dlz, dlz_e = dlogz[iT]
                     F += dlz
                     dF += dlz_e
-                    f.write(f"{kTs[iT+1]} {F} {dF} {dlz} {dlz_e}\n")
+                    f.write(f"{E2T*kTs[iT+1]} {F} {dF} {dlz} {dlz_e}\n")
     comm.Barrier()
diff --git a/abics/sampling/simple_parallel.py b/abics/sampling/simple_parallel.py
@@ -210,7 +210,8 @@ def __init__(
         model: Model,
         configs,
         kTs=None,
-        write_node=True,
+        write_node:bool=True,
+        T2E:float=1.0,
     ):
         """
 
@@ -234,7 +235,9 @@ def __init__(
             kTs = [0] * self.procs
         if isinstance(kTs, (int, float)):
             kTs = [kTs] * self.procs
-        self.kTs = kTs
+        self.T2E = T2E
+        self.E2T = 1.0 / T2E
+        self.kTs = [T2E * T for T in kTs]
         self.model = model
         self.nreplicas = len(configs)
         self.write_node = write_node
@@ -248,7 +251,7 @@ def __init__(
             sys.exit(1)
 
         myconfig = configs[self.rank]
-        mytemp = kTs[self.rank]
+        mytemp = self.kTs[self.rank]
         self.mycalc = MCalgo(model, mytemp, myconfig)
         self.obs_save: list[np.ndarray] = []
         self.kT_hist: list[float] = []
@@ -381,11 +384,11 @@ def run(
 
 
     def postproc(self, throw_out):
-        postproc(obs_save=np.array(self.obs_save), kTs=self.kTs, comm=self.comm, obsnames=self.obsnames, throw_out=throw_out)
+        postproc(obs_save=np.array(self.obs_save), kTs=self.kTs, comm=self.comm, obsnames=self.obsnames, throw_out=throw_out, E2T=self.E2T)
 
 class RandomSampling_MPI(ParallelMC):
     def __init__(
-        self, comm, MCalgo: type[MCAlgorithm], model: Model, configs, write_node=True
+        self, comm, MCalgo: type[MCAlgorithm], model: Model, configs, write_node=True, T2E:float=1.0,
     ):
         """
 
@@ -401,8 +404,8 @@ def __init__(
             Configuration
         """
 
-        super().__init__(comm, MCalgo, model, configs, kTs, write_node=write_node)
-        self.mycalc.kT = kTs[self.rank]
+        super().__init__(comm, MCalgo, model, configs, kTs, write_node=write_node, T2E=T2E)
+        self.mycalc.kT = self.kTs[self.rank]
         self.mycalc.config = configs[self.rank]
         self.betas = 1.0 / np.array(kTs)
         self.rank_to_T = np.arange(0, self.procs, 1, dtype=np.int64)
@@ -415,7 +418,9 @@ def __init__(
 
 
 def postproc(obs_save, kTs, comm,
-             obsnames, throw_out: int | float):
+             obsnames, throw_out: int | float,
+             E2T: float = 1.0,
+             ):
     assert throw_out >= 0
     rank = comm.Get_rank()
     nT = comm.Get_size()
@@ -456,7 +461,7 @@ def postproc(obs_save, kTs, comm,
                 f.write(f"# $6: <{oname}^2> - <{oname}>^2\n")
                 f.write(f"# $7: ERROR of <{oname}^2> - <{oname}>^2\n")
                 for iT in range(nT):
-                    f.write(f"{kTs[iT]}")
+                    f.write(f"{E2T*kTs[iT]}")
                     for itype in range(ntype):
                         f.write(f" {obs_all[iT, itype, iobs]}")
                     f.write("\n")

diff --git a/abics/scripts/main_dft_latgas.py b/abics/scripts/main_dft_latgas.py
@@ -89,7 +89,7 @@ def main_dft_latgas(params_root: MutableMapping):
 
         # RXMC parameters
         # specify temperatures for each replica, number of steps, etc.
-        kTs = kB * rxparams.kTs
+        kTs = rxparams.kTs
         kTstart = rxparams.kTs[0]
         kTend = rxparams.kTs[-1]
 
@@ -122,7 +122,7 @@ def main_dft_latgas(params_root: MutableMapping):
         kTend = pamcparams.kTs[-1]
         if kTstart < kTend:
             kTstart, kTend = kTend, kTstart
-        kTs = kB * pamcparams.kTs
+        kTs = pamcparams.kTs
 
         # Set Lreload to True when restarting
         Lreload = pamcparams.reload
@@ -168,7 +168,7 @@ def main_dft_latgas(params_root: MutableMapping):
         # specify temperatures for each replica, number of steps, etc.
         kTstart = rxparams.kTs[0]
         kTend = rxparams.kTs[-1]
-        kTs = kB * rxparams.kTs
+        kTs = rxparams.kTs
 
         # Set Lreload to True when restarting
         Lreload = rxparams.reload
@@ -350,7 +350,7 @@ def main_dft_latgas(params_root: MutableMapping):
     if sampler_type == "RXMC":
         # RXMC calculation
         RXcalc = TemperatureRX_MPI(
-            comm, mc_class, model, configs, kTs, write_node=write_node
+            comm, mc_class, model, configs, kTs, write_node=write_node, T2E=kB,
         )
         if Lreload:
             logger.info("-Reloading from previous calculation")
@@ -371,7 +371,7 @@ def main_dft_latgas(params_root: MutableMapping):
     elif sampler_type == "PAMC":
         # PAMC calculation
         PAcalc = PopulationAnnealing(
-            comm, mc_class, model, configs, kTs, write_node=write_node
+            comm, mc_class, model, configs, kTs, write_node=write_node, T2E=kB,
         )
         if Lreload:
             logger.info("-Reloading from previous calculation")
@@ -404,7 +404,7 @@ def main_dft_latgas(params_root: MutableMapping):
 
     elif sampler_type == "parallelMC":
         calc = EmbarrassinglyParallelSampling(
-            comm, mc_class, model, configs, kTs, write_node=write_node
+            comm, mc_class, model, configs, kTs, write_node=write_node, T2E=kB
         )
         if Lreload:
             calc.reload()

diff --git a/abics/scripts/postproc_dft_latgas.py b/abics/scripts/postproc_dft_latgas.py
@@ -387,6 +387,7 @@ def postproc_dft_latgas(params_root: MutableMapping):
             comm=comm,
             obsnames=obsnames,
             throw_out=throw_out,
+            E2T=1.0/kB,
         )
 
     elif sampler_type == "PAMC":
@@ -410,6 +411,7 @@ def postproc_dft_latgas(params_root: MutableMapping):
             obsnames=obsnames,
             isamples_offsets=isamples_offsets,
             comm=comm,
+            E2T=1.0/kB,
         )
 
     elif sampler_type == "parallelRand":
@@ -420,6 +422,7 @@ def postproc_dft_latgas(params_root: MutableMapping):
             comm=comm,
             obsnames=obsnames,
             throw_out=throw_out,
+            E2T=1.0/kB,
         )
 
     elif sampler_type == "parallelMC":
@@ -430,6 +433,7 @@ def postproc_dft_latgas(params_root: MutableMapping):
             comm=comm,
             obsnames=obsnames,
             throw_out=throw_out,
+            E2T=1.0/kB,
         )
     logger.info("--Sampling completed sucessfully.")
     logger.info("Exiting normally on {}\n".format(datetime.datetime.now()))
diff --git a/tests/integration/active_learn/input.toml b/tests/integration/active_learn/input.toml
@@ -1,4 +1,7 @@
 [sampling]
+sampler = "RXMC"
+# sampler = "PAMC"
+# kTnum = 10
 nreplicas = 2
 nprocs_per_replica = 1 
 kTstart = 1200.0