Merge pull request #165 from ICAMS/add_swap

Add swap

Author: srmnitc

calphy/alchemy.py

@@ -19,6 +19,10 @@
 
 For more information contact:
 [email protected]/[email protected]
+
+Notes
+-----
+- swapping is strictly only performed between types 1 and 2 at the moment; this needs to be refined further
 """
 
 import numpy as np
@@ -159,6 +163,9 @@ def run_integration(self, iteration=1):
         lmp.command(f'pair_coeff {self.calc.pair_coeff[0]}')
         lmp = ph.set_mass(lmp, self.calc)
 
+        #NEW ADDED
+        lmp.command("group g1 type 1")
+        lmp.command("group g2 type 2")
         #lmp = ph.set_double_hybrid_potential(lmp, self.options, self.calc._pressureair_style, self.calc._pressureair_coeff)
 
         #remap the box to get the correct pressure
@@ -173,9 +180,12 @@ def run_integration(self, iteration=1):
             lmp.command("fix             f1 all nvt temp %f %f %f"%(self.calc._temperature, self.calc._temperature, 
                 self.calc.md.thermostat_damping[1]))
 
+        
         lmp.command("thermo_style    custom step pe")
         lmp.command("thermo          1000")
         lmp.command("run             %d"%self.calc.n_equilibration_steps)
+        
+
         #equilibration run is over
 
         #---------------------------------------------------------------
@@ -226,29 +236,49 @@ def run_integration(self, iteration=1):
         lmp.command("variable        dU1 equal c_c1/atoms")             # Driving-force obtained from NEHI procedure.
         lmp.command("variable        dU2 equal c_c2/atoms")
 
+        #add swaps if n_swap is > 0
+        if self.calc.monte_carlo.n_swaps > 0:
+            self.logger.info(f'{self.calc.monte_carlo.n_swaps} swap moves are performed between {self.calc.monte_carlo.swap_types[0]} and {self.calc.monte_carlo.swap_types[1]} every {self.calc.monte_carlo.n_steps}')
+            lmp.command("fix  swap all atom/swap %d %d %d %f ke no types %d %d"%(self.calc.monte_carlo.n_steps,
+                                                                                self.calc.monte_carlo.n_swaps,
+                                                                                np.random.randint(1, 10000),
+                                                                                self.calc._temperature,
+                                                                                self.calc.monte_carlo.swap_types[0],
+                                                                                self.calc.monte_carlo.swap_types[1]))
+            lmp.command("variable a equal f_swap[1]")
+            lmp.command("variable b equal f_swap[2]")
+            lmp.command("fix             swap2 all print 1 \"${a} ${b} ${flambda}\" screen no file swap.forward_%d.dat"%iteration)
+
         # Thermo output.
-        lmp.command("thermo_style    custom step v_dU1 v_dU2")
+        if self.calc.monte_carlo.n_swaps > 0:
+            lmp.command("thermo_style    custom step v_dU1 v_dU2 v_a v_b")
+        else:
+            lmp.command("thermo_style    custom step v_dU1 v_dU2")
         lmp.command("thermo          1000")
 
-
+        
         #save the necessary items to a file: first step
         lmp.command("fix             f2 all print 1 \"${dU1} ${dU2} ${flambda}\" screen no file forward_%d.dat"%iteration)
-        lmp.command("run             %d"%self.calc._n_switching_steps)
-
+        lmp.command("run             %d"%self.calc._n_switching_steps)            
 
         #now equilibrate at the second potential
         lmp.command("unfix           f2")
         lmp.command("uncompute       c1")
         lmp.command("uncompute       c2")
 
+        #NEW SWAP
+        if self.calc.monte_carlo.n_swaps > 0:
+            lmp.command("unfix swap")
+            lmp.command("unfix swap2")
+
 
         lmp.command("pair_style      %s"%self.calc._pair_style_with_options[1])
         lmp.command("pair_coeff      %s"%self.calc.pair_coeff[1])
 
         # Thermo output.
         lmp.command("thermo_style    custom step pe")
         lmp.command("thermo          1000")
-
+        
         #run eqbrm run
         lmp.command("run             %d"%self.calc.n_equilibration_steps)
 
@@ -278,11 +308,39 @@ def run_integration(self, iteration=1):
         lmp.command("variable        dU1 equal c_c1/atoms")             # Driving-force obtained from NEHI procedure.
         lmp.command("variable        dU2 equal c_c2/atoms")
 
+        #add swaps if n_swap is > 0
+        if self.calc.monte_carlo.n_swaps > 0:
+            if self.calc.monte_carlo.reverse_swap:
+                self.logger.info(f'{self.calc.monte_carlo.n_swaps} swap moves are performed between {self.calc.monte_carlo.swap_types[1]} and {self.calc.monte_carlo.swap_types[0]} every {self.calc.monte_carlo.n_steps}')
+                lmp.command("fix  swap all atom/swap %d %d %d %f ke no types %d %d"%(self.calc.monte_carlo.n_steps,
+                                                                                    self.calc.monte_carlo.n_swaps,
+                                                                                    np.random.randint(1, 10000),
+                                                                                    self.calc._temperature,
+                                                                                    self.calc.monte_carlo.swap_types[1],
+                                                                                    self.calc.monte_carlo.swap_types[0]))
+            else:
+                self.logger.info(f'{self.calc.monte_carlo.n_swaps} swap moves are performed between {self.calc.monte_carlo.swap_types[0]} and {self.calc.monte_carlo.swap_types[1]} every {self.calc.monte_carlo.n_steps}')
+                self.logger.info('note that swaps are not reversed')
+                lmp.command("fix  swap all atom/swap %d %d %d %f ke no types %d %d"%(self.calc.monte_carlo.n_steps,
+                                                                                    self.calc.monte_carlo.n_swaps,
+                                                                                    np.random.randint(1, 10000),
+                                                                                    self.calc._temperature,
+                                                                                    self.calc.monte_carlo.swap_types[0],
+                                                                                    self.calc.monte_carlo.swap_types[1]))
+
+            lmp.command("variable a equal f_swap[1]")
+            lmp.command("variable b equal f_swap[2]")
+            lmp.command("fix             swap2 all print 1 \"${a} ${b} ${blambda}\" screen no file swap.backward_%d.dat"%iteration)
+
         # Thermo output.
-        lmp.command("thermo_style    custom step v_dU1 v_dU2")
+        if self.calc.monte_carlo.n_swaps > 0:
+            lmp.command("thermo_style    custom step v_dU1 v_dU2 v_a v_b")
+        else:        
+            lmp.command("thermo_style    custom step v_dU1 v_dU2")
         lmp.command("thermo          1000")
 
 
+        
         #save the necessary items to a file: first step
         lmp.command("fix             f2 all print 1 \"${dU1} ${dU2} ${flambda}\" screen no file backward_%d.dat"%iteration)
         lmp.command("run             %d"%self.calc._n_switching_steps)
@@ -293,6 +351,8 @@ def run_integration(self, iteration=1):
         lmp.command("uncompute       c1")
         lmp.command("uncompute       c2")
 
+        if self.calc.monte_carlo.n_swaps > 0:
+            lmp.command("unfix  swap")
         lmp.close()
 
 
@@ -329,8 +389,6 @@ def thermodynamic_integration(self):
             return flambda_arr, w_arr, q_arr, qerr_arr
 
 
-
-
     def mass_integration(self, flambda, ref_mass, target_masses, target_counts):
         mcorarr, mcorsum = integrate_mass(flambda, ref_mass, target_masses, target_counts,
     self.calc._temperature, self.natoms)

calphy/clitools.py

@@ -11,6 +11,7 @@
 from calphy.liquid import Liquid
 from calphy.solid import Solid
 from calphy.alchemy import Alchemy
+from calphy.phase_diagram import prepare_inputs_for_phase_diagram
 
 def _generate_job(calc, simfolder):
     if calc.mode == "alchemy" or calc.mode == "composition_scaling":
@@ -120,3 +121,8 @@ def convert_legacy_inputfile():
             yaml.safe_dump(data, fout)
 
 
+def phase_diagram():
+    arg = ap.ArgumentParser()
+    arg.add_argument("-i", "--input", required=True, type=str,
+    help="name of the input file")
+    prepare_inputs_for_phase_diagram(args['input'])

calphy/composition_transformation.py

@@ -30,6 +30,7 @@
 from ase.io import read, write
 from ase.atoms import Atoms
 from pyscal3.core import element_dict
+from calphy.integrators import kb
 
 class CompositionTransformation:
     """
@@ -124,6 +125,7 @@ def __init__(self, calc):
         self.atom_species = None
         self.mappings = None
         self.unique_mappings = None
+        self.mappingdict = None
         self.prepare_mappings()
 
     def dict_to_string(self, inputdict):
@@ -132,8 +134,32 @@ def dict_to_string(self, inputdict):
             strlst.append(str(key))
             strlst.append(str(val))
         return "".join(strlst)
-    
-    
+
+    @property
+    def entropy_contribution(self):
+        """
+        Find the entropy entribution of the transformation. To get
+        free energies, multiply by -T.
+        """
+        def _log(val):
+            if val == 0:
+                return 0
+            else:
+                return np.log(val)
+        ents = []
+        for key, val in self.output_chemical_composition.items():
+            if key in self.input_chemical_composition.keys():
+                t1 = self.input_chemical_composition[key]/self.natoms
+                t2 = self.output_chemical_composition[key]/self.natoms
+                cont = t2*_log(t2) - t1*_log(t1)
+            else:
+                t1 = 0
+                t2 = self.output_chemical_composition[key]/self.natoms
+                cont =  t2*_log(t2) - 0
+            ents.append(cont)
+        entropy_term = kb*np.sum(ents)
+        return entropy_term
+        
     def convert_to_pyscal(self):
         """
         Convert a given system to pyscal and give a dict of type mappings
@@ -159,7 +185,6 @@ def convert_to_pyscal(self):
         self.maxtype = self.actual_species + 1 #+ self.new_species
         #print(self.typedict)            
 
-
     def get_composition_transformation(self):
         """
         From the two given composition transformation, find the transformation dict
@@ -192,9 +217,10 @@ 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 = self.pyscal_structure.atoms.types
-            self.mappings = [f"{x}-{x}" for x in self.atom_type]
-            
+        
+        self.atom_type = self.pyscal_structure.atoms.types
+        self.mappings = [f"{x}-{x}" for x in self.atom_type]
+     
     def update_mark_atoms(self):
         self.marked_atoms = []
         for key, val in self.to_remove.items():
@@ -281,7 +307,7 @@ def get_mappings(self):
         self.update_typedicts()
         self.compute_possible_mappings()
         self.update_mappings()
-        
+    
     def prepare_pair_lists(self):
         self.pair_list_old = []
         self.pair_list_new = []
@@ -300,8 +326,12 @@ 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]]
-        for count in range(len(self.pyscal_structure.atoms.types)):
-            self.pyscal_structure.atoms.types[count] = self.atom_type[count]
+        
+        #smartify these loops
+        #npyscal = len(self.pyscal_structure.atoms.types)
+        self.pyscal_structure.atoms.types = self.atom_type
+        #for count in range(npyscal)):
+        #    self.pyscal_structure.atoms.types[count] = self.atom_type[count]
 
     def iselement(self, symbol):
         try:
@@ -331,6 +361,30 @@ def update_pair_coeff(self, pair_coeff):
         pc_old = " ".join([*pc_before, *self.pair_list_old, *pc_after])
         pc_new = " ".join([*pc_before, *self.pair_list_new, *pc_after])
         return pc_old, pc_new
+
+    def get_swap_types(self):
+        """
+        Get swapping types
+        """
+        swap_list = []
+        for mapping in self.unique_mappings:
+            map_split = mapping.split("-")
+            #conserved atom
+            if (map_split[0]==map_split[1]):
+                pass
+            else:
+                first_type = map_split[0]
+                second_type = map_split[1]
+                first_map = f'{first_type}-{first_type}'
+                second_map = mapping
+
+                #get the numbers from dict
+                first_swap_type = self.mappingdict[first_map]
+                second_swap_type = self.mappingdict[second_map]
+
+                swap_list.append([first_swap_type, second_swap_type])
+        return swap_list[0]
+            
 
     def write_structure(self, outfilename):
         #create some species dict

calphy/helpers.py

@@ -121,8 +121,12 @@ def create_structure(lmp, calc):
 
 
 def set_mass(lmp, options):
-    for i in range(options.n_elements):
-        lmp.command(f'mass {i+1} {options.mass[i]}')
+    if options.mode == 'composition_scaling':
+        lmp.command(f'mass * {options.mass[-1]}')
+
+    else:
+        for i in range(options.n_elements):
+            lmp.command(f'mass {i+1} {options.mass[i]}')
     return lmp
 
 

calphy/input.py

@@ -27,6 +27,7 @@
 from pydantic.functional_validators import AfterValidator, BeforeValidator
 from annotated_types import Len
 import mendeleev
+from tqdm import tqdm
 
 import yaml
 import numpy as np
@@ -72,7 +73,14 @@ def _to_float(val):
         return float(val)
     else:
         return [float(x) for x in val] 
-        
+
+class MonteCarlo(BaseModel, title='Options for Monte Carlo moves during particle swap moves'):
+    n_steps: Annotated[int, Field(default=1, gt=0, description='perform swap moves every n_step')]
+    n_swaps: Annotated[int, Field(default=0, ge=0, description='number of swap moves to perform')]
+    reverse_swap: Annotated[ bool, Field(default=True)]
+    swap_types: Annotated[conlist(int, min_length=2, max_length=2), Field(default=[1,2], 
+                description='which atoms to swap between')]
+
 class CompositionScaling(BaseModel, title='Composition scaling input options'):
     _input_chemical_composition: PrivateAttr(default=None)
     output_chemical_composition: Annotated[dict, Field(default={}, required=False)]
@@ -125,6 +133,7 @@ class MeltingTemperature(BaseModel, title='Input options for melting temperature
     attempts: Annotated[int, Field(default=5, ge=1)]
 
 class Calculation(BaseModel, title='Main input class'):
+    monte_carlo: Optional[MonteCarlo] = MonteCarlo()
     composition_scaling: Optional[CompositionScaling] = CompositionScaling()
     md: Optional[MD] = MD()
     nose_hoover: Optional[NoseHoover] = NoseHoover()
@@ -197,9 +206,16 @@ class Calculation(BaseModel, title='Main input class'):
     #add second level options; for example spring constants
     spring_constants: Annotated[Union[List[float],None], Field(default = None)]
 
+    #some input keywords that will be used for the phase diagram mode
+    phase_name: Annotated[str, Field(default="")]
+    reference_composition: Annotated[float, Field(default=0.00)]
+
     #structure items
     _structure: Any = PrivateAttr(default=None)
 
+    #just check for nlements in compscale
+    _totalelements = PrivateAttr(default=0)
+
     @model_validator(mode='after')
     def _validate_all(self) -> 'Input':
 
@@ -326,6 +342,7 @@ def _validate_all(self) -> 'Input':
 
         #generate temporary filename if needed
         write_structure_file = False
+        rename_structure_file = False
 
         if self.lattice == "":
             #fetch from dict
@@ -399,6 +416,7 @@ def _validate_all(self) -> 'Input':
                 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')
+                rename_structure_file = True
             else:
                 raise TypeError('Only lammps-data files are supported!')                
 
@@ -420,6 +438,12 @@ def _validate_all(self) -> 'Input':
             write(structure_filename, structure, format='lammps-data')
             self.lattice = structure_filename
 
+        if rename_structure_file:
+            structure_filename = ".".join([self.create_identifier(), str(self.kernel), "data"])
+            structure_filename = os.path.join(os.getcwd(), structure_filename)
+            shutil.copy(self.lattice, structure_filename)
+            self.lattice = structure_filename
+        
         if self.mode == 'composition_scaling':
             #we also should check if actual contents are present
             input_chem_comp = {}
@@ -548,7 +572,7 @@ def _read_inputfile(file):
     with open(file, 'r') as fin:
         data = yaml.safe_load(fin)
     calculations = []
-    for count, calc in enumerate(data['calculations']):
+    for count, calc in enumerate(tqdm(data['calculations'])):
         calc['kernel'] = count
         calc['inputfile'] = file
         if 'pressure' in calc.keys():