Merge pull request #45 from ICAMS/update/comp_integration

Update/comp integration

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.2.1
+current_version = 1.2.2
 commit = True
 tag = True
 

calphy/__init__.py

@@ -4,7 +4,7 @@
 from calphy.alchemy import Alchemy
 from calphy.routines import MeltingTemp
 
-__version__ = "1.2.1"
+__version__ = "1.2.2"
 
 def addtest(a,b):
     return a+b

calphy/alchemy.py

@@ -155,6 +155,7 @@ def run_integration(self, iteration=1):
         #remap the box to get the correct pressure
         lmp = ph.remap_box(lmp, self.lx, self.ly, self.lz)
 
+        lmp.command("velocity          all create %f %d mom yes rot yes dist gaussian"%(self.calc._temperature, np.random.randint(0, 10000)))
         # Integrator & thermostat.
         if self.calc._npt:
             lmp.command("fix             f1 all npt temp %f %f %f %s %f %f %f"%(self.calc._temperature, self.calc._temperature, 
@@ -321,8 +322,8 @@ def thermodynamic_integration(self):
 
 
     def mass_integration(self, flambda, ref_mass, target_masses, target_counts):
-        mcorarr = integrate_mass(flambda, ref_mass, target_masses, target_counts,
+        mcorarr, mcorsum = integrate_mass(flambda, ref_mass, target_masses, target_counts,
     self.calc._temperature, self.natoms)
-        return mcorarr      
+        return mcorarr, mcorsum      
 
 

calphy/input.py

@@ -30,6 +30,18 @@
 import shutil
 import numpy as np
 
+def read_report(folder):
+    """
+    Read the finished calculation report
+    """
+    repfile = os.path.join(folder, "report.yaml")
+    if not os.path.exists(repfile):
+        raise FileNotFoundError(f"file {repfile} not found")
+
+    with open(repfile, 'r') as fin:
+        data = yaml.safe_load(fin)
+    return data
+
 class InputTemplate:
     def __init__(self):
         pass
@@ -93,7 +105,6 @@ def merge_dicts(self, dicts):
                 merged_dict[key] = val
         return merged_dict
 
-
 class CompositionScaling(InputTemplate):
     def __init__(self):
         self._input_chemical_composition = None
@@ -157,6 +168,7 @@ def __init__(self):
         self._n_print_steps = 0
         self._n_iterations = 1
         self._equilibration_control = None
+        self._folder_prefix = None
 
         #add second level options; for example spring constants
         self._spring_constants = None
@@ -549,6 +561,14 @@ def spring_constants(self, val):
         val = self.check_and_convert_to_list(val, check_none=True)
         self._spring_constants = val
 
+    @property
+    def folder_prefix(self):
+        return self._folder_prefix
+
+    @folder_prefix.setter
+    def folder_prefix(self, val):
+        self._folder_prefix = val
+
     def fix_paths(self, potlist): 
         """
         Fix paths for potential files to complete ones
@@ -594,7 +614,11 @@ def create_identifier(self):
             l = self.lattice
             l = l.split('/')
             l = l[-1]
-        identistring = "-".join([prefix, l, str(ts), str(ps)])
+
+        if self.folder_prefix is None:
+            identistring = "-".join([prefix, l, str(ts), str(ps)])
+        else:
+            identistring = "-".join([self.folder_prefix, prefix, l, str(ts), str(ps)])
         return identistring
 
     def create_folders(self, prefix=None):
@@ -688,7 +712,8 @@ def read_inputfile(file):
         if mode == "melting_temperature":
             calc = Calculation.generate(indata)
             calc.add_from_dict(ci, keys=["mode", "pair_style", "pair_coeff", "repeat", "n_equilibration_steps",
-                                "n_switching_steps", "n_print_steps", "n_iterations", "spring_constants", "equilibration_control"])
+                                "n_switching_steps", "n_print_steps", "n_iterations", "spring_constants", "equilibration_control",
+                                "folder_prefix"])
             calc.pressure = Calculation.convert_to_list(ci["pressure"], check_none=True) if "pressure" in ci.keys() else 0
             calc.temperature = Calculation.convert_to_list(ci["temperature"]) if "temperature" in ci.keys() else None
             calc.lattice = Calculation.convert_to_list(ci["lattice"]) if "lattice" in ci.keys() else None
@@ -724,7 +749,7 @@ def read_inputfile(file):
                 calc = Calculation.generate(indata)
                 calc.add_from_dict(ci, keys=["mode", "pair_style", "pair_coeff", "repeat", "n_equilibration_steps",
                                 "n_switching_steps", "n_print_steps", "n_iterations", "potential_file", "spring_constants",
-                                "melting_cycle", "equilibration_control"])
+                                "melting_cycle", "equilibration_control", "folder_prefix"])
                 calc.lattice = combo[0]["lattice"]
                 calc.lattice_constant = combo[0]["lattice_constant"]
                 calc.reference_phase = combo[0]["reference_phase"]

calphy/integrators.py

@@ -360,15 +360,15 @@ def find_w(mainfolder, nelements=1, concentration=[1,], nsims=5, full=False, use
     if composition_integration:
         wsmean = np.mean(ws, axis=0)
         qsmean = np.mean(qs, axis=0)
-        qsstd = np.mean(qs, axis=0)
-        return wsmean, qsmean, qsstd, flambda
+        wsstd = np.std(ws, axis=0)
+        return wsmean, qsmean, wsstd, flambda
 
     wsmean = np.mean(ws)
     qsmean = np.mean(qs)
-    qsstd = np.mean(qs)
+    wsstd = np.std(ws)
 
     if full:
-        return wsmean, qsmean, qsstd
+        return wsmean, qsmean, wsstd
     else:
         return wsmean
 
@@ -671,7 +671,10 @@ def integrate_mass(flambda, ref_mass, target_masses, target_counts,
     temperature, natoms):
 
     mcorarr = np.zeros(len(flambda))
+    mcorsum = 0
+
     for i in range(len(target_masses)):
         mcorarr += 1.5*kb*temperature*(flambda*(target_counts[i]/natoms)*np.log(ref_mass/target_masses[i]))
-
-    return mcorarr
+        mcorsum += 1.5*kb*temperature*(1.0*(target_counts[i]/natoms)*np.log(ref_mass/target_masses[i]))
+
+    return mcorarr, mcorsum

calphy/phase.py

@@ -178,6 +178,15 @@ def __repr__(self):
         data = self.calc.__repr__()
         return data
 
+    def _from_dict(self, org_dict, indict):
+        for key, val in indict.items():
+            if isinstance(val, dict):
+                if key not in org_dict.keys():
+                    org_dict[key] = {}
+                self._from_dict(org_dict[key], val)
+            else:
+                org_dict[key] = val        
+
     def prepare_lattice(self):
         """
         Prepare the lattice for the simulation
@@ -564,13 +573,14 @@ def process_traj(self, filename, outfilename):
             os.remove(file)
 
 
-    def submit_report(self):
+    def submit_report(self, extra_dict=None):
         """
         Submit final report containing results
 
         Parameters
         ----------
-        None
+        extra_dict: dict
+            extra information to be written out
 
         Returns
         -------
@@ -602,6 +612,10 @@ def submit_report(self):
         report["results"]["reference_system"] = float(self.fref)
         report["results"]["work"] = float(self.w)
         report["results"]["pv"] = float(self.pv)
+
+        if extra_dict is not None:
+            self._from_dict(report, extra_dict)
+
         self.report = report
 
         reportfile = os.path.join(self.simfolder, "report.yaml")

calphy/routines.py

@@ -514,9 +514,12 @@ def routine_composition_scaling(job):
     flambda_arr, w_arr, q_arr, qerr_arr = job.thermodynamic_integration()
 
     #read the file
-    mcorrarr = job.mass_integration(flambda_arr, ref_mass, target_masses, target_counts)
+    mcorrarr, mcorsum = job.mass_integration(flambda_arr, ref_mass, target_masses, target_counts)
     netfe = w_arr - mcorrarr
 
+    job.fe = job.fe - mcorsum
+    job.submit_report(extra_dict = {"results":{"mass_correction": float(mcorsum)}})
+
     outfile = os.path.join(job.simfolder, "composition_sweep.dat")
     np.savetxt(outfile, np.column_stack((flambda_arr, netfe, w_arr, mcorrarr)))
     return job

docs/source/documentation/inputfile.md

@@ -12,7 +12,7 @@ The inputfile is `yaml` formatted. In this section the possible keys in the inpu
 | [mode](#mode) | [lattice](#lattice) | [reference_phase](#reference_phase) | [temperature](#temperature) | [pressure](#pressure) |
 | [temperature_high](#temperature_high) | [lattice_constant](#lattice_constant) | [repeat](#repeat) | [n_iterations](#n_iterations) | [n_switching_steps](#n_switching_steps) | 
 | [n_equilibration_steps](#n_equilibration_steps) | [pair_style](#pair_style) | [pair_coeff](#pair_coeff) | [n_print_steps](#n_print_steps) | [potential_file](#potential_file) |  
-| [spring_constants](#spring_constants) | [equilibration_control](#equilibration_control) | [melting_cycle](#melting_cycle) | | |   
+| [spring_constants](#spring_constants) | [equilibration_control](#equilibration_control) | [melting_cycle](#melting_cycle) | [folder_prefix](#folder_prefix) | |   
 
 | `md` block  | | | | |
 | :-: | :-: | :-: | :-: | :-: |
@@ -32,6 +32,10 @@ The inputfile is `yaml` formatted. In this section the possible keys in the inpu
 | :-: | :-: | :-: | :-: | :-: |
 | [step](#step) | [attempts](#attempts) |
 
+| `composition_scaling` block | | | | |
+| :-: | :-: | :-: | :-: | :-: |
+| [input_chemical_composition](#input_chemical_composition) | [output_chemical_composition](#output_chemical_composition) |
+
 | `nose_hoover` block | | | | |
 | :-: | :-: | :-: | :-: | :-: |
 | [thermostat_damping](#nose_hoover_thermostat_damping) | [barostat_damping](#nose_hoover_barostat_damping) |
@@ -384,6 +388,20 @@ melting_cycle: False
 
 If True, a melting cycle is carried out to melt the given input structure. Only used if the `reference_phase` is `"liquid"`.
 
+
+---
+
+#### <a name="folder_prefix"></a>`folder_prefix`        
+
+_type_: string    
+_default_: None     
+_example_:
+```
+folder_prefix: set1
+```  
+
+Prefix string to be added to folder names for calculation. Folders for calculations in calphy are named as `mode-lattice-temperature-pressure`. Therefore, if more than one calculation is run with the same parameters, they will be overwritten. To prevent this, `folder_prefix` can be used. If `folder_prefix` is provided, the folders will be named as `folder_prefix-mode-lattice-temperature-pressure`.
+
 ---
 ---
 
@@ -731,6 +749,51 @@ The number of maximum attempts to try find the melting temperature in a automate
 ---
 ---
 
+## `composition_scaling` block
+
+This block contains keywords that are used only for the mode `composition_scaling`.
+
+```
+composition_scaling:
+  input_chemical_composition:
+     - Cu: 512
+     - Zr: 512
+  output_chemical_composition:
+     - Cu: 513
+     - Zr: 511
+```
+
+---
+
+#### <a name="input_chemical_composition"></a>`input_chemical_composition`
+
+_type_: list         
+_example_:
+```
+input_chemical_composition:
+   - Cu: 512
+   - Zr: 512
+```
+
+The input chemical composition in number of atoms. It should be identical to the input structure provided.
+
+---
+
+#### <a name="output_chemical_composition"></a>`output_chemical_composition`
+
+_type_: list         
+_example_:
+```
+output_chemical_composition:
+   - Cu: 513
+   - Zr: 511
+```
+
+The output chemical composition in number of atoms. The total number of atoms should be equal to the input provided.
+
+---
+---
+
 
 ## <a name="nose_hoover"></a>`nose_hoover` block
 

docs/source/examples.rst

@@ -12,3 +12,4 @@ Examples
     examples/example_07/analysis
     examples/example_08/analysis
     examples/example_09/analysis
+    examples/example_10/analysis

environment-nolammps.yml

@@ -15,3 +15,4 @@ dependencies:
   - pyscal
   - matplotlib
   - tqdm
+  - uncertainties

environment.yml

@@ -13,3 +13,4 @@ dependencies:
   - pyscal
   - matplotlib
   - tqdm
+  - uncertainties