From a72c060bdcf2914155066b3abf244b9154837242 Mon Sep 17 00:00:00 2001 From: antonboev22 Date: Wed, 15 Jan 2025 13:04:46 +0300 Subject: [PATCH 1/6] 15.01.25 --- siman/analysis.py | 1 + siman/calculators/vasp.py | 44 ++++++++++++++ siman/geo.py | 2 +- siman/inout.py | 120 +++++++++++++++++++++++++++++-------- siman/matproj_functions.py | 3 +- 5 files changed, 141 insertions(+), 29 deletions(-) diff --git a/siman/analysis.py b/siman/analysis.py index 8cbdab7b..3511e130 100644 --- a/siman/analysis.py +++ b/siman/analysis.py @@ -1749,6 +1749,7 @@ def wulff(st, miller_list = None, e_surf_list = None, show = 0): from pymatgen.analysis.wulff import WulffShape WS = WulffShape(lat, miller_list, e_surf_list) # print(dir(WS)) + print(WS.area_fraction_dict) anisotropy = WS.anisotropy weighted_surface_energy = WS.weighted_surface_energy if show: diff --git a/siman/calculators/vasp.py b/siman/calculators/vasp.py index fd07200d..9a29edef 100644 --- a/siman/calculators/vasp.py +++ b/siman/calculators/vasp.py @@ -938,6 +938,50 @@ 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 + + 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 ] diff --git a/siman/geo.py b/siman/geo.py index 7ac91a06..60ac27ed 100644 --- a/siman/geo.py +++ b/siman/geo.py @@ -2826,7 +2826,7 @@ def create_ads_molecule(st, molecule = ['O'], mol_xc = [[0,0,0]], conf_i = [0], adsorbate = Molecule(molecule, mol_xc) - ads_structs = asf_pm.generate_adsorption_structures(adsorbate,repeat=None, min_lw=8.0, find_args=find_args) + ads_structs = asf_pm.generate_adsorption_structures(adsorbate,repeat=None, min_lw=5.0, find_args=find_args) print('\nI\'ve found ',len(ads_structs), ' configurations for ', molecule, ' on the surface\n') st_ads_pack = [] diff --git a/siman/inout.py b/siman/inout.py index 9277af4c..52ad4e24 100644 --- a/siman/inout.py +++ b/siman/inout.py @@ -164,7 +164,7 @@ def cif2poscar(cif_file, poscar_file): #check if not os.path.exists(poscar_file): - printlog("Error! cif2cell failed") + printlog("Error! cif2cell failed. Check that cif2cell is installed and available in your system") else: printlog('Error! Support of cif files requires pymatgen or cif2cell; install it with "pip install pymatgen" or provide POSCAR or Abinit input file') @@ -194,15 +194,17 @@ def smart_structure_read(filename = None, curver = 1, calcul = None, input_folde if filename: input_geo_file = filename + # print(input_geo_file) + # sys.exit() if input_geo_file: - printlog("You provided the following geo file explicitly ", input_geo_file, '; Version of file does not matter, I use *curver*=',curver, 'as a new version' ) + # sys.exit() elif input_folder: - printlog("I am searching for geofiles in folder "+input_folder+"\n" ) + printlog("I am searching for geofiles in folder "+input_folder+" with format", input_geo_format, "\n" ) if input_geo_format: geofilelist = glob.glob(input_folder+'/'+search_templates[input_geo_format]) #Find input_geofile of specific format @@ -296,12 +298,17 @@ def read_xyz(st, filename, rprimd = None): INPUT: st (Structure) filename (str) - path to file with xyz data - rprimd (list of lists) - if None or [None, ] then Tv are read; if Tv does not exist then create automatically + rprimd (list of lists) - if None or [None, ] then Tv or VEC are read; if Tv does not exist then create automatically """ with open(filename,'r') as f: nlines = int(f.readline()) - st.name = f.readline().strip().split()[0] + name = f.readline().strip().replace(' ', '_') + if name: + st.name = name + else: + st.name = 'noname' + # try: if 'SG' in st.name: @@ -312,13 +319,14 @@ def read_xyz(st, filename, rprimd = None): elements = [] st.xcart = [] st.rprimd = [] - for i in range(nlines): + for i in range(nlines+3): xc = f.readline().split() if len(xc) == 0: - printlog('Warning! xyz file is broken, not enough lines') + if i < nlines: + printlog('Warning! xyz file is broken') break - - if 'Tv' in xc[0]: + # print(xc[0]) + if 'Tv' in xc[0] or 'VEC' in xc[0]: st.rprimd.append(np.asarray(xc[1:], dtype = float) ) else: @@ -380,7 +388,8 @@ def read_xyz(st, filename, rprimd = None): printlog('Final rprimd = \n', np.round(st.rprimd, 3), imp = 'y') - + for i in 0,1,2: + st.rprimd[i] = np.array(st.rprimd[i]) st.nznucl = st.get_nznucl() @@ -1325,6 +1334,8 @@ def read_vasp_out(cl, load = '', out_type = '', show = '', voronoi = '', path_to maxforce = []; average = []; gstress =[] # mforce = [] self.list_e_sigma0 = [] + self.list_etotal = [] # list of MD energies + self.list_ekin = [] # list of kinetic energies self.list_e_without_entr = [] self.list_e_conv = [] # convergence of energy - all steps # try: @@ -1629,6 +1640,13 @@ def read_vasp_out(cl, load = '', out_type = '', show = '', voronoi = '', path_to #pass + if 0 and "g(Stress)" in line: + # print( line) + if 'ortho' in line: + gstr = float(line.split('=')[2].split('ortho')[0]) + else: + gstr = float(line.split('=')[2]) + gstress.append( round( gstr*1000 *100, 3 ) ) if "g(Stress)" in line: #print line gstress.append( round( float(line.split('=')[2].split()[0])*1000 *100, 3 ) ) @@ -1664,8 +1682,23 @@ def read_vasp_out(cl, load = '', out_type = '', show = '', voronoi = '', path_to if "-V(xc)+E(xc) XCENC" in line: self.energy.xc = float(line.split()[-1]) # Kohn-Sham exchange-correlation energy if "PAW double counting" in line: - self.energy.pawdc1 = float(line.split()[-2]) # - self.energy.pawdc2 = float(line.split()[-1]) # + # print('inout:', line) + if '*' in line: + line = 'p=0 -0' + numbers = line.split('=')[1] + e = numbers.split() + if len(e) == 1: + e = numbers.split('-') + e[1] = '-'+e[1] + if len(e) == 1: + e = ['0', '0'] + + # print(e) + self.energy.pawdc1 = float(e[0].strip()) # + # print(self.energy.pawdc1 ) + + self.energy.pawdc2 = float(e[1].strip()) # + if "eigenvalues EBANDS" in line: try: self.energy.bands = float(line.split()[-1]) # - Kohn Sham eigenvalues - include kinetic energy , but not exactly @@ -1698,8 +1731,15 @@ def read_vasp_out(cl, load = '', out_type = '', show = '', voronoi = '', path_to if "free energy TOTEN =" in line: #self.energy = float(line.split()[4]) - self.energy_free = float(line.split()[4]) #F free energy + self.energy_free = float(line.split()[4]) # F free energy + if "energy ETOTAL =" in line: + self.etotal = float(line.split()[4]) # Total energy in MD simulation + self.list_etotal.append(self.etotal) + + if "kinetic energy EKIN =" in line: + self.ekin = float(line.split()[4]) # Kinetic energy in MD simulation + self.list_ekin.append(self.ekin) @@ -1808,14 +1848,14 @@ def read_vasp_out(cl, load = '', out_type = '', show = '', voronoi = '', path_to if 'magnetization (x)' in line: # print(line) - if ifmaglist is not None: - mags = [] - for j in range(self.end.natom): - mags.append( float(outcarlines[i_line+j+4].split()[-1]) ) - - tot_mag_by_atoms.append(np.array(mags))#[ifmaglist]) - # print(ifmaglist) - tot_mag_by_mag_atoms.append(np.array(mags)[ifmaglist]) + # if ifmaglist is not None: + mags = [] + for j in range(self.end.natom): + mags.append( float(outcarlines[i_line+j+4].split()[-1]) ) + + tot_mag_by_atoms.append(np.array(mags))#[ifmaglist]) + # print(ifmaglist) + tot_mag_by_mag_atoms.append(np.array(mags)[ifmaglist]) # print tot_mag_by_atoms # magnetic_elements # ifmaglist @@ -2466,7 +2506,7 @@ def butter_lowpass_filter(data, cutoff, fs, order=5): # if not hasattr(cl, 'name'): # cl.name = 'noname' - outst_simple = '|'.join([cl.name, etot, lens, strs, Nmd]) + outst_simple = '|'.join([etot, lens, strs, Nmd]) # print("Bi2Se3.static.1 | -20.1543 | 10.27;10.27;10.27 | -680,-680,-657 | 1,13, 13 | ") if header.show_head : printlog(" | energy(eV)| Vector lenghts (A) | Stresses (MPa) | N MD, N SCF ", end = '\n', imp = 'Y') @@ -2512,15 +2552,27 @@ def read_aims_out(cl, load = '', out_type = '', show = ''): return outstr -def read_atat_fit_out(filename, filter_names = None, i_energy = 1): +def read_atat_fit_out(filename, filter_names = None, i_energy = 1, xmax = 1, e0e1exmax = None, subs = None): """ read fit.out of atat - filter_names - do not read name numbers greater than filter_name - i_energy - 1 for fit.out, 2 for predstr.out - return - concentration list, energy list + INPUT: + - filter_names - do not read name numbers greater than filter_name + - i_energy - 1 for fit.out, 2 for predstr.out + - return - concentration list, energy list + - xmax (float) - maximum value of concentration + - e0e1exmax (list*3) - e0, e1, and e_xmax, needed to recalculate if xmax is not equal to 1 """ X, E, nam = [] ,[], [] count = 0 + if subs is None: + subs = {} + + if e0e1exmax is None: + e0e1exmax = [0,0,0] + e0 = e0e1exmax[0] + e1 = e0e1exmax[1] + e_xmax = e0e1exmax[2] + with open(filename, 'r') as f: for line in f: @@ -2540,6 +2592,21 @@ def read_atat_fit_out(filename, filter_names = None, i_energy = 1): e = round(float(vals[i_energy]),4) x = float(vals[0]) + if x > xmax: + continue + print(x) + if x in subs: + print(x, 'in subs') + if subs[x] == 'skip': + continue + if x in subs: + print(x, e) + e = subs[x] + print('substituting with', e) + + if xmax != 1: + # print (e_xmax/xmax - e1 + (xmax-1)/xmax * e0) + e = e - x * ( e_xmax/xmax - e1 + (xmax-1)/xmax * e0 ) # print(e) # if e in E: # if X[E.index(e)] == x: @@ -2551,6 +2618,7 @@ def read_atat_fit_out(filename, filter_names = None, i_energy = 1): # if count > 10: # count = 0 X.append(x) # concentrations + E.append(e) # formation energies nam.append(name) # names of structures count+=1 diff --git a/siman/matproj_functions.py b/siman/matproj_functions.py index 3c9a9fc4..0cfee210 100644 --- a/siman/matproj_functions.py +++ b/siman/matproj_functions.py @@ -421,8 +421,7 @@ def calc_bulk_list2(data_list, suf = '', ise = '8m', status = None, up =None, co """ - if suf != '' - : suf = f'.{suf}' + if suf != '' : suf = f'.{suf}' for i in data_list: From 0b5531c302e7aeedff86360df7210958c97d1313 Mon Sep 17 00:00:00 2001 From: antonboev22 Date: Mon, 1 Sep 2025 10:25:50 +0300 Subject: [PATCH 2/6] 01.09.2025 --- siman/calculators/vasp.py | 67 ++++++++++++++++--- siman/core/calculation.py | 11 ++++ siman/core/structure.py | 13 +++- siman/geo.py | 131 ++++++++++++++++++++++---------------- siman/set_functions.py | 12 +++- 5 files changed, 169 insertions(+), 65 deletions(-) diff --git a/siman/calculators/vasp.py b/siman/calculators/vasp.py index 9a29edef..b07eb0b3 100644 --- a/siman/calculators/vasp.py +++ b/siman/calculators/vasp.py @@ -259,7 +259,7 @@ def calculate_nbands(self, curset, path_to_potcar = None, params = None): def make_incar(self): """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---------------------------------------------" @@ -304,8 +304,16 @@ def make_incar(self): '' elif key == 'MAGMOM' and hasattr(self.init, 'magmom') and self.init.magmom and any(self.init.magmom): # mag = self.init.magmom - 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 + + 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] + 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 elif vp[key] == None: '' @@ -517,6 +525,7 @@ def actualize_set(self, curset = None, params = None): curset = self.set vp = curset.vasp_params + print(params) if params is None: params = {} @@ -659,8 +668,12 @@ def actualize_set(self, curset = None, params = None): # print(self.init.magmom) # print(None in self.init.magmom) # sys.exit() - if hasattr(self.init, 'magmom') and hasattr(self.init.magmom, '__iter__') and not None in self.init.magmom and bool(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: @@ -720,7 +733,7 @@ def actualize_set(self, curset = None, params = None): if len(spec_mom_is) % 2 == 0 and len(spec_mom_is) > 0: printlog('Number of elements is even! trying to find all antiferromagnetic orderings:', imp = 'y') ns = len(spec_mom_is); - number_of_ord = int(math.factorial(ns) / math.factorial(0.5 * ns)**2) + number_of_ord = int(math.factorial(ns) / math.factorial(int(0.5 * ns))**2) if number_of_ord > 10000: printlog('Attention! Too much orderings (1000), skipping ...') @@ -805,10 +818,10 @@ def spin(ls, i): elif 'MAGMOM' in vp and vp['MAGMOM']: #just add * to magmom tag if it is provided without it printlog('Magnetic moments from vasp_params["MAGMOM"] are used\n') - - # if "*" not in vp['MAGMOM']: - # vp['MAGMOM'] = str(natom) +"*"+ vp['MAGMOM'] + # if "*" not in vp['MAGMOM']: + # vp['MAGMOM'] = str(natom) +"*"+ vp['MAGMOM'] + # print (self.init.magmom, 'asdfaaaaaaaaaaaa') @@ -1218,6 +1231,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 diff --git a/siman/core/calculation.py b/siman/core/calculation.py index e6404f0e..a2015240 100644 --- a/siman/core/calculation.py +++ b/siman/core/calculation.py @@ -559,6 +559,8 @@ def gmt(self, *args, **kwargs): return self.end.get_mag_tran(*args, **kwargs) + + def isggau(self): #check if calculation is gga+u #TODO - please finish @@ -970,6 +972,15 @@ def full(self, ise = None, up = 0, fit = 1, suf = '', add_loop_dic = None, up_r + 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) diff --git a/siman/core/structure.py b/siman/core/structure.py index 10fb6b71..f8438a4e 100644 --- a/siman/core/structure.py +++ b/siman/core/structure.py @@ -292,6 +292,10 @@ def exchange_axes_with_atoms(self, i1_r, i2_r): xr = copy.deepcopy(self.xred) r = copy.deepcopy(self.rprimd) + st.rprimd[i1_r] = r[i2_r] + st.rprimd[i2_r] = r[i1_r] + + for i in range(0,st.natom): st.xcart[i][i1_r] = xc[i][i2_r] st.xcart[i][i2_r] = xc[i][i1_r] @@ -550,7 +554,10 @@ def get_mag(self, silent = 0): return mag_grouped - + + + + def set_magnetic_config(self, element, moments): #set magnetic configuration based on symmetry non-equivalent positions @@ -947,6 +954,8 @@ def rprimd_len(self): + + def pvec(self): """ print primitive vectors in formatted way @@ -1873,6 +1882,8 @@ def del_atom(self, iat): '' # print ('No!') + if hasattr(st, 'select') and any(st.select): + del st.select[i] st.natom-=1 diff --git a/siman/geo.py b/siman/geo.py index 60ac27ed..6099b033 100644 --- a/siman/geo.py +++ b/siman/geo.py @@ -3129,43 +3129,71 @@ def triangle_area_points(v1, v2, v3): return a -def sl_misfit(st1,st2, silent = 0): +def slab_misfit(st1,st2, silent = 0): + ''' + Calculates difference between two slabs x,y sizes (Misfit) + ''' size1 = st1.rprimd_len() size2 = st2.rprimd_len() - misfit = [(j-i)*100/j for i,j in zip(size1,size2)] + misfit1 = [(j-i)*100/j for i,j in zip(size1,size2)] + # misfit2 = [(j-i)*100/j for i,j in zip(size1,[size2[1],size2[0],size2[2]])] # print('\n\nSize 1: {},{},{} A'.format(round(size1[0],2),round(size1[1],2),round(size1[2],2))) # print('Size 2: {},{},{} A'.format(round(size2[0],2),round(size2[1],2),round(size2[2],2))) if silent == 0: - print('Misfit: {},{} % \n\n'.format(round(misfit[0],2),round(misfit[1],2))) - return misfit + print('Misfit: {},{} % \n\n'.format(round(misfit1[0],2),round(misfit1[1],2))) + # print('Misfit: {},{} % \n\n'.format(round(misfit2[0],2),round(misfit2[1],2))) + return misfit1 -def fit2host(st_host, st_oxide): +def fit2host(st_host, st_film): replic = [1,1,1] - misf = sl_misfit(st_host,st_oxide, silent = 1) - for m in (0,1): - if 60 < abs(misf[m]) < 150: - replic[m] +=1 - elif 150 < abs(misf[m]) < 250: - replic[m] +=2 - elif 250 < abs(misf[m]) < 350: - replic[m] +=3 - elif 350 < abs(misf[m]) < 450: - replic[m] +=3 - st_oxide = st_oxide.replic(replic) + misf1 = slab_misfit(st_host,st_film, silent = 1) + st_film_rot = st_film.copy().exchange_axes_with_atoms(0,1) + # st_film_rot.update_xred() + misf2 = slab_misfit(st_host,st_film_rot, silent = 1) + st_films = [] + print(st_film.rprimd_len()) + print(st_film_rot.rprimd_len()) + + + for misf in [misf1,misf2]: + + for m in (0,1): + if 60 < abs(misf[m]) < 150: + replic[m] +=1 + elif 150 < abs(misf[m]) < 250: + replic[m] +=2 + elif 250 < abs(misf[m]) < 350: + replic[m] +=3 + elif 350 < abs(misf[m]) < 450: + replic[m] +=3 + st_films.append(st_film.copy().replic(replic)) + + misf1 = slab_misfit(st_host,st_films[0], silent = 1) + misf2 = slab_misfit(st_host,st_films[1], silent = 1) + + print(f'Misf1: {misf1[0]:f} {misf1[1]:f}\nMisf2: {misf2[0]:f} {misf2[1]:f}') + if sum(abs(i) for i in misf1) < sum(abs(i) for i in misf2): + st_film_new = st_films[0] + print(f'Final misfit: {misf1[0]:.1f},{misf1[1]:.1f} % \n') + else: + st_film_new = st_films[1] + print(f'Final misfit: {misf2[0]:.1f},{misf2[1]:.1f} % \n') - return st_oxide + return st_film_new -def hkl_slab(st, st_host, hkl, i_suf = None): - os.remove('/home/anton/media/vasp/log_best1') - f = open('/home/anton/media/vasp/log_best1', 'a') +def hkl_slab(st, st_host, hkl, i_suf = None, path = 'filename'): + + + os.remove(path) + f = open(path, 'a') if slabs2: if not i_suf: for sl_i in range(0,len(slabs2)): # print(hkl) st2_new = st.update_from_pymatgen(slabs2[sl_i]) - misf = sl_misfit(st_host,st2_new, silent = 0) + misf = slab_misfit(st_host,st2_new, silent = 0) replic = [1,1,1] for m in (0,1): @@ -3173,7 +3201,7 @@ def hkl_slab(st, st_host, hkl, i_suf = None): replic[m] +=1 st2_new = st2_new.replic(replic) # print(replic) - misf = sl_misfit(st_host,st2_new, silent = 1) + misf = slab_misfit(st_host,st2_new, silent = 1) print(hkl, sl_i) string = str(hkl) + ' ' + str(sl_i) + ' Misfit: {},{} % \n\n'.format(round(misf[0],2),round(misf[1],2)) if abs(misf[0]) < 20 and abs(misf[1])<20: @@ -3183,25 +3211,22 @@ def hkl_slab(st, st_host, hkl, i_suf = None): return misf, slabs2 -def create_interface_solid(st_host, st_oxide, suf_host, i_suf_host = 0, +def create_interface_solid(st_host, st_film, hkl_host, i_suf_host = 0, seek_mode = 0, seek_range = [0,2], check_shift = None, - hkl_lio = None, i_suf_lio = None, size = [5,5], ads_pos = None, z_shift = 1.5, lio_thick = 8): + hkl_film = None, i_suf_film = None, size = [5,5], ads_pos = None, z_shift = 1.5, film_thickness = 8): st1_init = st_host.copy() - st2_init = st_oxide.copy() + st2_init = st_film.copy() - if suf_host: + if hkl_host: sc1 = st1_init.get_conventional_cell()#.replic([2,2,1]) - slabs1 = create_surface2(sc1, suf_host, shift = None, min_slab_size = 10, min_vacuum_size = 15, - surface_i = 0, oxidation = None, suf = '', - symmetrize = 1, cut_thickness = 0, return_one = 0, lll_reduce = 1, primitive = 1) + slabs1 = create_surface2(sc1, hkl_host, surface_i = 0, min_slab_size =10, min_vacuum_size = 25, symmetrize = True, lll_reduce = 1, primitive = 1) slab1 = sc1.update_from_pymatgen(slabs1[i_suf_host]) - # slab1 = slabs1 mul_matrix = ortho_vec(slab1.rprimd, [5,5,15], silent = 1) # matrix which allows to obtain supercell close to 10x10x10 A cube - st1 = create_supercell(slab1, mul_matrix, silent = 1) + slab1 = create_supercell(slab1, mul_matrix, silent = 1) else: - st1 = st_host + slab1 = st_host @@ -3212,33 +3237,28 @@ def create_interface_solid(st_host, st_oxide, suf_host, i_suf_host = 0, hkl = [h,k,l] if hkl != [0,0,0]: - slabs2 = create_surface2(st2_init, hkl, shift = None, min_slab_size = 10, min_vacuum_size = 25, surface_i = 0, oxidation = None, suf = '', - primitive = 1, symmetrize = 0, cut_thickness = None, return_one = 0, write_poscar = 0, lll_reduce = 1) + slabs2 = create_surface2(st2_init, hkl, surface_i = 0, min_slab_size =10, min_vacuum_size = 25, symmetrize = True, lll_reduce = 1, primitive = 1) for sl in range(0,len(slabs2)): print(hkl, sl) - st2 = st2_init.update_from_pymatgen(slabs2[sl]) - slab_2 = fit2host(st1, st2) - sl_misfit(st1, slab_2, silent = 0) + slab2 = st2_init.update_from_pymatgen(slabs2[sl]) + slab2 = fit2host(slab1, slab2) + # slab_misfit(slab1, slab2, silent = 0) else: - slabs2 = create_surface2(st2_init, hkl_lio, shift = 0, min_slab_size = lio_thick, min_vacuum_size = 15, surface_i = i_suf_lio, oxidation = None, suf = '', - primitive = 0, symmetrize = 0, cut_thickness = 0, return_one = 1, write_poscar = 0, lll_reduce = 1) - # st2 = st2_init.update_from_pymatgen(slabs2[i_suf_lio]) - st2 = slabs2 - # st2.jmol() - slab_2 = fit2host(st1, st2) - st2_init = slab_2 + slabs2 = create_surface2(st2_init, hkl_film, surface_i = i_suf_film, min_slab_size =film_thickness, min_vacuum_size = 25, symmetrize = True, lll_reduce = 1, primitive = 1) + slab2 = st2_init.update_from_pymatgen(slabs2[i_suf_film]) + slab2 = fit2host(slab1, slab2) from siman.impurity import make_interface z_max1 = 0 - z_max2 = 50 - st1 = st1.add_vacuum(2,40) - for r in st1.xcart: + z_min2 = 50 + slab1 = slab1.add_vacuum(2,40) + for r in slab1.xcart: if r[2] > z_max1: z_max1 = r[2] - for r in slab_2.xcart: - if r[2] < z_max2: - z_max2 = r[2] + for r in slab2.xcart: + if r[2] < z_min2: + z_min2 = r[2] if 0: @@ -3316,15 +3336,16 @@ def create_interface_solid(st_host, st_oxide, suf_host, i_suf_host = 0, print('\nGood interface!\n\n') - print('\n\n\n', st1.rprimd, slab_2.rprimd, '\n\n\n',st1.get_angles(), slab_2.get_angles(),) + print('\n\n\n', slab1.rprimd, slab2.rprimd, '\n\n\n',slab1.get_angles(), slab2.get_angles(),) + + if ads_pos: interface_list = [] suf_ats2 = (slab_2.get_surface_atoms('OLiNa', surface = 0, surface_width=0.5)) - # print(suf_ats2) xc1 = ads_pos for at in suf_ats2[0:]: - xc2 = slab_2.xcart[at] - mat0, pas_scaled2 = make_interface(st1, xc1, slab_2, xc2) + xc2 = slab2.xcart[at] + mat0, pas_scaled2 = make_interface(slab1, xc1, slab2, xc2) # mat0.jmol(rep=[3,3,1]) interface_list.append(mat0) interface_list.append(pas_scaled2) @@ -3336,7 +3357,7 @@ def create_interface_solid(st_host, st_oxide, suf_host, i_suf_host = 0, - mat2, pas_scaled2 = make_interface(st1, [0, 4, z_max1+z_shift], slab_2, [0, 0.626, z_max2],) + mat2, pas_scaled2 = make_interface(slab1, [0, 4, z_max1+z_shift], slab2, [0, 0.626, z_min2],) # st_wide = st1.add_vacuum(2,15) # mat3, pas_scaled3 = make_interface(st_wide, [0, 4, z_max1+10.5], slab_2, [0, 0.626, z_max2],) diff --git a/siman/set_functions.py b/siman/set_functions.py index 052184b0..5812ea1d 100644 --- a/siman/set_functions.py +++ b/siman/set_functions.py @@ -188,7 +188,17 @@ 'MIXPRE', 'LANGEVIN_GAMMA_L', 'LANGEVIN_GAMMA', -'PMASS' +'PMASS', +'LBLUEOUT', +'ANDERSEN_PROB', +'LSPIRAL', +'LNONCOLLINEAR', +'QSPIRAL', +'ENMAX', +'ENINI', +'KPAR', +'LZEROZ' + ] vasp_keys = vasp_electronic_keys+vasp_ionic_keys+vasp_other_keys From 37fb90251639d622baeae18b429c1a7f127ef477 Mon Sep 17 00:00:00 2001 From: antonboev22 Date: Tue, 2 Dec 2025 12:38:33 +0300 Subject: [PATCH 3/6] 2_12_25 --- siman/calc_manage.py | 19 +++++++++++++++ siman/calculators/vasp.py | 10 ++++---- siman/core/calculation.py | 1 - siman/core/structure.py | 49 +++++++++++++++++++++++++++++++++++++-- siman/set_functions.py | 1 + 5 files changed, 73 insertions(+), 7 deletions(-) diff --git a/siman/calc_manage.py b/siman/calc_manage.py index b403b115..b5bbded0 100644 --- a/siman/calc_manage.py +++ b/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 @@ -455,6 +456,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 +1728,23 @@ 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 tempfile.NamedTemporaryFile('w', delete=False) as tmp: + # tmp.writelines(params['iconst_str']) + # tmp_path = tmp.name + + # shutil.copyfile(tmp_path, cl.dir+'/ICONST' ) # file is provided explicitly + # os.remove(tmp_path) + with open(cl.dir + '/ICONST', 'w') as f: + f.writelines(params['iconst_str']) + + + else: + # try: + # shutil.copyfile(dir_1+'/ICONST', cl.dir+'/ICONST' ) + # except: + printlog('Attention! no ICONST file was found!!!') # if cl.des: cl.des = ' '+struct_des[id[0]].des + '; ' + varset[id[1]].des diff --git a/siman/calculators/vasp.py b/siman/calculators/vasp.py index 8ec52b6c..af55c73a 100644 --- a/siman/calculators/vasp.py +++ b/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') @@ -334,6 +334,7 @@ def make_incar(self, filename = None): 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 @@ -868,7 +869,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: @@ -928,6 +929,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) @@ -1039,6 +1042,7 @@ def plot_energy_total_step2(self,): 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 @@ -1558,8 +1562,6 @@ def write_occmatrix(self): # print(self.get_path()) # sys.exit() - - return write_occmatrix(self.occ_matrices, self.get_path()) diff --git a/siman/core/calculation.py b/siman/core/calculation.py index 70e4a4b4..da1bad03 100644 --- a/siman/core/calculation.py +++ b/siman/core/calculation.py @@ -998,7 +998,6 @@ def full(self, ise=None, up=0, fit=1, suf='', add_loop_dic=None, up_res='up1'): self.children.append(child) return child -<<<<<<< HEAD diff --git a/siman/core/structure.py b/siman/core/structure.py index ae758fb6..790d42d5 100644 --- a/siman/core/structure.py +++ b/siman/core/structure.py @@ -2088,6 +2088,35 @@ def remove_atoms(self, atoms_to_remove, from_one = 0, clear_magmom = 1 ): # print(st.get_elements()) return st + def remove_atoms_randomly(self, el, n_at2remove, clear_magmom = 1): + """ + remove N atoms of element type *el* + st (Structure) + el (str) - element name + n_at2remove (int) - total number of atoms, which must be removed + clear_magmom - by default magmom is cleared + + """ + + from scipy.stats import qmc + + el_nums_list = self.get_elements_by_el_name(el) + + sampler = qmc.Sobol(d=1, scramble=True) + sobol_vals = sampler.random(n_at2remove).flatten() + indices = (sobol_vals * len(el_nums_list)).astype(int) + + atoms_to_remove = [] + for i in indices: + if i not in atoms_to_remove: + atoms_to_remove.append(i) + if len(atoms_to_remove) == n_at2remove: + break + st_new = self.remove_atoms(atoms_to_remove, clear_magmom = clear_magmom) + printlog('remove_atoms_randomly(): Atoms', atoms_to_remove, 'were removed') + + return st_new + def del_layers(self, xred_range = None, xcart_range = None): """ @@ -2294,24 +2323,40 @@ def remove_part(self, element, new_conc): - def add_vacuum(self, vector_i, thickness): + def add_vacuum(self, vector_i, thickness, direction = '+'): """ Allows to add or remove vacuum along one of the rprimd vectors vector_i (int) - index of vector along which vacuum should be added 0, 1, 2 thickness (float) - thickness of added (positive) or removed (negative) vacuum + direction (str, optional) - Direction where vacuum is added: + '+' (default) -> vacuum added "above" (positive direction of lattice vector) + '-' -> vacuum added "below" (negative direction of lattice vector) + TODO: add capability to add vacuum normal to surface in case of non-orthogonal cells """ st = copy.deepcopy(self) - v = st.rprimd[vector_i] + v = st.rprimd[vector_i].copy() v_l = np.linalg.norm(v) new_len = v_l+thickness + # scale the vector st.rprimd[vector_i]*=new_len/v_l + if direction == '+': + # default: do nothing, vacuum above + pass + elif direction == '-': + # shift atoms up by thickness + shift_vec = st.rprimd[vector_i] - v + st = st.shift_atoms(vector_cart=shift_vec, return2cell=0) + else: + raise ValueError("direction must be '+' or '-'") + + st.update_xred() st.name+='_vac' # st.write_xyz() diff --git a/siman/set_functions.py b/siman/set_functions.py index 1c4ab76a..3e7c03be 100644 --- a/siman/set_functions.py +++ b/siman/set_functions.py @@ -201,6 +201,7 @@ 'LMONO', 'LWANNIER90', 'LWRITE_MMN_AMN', 'LWRITE_UNK', 'NUM_WANN', 'LWRITE_SPN', 'WANNIER90_WIN', 'LWANNIER90_RUN', 'QSPIRAL', 'LSPIRAL', 'ENINI', 'ENMAX', 'LZEROZ', +'VCA' ] vasp_keys = vasp_electronic_keys+vasp_ionic_keys+vasp_other_keys From 2aa56382bac4f99f7d37240a58bc946e5c45029d Mon Sep 17 00:00:00 2001 From: antonboev22 Date: Tue, 2 Dec 2025 12:55:50 +0300 Subject: [PATCH 4/6] 2_12_25 --- siman/calculators/vasp.py | 1 - siman/core/calculation.py | 20 ++++++-------------- 2 files changed, 6 insertions(+), 15 deletions(-) diff --git a/siman/calculators/vasp.py b/siman/calculators/vasp.py index b63b9ac2..e1879042 100644 --- a/siman/calculators/vasp.py +++ b/siman/calculators/vasp.py @@ -568,7 +568,6 @@ def actualize_set(self, curset = None, params = None): curset = self.set vp = curset.vasp_params - print(params) if params is None: params = {} diff --git a/siman/core/calculation.py b/siman/core/calculation.py index da1bad03..9655cdc5 100644 --- a/siman/core/calculation.py +++ b/siman/core/calculation.py @@ -1000,23 +1000,15 @@ def full(self, ise=None, up=0, fit=1, suf='', add_loop_dic=None, up_res='up1'): 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 - - 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) - - - - + # total_mags, atom_mags_per_step = read_oszicar(path_to_file) + # print(total_mags, atom_mags_per_step) From 7947d23547894b992643e7e83ffa500b581e484b Mon Sep 17 00:00:00 2001 From: antonboev22 Date: Fri, 6 Feb 2026 11:37:59 +0300 Subject: [PATCH 5/6] feb 6 2026 --- siman/analysis.py | 12 +++- siman/calc_manage.py | 3 +- siman/calculators/vasp.py | 1 - siman/core/cluster_batch_script.py | 8 ++- siman/geo.py | 93 +++++++++++++++++------------- siman/inout.py | 24 ++++---- siman/neb.py | 4 +- 7 files changed, 83 insertions(+), 62 deletions(-) diff --git a/siman/analysis.py b/siman/analysis.py index 6ddff416..16955591 100644 --- a/siman/analysis.py +++ b/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,12 +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)) print(WS.area_fraction_dict) 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 diff --git a/siman/calc_manage.py b/siman/calc_manage.py index b5bbded0..69bcb3c0 100644 --- a/siman/calc_manage.py +++ b/siman/calc_manage.py @@ -366,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. @@ -3538,7 +3539,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 diff --git a/siman/calculators/vasp.py b/siman/calculators/vasp.py index e1879042..41bbde27 100644 --- a/siman/calculators/vasp.py +++ b/siman/calculators/vasp.py @@ -866,7 +866,6 @@ def spin(ls, i): elif 'MAGMOM' in vp and vp['MAGMOM']: #just add * to magmom tag if it is provided without it printlog('Magnetic moments from vasp_params["MAGMOM"] are used\n') - magmom = vp['MAGMOM'] # if "*" not in vp['MAGMOM']: # vp['MAGMOM'] = str(natom) +"*"+ vp['MAGMOM'] diff --git a/siman/core/cluster_batch_script.py b/siman/core/cluster_batch_script.py index a0fe33f0..b16cbb71 100644 --- a/siman/core/cluster_batch_script.py +++ b/siman/core/cluster_batch_script.py @@ -171,6 +171,8 @@ def write_batch_header(cl, batch_script_filename = None, if 'partition' in hc: f.write('#SBATCH -p '+hc['partition']+'\n') + if 'exclude' in hc: + f.write('#SBATCH --exclude='+hc['exclude']+'\n') if 'any_commands' in header.cluster: lines = header.cluster['any_commands'] @@ -182,9 +184,9 @@ def write_batch_header(cl, batch_script_filename = None, # f.write("#SBATCH -I other=avx\n") # AVX2 instructions for new node to improve speed by 18% # f.write("#SBATCH --nodelist=node-amg03\n") - if header.siman_run: #only for me - if header.EXCLUDE_NODES: - f.write("#SBATCH --exclude=node-amg11\n") + # if header.siman_run: #only for me + # if header.EXCLUDE_NODES: + # f.write("#SBATCH --exclude=node-amg11\n") # f.write("#SBATCH --mail-user=\n") # f.write("#SBATCH --mail-type=END\n") diff --git a/siman/geo.py b/siman/geo.py index 16f93ed8..fbe08fd9 100644 --- a/siman/geo.py +++ b/siman/geo.py @@ -2446,52 +2446,63 @@ def create_surface(st, miller_index, min_slab_size = 10, min_vacuum_size = 10, s return slabs[surface_i] -def stoichiometry_criteria(st1,st2): - - natom1 = st1.get_natom() - natom2 = st2.get_natom() - - tra1 = st1.get_transition_elements() - tra2 = st2.get_transition_elements() - ntra1 = len(tra1) - if ntra1 == 0: - ntra1 = natom1 - ntra2 = len(tra2) - if ntra2 == 0: - ntra2 = natom2 - rat1 = natom1/ntra1 - rat2 = natom2/ntra2 - mul = ntra1/ntra2 - - if rat1 == rat2: - return 1 - else: - return 0 +def stoichiometry_criteria(st_slab, st_bulk, silent=True, tol=1e-6): + from collections import Counter + """ + Check whether slab stoichiometry is preserved with respect to the bulk structure. -def stoichiometry_criteria2(st1,st2, silent = 1): - atoms1 = st1.get_elements() - atoms2 = st2.get_elements() + Stoichiometry is preserved if there exists a single scaling factor k such that: + N_slab(element) = k * N_bulk(element) + for all elements. - from collections import Counter - el_dict1 = Counter(atoms1) - el_dict2 = Counter(atoms2) - el1 = list(el_dict1.keys())[0] - el2 = list(el_dict1.keys())[1] - # print(el_dict1) - # print(el_dict2) - ratio1 = el_dict1[el1]/el_dict1[el2] - ratio2 = el_dict2[el1]/el_dict2[el2] - - if ratio1 == ratio2: - if not silent: - print('Stoichiometric') - return 1 - else: + Parameters + ---------- + st_bulk : Structure() - Reference bulk structure + st_slab : Structure() - Slab or supercell structure + silent : bool + If False, print diagnostic messages + tol : float + Numerical tolerance for ratio comparison + """ + + # Count elements in bulk and slab + c_bulk = Counter(st_bulk.get_elements()) + c_slab = Counter(st_slab.get_elements()) + + # 1. Check that both structures contain the same set of elements + if set(c_bulk) != set(c_slab): if not silent: - print('Non-stoichiometric') - print(round(ratio1,2), round(ratio2,2)) + print("Warning! Different set of elements") + print("bulk:", dict(c_bulk)) + print("slab:", dict(c_slab)) return 0 + # 2. Compute scaling ratios for each element + ratios = [] + for el in c_bulk: + if c_bulk[el] == 0: + return 0 + ratios.append(c_slab[el] / c_bulk[el]) + + # 3. Check that all ratios are equal within tolerance + ref = ratios[0] + for r in ratios[1:]: + if abs(r - ref) > tol: + if not silent: + print("❌ Non-stoichiometric") + for el in c_bulk: + print(f"{el}: bulk={c_bulk[el]}, slab={c_slab[el]}, ratio={c_slab[el]/c_bulk[el]:.3f}") + return 0 + + if not silent: + print(f"✅ Stoichiometric (scale factor ≈ {ref:.3f})") + # for el in c_bulk: + # print(f"{el}: bulk={c_bulk[el]}, slab={c_slab[el]}, ratio={c_slab[el]/c_bulk[el]:.3f}") + + + + return 1 + def symmetry_criteria(st): from pymatgen.symmetry.analyzer import SpacegroupAnalyzer st = st.convert2pymatgen() diff --git a/siman/inout.py b/siman/inout.py index 81d832d4..d91d6375 100644 --- a/siman/inout.py +++ b/siman/inout.py @@ -1756,10 +1756,7 @@ def read_vasp_out(cl, load = '', out_type = '', show = '', voronoi = '', path_to de_each_md_list.append(de_each_md) - if "kinetic energy EKIN =" in line: - self.ekin = float(line.split()[4]) - self.list_ekin.append(self.ekin) - + if "(temperature" in line: self.list_temp.append(float(line.split()[-2])) @@ -1857,14 +1854,17 @@ def read_vasp_out(cl, load = '', out_type = '', show = '', voronoi = '', path_to if "Iteration" in line: - self.mdstep = int(line.split('(')[0].split()[2].strip()) - iterat +=1 - # print self.mdstep - # print line - if mdstep_old != self.mdstep: - nscflist.append( niter ) # add to list number of scf iterations during mdstep_old - niter = int(line.split(')')[0].split('(')[-1].strip()) #number of scf iterations - mdstep_old = self.mdstep + try: + self.mdstep = int(line.split('(')[0].split()[2].strip()) + iterat +=1 + # print self.mdstep + # print line + if mdstep_old != self.mdstep: + nscflist.append( niter ) # add to list number of scf iterations during mdstep_old + niter = int(line.split(')')[0].split('(')[-1].strip()) #number of scf iterations + mdstep_old = self.mdstep + except ValueError: + print(line) if 'number of electron ' in line: diff --git a/siman/neb.py b/siman/neb.py index 63dcbad4..a82134ea 100644 --- a/siman/neb.py +++ b/siman/neb.py @@ -370,7 +370,7 @@ def add_neb(starting_calc = None, st = None, st_end = None, # printlog('atom ', type_atom_to_move, 'will be moved', imp ='y') - if i_atom_to_move: + if i_atom_to_move != None: printlog('add_neb(): *i_atom_to_move* = ', i_atom_to_move, 'is used', imp ='y') numbers = [[i_atom_to_move]] i_void_start = 1 @@ -379,7 +379,7 @@ def add_neb(starting_calc = None, st = None, st_end = None, numbers = determine_symmetry_positions(st, atom_to_move) - # print(numbers) + print('+++++++++++++++++++++++++++',numbers) # sys.exit() if len(numbers)>0: printlog('Please choose position using *i_void_start* :', [i+1 for i in range(len(numbers))],imp = 'y' ) From 1f1ede36dc0507266f073a8c2d95b34d16a32141 Mon Sep 17 00:00:00 2001 From: antonboev22 Date: Thu, 12 Feb 2026 13:50:34 +0300 Subject: [PATCH 6/6] 12 02 26 --- siman/external_tools/nebmaker.py | 2 +- siman/neb.py | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/siman/external_tools/nebmaker.py b/siman/external_tools/nebmaker.py index 24aec1d7..14866a45 100644 --- a/siman/external_tools/nebmaker.py +++ b/siman/external_tools/nebmaker.py @@ -52,7 +52,7 @@ def nebmaker(init, fin, n, idpp = True): neb.interpolate(mic=True) else: neb.interpolate('idpp', mic=True) - dir_names = ['0'+str(i) if i < 10 else str(i) for i in range(len(images))] + dir_names = ['0'+str(i) if i < 10 else str(i) for i in range(len(images))] for i, image in zip(dir_names,neb.images): if not os.path.isdir(i): os.mkdir(i) diff --git a/siman/neb.py b/siman/neb.py index a82134ea..08c518f5 100644 --- a/siman/neb.py +++ b/siman/neb.py @@ -99,7 +99,8 @@ def add_neb(starting_calc = None, st = None, st_end = None, x_final = None, xr_final = None, upload_vts = False, center_on_moving = True, - run = False, add_loop_dic = None, old_behaviour = None, params = None + run = False, add_loop_dic = None, old_behaviour = None, params = None, + idpp = 0 ): @@ -867,7 +868,8 @@ def unique_by_key(elements, key=None): else: # a = runBash(header.PATH2NEBMAKE+' POSCAR1 POSCAR2 3') # print(a) - nebmaker('POSCAR1','POSCAR2', images) + printlog('Starting nebmaker ...', imp = 'y') + nebmaker('POSCAR1','POSCAR2', images, idpp = False) dst = it_new+'_all' makedir(dst+'/any') for f in ['00', '01', '02', '03', '04']: