Merge pull request #86 from ICAMS/switch_to_pydantic

Switch to pydantic

calphy/alchemy.py

@@ -47,11 +47,12 @@ class Alchemy(cph.Phase):
         base folder for running calculations
 
     """
-    def __init__(self, calculation=None, simfolder=None):
+    def __init__(self, calculation=None, simfolder=None, log_to_screen=False):
 
         #call base class
         super().__init__(calculation=calculation,
-        simfolder=simfolder)
+        simfolder=simfolder,
+        log_to_screen=log_to_screen)
 
 
     def run_averaging(self):
@@ -160,7 +161,7 @@ def run_integration(self, iteration=1):
 
         lmp.command("velocity          all create %f %d mom yes rot yes dist gaussian"%(self.calc._temperature, np.random.randint(1, 10000)))
         # Integrator & thermostat.
-        if self.calc._npt:
+        if self.calc.npt:
             lmp.command("fix             f1 all npt temp %f %f %f %s %f %f %f"%(self.calc._temperature, self.calc._temperature, 
                 self.calc.md.thermostat_damping[1], self.iso, self.calc._pressure, self.calc._pressure, self.calc.md.barostat_damping[1]))        
         else:
@@ -198,8 +199,8 @@ def run_integration(self, iteration=1):
 
 
         lmp.command("pair_style       hybrid/scaled v_flambda %s v_blambda %s"%(
-            self.calc.pair_style_with_options[0],
-            self.calc.pair_style_with_options[1]
+            self.calc._pair_style_with_options[0],
+            self.calc._pair_style_with_options[1]
             )
         )
         lmp.command("pair_coeff       %s"%pc1)
@@ -236,7 +237,7 @@ def run_integration(self, iteration=1):
         lmp.command("uncompute       c2")
 
 
-        lmp.command("pair_style      %s"%self.calc.pair_style_with_options[1])
+        lmp.command("pair_style      %s"%self.calc._pair_style_with_options[1])
         lmp.command("pair_coeff      %s"%self.calc.pair_coeff[1])
 
         # Thermo output.
@@ -252,8 +253,8 @@ def run_integration(self, iteration=1):
         lmp.command("variable         blambda equal ramp(${li},${lf})")
 
 
-        lmp.command("pair_style       hybrid/scaled v_flambda %s v_blambda %s"%(self.calc.pair_style_with_options[0], 
-            self.calc.pair_style_with_options[1]))
+        lmp.command("pair_style       hybrid/scaled v_flambda %s v_blambda %s"%(self.calc._pair_style_with_options[0], 
+            self.calc._pair_style_with_options[1]))
         lmp.command("pair_coeff       %s"%pc1)
         lmp.command("pair_coeff       %s"%pc2)
 
@@ -309,7 +310,8 @@ def thermodynamic_integration(self):
         the calculated free energy is the same as the work.
         """
         w, q, qerr = find_w(self.simfolder, nelements=self.calc.n_elements, 
-            concentration=self.concentration, nsims=self.calc.n_iterations, 
+            concentration=[val['composition'] for key, val in self.calc._element_dict.items()], 
+            nsims=self.calc.n_iterations, 
             full=True, solid=False, alchemy=True)
 
         self.w = w
@@ -318,7 +320,8 @@ def thermodynamic_integration(self):
 
         if self.calc.mode == "composition_scaling":
             w_arr, q_arr, qerr_arr, flambda_arr = find_w(self.simfolder, nelements=self.calc.n_elements, 
-                concentration=self.concentration, nsims=self.calc.n_iterations, 
+                concentration=[val['composition'] for key, val in self.calc._element_dict.items()], 
+                nsims=self.calc.n_iterations, 
                 full=True, solid=False, alchemy=True, composition_integration=True)
 
             #now we need to process the comp scaling

calphy/clitools.py

@@ -7,7 +7,7 @@
 import time
 import datetime
 
-from calphy.input import read_inputfile, load_job, save_job
+from calphy.input import read_inputfile, load_job, save_job, _convert_legacy_inputfile
 from calphy.liquid import Liquid
 from calphy.solid import Solid
 from calphy.alchemy import Alchemy
@@ -90,3 +90,33 @@ def process_integration():
     job.thermodynamic_integration()
     job.submit_report()
     save_job(job)
+
+def convert_legacy_inputfile():
+    arg = ap.ArgumentParser()
+    arg.add_argument("-i", "--input", required=True, type=str,
+    help="name of the input file")
+    arg.add_argument("-s", "--split", required=False, type=bool, 
+    help="split each calculation into new file.", default=True)
+    arg.add_argument("-o", "--output", required=False, type=str, 
+    help="output file string, calculations will be named <outputstring>.*.yaml", 
+        default='input')    
+    args = vars(arg.parse_args())
+    calculations = _convert_legacy_inputfile(args['input'], return_calcs=True)
+    outputstr = args['output']
+
+    if args['split']:
+        #now we have to write this out to file
+        for count, calc in enumerate(calculations):
+            data = {}
+            data['calculations'] = [calc]
+            outfile = ".".join([outputstr, str(count+1), 'yaml'])
+            with open(outfile, 'w') as fout:
+                yaml.safe_dump(data, fout)
+    else:
+        data = {}        
+        data['calculations'] = calculations
+        outfile = ".".join([outputstr, 'yaml'])
+        with open(outfile, 'w') as fout:
+            yaml.safe_dump(data, fout)
+
+

calphy/composition_transformation.py

@@ -25,10 +25,11 @@
 import numpy as np
 import os
 import random
-import pyscal.core as pc
+import pyscal3.core as pc
 from mendeleev import element
 from ase.io import read, write
 from ase.atoms import Atoms
+from pyscal3.core import element_dict
 
 class CompositionTransformation:
     """
@@ -110,16 +111,12 @@ class CompositionTransformation:
     ```
     The output is written in LAMMPS dump format.
     """
-    def __init__(self, input_structure, input_chemical_composition, 
-        output_chemical_composition, restrictions=None):
+    def __init__(self, calc):
 
-        self.structure = self.prepare_structure(input_structure)
-        self.input_chemical_composition = input_chemical_composition
-        self.output_chemical_composition = output_chemical_composition
-        if restrictions is None:
-            self.restrictions = []
-        else:
-            self.restrictions = restrictions
+        self.input_chemical_composition = calc.composition_scaling._input_chemical_composition
+        self.output_chemical_composition = calc.composition_scaling.output_chemical_composition
+        self.restrictions = calc.composition_scaling.restrictions
+        self.calc = calc
         self.actual_species = None
         self.new_species = None
         self.maxtype = None        
@@ -136,55 +133,21 @@ def dict_to_string(self, inputdict):
             strlst.append(str(val))
         return "".join(strlst)
 
-    def is_data_file(self, filename):
-        try:
-            atoms = read(filename, format="lammps-data", style="atomic")
-            return atoms
-        except:
-            return None
-
-    def is_dump_file(self, filename):
-        try:
-            atoms = read(filename, format="lammps-dump-text")
-            return atoms
-        except:
-            return None
-
-    def prepare_structure(self, input_structure):
-        """
-        Check the format of a given input file and validate it.
-        """
-        if isinstance(input_structure, Atoms):
-            return input_structure
-        elif os.path.exists(input_structure):
-            atoms = self.is_data_file(input_structure)
-            if atoms is None:
-                atoms = self.is_dump_file(input_structure)
-            return atoms
 
     def convert_to_pyscal(self):
         """
         Convert a given system to pyscal and give a dict of type mappings
         """
-        pstruct = pc.System()
-        pstruct.read_inputfile(self.structure, format="ase")
+        aseobj = read(self.calc.lattice, format='lammps-data', style='atomic')
+        pstruct = pc.System(aseobj, format='ase')
 
         #here we have to validate the input composition dict; and map it
-        composition = pstruct.get_concentration()
-        atomsymbols = []
-        atomtypes = []
-        for key, val in self.input_chemical_composition.items():
-            if not val==0:
-                found = False
-                for key1, val1 in composition.items():
-                    if val1==val:
-                        found = True
-                        atomsymbols.append(key)
-                        atomtypes.append(int(key1))
-                        del composition[key1]
-                        break
-                if not found:
-                    raise ValueError("Input structure and composition do not match!")
+        typelist = pstruct.atoms.species
+        types, typecounts = np.unique(typelist, return_counts=True)
+        composition = {types[x]: typecounts[x] for x in range(len(types))}
+
+        atomsymbols = self.calc.element
+        atomtypes = [x+1 for x in range(len(self.calc.element))]
 
         self.pyscal_structure = pstruct
         self.typedict = dict(zip(atomsymbols, atomtypes))
@@ -196,14 +159,6 @@ def convert_to_pyscal(self):
         self.maxtype = self.actual_species + 1 #+ self.new_species
         #print(self.typedict)            
 
-    def check_if_atoms_conserved(self):
-        """
-        Check if a given transformation is possible by checking if number of atoms are conserved
-        """
-        natoms1 = np.sum([val for key, val in self.input_chemical_composition.items()])
-        natoms2 = np.sum([val for key, val in self.output_chemical_composition.items()])
-        if not (natoms1==natoms2==self.natoms):
-            raise ValueError(f"Input and output number of atoms are not conserved! Input {self.dict_to_string(self.input_chemical_composition)}, output {self.dict_to_string(self.output_chemical_composition)}, total atoms in structure {self.natoms}")
 
     def get_composition_transformation(self):
         """
@@ -237,7 +192,7 @@ def get_random_index_of_species(self, species):
     def mark_atoms(self):
         for i in range(self.natoms):
             self.atom_mark.append(False)
-            self.atom_type = [atom.type for atom in self.pyscal_structure.iter_atoms()]
+            self.atom_type = self.pyscal_structure.atoms.types
             self.mappings = [f"{x}-{x}" for x in self.atom_type]
 
     def update_mark_atoms(self):
@@ -345,10 +300,8 @@ def prepare_pair_lists(self):
     def update_types(self):
         for x in range(len(self.atom_type)):
             self.atom_type[x] = self.mappingdict[self.mappings[x]]
-        atoms = self.pyscal_structure.atoms
-        for count, atom in enumerate(atoms):
-            atom.type = self.atom_type[count]
-        self.pyscal_structure.atoms = atoms
+        for count in range(len(self.pyscal_structure.atoms.types)):
+            self.pyscal_structure.atoms.types[count] = self.atom_type[count]
 
     def iselement(self, symbol):
         try:
@@ -380,7 +333,17 @@ def update_pair_coeff(self, pair_coeff):
         return pc_old, pc_new
 
     def write_structure(self, outfilename):
-        self.pyscal_structure.to_file(outfilename)
+        #create some species dict
+        #just to trick ase to write
+        utypes = np.unique(self.pyscal_structure.atoms["types"])
+        element_list = list(element_dict.keys())       
+        element_type_dict = {str(u):element_list[count] for count, u in enumerate(utypes)}
+        species = [element_type_dict[str(x)] for x in self.pyscal_structure.atoms["types"]]
+        self.pyscal_structure.atoms["species"] = species
+        self.pyscal_structure.write.file(outfilename, format='lammps-data')
+
+
+
 
     def prepare_mappings(self):
         self.atom_mark = []
@@ -390,7 +353,6 @@ def prepare_mappings(self):
 
         self.get_composition_transformation()
         self.convert_to_pyscal()
-        self.check_if_atoms_conserved()
 
         self.mark_atoms()
         self.update_mark_atoms()

calphy/helpers.py

@@ -98,7 +98,7 @@ def create_object(cores, directory, timestep, cmdargs=None,
     return lmp
 
 
-def create_structure(lmp, calc, species=None):
+def create_structure(lmp, calc):
     """
     Create structure using LAMMPS
 
@@ -113,33 +113,7 @@ def create_structure(lmp, calc, species=None):
     -------
     lmp : LammpsLibrary object
     """
-    l, alat, apc, conc, dumpfile = pl.prepare_lattice(calc)
-
-    if species is None:
-        species = len(conc)
-
-    if l == "file":
-        if dumpfile:
-            # reset_timestep(calc.lattice, calc.lattice, keys=None)
-            lmp.command("lattice          fcc 4.0")
-            lmp.command("region           box block 0 2 0 2 0 2")
-
-            #comment this out
-            warnings.warn("If the box is triclinic, please provide a data file instead")
-            #lmp.command("box tilt large")
-
-            lmp.command("create_box       %d box" % species)
-            lmp.command(
-                "read_dump        %s 0 x y z scaled no box yes add keep" % calc.lattice
-            )
-            #lmp.command("change_box       all triclinic")
-        else:
-            lmp.command("read_data      %s" % calc.lattice)
-    else:
-        lmp.command(f'lattice {l} {alat}')
-        lmp.command(f'region box block 0 {calc.repeat[0]} 0 {calc.repeat[1]} 0 {calc.repeat[2]}')
-        lmp.command(f'create_box 1 box')
-        lmp.command(f'create_atoms 1 box')
+    lmp.command("read_data      %s" % calc.lattice)
     return lmp
 
 
@@ -171,7 +145,7 @@ def set_potential(lmp, options, ghost_elements=0):
     """
     #lmp.pair_style(options.pair_style_with_options[0])
     #lmp.pair_coeff(options.pair_coeff[0])
-    lmp.command(f'pair_style {options.pair_style_with_options[0]}')
+    lmp.command(f'pair_style {options._pair_style_with_options[0]}')
     lmp.command(f'pair_coeff {options.pair_coeff[0]}')
 
     lmp = set_mass(lmp, options, ghost_elements=ghost_elements)
@@ -234,15 +208,15 @@ def set_hybrid_potential(lmp, options, eps, ghost_elements=0):
         [
             *pcraw[:2],
             *[
-                options.pair_style_with_options[0],
+                options._pair_style_with_options[0],
             ],
             *pcraw[2:],
         ]
     )
 
     lmp.command(
         "pair_style       hybrid/overlay %s ufm 7.5"
-        % options.pair_style_with_options[0]
+        % options._pair_style_with_options[0]
     )
     lmp.command("pair_coeff       %s" % pcnew)
     lmp.command("pair_coeff       * * ufm %f 1.5" % eps)
@@ -303,7 +277,7 @@ def set_double_hybrid_potential(lmp, options, ghost_elements=0):
 
     lmp.command(
         "pair_style       hybrid/overlay %s %s"
-        % (options.pair_style_with_options[0], options.pair_style_with_options[1])
+        % (options._pair_style_with_options[0], options._pair_style_with_options[1])
     )
 
     lmp.command("pair_coeff       %s" % pcnew1)
@@ -414,16 +388,21 @@ def reset_timestep(conf, file="current.data", init_commands=None):
 """
 
 
-def prepare_log(file):
+def prepare_log(file, screen=False):
     logger = logging.getLogger(__name__)
     handler = logging.FileHandler(file)
     formatter = logging.Formatter("%(asctime)s %(name)-12s %(levelname)-8s %(message)s")
     handler.setFormatter(formatter)
     logger.addHandler(handler)
     logger.setLevel(logging.DEBUG)
     logger.propagate = False
-    return logger
 
+    if screen:
+        scr = logging.StreamHandler()
+        scr.setLevel(logging.INFO)
+        scr.setFormatter(formatter)
+        logger.addHandler(scr)
+    return logger
 
 def check_if_any_is_none(data):
     """