ScalableWorkflow_VASP_Calculations/generate_vasp_atomic_configurations_binary.py at main · ORNL/ScalableWorkflow_VASP_Calculations · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 25 11:56:29 2022

@author: 7ml
"""

import os
import numpy as np
import matplotlib.pyplot as plt
import shutil

from ase.io.vasp import read_vasp, write_vasp

from itertools import combinations, permutations

plt.rcParams.update({"font.size": 20})

from math import factorial

# Fixing the seed
np.random.seed(42)

def nPr(n, r):
    return int(factorial(n) / factorial(n - r))


def nCr(n, r):
    return int(factorial(n) / (factorial(n - r) * factorial(r)))


def create_vasp_case(case_path, composition_ase_object):
    os.makedirs(case_path, exist_ok=False)
    shutil.copy(source_path + "/START", case_path)
    shutil.copy(source_path + "/POTCAR", case_path)
    shutil.copy(source_path + "/KPOINTS", case_path)
    write_vasp(case_path + "/POSCAR", direct=True, sort=True, atoms=composition_ase_object)


def generate_VASP_randomized_binary_configurations(source_path, destination_path, atomic_increment, symmetry_savings,
                                                   num_atomic_configurations=1):
    " Read atomic configuration from POSCAR file "

    os.makedirs(destination_path, exist_ok=False)

    assert os.path.exists(source_path + "/START")
    assert os.path.exists(source_path + "/POTCAR")
    assert os.path.exists(source_path + "/KPOINTS")
    assert os.path.exists(source_path + "/POSCAR")

    ase_object = read_vasp(source_path + "/" + "POSCAR")

    # Extract set of atom elements
    elements_list = list(set(list(ase_object.numbers)))
    assert 2 == len(elements_list), "This function can only be used for binary alloys"
    num_atoms = ase_object.numbers.shape[0]

    chemical_composition = np.zeros_like(ase_object.numbers)

    for num_element_atoms in range(0, num_atoms + 1, atomic_increment):

        for composition_idx in range(0, num_element_atoms):
            chemical_composition[composition_idx] = elements_list[0]

        for composition_idx in range(num_element_atoms, num_atoms):
            chemical_composition[composition_idx] = elements_list[1]

        composition_ase_object = ase_object
        composition_ase_object.numbers = chemical_composition
        chemical_composition_dir = destination_path + "/" + str(composition_ase_object.symbols)

        os.makedirs(chemical_composition_dir, exist_ok=False)

        if nCr(num_atoms, num_element_atoms) > num_atomic_configurations:
            for randomized_configuration_id in range(0, num_atomic_configurations):
                case_path = chemical_composition_dir + "/case-" + str(randomized_configuration_id + 1)
                # When the histogram of the distribution is cut,
                # a check is performed to see if the total number of configurations is lower or higher than the number of atomic configurations needed
                np.random.shuffle(composition_ase_object.positions)
                create_vasp_case(case_path, composition_ase_object)
        elif symmetry_savings:
            # FIXME not implemented yet
            raise TypeError("Error: symmetry savings not implemented yet")
        elif not symmetry_savings:
            indices_list = [index for index in range(0, num_atoms)]
            if num_element_atoms < int(num_atoms / 2):
                element_original_indices = np.where(chemical_composition == elements_list[0])[0]
                num_swapping = num_element_atoms
            else:
                element_original_indices = np.where(chemical_composition == elements_list[1])[0]
                num_swapping = num_atoms - num_element_atoms
            permutation_count = 0
            # if the total number of permutations is smaller than the total number of configurations needed,
            # randomization is not needed because we need to include all the permutations in the dataset
            for permutation in permutations(indices_list, num_swapping):
                case_path = chemical_composition_dir + "/case-" + str(permutation_count + 1)
                new_ase_object = composition_ase_object
                # we need to exchange the rows of the position matrix associated with the atom types that need to be swapped
                for original_index, new_index in zip(element_original_indices, list(permutation)):
                    new_ase_object.positions[[original_index, new_index]] = composition_ase_object.positions[
                        [new_index, original_index]]
                create_vasp_case(case_path, new_ase_object)
                permutation_count += 1
            assert nCr(num_atoms, num_element_atoms) == permutation_count


if __name__ == '__main__':
    current_directory = os.getcwd()
    source_path = current_directory + '/template_Ti-Cr'
    destination_path = current_directory + '/bcc/Ti-Cr'
    atomic_increment = 4
    num_atomic_configurations_per_composition = 100
    symmetry_savings = False
    generate_VASP_randomized_binary_configurations(source_path, destination_path, atomic_increment, symmetry_savings,
                                                   num_atomic_configurations_per_composition)