Merge pull request #100 from ICAMS/fix_automated_tm

Fix automated tm

Author: srmnitc

calphy/input.py

@@ -26,6 +26,7 @@
 from pydantic import BaseModel, Field, ValidationError, model_validator, conlist, PrivateAttr
 from pydantic.functional_validators import AfterValidator, BeforeValidator
 from annotated_types import Len
+import mendeleev
 
 import yaml
 import numpy as np
@@ -127,6 +128,7 @@ class Calculation(BaseModel, title='Main input class'):
                                       Field(default=[])]
     _element_dict: dict = PrivateAttr(default={})
     kernel: Annotated[int, Field(default=0)]
+    inputfile: Annotated[str, Field(default='')]
 
     mode: Annotated[str, Field(default=None)]
     lattice: Annotated[str, Field(default="")]
@@ -235,14 +237,38 @@ def _validate_all(self) -> 'Input':
             raise ValueError('Unknown format for pressure')
 
         self._temperature_input = copy.copy(self.temperature)
-        if self.temperature is None:
-            self._temperature = None
+        #guess a melting temp of the system, this will be mostly ignored
+        #chem = mendeleev.element(self.element[0])
+        #self._melting_temperature = chem.melting_point
+        try:
+            chem = mendeleev.element(self.element[0])
+            self._melting_temperature = chem.melting_point
+        except:
+            self._melting_temperature = None
+
+        if self.temperature == 0:
+            #the only situation in which it can be None is if mode is melting temp
+            if len(self.element) > 1:
+                raise ValueError("Cannot guess start temperature for more than one species, please specify")
+            #now try to guess
+            if self._melting_temperature is None:
+                raise ValueError("Could not guess start temperature for more than one species, please specify")
+            self._temperature = self._melting_temperature
+            self._temperature_stop = self._melting_temperature
+            if self.temperature_high == 0:
+                self._temperature_high = 2*self._melting_temperature
+            else:
+                self._temperature_high = self.temperature_high
+
         elif np.shape(np.atleast_1d(self.temperature)) == (1,):
             temp = np.atleast_1d(self.temperature)
             self._temperature = temp[0]
             self._temperature_stop = temp[0]
-            if self._temperature_high is None:
+            if self.temperature_high == 0:
                 self._temperature_high = 2*temp[0]
+            else:
+                self._temperature_high = self.temperature_high
+        
         elif np.shape(self.temperature) == (2,):
             temp = self.temperature
             self._temperature = temp[0]
@@ -293,16 +319,16 @@ def _validate_all(self) -> 'Input':
         if self.lattice == "":
             #fetch from dict
             if len(self.element) > 1:
-                raise ValueError("MeltingTemperature can be used only with one element")
+                raise ValueError("Cannot create lattice for more than one element")
             if self.element[0] in element_dict.keys():
-                self.lattice = element_dict[self.element[0]]['structure']           
+                self.lattice = element_dict[self.element[0]]['structure']
                 self.lattice_constant = element_dict[self.element[0]]['lattice_constant']
             else:
                 raise ValueError("Could not find structure, please provide lattice and lattice_constant explicitely")                
 
             if self.repeat == [1,1,1]:
                 self.repeat = [5,5,5]
-                
+
             structure = _make_crystal(self.lattice.lower(),
                 lattice_constant=self.lattice_constant,
                 repetitions=self.repeat,
@@ -317,7 +343,7 @@ def _validate_all(self) -> 'Input':
                 self._element_dict[t]['composition'] = typecounts[c]/np.sum(typecounts)
 
             self._natoms = structure.natoms
-            self._original_lattice = self.lattice
+            self._original_lattice = self.lattice.lower()
             write_structure_file = True
 
         elif self.lattice.lower() in structure_dict.keys():
@@ -347,7 +373,7 @@ def _validate_all(self) -> 'Input':
             #self._composition = concdict_frac
             #self._composition_counts = concdict_counts
             self._natoms = structure.natoms
-            self._original_lattice = self.lattice
+            self._original_lattice = self.lattice.lower()
             write_structure_file = True
 
         else:
@@ -369,7 +395,7 @@ def _validate_all(self) -> 'Input':
                 self._element_dict[t]['count'] = typecounts[c]
                 self._element_dict[t]['composition'] = typecounts[c]/np.sum(typecounts)
             self._natoms = structure.natoms
-            self._original_lattice = self.lattice
+            self._original_lattice = os.path.basename(self.lattice)
             self.lattice = os.path.abspath(self.lattice)
 
         #if needed, write the file out
@@ -435,19 +461,14 @@ def create_identifier(self):
         """
         #lattice processed
         prefix = self.mode
-        if prefix == 'melting_temperature':
-            ts = int(0)
-            ps = int(0)
-            l = 'tm'
+        ts = int(self._temperature)
+        if self._pressure is None:
+            ps = "None"
         else:
-            ts = int(self._temperature)
-            if self._pressure is None:
-                ps = "None"
-            else:
-                ps = "%d"%(int(self._pressure))
-            l = self._original_lattice
-            l = l.split('/')
-            l = l[-1]
+            ps = "%d"%(int(self._pressure))
+        l = self._original_lattice
+        l = l.split('/')
+        l = l[-1]
 
         if self.folder_prefix is None:
             identistring = "-".join([prefix, l.lower(), self.reference_phase, str(ts), str(ps)])
@@ -524,6 +545,7 @@ def _read_inputfile(file):
     calculations = []
     for count, calc in enumerate(data['calculations']):
         calc['kernel'] = count
+        calc['inputfile'] = file
         calculations.append(Calculation(**calc))
     return calculations
 
@@ -559,6 +581,7 @@ def _convert_legacy_inputfile(file, return_calcs=False):
             #calc['reference_phase'] = str(ci["reference_phase"]) if "reference_phase" in ci.keys() else 'none'
             #calc['lattice_constant'] = float(ci["lattice_constant"]) if "lattice_constant" in ci.keys() else 0
             calc['kernel'] = cc
+            calc['inputfile'] = file
             calculations.append(calc)
 
         else:
@@ -608,6 +631,7 @@ def _convert_legacy_inputfile(file, return_calcs=False):
                 calc["pressure"] = _to_float(combo[1])
                 calc["temperature"] = _to_float(combo[2])
                 calc['kernel'] = cc
+                calc['inputfile'] = file
                 cc += 1
                 calculations.append(calc)
 

calphy/routines.py

@@ -26,6 +26,7 @@
 import os
 import time
 from mendeleev import element
+import yaml
 
 from calphy.input import read_inputfile
 #import calphy.queuekernel as cq
@@ -56,19 +57,16 @@ def __init__(self, calculation=None, simfolder=None, log_to_screen=False):
         self.calc = calculation
         self.simfolder = simfolder
         self.log_to_screen = log_to_screen
-        self.org_tm = 0
         self.dtemp = self.calc.melting_temperature.step
         self.maxattempts = self.calc.melting_temperature.attempts
         self.attempts = 0
-        self.exp_tm = self.calc._temperature
         self.calculations = []
 
-        self.get_props(self.calc.element[0])
         self.get_trange()
         self.arg = None
 
 
-        logfile = os.path.join(os.getcwd(), "calphy.log")
+        logfile = os.path.join(os.getcwd(), f'{self.calc.create_identifier()}.log')
         self.logger = ph.prepare_log(logfile, screen=log_to_screen)
 
     def prepare_calcs(self):
@@ -83,53 +81,39 @@ def prepare_calcs(self):
         -------
         None 
         """
-        self.calc.temperature = [int(self.tmin), int(self.tmax)]
-        self.calc._temperature = int(self.tmin)
-        self.calc._temperature_stop = int(self.tmax)
-        csol = copy.deepcopy(self.calc)
-        clqd = copy.deepcopy(self.calc)
-        
-        #csol.lattice = self.lattice.upper()
-        #clqd.lattice = 'LQD'
-        csol.reference_phase = 'solid'
-        clqd.reference_phase = 'liquid'
-        #csol.lattice_constant = self.lattice_constant
-        #clqd.lattice_constant = self.lattice_constant
-        csol._temperature_high = self.tmin
-        clqd._temperature_high = 1.5*self.tmax
-        csol.mode = 'ts'
-        clqd.mode = 'ts'
-        
-        #csol['directory'] = create_identifier(csol)
-        #clqd['directory'] = create_identifier(clqd)
-        self.calculations = [csol, clqd]
-        
-        
-    def get_props(self, elem):
-        """
-        Get properties from mendeleev
-
-        Parameters
-        ----------
-        elem : string
-            Chemical symbol of the element
 
-        Returns
-        -------
-        None
-        """
-        chem = element(elem)
-        #lattice = chem.lattice_structure
-        #self.lattice_constant = chem.lattice_constant
-        self.org_tm = chem.melting_point
+        #here, we need to prepare a new calculation
+        #protocol, read in, modify, write a output
+        #read input again
+        calculations = {"calculations": []}
 
-        if self.exp_tm == 0:
-            self.exp_tm = chem.melting_point
-
-        #if lattice == "HEX":
-        #    lattice = "HCP"
-        #self.lattice = lattice.lower()
+        with open(self.calc.inputfile, 'r') as fin:
+            data = yaml.safe_load(fin)
+        calc = data["calculations"][int(self.calc.kernel)]
+
+        calc["mode"] = "ts"
+        calc["temperature"] = [int(self.tmin), int(self.tmax)]
+        calc["reference_phase"] = 'solid'
+        calculations["calculations"].append(calc)
+
+        with open(self.calc.inputfile, 'r') as fin:
+            data = yaml.safe_load(fin)
+        calc = data["calculations"][int(self.calc.kernel)]
+        
+        calc["mode"] = "ts"
+        calc["temperature"] = [int(self.tmin), int(self.tmax)]
+        calc["reference_phase"] = 'liquid'
+        calculations["calculations"].append(calc)
+
+        outfile = f'{self.calc.create_identifier()}.{self.attempts}.yaml'
+        with open(outfile, "w") as fout:
+            yaml.safe_dump(calculations, fout)
+
+        #now read in again, which would allow for checking and so on
+        #one could do this smartly, and simply create from here.
+        self.calculations = read_inputfile(outfile)
 
+                
     def get_trange(self):
         """
         Get temperature range for calculations
@@ -142,10 +126,10 @@ def get_trange(self):
         -------
         None
         """
-        tmin = self.exp_tm - self.dtemp
+        tmin = self.calc._temperature - self.dtemp
         if tmin < 0:
             tmin = 10
-        tmax = self.exp_tm + self.dtemp
+        tmax = self.calc._temperature + self.dtemp
         self.tmax = tmax
         self.tmin = tmin
 
@@ -332,8 +316,9 @@ def calculate_tm(self):
         self.start_calculation()
         tm, tmerr = self.find_tm()
         self.logger.info('Found melting temperature = %.2f +/- %.2f K '%(tm, tmerr))
-        self.logger.info('Experimental melting temperature = %.2f K '%(self.org_tm))
-        self.logger.info('STATE: Tm = %.2f K +/- %.2f K, Exp. Tm = %.2f K'%(tm, tmerr, self.org_tm))
+        if self.calc._melting_temperature is not None:
+            self.logger.info('Experimental melting temperature = %.2f K '%(self.calc._melting_temperature))
+        self.logger.info('STATE: Tm = %.2f K +/- %.2f K'%(tm, tmerr))
 
 def routine_fe(job):
     """

examples/example_04/input.1.yaml

@@ -9,12 +9,9 @@ calculations:
   n_switching_steps: 15000
   pair_coeff: '* * ../potentials/Cu01.eam.alloy Cu'
   pair_style: eam/alloy
+  equilibration_control: berendsen
   queue:
     commands:
     - conda activate calphy
     cores: 4
     scheduler: local
-  repeat:
-  - 5
-  - 5
-  - 5