02 12 25

siman/analysis.py

@@ -4,6 +4,7 @@
 import os, copy, shutil, sys
 import numpy as np
 
+
 try:
     import scipy
     from scipy import interpolate
@@ -1733,7 +1734,14 @@ def ads_en(cl_slab_ads, cl_slab, ads_at = 'O'):
 
 def wulff(st, miller_list = None, e_surf_list = None, show = 0):
 
+    import os
+    import matplotlib.pyplot as plt
+    os.environ["QT_QPA_PLATFORM"] = "xcb"
     from pymatgen.core.structure import Structure
+    from pymatgen.analysis.wulff import WulffShape
+
+
+
     stpm = st.convert2pymatgen()
 
     lat = stpm.lattice
@@ -1746,11 +1754,12 @@ def wulff(st, miller_list = None, e_surf_list = None, show = 0):
     print(dire.get_space_group_info())
     print(recp.get_space_group_info())
     # print(lat)
-    from pymatgen.analysis.wulff import WulffShape
     WS = WulffShape(lat, miller_list, e_surf_list)
     # print(dir(WS))
+    printlog('Area fraction: ',WS.area_fraction_dict, imp = 'y')
     anisotropy = WS.anisotropy
     weighted_surface_energy = WS.weighted_surface_energy
     if show:
-        WS.show()
+        fig = WS.get_plot()
+        plt.show(block=True)
     return anisotropy, weighted_surface_energy

siman/calc_manage.py

@@ -4,6 +4,7 @@
 from operator import itemgetter
 import copy, traceback, datetime, sys, os, glob, shutil, re, io, json
 from itertools import product
+import tempfile
 
 import numpy as np
 
@@ -365,6 +366,7 @@ def add_loop(it = None, setlist = None, verlist = 1, calc = None, varset = None,
             't' - XDATCAR
             'z' - OSZICAR
             'w' - WAVECAR
+            'r' - REPORT
 
         - ifolder - explicit path to folder where to search for input geo file.
 
@@ -455,6 +457,7 @@ def add_loop(it = None, setlist = None, verlist = 1, calc = None, varset = None,
         - params (dic) - dictionary of additional parameters, please move here numerous arguments
 
             - 'occmatrix' - explicit path to occmatrix file
+            - 'iconst_str' - list of strings of ICONST file needed for constrained AIMD
             - 'update_set_dic' (dict) - additional parameters to override the existing set
             - 'monte' - dictionary with parameters for Monte-Carlo regime
 
@@ -1726,6 +1729,15 @@ def write_configuration_file_for_cluster(name, vasp_run_com, params):
                 except:
                     printlog('Attention! no OCCMATRIX file was found!!!')
 
+        if 'iconst_str' in params or ('LBLUEOUT' in cl.set.vasp_params and cl.set.vasp_params['LBLUEOUT'] == 1): #copy iconst file
+            if 'iconst_str' in params:
+
+                with open(cl.dir + '/ICONST', 'w') as f:
+                    f.writelines(params['iconst_str'])
+
+
+            else:
+                printlog('Attention! no ICONST file was found!!!', imp = 'Y')
 
         # if cl.des:
         cl.des = ' '+struct_des[id[0]].des + '; ' + varset[id[1]].des
@@ -3519,7 +3531,7 @@ def get_structure_from_matproj(it = None, it_folder = None, ver = None, mat_proj
             # print groundstate_st_id
             # print m.get_data(groundstate_st_id, data_type='vasp', prop='hubbards')
 
-        st_pmg =  m.get_structure_by_material_id(groundstate_st_id, final=True)
+        st_pmg =  m.get_structure_by_material_id(groundstate_st_id)
 
         if 'conv' in mat_proj_cell:
             from pymatgen.symmetry.analyzer import SpacegroupAnalyzer

siman/calculators/vasp.py

@@ -255,7 +255,7 @@ def calculate_nbands(self, curset, path_to_potcar = None, params = None):
                 elif vp['LNONCOLLINEAR']: 
                     printlog('LNONCOLLINEAR calculation detected; I double number of bands up to', vp['NBANDS'], imp = 'Y')
 
-                if 'LREAL' in vp and 'False' not in vp['LREAL']:
+                if 'LREAL' in vp and 'False' not in str(vp['LREAL']):
                     printlog('Warning! LREAL = .False. is suggested for SOC!', imp = 'Y')
 
 
@@ -277,7 +277,7 @@ def calculate_nbands(self, curset, path_to_potcar = None, params = None):
 
     def make_incar(self, filename = None):
         """Makes Incar file for current calculation and copy all
-        TO DO: there is no need to send all POSCAR files; It is enothg to send only one. However for rsync its not that crucial
+        TO DO: there is no need to send all POSCAR files; It is enough to send only one. However for rsync its not that crucial
         """
         #print "Begin make---------------------------------------------"
 
@@ -327,20 +327,18 @@ def make_incar(self, filename = None):
                         ''
                     elif key == 'MAGMOM' and hasattr(self.init, 'magmom') and self.init.magmom and any(self.init.magmom): #
                         mag = self.init.magmom
-
-                        if len(mag) == natom:
-                            magmom_aligned_with_poscar = [mag[i] for i in poscar_atom_order ]
-
-                        elif len(mag) == 3*natom:
-
-                            mag3 = [mag[i:i+3] for i in range(0, len(mag), 3)]
-
-                            mag3_aligned = [mag3[i] for i in poscar_atom_order ]
 
-                            magmom_aligned_with_poscar = [x for sublist in mag3_aligned for x in sublist]
+                        if isinstance(mag[0], (list, tuple)):  # vector magmom
+                            magmom_aligned_with_poscar = [mag[i] for i in poscar_atom_order]
+                            flat_list = [str(c) for vec in magmom_aligned_with_poscar for c in vec]
+                            f.write('MAGMOM = ' + ' '.join(flat_list) + '\n')
+                        else:  # scalar magmom
+                            magmom_aligned_with_poscar = [mag[i] for i in poscar_atom_order]
+                            print(magmom_aligned_with_poscar)
+                            f.write('MAGMOM = ' + list2string(magmom_aligned_with_poscar) + '\n')
+                        # magmom_aligned_with_poscar = [mag[i] for i in poscar_atom_order ]
+                        # f.write('MAGMOM = '+list2string(magmom_aligned_with_poscar)+"\n") #magmom from geo file has higher preference
 
-
-                        f.write('MAGMOM = '+list2string(magmom_aligned_with_poscar)+"\n") #magmom from geo file has higher preference
 
                     elif vp[key] == None:
                         ''
@@ -717,6 +715,11 @@ def actualize_set(self, curset = None, params = None):
         if hasattr(self.init, 'magmom') and hasattr(self.init.magmom, '__iter__') and not None in self.init.magmom and bool(self.init.magmom):
             magmom = self.init.magmom
 
+        # if params['update_set_dic']['MAGMOM']:
+        #     printlog('actualize_set(): Magnetic moments are determined from params:',params['update_set_dic']['MAGMOM'], imp = 'y')
+
+        # else:
+        if hasattr(self.init, 'magmom') and hasattr(self.init.magmom, '__iter__') and not None in self.init.magmom and bool(self.init.magmom):
             printlog('actualize_set(): Magnetic moments are determined from self.init.magmom:',self.init.magmom, imp = 'y')
 
         elif hasattr(curset, 'magnetic_moments') and curset.magnetic_moments:
@@ -870,7 +873,7 @@ def spin(ls, i):
         if ('LSORBIT' in vp and vp['LSORBIT']) or ('LNONCOLLINEAR' in vp and vp['LNONCOLLINEAR']):
             if len(magmom) == 0:
                 printlog('Warning! magmom is empty')
-            elif len(magmom) == self.init.natom:
+            elif len(magmom) == self.init.natom and len(magmom[0]) == 1:
                 printlog('magmom is converted to [0,0, m, 0, 0, m ...] required for noncollinear calculations', imp = 'y')
                 magmom = [x for m in magmom for x in (0, 0, m)]
             else:
@@ -930,6 +933,8 @@ def copy_to_cluster(self, list_to_copy, update):
         if 'OCCEXT' in self.set.vasp_params and self.set.vasp_params['OCCEXT'] == 1:
             list_to_copy.append(  os.path.join(d, 'OCCMATRIX')   )
 
+        if 'LBLUEOUT' in self.set.vasp_params and self.set.vasp_params['LBLUEOUT'] in [1, '.TRUE.']:
+            list_to_copy.append(  os.path.join(d, 'ICONST')   )
 
         if "up" in update: #Copy to server
             printlog('Files to copy:', list_to_copy)
@@ -1019,6 +1024,51 @@ def plot_energy_step(self,):
         plt.show()
         return
 
+    def plot_energy_total_step(self,):
+        # print(self.maxforce)
+        # maxf = [m[1] for m in self.maxforce_list ]
+        # print(maxf)
+        from siman.picture_functions import plt
+        steps = range(len(self.list_energy_total))
+        plt.plot(steps, 1000*(np.array(self.list_etotal)-self.list_etotal[-1]) , '-o')
+        # plt.xlabel('MD step')
+        # plt.ylabel('Energy per cell (eV')
+        plt.xlabel('Step')
+        plt.ylabel('Total energy rel. to min (meV)')
+        plt.tight_layout()
+        plt.show()
+        return
+
+    def plot_energy_total_step2(self,):
+        # print(self.maxforce)
+        # maxf = [m[1] for m in self.maxforce_list ]
+        # print(maxf)
+        from siman.picture_functions import plt
+        import statsmodels.api as sm  # module to build a LOWESS model
+        import pandas as pd
+        import numpy as np
+
+        lowess = sm.nonparametric.lowess
+
+        steps = range(len(self.list_e_sigma0))
+        e_total = np.array(self.list_e_sigma0) + np.array(self.list_ekin)
+        de = 1000*(np.array(e_total)-e_total[-1])
+
+        temp_lowess = lowess(endog=de, exog=steps)
+        temp_lowess = pd.Series(temp_lowess[:, 1], index=range(len(de)))
+
+
+        plt.plot(steps, de , '-o')
+
+        plt.plot(temp_lowess,    marker="o",    markersize=2,    linestyle="-")
+        # plt.xlabel('MD step')
+        # plt.ylabel('Energy per cell (eV')
+        plt.xlabel('Step')
+        plt.ylabel('Total energy rel. to min (meV)')
+        plt.tight_layout()
+        plt.show()
+        return
+
     def plot_energy_conv(self,):
         # print(self.maxforce)
         # maxf = [m[1] for m in self.maxforce_list ]
@@ -1259,6 +1309,44 @@ def read_results(self, load = '', out_type = '', voronoi = False, show = '', cho
 
 
 
+    def read_oszicar(path_to_file):
+
+        total_mags = []          # Список полной магнетизации на каждый шаг
+        atom_mags_per_step = []  # Список: каждый элемент — массив магнитных моментов атомов на шаге
+
+        with open(filename, 'r') as f:
+            lines = f.readlines()
+
+        current_atom_mags = []
+
+        for line in lines:
+            line = line.strip()
+
+            # Полная магнетизация
+            if "mag=" in line:
+                match = re.search(r"mag=\s*([-0-9.Ee+]+)", line)
+                if match:
+                    total_mags.append(float(match.group(1)))
+
+                if current_atom_mags:
+                    atom_mags_per_step.append(current_atom_mags)
+                    current_atom_mags = []
+
+            # Поатомные магнитные моменты (по колонке M_int)
+            elif re.match(r'^\s*\d+\s', line):  # строка начинается с номера атома
+                parts = line.split()
+                if len(parts) >= 8:
+                    # M_int: три последние компоненты
+                    m_int = list(map(float, parts[5:8]))
+                    current_atom_mags.append(m_int)
+
+        # На случай, если последний шаг не завершился mag=
+        if current_atom_mags:
+            atom_mags_per_step.append(current_atom_mags)
+
+        return total_mags, atom_mags_per_step
+
+
     def determine_filenames(self, nametype = 'asoutcar'):
         """
         try to determine correct filenames
@@ -1478,8 +1566,6 @@ def write_occmatrix(self):
         # print(self.get_path())
         # sys.exit()
 
-
-
         return write_occmatrix(self.occ_matrices, self.get_path())
 
 

siman/chg/chg_func.py

@@ -172,7 +172,7 @@ def cal_chg_diff_files(chg_file1, chg_file2, cl3=None, wcell=0, chg = 'CHGCAR'):
 
 def chg_at_z_direct(cl, k_p = 20, plot = None, filetype = 'CHGCAR'):
     """
-    Return the the value of charge density (filetype = CHGCAR) or electrostatic potential (filetype = LOCPOT) 
+    Return the the value of charge density (filetype = CHGCAR) 
     along z direction of slab; 
 
     chgfile - full path to the file with charge density
@@ -228,6 +228,88 @@ def chg_at_z_direct(cl, k_p = 20, plot = None, filetype = 'CHGCAR'):
     if plot:
         # print(z_coord, elst)
         fit_and_plot(a=(z_coord, elst, '-b'), xlabel = 'Z coordinate, $\AA$', 
-            ylabel = 'Potential, eV', filename = 'figs/'+cl.id[0]+'_pot')
+            ylabel = 'Potential, eV', filename = 'figs/'+cl.id[0]+'_pot', show=1)
 
-    return z_coord, elst
+    return z_coord, elst
+
+
+def locpot_z_profile(cl, axis=2,filename='1.LOCPOT', plot=False, return_numpy=False):
+    """
+    Calculate planar-averaged electrostatic potential from VASP LOCPOT.
+
+    INPUT
+    ----------
+    cl : Calculation object
+
+    axis : int
+        Averaging direction:
+            0 -> x
+            1 -> y
+            2 -> z (default)
+
+    filename : str
+        LOCPOT filename.
+
+    plot : bool
+        Plot potential profile.
+
+    return_numpy : bool
+        Return numpy arrays if True.
+
+    RETURN
+    -------
+    coord : list or np.ndarray
+        Coordinates in Angstrom.
+
+    potential : list or np.ndarray
+        Planar-averaged electrostatic potential in eV.
+
+    AUTHOR:
+
+    A.Boev
+    """
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    from pymatgen.io.vasp.outputs import Locpot
+
+    # Read LOCPOT
+    locpot_path = cl.get_file(filetype=filename)
+
+    locpot = Locpot.from_file(locpot_path)
+
+    # Planar average
+    potential = locpot.get_average_along_axis(axis)
+
+    potential = np.array(potential)
+
+    # Lattice vector length
+    lattice = locpot.structure.lattice
+
+    axis_length = lattice.abc[axis]
+
+    # Coordinate array
+    coord = np.linspace(
+        0,
+        axis_length,
+        len(potential),
+        endpoint=False
+    )
+
+    # Plot
+    if plot:
+
+        plt.figure(figsize=(7,4))
+        plt.plot(coord, potential)
+        xlabel = ['x (Å)', 'y (Å)', 'z (Å)'][axis]
+        plt.xlabel(xlabel)
+        plt.ylabel('Electrostatic potential (eV)')
+        plt.tight_layout()
+        plt.show()
+
+    # Return
+    if return_numpy:
+        return coord, potential
+
+    return coord.tolist(), potential.tolist()

siman/core/calculation.py

@@ -998,7 +998,20 @@ def full(self, ise=None, up=0, fit=1, suf='', add_loop_dic=None, up_res='up1'):
             self.children.append(child)
 
         return child
-
+
+
+    # def get_mag_oszi(self):
+    #     #get magnetic moments from OSZICAR
+
+    #     path_to_file = self.get_file('*OSZICAR')
+
+    #     from siman.calculators.vasp import read_oszicar
+
+    #     total_mags, atom_mags_per_step = read_oszicar(path_to_file)
+    #     print(total_mags, atom_mags_per_step)
+
+
+
     def plot_locpot(self, filename=None):
         """
         Plot local electrostatic potential