diff --git a/Makefile b/Makefile index 65ceda1..0a8db3d 100644 --- a/Makefile +++ b/Makefile @@ -38,10 +38,9 @@ light: venv @$(INSIDE_VENV) pip install src/ nogui: VENV_PATH = ./_venv -nogui: venv - @$(INSIDE_VENV) pip install --no-cache-dir --upgrade -r src/pip.freeze +nogui: venv pybinding @$(INSIDE_VENV) pip install -e src/ - @+$(INSIDE_VENV) $(MAKE) -C src pybinding + _attrdict: # Check if virtualenv python version > 3.3.0 diff --git a/src/msspec/iodata.py b/src/msspec/iodata.py index c60ccc3..61b6772 100644 --- a/src/msspec/iodata.py +++ b/src/msspec/iodata.py @@ -17,7 +17,7 @@ # along with this msspec. If not, see . # # Source file : src/msspec/iodata.py -# Last modified: Tue, 22 Oct 2024 12:39:54 +0200 +# Last modified: Thu, 27 Feb 2025 16:33:09 +0100 # Committed by : Sylvain Tricot diff --git a/src/msspec/looper.py b/src/msspec/looper.py index 0dbfdc4..bf91aa6 100644 --- a/src/msspec/looper.py +++ b/src/msspec/looper.py @@ -17,8 +17,8 @@ # along with this msspec. If not, see . # # Source file : src/msspec/looper.py -# Last modified: Mon, 27 Sep 2021 17:49:48 +0200 -# Committed by : sylvain tricot +# Last modified: Thu, 27 Feb 2025 16:33:09 +0100 +# Committed by : Sylvain Tricot from collections import OrderedDict from functools import partial @@ -92,9 +92,8 @@ class Sweep: class SweepRange: - def __init__(self, *sweeps, passindex=False): + def __init__(self, *sweeps): self.sweeps = sweeps - self.passindex = passindex self.index = 0 # First check that sweeps that are linked to another on are all included @@ -158,17 +157,15 @@ class SweepRange: for s in [sweep,] + children: key, value = s[idx] row[key] = value - if self.passindex: - row['sweep_index'] = i + row['sweep_index'] = i return row else: raise StopIteration @property def columns(self): - cols = [sweep.key for sweep in self.sweeps] - if self.passindex: - cols.append('sweep_index') + cols = ['sweep_index'] + cols += [sweep.key for sweep in self.sweeps] return cols @property @@ -202,31 +199,27 @@ class Looper: logger.debug("Pipeline called with {}".format(x)) return self.pipeline(**x) - def run(self, *sweeps, ncpu=1, passindex=False): + def run(self, *sweeps, ncpu=1, **kwargs): logger.info("Loop starts...") # prepare the list of inputs - sr = SweepRange(*sweeps, passindex=passindex) + sr = SweepRange(*sweeps) items = sr.items data = [] + t0 = time.time() + if ncpu == 1: # serial processing... logger.info("serial processing...") - t0 = time.time() - for i, values in enumerate(items): + values.update(kwargs) result = self._wrapper(values) data.append(result) - - t1 = time.time() - dt = t1 - t0 - logger.info("Processed {:d} sets of inputs in {:.3f} seconds".format( - len(sr), dt)) - else: # Parallel processing... chunksize = 1 #int(nsets/ncpu) + [values.update(kwargs) for values in items] logger.info(("Parallel processing over {:d} cpu (chunksize={:d})..." "").format(ncpu, chunksize)) t0 = time.time() @@ -236,21 +229,23 @@ class Looper: pool.close() pool.join() - t1 = time.time() - dt = t1 - t0 - logger.info(("Processed {:d} sets of inputs in {:.3f} seconds" - "").format(len(sr), dt)) + t1 = time.time() + dt = t1 - t0 + logger.info(("Processed {:d} sets of inputs in {:.3f} seconds" + "").format(len(sr), dt)) # Create the DataFrame dfdata = [] - columns = sr.columns + ['output',] + columns = sr.columns + list(kwargs.keys()) + ['output',] for i in range(len(sr)): row = list(items[i].values()) row.append(data[i]) dfdata.append(row) + df = pd.DataFrame(dfdata, columns=columns) + df = df.drop(columns=['sweep_index']) self.data = df @@ -259,14 +254,14 @@ class Looper: # of corresponding dict of parameters {'keyA': [val0,...valn], # 'keyB': [val0,...valn], ...} - all_xy = [] - for irow, row in df.iterrows(): - all_xy.append(row.output[0]) - all_xy.append(row.output[1]) - parameters = df.to_dict() - parameters.pop('output') + # all_xy = [] + # for irow, row in df.iterrows(): + # all_xy.append(row.output[0]) + # all_xy.append(row.output[1]) + # parameters = df.to_dict() + # parameters.pop('output') - return all_xy, parameters + return self.data #all_xy, parameters @@ -276,17 +271,16 @@ class Looper: if __name__ == "__main__": import numpy as np import time + import logging + + + logging.basicConfig(level=logging.DEBUG) - logger.setLevel(logging.DEBUG) def bar(**kwargs): - return 0 - - def post_process(data): - x = data.x.unique() - y = data.y.unique() - + i = kwargs.get('sweep_index') + return np.linspace(0,i,10) theta = Sweep(key='theta', comments="The polar angle", start=-70, stop=70, num=3, @@ -314,7 +308,16 @@ if __name__ == "__main__": looper = Looper() looper.pipeline = bar - data = looper.run(emitter, emitter_plane, uij, theta, levels, ncpu=4, - passindex=True) + other_kws = {'un':1, 'deux':2} + data = looper.run(emitter, emitter_plane, uij, theta, levels, ncpu=4, **other_kws) + + # Print the dataframe print(data) - #print(data[data.emitter_plane.eq(0)].theta.unique()) + + # Accessing the parameters and ouput values for a given sweep (e.g the last one) + print(looper.data.iloc[-1]) + + # Post-process the output values. For example here, the output is a 1D-array, + # make the sum of sweeps with 'Sr' emitter + X = np.array([ x for x in data[data.emitter == 'Sr'].output]).sum(axis=0) + print(X) diff --git a/src/msspec/utils.py b/src/msspec/utils.py index b5f4843..688d200 100644 --- a/src/msspec/utils.py +++ b/src/msspec/utils.py @@ -19,8 +19,8 @@ # along with this msspec. If not, see . # # Source file : src/msspec/utils.py -# Last modified: Thu, 06 Oct 2022 18:27:24 +0200 -# Committed by : Sylvain Tricot 1665073644 +0200 +# Last modified: Thu, 27 Feb 2025 16:33:09 +0100 +# Committed by : Sylvain Tricot """ @@ -468,8 +468,12 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0, # the symbol of your emitter symbol = cluster[np.where(cluster.get_tags() == emitter_tag)[0][0]].symbol - assert (diameter != 0 or planes != 0), \ - "At least one of diameter or planes parameter must be use." + if shape.lower() in ('spherical'): + assert (diameter != 0 or planes != 0), \ + "At least one of diameter or planes parameter must be use." + elif shape.lower() in ('cylindrical'): + assert (diameter != 0 and planes != 0), \ + "Diameter and planes parameters must be defined for cylindrical shape." if diameter == 0: # calculate the minimal diameter according to the number of planes @@ -479,6 +483,7 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0, # number of repetition in each direction rep = int(3*min_diameter/min(a, c)) + #print("rep = ", rep) # repeat the cluster cluster = cluster.repeat((rep, rep, rep)) @@ -542,7 +547,7 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0, xplan, yplan = get_xypos(cluster, zplan) radius = np.sqrt(xplan**2 + yplan**2 + zplan**2) - if diameter != 0: + if diameter != 0 and shape in ('spherical'): assert (radius <= diameter/2), ("The number of planes is too high " "compared to the diameter.") radius = max(radius, diameter/2) @@ -575,3 +580,90 @@ def hemispherical_cluster(cluster, emitter_tag=0, emitter_plane=0, diameter=0, Atoms.translate(cluster, [0, 0, -ze]) # put the emitter in (0,0,0) return cluster + +def shape_cluster(primitive, emitter_tag=0, emitter_plane=0, diameter=0, + planes=0, shape='spherical'): + + """Creates and returns a cluster based on an Atoms object and some + parameters. + + :param cluster: the Atoms object used to create the cluster + :type cluster: Atoms object + :param emitter_tag: the tag of your emitter + :type emitter_tag: integer + :param diameter: the diameter of your cluster in Angströms + :type diameter: float + :param planes: the number of planes of your cluster + :type planes: integer + :param emitter_plane: the plane where your emitter will be starting by 0 + for the first plane + :type emitter_plane: integer + + See :ref:`hemispherical_cluster_faq` for more informations. + """ + # We need the radius of the cluster and the number of planes + if shape.lower() in ('ispherical', 'cylindrical'): + assert (nplanes != 0 and diameter != 0), "nplanes and diameter cannot be zero for '{}' shape".format(shape) + elif shape.lower() in ('spherical'): + if diameter <= 0: + # find the diameter based on the number of planes + assert planes != 0, "planes should be > 0" + + + n = 3 + natoms = 0 + while True: + n += 2 + cluster = primitive.copy() + # Repeat the primitive cell + cluster = cluster.repeat((n, n, n)) + center_cluster(cluster) + + # Find the emitter closest to the origin + all_tags = cluster.get_tags() + are_emitters = all_tags == emitter_tag + _ie = np.linalg.norm(cluster[are_emitters].positions, axis=1).argmin() + ie = np.nonzero(are_emitters)[0][_ie] + # Translate the cluster to this emitter position + cluster.translate(-cluster[ie].position) + # cut plane at surface and at bottom + all_z = np.unique(cluster.positions[:,2]) + try: + zsurf = all_z[all_z >= 0][emitter_plane] + except IndexError: + # There are not enough planes above the emitter + zsurf = all_z.max() + try: + zbottom = all_z[all_z <= 0][::-1][planes - (emitter_plane+1)] + except IndexError: + # There are not enough planes below the emitter + zbottom = all_z.min() + cluster = cut_plane(cluster, z=[zbottom,zsurf]) + # spherical shape + if shape.lower() in ('spherical'): + cluster = cut_sphere(cluster, radius=diameter/2, center=(0,0,zsurf)) + if shape.lower() in ('ispherical'): + cluster = cut_sphere(cluster, radius=diameter/2, center=(0,0,0)) + elif shape.lower() in ('cylindrical'): + cluster = cut_cylinder(cluster, radius=diameter/2) + else: + raise NameError("Unknown shape") + cluster.set_cell(primitive.cell) + if len(cluster) <= natoms: + break + else: + natoms = len(cluster) + + + return cluster + + + + + + + + + + +