Merge pull request #108 from ICAMS/standardise_pair_styles

Standardise pair styles

calphy/alchemy.py

@@ -84,7 +84,7 @@ def run_averaging(self):
 
         #set up potential
         lmp.command(f'pair_coeff {self.calc.pair_coeff[0]}')
-        lmp = ph.set_mass(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
+        lmp = ph.set_mass(lmp, self.calc)
 
         #add some computes
         lmp.command("variable         mvol equal vol")
@@ -157,7 +157,7 @@ def run_integration(self, iteration=1):
         #set up hybrid potential
         #here we only need to set one potential
         lmp.command(f'pair_coeff {self.calc.pair_coeff[0]}')
-        lmp = ph.set_mass(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
+        lmp = ph.set_mass(lmp, self.calc)
 
         #lmp = ph.set_double_hybrid_potential(lmp, self.options, self.calc._pressureair_style, self.calc._pressureair_coeff)
 

calphy/helpers.py

@@ -118,19 +118,13 @@ def create_structure(lmp, calc):
     return lmp
 
 
-def set_mass(lmp, options, ghost_elements=0):
-    count = 1
+def set_mass(lmp, options):
     for i in range(options.n_elements):
         lmp.command(f'mass {i+1} {options.mass[i]}')
-        count += 1
-
-    for i in range(ghost_elements):
-        lmp.command(f'mass {count+i} 1.00')
-
     return lmp
 
 
-def set_potential(lmp, options, ghost_elements=0):
+def set_potential(lmp, options):
     """
     Set the interatomic potential
 
@@ -149,7 +143,7 @@ def set_potential(lmp, options, ghost_elements=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)
+    lmp = set_mass(lmp, options)
 
     return lmp
 

calphy/input.py

@@ -202,7 +202,6 @@ class Calculation(BaseModel, title='Main input class'):
 
     #add second level options; for example spring constants
     spring_constants: Annotated[List[float], Field(default = None)]
-    _ghost_element_count: int = PrivateAttr(default=0)
 
     #structure items
     _structure: Any = PrivateAttr(default=None)
@@ -329,14 +328,15 @@ def _validate_all(self) -> 'Input':
             self._element_dict[element]['mass'] = self.mass[count]
             self._element_dict[element]['count'] = 0
             self._element_dict[element]['composition'] = 0.0
+            self._element_dict[element]['atomic_number'] = mendeleev.element(element).atomic_number
 
         #generate temporary filename if needed
         write_structure_file = False
 
         if self.lattice == "":
             #fetch from dict
             if len(self.element) > 1:
-                raise ValueError("Cannot create lattice for more than one element")
+                raise ValueError("Cannot create lattice for more than one element, provide a lammps-data file explicitly")
             if self.element[0] in element_dict.keys():
                 self.lattice = element_dict[self.element[0]]['structure']
                 self.lattice_constant = element_dict[self.element[0]]['lattice_constant']
@@ -350,20 +350,23 @@ def _validate_all(self) -> 'Input':
                 lattice_constant=self.lattice_constant,
                 repetitions=self.repeat,
                 element=self.element)
+            structure = structure.write.ase()
 
             #extract composition
-            typelist = structure.atoms.species
-            types, typecounts = np.unique(typelist, return_counts=True)
+            types, typecounts = np.unique(structure.get_chemical_symbols(), return_counts=True)
 
             for c, t in enumerate(types):
                 self._element_dict[t]['count'] = typecounts[c]
                 self._element_dict[t]['composition'] = typecounts[c]/np.sum(typecounts)
 
-            self._natoms = structure.natoms
+            self._natoms = len(structure)
             self._original_lattice = self.lattice.lower()
             write_structure_file = True
 
         elif self.lattice.lower() in structure_dict.keys():
+            if len(self.element) > 1:
+                raise ValueError("Cannot create lattice for more than one element, provide a lammps-data file explicitly")
+
             #this is a valid structure
             if self.lattice_constant == 0:
                 #we try try to get lattice_constant
@@ -376,10 +379,10 @@ def _validate_all(self) -> 'Input':
                 lattice_constant=self.lattice_constant,
                 repetitions=self.repeat,
                 element=self.element)
+            structure = structure.write.ase()
 
             #extract composition
-            typelist = structure.atoms.species
-            types, typecounts = np.unique(typelist, return_counts=True)
+            types, typecounts = np.unique(structure.get_chemical_symbols(), return_counts=True)
 
             for c, t in enumerate(types):
                 self._element_dict[t]['count'] = typecounts[c]
@@ -389,7 +392,7 @@ def _validate_all(self) -> 'Input':
             #concdict_frac = {str(t): typecounts[c]/np.sum(typecounts) for c, t in enumerate(types)}
             #self._composition = concdict_frac
             #self._composition_counts = concdict_counts
-            self._natoms = structure.natoms
+            self._natoms = len(structure)
             self._original_lattice = self.lattice.lower()
             write_structure_file = True
 
@@ -398,34 +401,32 @@ def _validate_all(self) -> 'Input':
             if not os.path.exists(self.lattice):
                 raise ValueError(f'File {self.lattice} could not be found')
             if self.file_format == 'lammps-data':
-                aseobj = read(self.lattice, format='lammps-data', style='atomic')
-                structure = System(aseobj, format='ase')
+                #create atomic numbers for proper reading
+                Z_of_type = dict([(count+1, self._element_dict[element]['atomic_number']) for count, element in enumerate(self.element)])
+                structure = read(self.lattice, format='lammps-data', style='atomic', Z_of_type=Z_of_type)
+                #structure = System(aseobj, format='ase')
             else:
                 raise TypeError('Only lammps-data files are supported!')                
 
             #extract composition
-            typelist = structure.atoms.types
-            #convert to species
-            typelist = [self.element[x-1] for x in typelist]
-            types, typecounts = np.unique(typelist, return_counts=True)
+            #this is the types read in from the file
+            types, typecounts = np.unique(structure.get_chemical_symbols(), return_counts=True)
             for c, t in enumerate(types):
                 self._element_dict[t]['count'] = typecounts[c]
                 self._element_dict[t]['composition'] = typecounts[c]/np.sum(typecounts)
-            self._natoms = structure.natoms
+
+            self._natoms = len(structure)
             self._original_lattice = os.path.basename(self.lattice)
             self.lattice = os.path.abspath(self.lattice)
 
         #if needed, write the file out
         if write_structure_file:
             structure_filename = ".".join([self.create_identifier(), str(self.kernel), "data"])
             structure_filename = os.path.join(os.getcwd(), structure_filename)
-            structure.write.file(structure_filename, format='lammps-data')
+            write(structure_filename, structure, format='lammps-data')
             self.lattice = structure_filename
 
         if self.mode == 'composition_scaling':
-            aseobj = read(self.lattice, format='lammps-data', style='atomic')
-            structure = System(aseobj, format='ase')
-
             #we also should check if actual contents are present
             input_chem_comp = {}
             for key, val in self._element_dict.items():
@@ -439,7 +440,7 @@ def _validate_all(self) -> 'Input':
 
             natoms1 = np.sum([val for key, val in self.composition_scaling._input_chemical_composition.items()])
             natoms2 = np.sum([val for key, val in self.composition_scaling.output_chemical_composition.items()])
-            if not (natoms1==natoms2==structure.natoms):
+            if not (natoms1==natoms2):
                 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 {structure.natoms}")
         return self
 

calphy/integrators.py

@@ -46,10 +46,12 @@
 #             TI PATH INTEGRATION ROUTINES
 #--------------------------------------------------------------------
 
-def integrate_path(fwdfilename, bkdfilename, 
-    nelements=1, concentration=[1,], 
-    usecols=(0, 1, 2), solid=True,
-    alchemy=False, composition_integration=False):
+def integrate_path(calc,
+    fwdfilename, 
+    bkdfilename,  
+    solid=True,
+    alchemy=False, 
+    composition_integration=False):
     """
     Get a filename with columns du and dlambda and integrate
 
@@ -72,23 +74,36 @@ def integrate_path(fwdfilename, bkdfilename,
     q : float
         heat dissipation during switching of system
     """
+    natoms = np.array([calc._element_dict[x]['count'] for x in calc.element])
+    concentration = np.array([calc._element_dict[x]['composition'] for x in calc.element])
+    fdata = np.loadtxt(fwdfilename, unpack=True, comments="#")
+    bdata = np.loadtxt(bkdfilename, unpack=True, comments="#")
+
     if solid:
-        fdui = np.loadtxt(fwdfilename, unpack=True, comments="#", usecols=(0,))
-        bdui = np.loadtxt(bkdfilename, unpack=True, comments="#", usecols=(0,))
+        fdui = fdata[0]
+        bdui = bdata[0]
 
         fdur = np.zeros(len(fdui))
         bdur = np.zeros(len(bdui))
 
-        for i in range(nelements):
-            fdur += concentration[i]*np.loadtxt(fwdfilename, unpack=True, comments="#", usecols=(i+1,))
-            bdur += concentration[i]*np.loadtxt(bkdfilename, unpack=True, comments="#", usecols=(i+1,))
+        for i in range(calc.n_elements):
+            if natoms[i] > 0:
+                fdur += concentration[i]*fdata[i+1]/natoms[i]
+                bdur += concentration[i]*bdata[i+1]/natoms[i]
 
-        flambda = np.loadtxt(fwdfilename, unpack=True, comments="#", usecols=(nelements+1,))
-        blambda = np.loadtxt(bkdfilename, unpack=True, comments="#", usecols=(nelements+1,))
+        flambda = fdata[calc.n_elements+1]
+        blambda = bdata[calc.n_elements+1]
 
     else:
-        fdui, fdur, flambda = np.loadtxt(fwdfilename, unpack=True, comments="#", usecols=usecols)
-        bdui, bdur, blambda = np.loadtxt(bkdfilename, unpack=True, comments="#", usecols=usecols)
+        fdui = fdata[0]
+        bdui = bdata[0]
+
+        fdur = fdata[1]
+        bdur = bdata[1]
+
+        flambda = fdata[2]
+        blambda = bdata[2]
+
 
     fdu = fdui - fdur
     bdu = bdui - bdur
@@ -106,11 +121,12 @@ def integrate_path(fwdfilename, bkdfilename,
     return w, q, flambda
 
 
-def find_w(mainfolder, nelements=1, 
-    concentration=[1,], nsims=5, 
-    full=False, usecols=(0,1,2), 
+def find_w(mainfolder,
+    calc, 
+    full=False, 
     solid=True,
-    alchemy=False, composition_integration=False):
+    alchemy=False, 
+    composition_integration=False):
     """
     Integrate the irreversible work and dissipation for independent simulations
 
@@ -142,15 +158,20 @@ def find_w(mainfolder, nelements=1,
     ws = []
     qs = []
 
-    for i in range(nsims):
+    for i in range(calc.n_iterations):
         fwdfilestring = 'forward_%d.dat' % (i+1)
         fwdfilename = os.path.join(mainfolder,fwdfilestring)
+
         bkdfilestring = 'backward_%d.dat' % (i+1)
         bkdfilename = os.path.join(mainfolder,bkdfilestring)
-        w, q, flambda = integrate_path(fwdfilename, bkdfilename, 
-            nelements=nelements, concentration=concentration, 
-            usecols=usecols, solid=solid,
-            alchemy=alchemy, composition_integration=composition_integration)
+
+        w, q, flambda = integrate_path(calc,
+            fwdfilename, 
+            bkdfilename, 
+            solid=solid,
+            alchemy=alchemy, 
+            composition_integration=composition_integration)
+
         ws.append(w)
         qs.append(q)
 
@@ -444,7 +465,11 @@ def integrate_dcc(folder1, folder2, nsims=1, scale_energy=True,
 #             REF. STATE ROUTINES: SOLID
 #--------------------------------------------------------------------
 
-def get_einstein_crystal_fe(temp, natoms, mass, vol, k, concentration, cm_correction=True):
+def get_einstein_crystal_fe(
+    calc,
+    vol, 
+    k, 
+    cm_correction=True):
     """
     Get the free energy of einstein crystal
 
@@ -475,7 +500,10 @@ def get_einstein_crystal_fe(temp, natoms, mass, vol, k, concentration, cm_correc
 
     """
     #convert mass first for single particle in kg
-    mass = (np.array(mass)/Na)*1E-3
+    mass = np.array([calc._element_dict[x]['mass'] for x in calc.element])
+    mass = (mass/Na)*1E-3
+    natoms = np.sum([calc._element_dict[x]['count'] for x in calc.element])
+    concentration = np.array([calc._element_dict[x]['composition'] for x in calc.element])
 
     #convert k from ev/A2 to J/m2
     k = np.array(k)*(eV2J/1E-20)
@@ -488,12 +516,13 @@ def get_einstein_crystal_fe(temp, natoms, mass, vol, k, concentration, cm_correc
     F_cm = 0
 
     for count, om in enumerate(omega):
-        F_harm += concentration[count]*np.log((hbar*om)/(kb*temp))
-        if cm_correction:
-            F_cm += np.log((natoms*concentration[count]/vol)*(2*np.pi*kbJ*temp/(natoms*concentration[count]*k[count]))**1.5)
+        if concentration[count] > 0:
+            F_harm += concentration[count]*np.log((hbar*om)/(kb*calc._temperature))
+            if cm_correction:
+                F_cm += np.log((natoms*concentration[count]/vol)*(2*np.pi*kbJ*calc._temperature/(natoms*concentration[count]*k[count]))**1.5)
             #F_cm = 0
-    F_harm = 3*kb*temp*F_harm
-    F_cm = (kb*temp/natoms)*F_cm
+    F_harm = 3*kb*calc._temperature*F_harm
+    F_cm = (kb*calc._temperature/natoms)*F_cm
 
     F_harm = F_harm + F_cm
 

calphy/liquid.py

@@ -120,7 +120,7 @@ def run_averaging(self):
 
         #set up potential
         lmp.command(f'pair_coeff {self.calc.pair_coeff[0]}')
-        lmp = ph.set_mass(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
+        lmp = ph.set_mass(lmp, self.calc)
 
         #Melt regime for the liquid
         lmp.velocity("all create", self.calc._temperature_high, np.random.randint(1, 10000))
@@ -189,7 +189,7 @@ def run_integration(self, iteration=1):
         #set hybrid ufm and normal potential
         #lmp = ph.set_hybrid_potential(lmp, self.options, self.eps)
         lmp.command(f'pair_coeff {self.calc.pair_coeff[0]}')
-        lmp = ph.set_mass(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
+        lmp = ph.set_mass(lmp, self.calc)
 
         #remap the box to get the correct pressure
         lmp = ph.remap_box(lmp, self.lx, self.ly, self.lz)
@@ -331,7 +331,8 @@ def thermodynamic_integration(self):
             self.rho, 50, 1.5)
 
         #Get ideal gas fe
-        f2 = get_ideal_gas_fe(self.calc._temperature, self.rho, 
+        f2 = get_ideal_gas_fe(self.calc._temperature, 
+            self.rho, 
             self.natoms, 
             [val['mass'] for key, val in self.calc._element_dict.items()], 
             [val['composition'] for key, val in self.calc._element_dict.items()])