grand canonical sampling

abics/applications/latgas_abinitio_interface/default_observer.py

@@ -26,8 +26,8 @@
 from pymatgen.analysis.structure_matcher import StructureMatcher, FrameworkComparator
 
 from abics import __version__
+from abics.util import expand_path, expand_cmd_path
 from abics.exception import InputError
-from abics.util import expand_path
 from abics.mc import ObserverBase, MCAlgorithm
 from abics.applications.latgas_abinitio_interface.base_solver import create_solver
 from abics.applications.latgas_abinitio_interface.run_base_mpi import (
@@ -64,6 +64,8 @@ def similarity(
     nsites += len(sites)
     return matched / nsites
 
+import logging
+logger = logging.getLogger(__name__)
 
 class DefaultObserver(ObserverBase):
     """
@@ -90,15 +92,15 @@ def __init__(self, comm, Lreload=False, params={}):
             Reload or not
         """
         super().__init__(comm, Lreload, params)
-        self.minE = 100000.0
+        self.minE = None
         myrank = comm.Get_rank()
         if Lreload:
             with open(os.path.join(str(myrank), "minEfi.dat"), "r") as f:
                 self.minE = float(f.readlines()[-1])
             with open(os.path.join(str(myrank), "obs.dat"), "r") as f:
                 self.lprintcount = int(f.readlines()[-1].split()[0]) + 1
 
-        self.names = ["energy"]
+        self.names = ["energy_internal", "energy"]
 
         params_solvers = params.get("solver", [])
         self.calculators = []
@@ -159,13 +161,18 @@ def __init__(self, comm, Lreload=False, params={}):
     def logfunc(self, calc_state: MCAlgorithm) -> tuple[float, ...]:
         assert calc_state.config is not None
         structure: Structure = calc_state.config.structure_norel
-        energy = calc_state.energy
-        result = [energy]
-        if energy < self.minE:
+
+        energy_internal = calc_state.model.internal_energy(calc_state.config)
+        energy = calc_state.model.energy(calc_state.config)
+
+        result = [energy_internal, energy]
+
+        if self.minE is None or energy < self.minE:
             self.minE = energy
             with open("minEfi.dat", "a") as f:
                 f.write(str(self.minE) + "\n")
             structure.to(fmt="POSCAR", filename="minE.vasp")
+
         for calculator in self.calculators:
             name = calculator["name"]
             runners: list[Runner] = calculator["runners"]
@@ -212,12 +219,16 @@ def __init__(self, comm, energy_calculators, Lreload=False):
         self.comm = comm
 
     def logfunc(self, calc_state: MCAlgorithm):
-        if calc_state.energy < self.minE:
-            self.minE = calc_state.energy
+        energy_internal = calc_state.model.internal_energy(calc_state.config)
+        energy = calc_state.model.energy(calc_state.config)
+
+        if self.minE is None or energy_internal < self.minE:
+            self.minE = energy_internal
             with open("minEfi.dat", "a") as f:
                 f.write(str(self.minE) + "\n")
             calc_state.config.structure.to(fmt="POSCAR", filename="minE.vasp")
-        energies = [calc_state.energy]
+
+        energies = [energy_internal]
         npar = self.comm.Get_size()
         if npar > 1:
             assert npar == len(self.calculators)
@@ -240,6 +251,7 @@ def logfunc(self, calc_state: MCAlgorithm):
             std = np.std(energies_tmp, ddof=1)
         energies.extend(energies_tmp)
         energies.append(std)
+
         return np.asarray(energies)
 
 
@@ -268,7 +280,8 @@ def from_dict(cls, d):
             map(lambda x: expand_path(x, os.getcwd()), base_input_dirs)
         )
         params.solver = d["type"]
-        params.path = expand_path(d["path"], os.getcwd())
+        #params.path = expand_path(d["path"], os.getcwd())
+        params.path = expand_cmd_path(d["path"])
         params.perturb = d.get("perturb", 0.1)
         params.solver_run_scheme = d.get("run_scheme", "mpi_spawn_ready")
         params.ignore_species = d.get("ignore_species", None)

abics/applications/latgas_abinitio_interface/defect.py

@@ -18,6 +18,7 @@
 
 from typing import Any, MutableMapping
 import os, sys
+import numpy as np
 
 from pymatgen.core import Lattice, Structure
 
@@ -26,6 +27,9 @@
 from abics.exception import InputError
 from abics.util import read_matrix, expand_path
 
+import logging
+logger = logging.getLogger("main")
+
 
 class DFTConfigParams:
 
@@ -108,6 +112,50 @@ def __init__(self, dconfig: MutableMapping[str, Any]) -> None:
             for ds in dconfig["defect_structure"]
         ]
 
+        if "chemical_potential" in dconfig:
+            self.chemical_potential = []
+            for entry in dconfig["chemical_potential"]:
+                species = entry.get("species", None)
+                mu = entry.get("mu", None)
+                if type is None or mu is None:
+                    raise InputError("invalid input for chemical potentials")
+                species = [species] if type(species) is str else species
+                self.chemical_potential.append({'species': species, 'mu': mu})
+        else:
+            self.chemical_potential = None
+        logger.debug("chemical_potential = {}".format(self.chemical_potential))
+
+        def _uniq_count(xs):
+            return [(i,j) for i,j in zip(*np.unique(xs,return_counts=True))]
+
+        if "grandcanonical_move" in dconfig:
+            self.gc_move = []
+            for entry in dconfig["grandcanonical_move"]:
+                if "species" in entry.keys():
+                    species = entry["species"]
+                    species = [ species ] if type(species) is not list else species
+                    self.gc_move.append({ 'from': _uniq_count(species), 'to': [], 'reverse': True })
+                else:
+                    sp_from = entry.get("from", [])
+                    sp_to   = entry.get("to", [])
+                    #sp_rev  = entry.get("reverse", True)
+                    sp_rev  = True  # consider only reversible case
+                    if sp_from is [] and sp_to is []:
+                        raise("invalid input for grandcanonical_move")
+                    sp_from = [ sp_from ] if type(sp_from) is not list else sp_from
+                    sp_to   = [ sp_to ]   if type(sp_to)   is not list else sp_to
+                    self.gc_move.append({ 'from': _uniq_count(sp_from), 'to': _uniq_count(sp_to), 'reverse': sp_rev })
+        else:
+            if self.chemical_potential is None:
+                self.gc_move = None
+            else:
+                # taken from chemical_potential table: assume add/remove of species
+                self.gc_move = []
+                for entry in self.chemical_potential:
+                    species = entry["species"]
+                    self.gc_move.append({ 'from': _uniq_count(species), 'to': [], 'reverse': True })
+        logger.debug("grandcanonical_move = {}".format(self.gc_move))
+
     @classmethod
     def from_dict(cls, d: MutableMapping) -> "DFTConfigParams":
         return cls(d)
@@ -149,12 +197,14 @@ def defect_config(cparam: DFTConfigParams) -> Config:
         spinel configure object
     """
     spinel_config = Config(
-        cparam.base_structure,
-        cparam.defect_sublattices,
-        cparam.num_defects,
-        cparam.supercell,
-        cparam.constraint_func,
-        cparam.constraint_energy
+        base_structure = cparam.base_structure,
+        defect_sublattices = cparam.defect_sublattices,
+        num_defects = cparam.num_defects,
+        cellsize = cparam.supercell,
+        constraint_func = cparam.constraint_func,
+        constraint_energy = cparam.constraint_energy,
+        chemical_potential = cparam.chemical_potential,
+        gc_move = cparam.gc_move,
     )
     spinel_config.shuffle()
     return spinel_config