Merge branch 'master' of https://github.com/lucianogsilvestri/sarkas

Merging QMC changes with magnetic simchanges.

sarkas/particles.py

@@ -24,6 +24,7 @@
 from os.path import join
 from scipy.linalg import norm
 from scipy.spatial.distance import pdist
+from scipy.stats import moment, qmc
 from warnings import warn
 
 from .utilities.exceptions import ParticlesError, ParticlesWarning
@@ -158,6 +159,7 @@ def __init__(self):
         self.species_thermostat_temperatures = None
         self.species_masses = None
         self.species_charges = None
+        self.species_velocity_moments = None
 
         self.no_grs = None
         self.rdf_hist = None
@@ -170,6 +172,7 @@ def __init__(self):
             "Kinetic Energy": self.calculate_species_kinetic_temperature,
             "Potential Energy": self.calculate_species_potential_energy,
             "Momentum": self.calculate_species_momentum,
+            "Velocity Moments": self.calculate_species_velocity_moments,
             "Electric Current": self.calculate_species_electric_current,
             "Pressure Tensor": self.calculate_species_pressure_tensor,
             "Enthalpy": self.calculate_species_enthalpy,
@@ -185,6 +188,9 @@ def __init__(self):
             "Enthalpy": self.calculate_species_enthalpy,
             "Heat Flux": self.calculate_species_heat_flux,
         }
+        self.qmc_sequence = None
+        self.available_qmc_sequences = ["halton", "sobol", "poissondisk", "latinhypercube"]
+        self.max_velocity_distribution_moment = 6
 
     def __repr__(self):
         sortedDict = dict(sorted(self.__dict__.items(), key=lambda x: x[0].lower()))
@@ -297,6 +303,14 @@ def copy_params(self, params):
         self.box_volume = params.box_volume
         self.pbox_volume = params.pbox_volume
         self.load_method = params.load_method
+
+        if hasattr(params, "qmc_sequence"):
+            self.qmc_sequence = params.qmc_sequence
+
+        if hasattr(params, "max_velocity_distribution_moment"):
+            self.max_velocity_distribution_moment = params.max_velocity_distribution_moment
+            self.species_velocity_moments = zeros((self.num_species, self.max_velocity_distribution_moment))
+
         self.restart_step = params.restart_step
         self.particles_input_file = params.particles_input_file
         self.load_perturb = params.load_perturb
@@ -590,9 +604,82 @@ def initialize_positions(self):
                     self.box_lengths[0] / 2.0, self.load_gauss_sigma[sp], (sp_num, 3)
                 )
                 sp_start += sp_num
+        elif self.load_method in ["qmc", "quasi_monte_carlo"]:
+            # Ensure that we are not making silly typos
+            self.qmc_sequence = self.qmc_sequence.lower()
+
+            if self.qmc_sequence not in self.available_qmc_sequences:
+                raise AttributeError(
+                    f"Quasi Monte Carlo sequence not recognized. Please choose from {self.available_qmc_sequences}"
+                )
+
+            self.pos = self.quasi_monte_carlo(self.qmc_sequence, self.total_num_ptcls, self.dimensions, self.box_lengths)
         else:
             raise AttributeError("Incorrect particle placement scheme specified.")
 
+    @staticmethod
+    def quasi_monte_carlo(
+        sequence: str = "sobol", num_ptcls: int = 2**10, dimensions: int = 3, box_lengths: ndarray = None
+    ):
+        """
+        Place particles according to a quasi-Monte Carlo sequence.
+
+        This method uses the classes in the Quasi Monte Carlo module of `scipy` to place particles in a simulation box
+        according to a quasi-Monte Carlo sequence. The sequence can be chosen from the following options: Halton, Sobol,
+        PoissonDisk, LatinHypercube.
+
+        Parameters
+        ----------
+        sequence : str, optional
+            Name of the sequence to use. Options: Halton, Sobol, PoissonDisk, LatinHypercube. Default is "Sobol".
+
+        num_ptcls : int, optional
+            Number of particles to place. Default is 2**10.
+
+        dimensions : int, optional
+            Number of dimensions. Default is 3.
+
+        box_lengths : array-like, optional
+            Lengths of the simulation box in each dimension. Default is None.
+
+        Returns
+        -------
+        array-like
+            Randomly generated positions of the particles, scaled by the box lengths.
+
+        Raises
+        ------
+        ValueError
+            If an invalid sequence name is provided.
+
+        Notes
+        -----
+        The `quasi_monte_carlo` method initializes a quasi-Monte Carlo sequence generator based on the provided sequence
+        name. It then generates random positions for the specified number of particles in the specified number of
+        dimensions. The generated positions are scaled by the box lengths of the simulation.
+
+        Examples
+        --------
+        >>> positions = quasi_monte_carlo("Halton", 100, 3, [10, 10, 10])
+
+        """
+        if sequence == "sobol":
+            qmc_sampler = qmc.Sobol(dimensions)
+        elif sequence == "latinhypercube":
+            qmc_sampler = qmc.LatinHypercube(dimensions)
+        elif sequence == "poissondisk":
+            qmc_sampler = qmc.PoissonDisk(dimensions)
+        elif sequence == "halton":
+            qmc_sampler = qmc.Halton(dimensions)
+        else:
+            raise ValueError("Invalid sequence name.")
+
+        positions = qmc_sampler.random(num_ptcls)
+        if box_lengths is not None:
+            positions *= box_lengths
+
+        return positions
+
     def initialize_velocities(self, species):
         """
         Initialize particles' velocities based on the species input values. The velocities can be initialized from a
@@ -1036,6 +1123,19 @@ def calculate_species_momentum(self):
         velocity = vector_species_loop(self.vel.transpose(), self.species_num)
         self.species_momentum = self.species_masses * velocity
 
+    def calculate_species_velocity_moments(self):
+        """Calculate the moments of the velocity distribution using the velocity of each species and stores them into :attr:`species_velocity_moments`."""
+        species_start = 0
+        species_end = 0
+
+        for i, num in enumerate(self.species_num):
+            species_end += num
+            for mom in range(self.max_velocity_distribution_moment):
+                self.species_velocity_moments[i, mom] = moment(
+                    self.vel[species_start:species_end, :], moment=mom + 1, axis=0
+                )
+            species_start += num
+
     def calculate_species_potential_energy(self):
         """Calculate the potential energy of each species from :attr:`potential_energy`, calculated in the force loop, and stores it into :attr:`species_potential_energy`."""
         self.species_potential_energy = scalar_species_loop(self.potential_energy, self.species_num)
@@ -1351,7 +1451,7 @@ def setup(self, params, species):
             self.initialize_velocities(species)
             self.initialize_accelerations()
 
-        if len(self.observables_list) > 2:
+        if len(self.observables_list) > 3:
             self.calculate_thermodynamic_quantities = self.calculate_thermodynamic_quantities_full
             self.make_thermodynamics_dictionary = self.make_thermodynamics_dictionary_full
         else:

sarkas/processes.py

@@ -139,7 +139,7 @@ def update_subclasses_from_dict(self, nested_dict: dict):
                 # args = {"Particles" : [ { "Species" : { "name": "O" } } ] }
 
                 # Check if you already have a non-empty dict of species
-                if len(self.species) == 0:
+                if len(self.species) > 0:
                     # If so do you want to replace or update?
                     # Update species attributes
 
@@ -149,11 +149,11 @@ def update_subclasses_from_dict(self, nested_dict: dict):
                             self.species[sp].__dict__.update(spec.__dict__)
                         else:
                             self.species.append(spec)
-                    else:
-                        # Append new species
-                        for sp, species in enumerate(vals):
-                            spec = Species(species["Species"])
-                            self.species.append(spec)
+                else:
+                    # Append new species
+                    for sp, species in enumerate(vals):
+                        spec = Species(species["Species"])
+                        self.species.append(spec)
 
             elif lkey in ["Observables"]:
                 for obs_dict in vals:
@@ -491,7 +491,6 @@ def setup(self, read_yaml: bool = True, input_file: str = None, other_inputs: di
             self.update_subclasses_from_dict(other_inputs)
 
         if self.__name__ == "postprocessing":
-
             # Create the file paths without creating directories and redefining io attributes
             self.io.make_directory_tree()
             self.io.make_directories()
@@ -547,7 +546,6 @@ def setup_from_dict(self, input_dict: dict):
         self.instantiate_subclasses_from_dict(input_dict)
 
         if self.__name__ == "postprocessing":
-
             # Create the file paths without creating directories and redefining io attributes
             self.io.make_directory_tree()
             self.io.make_directories()
@@ -601,7 +599,6 @@ class PostProcess(Process):
     """
 
     def __init__(self, input_file: str = None, grab_last_step: bool = False):
-
         self.__name__ = "postprocessing"
         self.grab_last_step = grab_last_step
 
@@ -614,7 +611,6 @@ def run(self):
             print("No observables found in observables_dict")
         else:
             for obs_key, obs_class in self.observables_dict.items():
-
                 obs_class.setup(self.parameters)
                 msg = obs_class.pretty_print_msg()
                 self.io.write_to_logger(msg)
@@ -629,7 +625,6 @@ def run(self):
             print("No transport coefficients found in tranport_dict")
         else:
             for obs_key, obs_class in self.transport_dict.items():
-
                 obs_class.setup(self.parameters)
                 msg = obs_class.pretty_print_msg()
                 self.io.write_to_logger(msg)
@@ -1124,7 +1119,6 @@ def make_force_v_timing_plot(self, data_df: DataFrame = None):
             fig.savefig(join(fig_path, f"ForceErrorMap_v_Timing_cao_{cao}_{self.io.job_id}.png"))
 
     def postproc_estimates(self):
-
         # POST- PROCESSING
         self.io.postprocess_info(self, observable="header")
         # Header of process
@@ -1406,7 +1400,6 @@ def timing_study_calculation(self):
         for _, m in enumerate(
             tqdm(self.pm_meshes, desc="Looping over the PM meshes", disable=not self.parameters.verbose)
         ):
-
             # Setup PM params
             self.potential.pppm_mesh = m * array([1, 1, 1], dtype=int)
             self.potential.pppm_alpha_ewald = 0.3 * m / self.potential.box_lengths.min()
@@ -1708,7 +1701,6 @@ def magnetize(self):
         self.io.time_stamp("Magnetization", self.timer.time_division(time_eq))
 
     def produce(self):
-
         it_start = self.check_restart(phase="production")
 
         self.integrator.update = self.integrator.type_setup(self.integrator.production_type)