Master

dfttk/analysis/debye.py

@@ -202,8 +202,10 @@ def debye_temperature(self, volume):
             # Low temperature limit: 1
             # take 0 K E-V curve properties
             dBdP = self.ev_eos_fit.b1  # bulk modulus/pressure derivative
+            debye = s*A * (self.ev_eos_fit.v0*1.e-30/self.natoms) ** (1. / 6.) * np.sqrt(self.bulk_modulus*1e9/self.avg_mass)
             gamma = (1+dBdP)/2 - self.bp2gru  # 0K equilibrium Gruneisen parameter
-            return debye * (self.ev_eos_fit.v0 / volume) ** (gamma)
+            debye = debye*(self.ev_eos_fit.v0 / volume) ** (gamma)
+            return debye 
         else:
             return debye
 

dfttk/analysis/phonon.py

@@ -42,7 +42,7 @@ def get_f_vib_phonopy(structure, supercell_matrix, vasprun_path,
     # get codename and version from vasprun.xml file
     code_name, code_version = get_code_version(xml=vasprun_path)
     force_constant_factor = 1.0
-    if code_version[0:1] >= '6':
+    if code_version[0:3] >= '6.2':
         force_constant_factor = 0.004091649655126895
 
     # get the force constants from a vasprun.xml file

dfttk/analysis/ywplot.py

@@ -617,7 +617,7 @@ def __init__(self, folder,thermodynamicproperty,x,y,reflin=None, yzero=None,fitt
         elif self.thermodynamicproperty.lower()!="heat capacities (J/mol-atom/K)".lower(): self.plot_default()
         else: self.plot_Heat_Capacity()
 
-        if self.plottitle!=None: plt.title(self.plottitle)
+        #if self.plottitle!=None: plt.title(self.plottitle)
         plt.xlabel(self._xlabel)
         plt.ylabel(self._ylabel)
         self.ax.legend(loc=0, prop={'size': 24})
@@ -642,10 +642,11 @@ def plot_EV(self):
             self.ax.plot(self.x, self.y, fillstyle='none', marker='o', markersize=12,
                 color='k', linestyle='None', label=self._label)
         xnew = np.linspace(min(self.x)*0.95,max(self.x)*1.05, 300)
-        from dfttk.pythelec import BMvol4, BMvol, alt_curve_fit
-        f2, pcov = alt_curve_fit(BMvol4, self.x, self.y)
-        ynew = BMvol(xnew, f2)
-        self.ax.plot(xnew,ynew,'-',linewidth=1,color='b', label="BMvol4")
+        if len(self.x)>=4:
+            from dfttk.pythelec import BMvol4, BMvol, alt_curve_fit
+            f2, pcov = alt_curve_fit(BMvol4, self.x, self.y)
+            ynew = BMvol(xnew, f2)
+            self.ax.plot(xnew,ynew,'-',linewidth=1,color='b', label="BMvol4")
 
 
     def plot_Helmholtz_energy_v0(self):
@@ -1432,7 +1433,7 @@ def mkDict(line):
   idx = 1
   while True:
     if not os.path.exists(dir0): break
-    recordfile = dir0,"record.json"
+    recordfile = os.path.join(dir0,"record.json")
     newdir = False
     try:
       if os.path.exists(recordfile):
@@ -1743,7 +1744,6 @@ def Phonon298(dir0, pvdos=False):
     copyfile(os.path.join(plotdatabase,dfile),os.path.join(phdir298,dfile))
     cwd = os.getcwd()
     os.chdir( phdir298 )
-    import platform
     if platform.system()=="Linux":
       cmd = 'timeout 6 pos2s Symmetry.pos -THR 3.e-4 >&symmetry.out'
       output = subprocess.run(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT,
@@ -2366,19 +2366,29 @@ def Plot298(folder, V298, volumes, debug=False, plottitle=None, local=None):
     ydir = os.path.join(folder,'..')
 
   vdict = {}
+  emin = 1.0e36
   for root, dirs, files in os.walk(ydir):
     for dir in dirs:
       poscar = os.path.join(ydir,dir,'POSCAR')
+
       if os.path.exists(poscar):
         structure = Structure.from_file(poscar)
         vol = 'V{:010.6f}'.format(structure.volume)
         vdict[vol]=dir
-  try:
-    natom = len(structure.sites)
-    sa = SpacegroupAnalyzer(structure)
-    ngroup = sa.get_space_group_number()
-  except:
-    return
+        oszicar = os.path.join(ydir,dir, 'OSZICAR')
+        if not os.path.exists(poscar): continue
+        with open(oszicar,"r") as fp:
+          lines = fp.readlines()
+          for line in lines:
+            dat = [s for s in line.split() if s!=""]
+            if len(dat) < 5: continue
+            if dat[1]!="F=" or dat[3]!="E0=": continue
+            e = float(dat[4])
+            if e < emin:
+              emin = e
+              natom = len(structure.sites)
+              sa = SpacegroupAnalyzer(structure)
+              ngroup = sa.get_space_group_number()
 
   #print(natom,ngroup)
   i1 = 0
@@ -2435,7 +2445,6 @@ def Plot298(folder, V298, volumes, debug=False, plottitle=None, local=None):
       move("vdos.png", os.path.join(cwd,folder,'vdos298.15.png'))
 
   if not os.path.exists('symmetry.mode'):
-    import platform
     if platform.system()=="Linux":
       cmd = "pos2s Symmetry.pos -THR 0.001"
       output = subprocess.run(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT,

dfttk/analysis/ywutils.py

@@ -227,10 +227,19 @@ def get_Magnetic_State(calc):
     magmoms = [k for k in magmoms]
     magmoms = np.array(magmoms)
     fmax = max(magmoms)
+    sdw = magmoms[magmoms>0.1]
+    if len(sdw)>0:
+        smin = min(sdw)
+        smax = max(sdw)
+    else:
+        smin = 0.0
+        smax = 0.0
+    fmax = max(magmoms)
     fmin = min(magmoms)
     fsum = sum(magmoms)
     if fmax>0.1 and fmin<=-0.1:
         if abs(fsum) > 0.1: return "_FIM"
+        elif smax-smin > 0.1: return "_SDW"
         else: return "_AFM"
     elif fmax>0.1 or fmin<=-0.1: return "_FM"
     else: return ""
@@ -262,7 +271,7 @@ def get_rec_from_metatag(vasp_db,m, test=False):
     for calc in static_calculations:
         vol = calc['output']['structure']['lattice']['volume']
         if kpoints is None: kpoints = calc['orig_inputs']['kpoints']['kpoints']
-        if vol_within(vol, volumes): continue
+        if vol_within(vol, volumes, thr=1.e-6): continue
         natoms = len(calc['output']['structure']['sites'])
         try:
             sites = calc['output']['structure']['sites']
@@ -288,16 +297,16 @@ def get_rec_from_metatag(vasp_db,m, test=False):
         else: pressures.append(None)
         if not gapfound: gapfound = float(gap) > 0.0
     tvolumes = np.array(sorted(volumes))
-    if len(tvolumes)>1:
-        dvolumes = tvolumes[1:-1] - tvolumes[0:-2]
+    if len(tvolumes)>=3:
+        dvolumes = tvolumes[1:] - tvolumes[0:-1]
         dvolumes = sorted(dvolumes)
         if abs(dvolumes[-1]-dvolumes[-2]) > 0.01*dvolumes[-1]:
             all_static_calculations = vasp_db.collection.\
                 find({'$and':[ {'metadata.tag': m}, {'adopted': True} ]})
             for calc in all_static_calculations:
                 if len(calc['metadata'])<=1:continue # only check constrained calculation
                 vol = calc['output']['structure']['lattice']['volume']
-                if vol_within(vol, volumes): continue
+                if vol_within(vol, volumes, thr=1.e-6): continue
                 natoms = len(calc['output']['structure']['sites'])
                 try:
                     sites = calc['output']['structure']['sites']

dfttk/elasticity/parse_outputs.py

@@ -305,9 +305,14 @@ def run_task(self, fw_spec):
         }
         vasp_input_set = env_chk(self.get('vasp_input_set'), False)
         if vasp_input_set:
+            from pymatgen.io.vasp.inputs import Kpoints
             d['vasp_input_set_all'] = vasp_input_set.config
-            d['vasp_input_set'] = {'INCAR':vasp_input_set.config['INCAR'], \
-            'KPOINTS':vasp_input_set.config['KPOINTS'].as_dict()}
+            if isinstance(vasp_input_set.config['KPOINTS'], Kpoints):
+                d['vasp_input_set'] = {'INCAR':vasp_input_set.config['INCAR'], \
+                'KPOINTS':vasp_input_set.config['KPOINTS'].as_dict()}
+            else:
+                d['vasp_input_set'] = {'INCAR':vasp_input_set.config['INCAR'], \
+                'KPOINTS':vasp_input_set.config['KPOINTS']}
 
         # Get optimized structure
         calc_locs_opt = [cl for cl in fw_spec.get('calc_locs', []) if 'optimiz' in cl['name']]

dfttk/ftasks.py

@@ -32,8 +32,8 @@
 from atomate import __version__ as atomate_ver
 from dfttk import __version__ as dfttk_ver
 from pymatgen.core import __version__ as pymatgen_ver
-from dfttk.pythelec import get_static_calculations
-from dfttk.scripts.assign_fworker_name import Customizing_Workflows, get_powerups_options
+from dfttk.pythelec import get_static_calculations, vol_within
+from dfttk.scripts.assign_fworker_name import Customizing_Workflows
 
 def extend_calc_locs(name, fw_spec):
     """
@@ -340,9 +340,13 @@ def run_task(self, fw_spec):
         #always perform phonon calculations when when enough phonon calculations found
         #to perform a quasiharmonic phonon calculations, one needs at least phonon results five volumes
         #phonon_calculations= list(vasp_db.db['phonon'].find({'$and':[ {'metadata.tag': tag}, {'adopted': True} ]}))     
-        phonon_calculations= list(vasp_db.db['phonon'].find({'$and':[ {'metadata': {'tag':tag}}, {'adopted': True} ]}))     
+        phonon_calculations= list(vasp_db.db['phonon'].find({'$and':[ {'metadata': {'tag':tag}}, {'adopted': True}, {'S_vib': {'$exists': True}} ]}))     
         num_phonon_finished = len(phonon_calculations)       
         qha_result['has_phonon'] = num_phonon_finished >= 5
+        if not qha_result['has_phonon']:
+            phonon_calculations= list(vasp_db.db['phonon'].find({'$and':[ {'metadata': {'tag':tag}}, {'S_vib': {'$exists': True}} ]}))     
+            num_phonon_finished = len(phonon_calculations)       
+            qha_result['has_phonon'] = num_phonon_finished >= 5
         #if self['phonon']:
         if qha_result['has_phonon']:
             # get the vibrational properties from the FW spec
@@ -355,12 +359,15 @@ def run_task(self, fw_spec):
             vol_s_vib = []
             vol_c_vib = []
             for calc in phonon_calculations:
-                if calc['volume'] in vol_vol: continue
-                if calc['volume'] not in volumes: continue
-                vol_vol.append(calc['volume'])
+                #if calc['volume'] in vol_vol: continue
+                #if calc['volume'] not in volumes: continue
+                if vol_within(calc['volume'], vol_vol, thr=1.e-6): continue
+                if not vol_within(calc['volume'],volumes, thr=1.e-6): continue
                 vol_f_vib.append(calc['F_vib'][::everyT])
                 vol_s_vib.append(calc['S_vib'])
                 vol_c_vib.append(calc['CV_vib'])
+                vol_vol.append(calc['volume'])
+
             # sort them order of the unit cell volumes
             vol_f_vib = sort_x_by_y(vol_f_vib, vol_vol)
             vol_s_vib = sort_x_by_y(vol_s_vib, vol_vol)
@@ -371,7 +378,8 @@ def run_task(self, fw_spec):
             _energies = []
             _dos_objs = []
             for iv,vol in enumerate(volumes):
-                if vol not in vol_vol:  continue
+                #if vol not in vol_vol:  continue
+                if not vol_within(vol,vol_vol, thr=1.e-6): continue
                 _volumes.append(vol)
                 _energies.append(energies[iv])
                 _dos_objs.append(dos_objs[iv])
@@ -740,6 +748,88 @@ def run_task(self, fw_spec):
         kpoint.write_file('KPOINTS')
 
 
+@explicit_serialize
+class Crosscom_Calculation(FiretaskBase):
+    """Continue Static/Phonon calculations
+    """
+
+    required_params = []
+    optional_params = ['db_file', 'vasp_cmd', 'a_kwargs',
+                       'db_insert', 'tag', 'metadata', 'name', 'vasp_input_set',
+                       't_min', 't_max', 't_step', 
+                       'verbose', 'modify_incar_params', 'modify_kpoints_params', 
+                       'override_default_vasp_params', 
+                       'store_volumetric_data', 'static', 'defo']
+
+    def run_task(self, fw_spec):
+        db_file = self.get('db_file') or DB_FILE
+        vasp_cmd = self.get('vasp_cmd') or VASP_CMD
+        db_insert = self.get('db_insert', None)
+        tag = self.get('tag')
+        metadata = self.get('metadata')
+        name = self.get('name', "Crosscom_Calculation")
+        t_min = self.get('t_min', None)
+        t_max = self.get('t_max', None)
+        t_step = self.get('t_step', None)
+        modify_incar_params = self.get('modify_incar_params', {})
+        modify_kpoints_params = self.get('modify_kpoints_params', {})
+        override_default_vasp_params = self.get('override_default_vasp_params', {})
+        store_volumetric_data = self.get('store_volumetric_data', False)
+
+        return FWAction(detours=self.get_detour_workflow(
+            db_file, vasp_cmd, db_insert, tag, metadata, name, 
+            t_min, t_max, t_step, 
+            modify_incar_params, modify_kpoints_params,
+            override_default_vasp_params,
+            store_volumetric_data)
+            )
+
+    def get_detour_workflow(self,
+        db_file, vasp_cmd, db_insert, tag, metadata, name, 
+        t_min, t_max, t_step, 
+        modify_incar_params, modify_kpoints_params, 
+        override_default_vasp_params, 
+        store_volumetric_data):
+        from fireworks import Workflow
+        from .fworks import PhononFW, StaticFW
+        a_kwargs = self.get('a_kwargs',{})
+        settings = a_kwargs.get('settings', {})
+        stable_tor = settings.get('stable_tor', 0.01)
+
+        phonon = settings.get('phonon', False)
+        phonon_supercell_matrix = a_kwargs.get('phonon_supercell_matrix', None)
+        structure=a_kwargs.get('structure', None)
+        site_properties = structure.site_properties
+
+        detour_fws = []
+        inp_structure = Structure.from_file('CONTCAR')
+        if len(site_properties)>0:
+            for prop, vals in site_properties.items():
+                inp_structure.add_site_property(prop, vals)
+
+        detour_fws.append(StaticFW(inp_structure, name="crosscom-static"+'-defo={:5.3f}'.format(self.get('defo',1.0)), 
+                 vasp_cmd=vasp_cmd, metadata=metadata, prev_calc_loc=False, modify_incar=modify_incar_params, 
+                 db_file=db_file, tag=tag, 
+                 a_kwargs=a_kwargs,
+                 override_default_vasp_params=override_default_vasp_params,
+                 store_volumetric_data=store_volumetric_data))
+
+        if phonon:
+            t_kwargs = {'t_min': t_min, 't_max': t_max, 't_step': t_step}
+            common_kwargs = {'vasp_cmd': vasp_cmd, 'db_file': db_file, "metadata": metadata, "tag": tag,
+                'override_default_vasp_params': override_default_vasp_params}
+            detour_fws.append(PhononFW(inp_structure, phonon_supercell_matrix, 
+                name='crosscom-phonon'+'-defo={:5.3f}'.format(self.get('defo',1.0)), 
+                prev_calc_loc=False, stable_tor=stable_tor,
+                a_kwargs=a_kwargs,
+                **t_kwargs, **common_kwargs))
+        """
+        override_default_vasp_params = override_default_vasp_params or {}
+        user_incar_settings = override_default_vasp_params.get('user_incar_settings',{})
+        """
+        return Customizing_Workflows(detour_fws, powerups_options=settings.get('powerups', None))
+
+
 @explicit_serialize
 class CheckRelaxation(FiretaskBase):
     """Run VASP calculations to get symmetry conserved and symmetry broken structures.
@@ -806,7 +896,8 @@ class CheckRelaxation(FiretaskBase):
 
     required_params = ["db_file", "tag", "common_kwargs"]
     optional_params = ["metadata", "tol_energy", "tol_strain", "tol_bond", 'level', 'isif4',  "energy_with_isif",
-                       "static_kwargs", "relax_kwargs", 'store_volumetric_data', 'site_properties']
+                       "static_kwargs", "relax_kwargs", 'store_volumetric_data', 
+                       'a_kwargs', 'site_properties']
 
     def run_task(self, fw_spec):
         self.db_file = env_chk(self.get("db_file"), fw_spec)
@@ -930,6 +1021,8 @@ def get_detour_workflow(self, next_steps, final_energy, energy_with_isif={}):
 
         symmetry_options = self.symmetry_options
         static_kwargs = self.get('static_kwargs', {})
+        a_kwargs = self.get('a_kwargs', {})
+        settings = a_kwargs.get('settings', {})
 
         # Assume the data for the current step is already in the database
         db = VaspCalcDb.from_db_file(self.db_file, admin=True).db['relaxations']
@@ -949,17 +1042,21 @@ def _get_input_structure_for_step(step):
                 md['symmetry_type'] = step["symmetry_type"]
                 common_copy["metadata"] = md
                 detour_fws.append(StaticFW(inp_structure, isif=step['isif'], store_volumetric_data=self.store_volumetric_data,
-                                           **static_kwargs, **common_copy))
+                        a_kwargs=a_kwargs,
+                        **static_kwargs, **common_copy))
             elif job_type == "relax":
                 detour_fws.append(RobustOptimizeFW(inp_structure, isif=step["isif"], energy_with_isif=energy_with_isif,
-                                override_symmetry_tolerances=symmetry_options, store_volumetric_data=self.store_volumetric_data, **self["common_kwargs"]))
+                                override_symmetry_tolerances=symmetry_options, store_volumetric_data=self.store_volumetric_data, 
+                                a_kwargs=a_kwargs, **self["common_kwargs"]))
             else:
                 raise ValueError(f"Unknown job_type {job_type} for step {step}.")
+        """
         common_kwargs = self.get('common_kwargs',{})
         override_default_vasp_params = common_kwargs.get('override_default_vasp_params',{})
         user_incar_settings = override_default_vasp_params.get('user_incar_settings',{})
         #user_incar_settings = common_kwargs.get('modify_incar_params',{})
-        return Customizing_Workflows(detour_fws, powerups_options=user_incar_settings.get('powerups', None))
+        """
+        return Customizing_Workflows(detour_fws, powerups_options=settings.get('powerups', None))
 
 
 @explicit_serialize