utils.py 3.05 KB
Newer Older
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
import numpy as np
from abc import ABC, abstractmethod
from ase.data import covalent_radii
from ase.geometry import get_distances
from ase import Atoms

from ase.visualize import view

def check_valid_bondlengths(a, blmin=None, blmax=None, indices=None, check_too_close=True, check_isolated=True):
    """Calculates if the bondlengths between atoms with indices
    in 'indices' and all other atoms are valid. The validity is
    determined by blmin and blmax.

    Parameters:

    a: Atoms object

    blmin: The minimum allowed distance between atoms in units of
    the covalent distance between atoms, where d_cov=r_cov_i+r_cov_j.
    
    blmax: The maximum allowed distance, in units of the covalent 
    distance, from a single isolated atom to the closest atom. If
    blmax=None, no constraint is enforced on isolated atoms.

    indices: The indices of the atoms of which the bondlengths
    with all other atoms is checked. if indices=None all bondlengths
    are checked.
    """
    bl = get_distances_as_fraction_of_covalent(a,indices)
    
    # Filter away self interactions.
    bl = bl[bl > 1e-6].reshape(bl.shape[0],bl.shape[1]-1)
    
    # Check if atoms are too close
    if blmin is not None and check_too_close:
        tc = np.any(bl < blmin)
    else:
        tc = False
        
    # Check if there are isolated atoms
    if blmax is not None and check_isolated:
        isolated = np.any(np.all(bl > blmax, axis=1))
    else:
        isolated = False
        
    is_valid = not tc and not isolated
    return is_valid
    
def get_distances_as_fraction_of_covalent(a, indices=None, covalent_distances=None):
    if indices is None:
        indices = np.arange(len(a))
        
    if covalent_distances is None:
        cd = get_covalent_distance_from_atom_numbers(a, indices=indices)
    else:
        cd = covalent_distances[indices,:]
    _, d = get_distances(a[indices].positions,
                         a.positions,
                         cell=a.get_cell(),
                         pbc=a.get_pbc())
    bl = d/cd
    return bl

def get_covalent_distance_from_atom_numbers(a, indices=None):
    r_cov = np.array([covalent_radii[n] for n in a.get_atomic_numbers()])
    if indices is None:
        r_cov_sub = r_cov
    else:
        r_cov_sub = r_cov[indices]
    cd_mat = r_cov_sub.reshape(-1,1) + r_cov.reshape(1,-1)
    return cd_mat
72
73
74
75
76
77
78

def get_min_distances_as_fraction_of_covalent(a, indices=None, covalent_distances=None):
    bl = get_distances_as_fraction_of_covalent(a,indices)
    
    # Filter away self interactions.
    bl = bl[bl > 1e-6].reshape(bl.shape[0],bl.shape[1]-1)
    
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
79
80
81
82
83
84
85
86
87
88
89
90
91
92
    return np.min(bl), bl.min(axis=1).argmin()

def array_to_string(arr, unit='', format='0.4f', max_line_length=80):
    msg = ''
    line_length_counter = 0
    for i, x in enumerate(arr):
        string = f'{i} = {x:{format}}{unit},  '
        #string = f"{f'{i}={x:{format}}{unit},':15}"
        line_length_counter += len(string)
        if line_length_counter >= max_line_length:
            msg += '\n'
            line_length_counter = len(string)
        msg += string
    return msg