Storion Wiki

This is an example of code which allows to calculate percolation barrier. The module is available under request from the CMG group by Aksyonov Dmitry

import os
import pandas as pd
from ase.io import read, write
from pymatgen.io.ase import AseAtomsAdaptor
from siman.header import db
from siman.calc_manage import res_loop
from siman.neb import add_neb
from siman import header
from siman.inout import smart_structure_read
from siman.set_functions import read_vasp_sets
from project_sets import user_vasp_sets
 
 
varset = read_vasp_sets(user_vasp_sets, override_global = 1) #read user sets
 
def _counter(folder):
    files = os.listdir(folder)
    init_traj_files = [f for f in files if 'trajectory_init_edge' in f]
    relaxed_traj_files = [f for f in files if 'trajectory_init_edge' in f]
    if len(init_traj_files) != len(relaxed_traj_files):
        print('n_relaxed < n_init trajectories')
        raise
    return len(init_traj_files)
 
# select materials with E_1D < 0.6
barrier_file = '/Users/artemdembitskiy/Desktop/projects/rsf_solid_state/src/data/Li_percolation_barriers_MACE.csv'
df = pd.read_csv(barrier_file)
df = df[df.e3d > 0]
df = df[df.e3d < .6].sort_values(by = 'e3d').reset_index()
 
# base folder
geo_path = '/Users/artemdembitskiy/Desktop/projects/rsf_solid_state/MACE_Li_percolation_barriers'
 
# calc setup
ise_new = 'mprelax_neb'
 
print(len(df))
 
batch_start, batch_stop = 0, 500
df = df[batch_start: batch_stop]
 
 
read_results = False
 
 
for material_id, em_mlip in zip(df.material_id.values, df.e1d.values):
    material_folder = f"{material_id}.neb"
    n_edges = _counter(f"{geo_path}/{material_folder}")
 
    for edge_id in range(n_edges):
 
        init_neb_geo_fld = f'{material_folder}/edge_id_{edge_id}'
        init_neb_geo_full_path = os.path.join(geo_path, init_neb_geo_fld)
 
        it_new = init_neb_geo_fld.replace('/', '.')
        it_new_folder = 'HT_neb_coatings/'+os.path.dirname(init_neb_geo_fld)
 
        print(it_new, it_new_folder, round(em_mlip, 3))
        add_loop_dic = {'coord':'car'}
 
        if not read_results:
            st = smart_structure_read(f'{init_neb_geo_full_path}/00/POSCAR')
            nelectrons = st.get_total_number_electrons(header.varset['mprelax_neb'])
            print(f"NELECT: {nelectrons}")
 
            if nelectrons > 800:
                continue
 
            save_dir = f'/Users/artemdembitskiy/dft_env/src/{it_new_folder}.edge_id_{edge_id}'
            print('\n\n\n\n\n',save_dir, os.path.exists(save_dir), 'n\n\n\n\n\n')
            if os.path.exists(save_dir):
                print('Calculation folder already exists. Go to the next step.')
                continue
            add_neb(
                it_new = it_new,                     # structure name
                ise_new = ise_new,                   # setup      ->>>>> calc_name.setup.image_number - name in db 
                it_new_folder = it_new_folder,       # parent folder
                init_neb_geo_fld = init_neb_geo_full_path, # trajectory
                add_loop_dic = add_loop_dic,         # optional
                up = 'up2'
                )
 
        else:
            traj = read(f"{geo_path}/{material_folder}/edge_id_{edge_id}/traj_init.xyz", index = ':')
            n_images = len(traj)
            ids = [i for i in range(n_images + 1)]
            res_loop(it_new, ise_new, ids,
                      #show =  'fomepp.neb_geo2',
                    analys_type = 'neb')
 
            try:
                calcs = []
                for i in range(1, n_images + 1):
                    calc = db[it_new, ise_new, i]
                    calcs.append(calc)
 
                calcs_sorted = []
                calcs_sorted.append(calcs[0])
                for calc in calcs[2:]:
                    calcs_sorted.append(calc)
                calcs_sorted.append(calcs[1])
 
                trajectory_relaxed = []
                for calc in calcs_sorted:
                    st = calc.end.convert2pymatgen()
                    image = AseAtomsAdaptor.get_atoms(st)
                    image.info.update({'potential_energy': calc.e0})
                    image.info.update({'max_force': calc.maxforce})
                    trajectory_relaxed.append(image)
                save_pth = f"{geo_path}/{material_folder}/edge_id_{edge_id}/traj_relaxed_PBE.xyz"
                write(save_pth, trajectory_relaxed)
                print('Save relaxed trajectory to:', save_pth)
            except:
                print('Bad calc')