diff --git a/siman/analysis.py b/siman/analysis.py index 1dec537..1695559 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,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)) + 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 b403b11..69bcb3c 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 @@ -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,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 @@ -3519,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 1679e5c..41bbde2 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') @@ -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()) diff --git a/siman/core/calculation.py b/siman/core/calculation.py index eea94af..9655cdc 100644 --- a/siman/core/calculation.py +++ b/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 diff --git a/siman/core/cluster_batch_script.py b/siman/core/cluster_batch_script.py index a0fe33f..b16cbb7 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/core/structure.py b/siman/core/structure.py index 1df5c87..f8e8c49 100644 --- a/siman/core/structure.py +++ b/siman/core/structure.py @@ -294,6 +294,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] @@ -587,7 +591,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 @@ -984,6 +991,8 @@ def rprimd_len(self): + + def pvec(self): """ print primitive vectors in formatted way @@ -1911,6 +1920,8 @@ def del_atom(self, iat): '' # print ('No!') + if hasattr(st, 'select') and any(st.select): + del st.select[i] st.natom-=1 @@ -2077,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): """ @@ -2283,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/external_tools/nebmaker.py b/siman/external_tools/nebmaker.py index 24aec1d..14866a4 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/geo.py b/siman/geo.py index dd10188..fbe08fd 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() @@ -2826,7 +2837,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 = [] @@ -3195,43 +3206,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()) - return st_oxide + + 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_film_new -def hkl_slab(st, st_host, hkl, i_suf = None): +def hkl_slab(st, st_host, hkl, i_suf = None, path = 'filename'): - os.remove('/home/anton/media/vasp/log_best1') - f = open('/home/anton/media/vasp/log_best1', 'a') + + 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): @@ -3239,7 +3278,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: @@ -3249,25 +3288,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 @@ -3278,33 +3314,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: @@ -3382,15 +3413,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) @@ -3402,7 +3434,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/inout.py b/siman/inout.py index cbfcce6..d91d637 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])) @@ -1774,12 +1771,16 @@ 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.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) + @@ -1853,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: @@ -1886,14 +1890,13 @@ def read_vasp_out(cl, load = '', out_type = '', show = '', voronoi = '', path_to if 'magnetization (x)' in line: - # print(line) + mags = [] if ifmaglist is not None: 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 diff --git a/siman/matproj_functions.py b/siman/matproj_functions.py index 3c9a9fc..0cfee21 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: diff --git a/siman/neb.py b/siman/neb.py index 63dcbad..08c518f 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 ): @@ -370,7 +371,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 +380,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' ) @@ -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']: diff --git a/siman/set_functions.py b/siman/set_functions.py index cbe6ba4..e5759cc 100644 --- a/siman/set_functions.py +++ b/siman/set_functions.py @@ -187,6 +187,10 @@ 'LANGEVIN_GAMMA_L', 'LANGEVIN_GAMMA', 'PMASS', +'LBLUEOUT', +'ANDERSEN_PROB', +'LNONCOLLINEAR', +'KPAR', 'ZAB_VDW', 'IWAVPR', 'NBLOCK', @@ -194,13 +198,10 @@ 'KPAR', 'LMIXTAU', 'LORBMOM', -'LNONCOLLINEAR', '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