Crystal-MET/src/utils.py

import numpy as np
import operator
import sys

def distance(a,b):
    # Returns the 3D distance between a and b where 
    # a and b are array where x, y and z are at the
    # position 0, 1 and 2

    x = a[0]-b[0]
    y = a[1]-b[1]
    z = a[2]-b[2]

    return np.sqrt(x**2+y**2+z**2)

def get_cell_matrix(a,b,c,alpha,beta,gamma):
    # Computing the volume of the primitive cell
    omega = a*b*c*np.sqrt(1-np.cos(alpha)**2-np.cos(beta)**2-np.cos(gamma)**2+2*np.cos(alpha)*np.cos(beta)*np.cos(gamma))

    # Computing the matrix 
    M = [
            [a,b*np.cos(gamma),c*np.cos(beta)],
            [0,b*np.sin(gamma),c*(np.cos(alpha)-np.cos(beta)*np.cos(gamma))/(np.sin(gamma))],
            [0,0,omega/(a*b*np.sin(gamma))]
        ]
    return M

def big_cell(generator,symGenerator,a,b,c,alpha,beta,gamma,nA,nB,nC):
    coords = []
    
    # Computing the matrix converting fractional to cartesian
    fracToCart = get_cell_matrix(a,b,c,alpha,beta,gamma)

    for gen in generator:
        x = gen[1]
        y = gen[2]
        z = gen[3]

        for sym in symGenerator:
            u = eval(sym[0])
            v = eval(sym[1])
            w = eval(sym[2])

            # Making sure the value is within the range [0,1]
            u = u + 1*(u<0) - 1*(u>1)
            v = v + 1*(v<0) - 1*(v>1)
            w = w + 1*(w<0) - 1*(w>1)
            
            coords.append([u,v,w,gen[0]])

    # Deleting the redundant atoms
    toDel = []
    for i in range(len(coords)-1):
            for j in range(i+1,len(coords)):
                # Computing the distance using the minimum image convention
                # as described in Appendix B equation 9 of 
                # "Statistical Mechanics : Theory and Molecular Simulations
                # Mark E. Tuckerman"

                r1 = np.array(coords[i][:3])
                r2 = np.array(coords[j][:3])
                r12 = r1-r2
                da = np.sqrt(r12[0]**2+r12[1]**2+r12[2]**2)
                r12 = r12 - np.round(r12)
                db = da - np.sqrt(r12[0]**2+r12[1]**2+r12[2]**2)
                r12 = np.matmul(fracToCart,r12)
                d = np.sqrt(r12[0]**2+r12[1]**2+r12[2]**2)

                if(d<1e-2):
                    # We check if we don't already want to delete this atom
                    if j not in toDel:
                        toDel.append(j)

    toDel = sorted(toDel)

    # We delete the atoms in the list
    for i in range(len(toDel)):
        coords.pop(toDel[i]-i)

    newCoords = []

    # We replicate the cell nA, nB, nC times 
    for at in coords:
        newCoords.append([at[0],at[1],at[2],at[3]])
        for a in range(1,nA):
            newCoords.append([at[0]+a,at[1],at[2],at[3]])
            for b in range(1,nB):
                newCoords.append([at[0]+a,at[1]+b,at[2],at[3]])
                for c in range(1,nC):
                    newCoords.append([at[0]+a,at[1]+b,at[2]+c,at[3]])
            for c in range(1,nC):
                newCoords.append([at[0]+a,at[1],at[2]+c,at[3]])
        for b in range(1,nB):
            newCoords.append([at[0],at[1]+b,at[2],at[3]])
            for c in range(1,nC):
                newCoords.append([at[0],at[1]+b,at[2]+c,at[3]])
        for c in range(1,nC):
            newCoords.append([at[0],at[1],at[2]+c,at[3]])

    # Now we convert the fractionnal coordinates to cartesian coordinates
    coords = []

    for at in newCoords:
        r = [at[0],at[1],at[2]]
        rxyz = np.matmul(fracToCart,r)
        coords.append([rxyz[0],rxyz[1],rxyz[2],at[3],'C'])


    # Returns the list of the atoms [x,y,z,label,second_label]
    return coords

# Translates all the coordinates with the vector v
def translate(v,coordinates):
    for c in coordinates:
        c[0] += v[0]
        c[1] += v[1]
        c[2] += v[2]
    return coordinates

# Finds the point at the center of the given atoms that are the 
# closest to the origin
def find_center(centerList, coordinates, without=[]):

    centers = []
    for i in range(len(centerList)):
        centers.append([100,100,100])  # Setting a large value for each center

    for c in centers:
        c.append(distance(c,[0,0,0]))  # Computing the distance to the origin

    for at in coordinates:
        w = True
        for a in without:
            if distance(at, a) < 1e-6:
                w = False
                break
        if not w:
            continue
        if at[3] in centerList:
            centers = sorted(centers, key=operator.itemgetter(3)) # Sorting the list with respect to the distance to the origin
            d = distance(at,[0,0,0])
            if d <= centers[-1][-1] and d > 0.0:
                centers[-1] = [at[0],at[1],at[2],d]

    center = np.mean(centers,axis=0)[:3] # Computing the barycenter

    return center

# Defines a rotation matrix that will put r1 at the position r2
def rotation_matrix(r1,r2):

    r1 = np.array(r1)/np.linalg.norm(r1)
    r2 = np.array(r2)/np.linalg.norm(r2)

    # Computing the cross product which is the vector around which
    # the rotation is done
    crossProduct = np.cross(r1,r2)
    crossProduct = crossProduct/np.linalg.norm(crossProduct)

    # Computing the angle of the rotation
    a = np.arccos(np.dot(r1,r2))

    c = np.cos(a)
    s = np.sin(a)
    x = crossProduct[0]
    y = crossProduct[1]
    z = crossProduct[2]

    M = [
            [x**2*(1-c)+c,x*y*(1-c)-z*s,x*z*(1-c)+y*s],
            [x*y*(1-c)+z*s,y**2*(1-c)+c,y*z*(1-c)-x*s],
            [x*z*(1-c)-y*s,y*z*(1-c)+x*s,z**2*(1-c)+c]
            ]

    return M

# Rotates all the coordinates using the rotation matric M
def rotate(M,coordinates):
    for i in range(len(coordinates)):
        r = [coordinates[i][0],coordinates[i][1],coordinates[i][2]]
        rV = np.matmul(M,r)
        coordinates[i][0] = rV[0]
        coordinates[i][1] = rV[1]
        coordinates[i][2] = rV[2]

    return coordinates

# Cuts a sphere centered on the origin in the coordinates
def cut_sphere(coordinates,r):
    sphere = []
    for i in range(len(coordinates)):
        if distance(coordinates[i],[0,0,0]) <= r:
            sphere.append(coordinates[i])

    return sphere

# Finds the fragment in the coordinates
def find_fragment(coordinates, patterns, npatterns,notInFrag):

    inFrag = []

    for n in range(len(patterns)):
        pattern = patterns[n]
        npattern = npatterns[n]
        for i in range(npattern):
            c = [100,100,100]
            dc = distance([0,0,0],c)

            inPattern = []
            # Finding the closest atom of the first type in the pattern
            for at in coordinates:
                if at[3] == pattern[1]:
                    d = distance([0,0,0],at)
                    if d > 10:
                        break
                    if d < dc :
                        accept = True
                        for exc in notInFrag:
                            d = distance(exc,at)
                            if d < 1e-5:
                                accept = False
                        if accept and coordinates.index(at) not in inFrag:
                            c = [at[0],at[1],at[2],0.0, coordinates.index(at)]
                            dc = distance([0,0,0],c)
            # Finding the rest of the pattern around the atom previously found
            atIn = []
            for j in range(0,len(pattern),2):
                d = distance(c,[100,100,100])
                # Initializing the atoms 
                for k in range(pattern[j]):
                    atIn.append([100,100,100,d])

                for at in coordinates:
                    if distance(at,[0,0,0]) > 10:
                        break
                    if at[3] == pattern[j+1]:
                        atIn = sorted(atIn,key=operator.itemgetter(3))
                        d = distance(at,c)
                        trial = [at[0],at[1],at[2],d,coordinates.index(at)]
                        if d < atIn[-1][3] and trial not in atIn:
                            accept = True
                            for exc in notInFrag:
                                d = distance(exc,trial)
                                if d < 1e-5:
                                    accept = False
                            if accept:
                                atIn[-1] = trial
            for at in atIn:
                inPattern.append(at[4])

            for at in inPattern:
                if at not in inFrag:
                    inFrag.append(at)

    for at in inFrag:
        coordinates[at][4] = 'O'

    return len(inFrag), coordinates


# Finds the pseudopotential layer around
# the fragment
def find_pseudo(coordinates, rPP, notInPseudo):

    for at in coordinates:
        if at[4] != 'O':
            continue
        for i in range(len(coordinates)):
            if coordinates[i][4] != 'C':
                continue
            d = distance(at,coordinates[i])
            if d < rPP:
                coordinates[i][4] = 'Cl'

    return coordinates

# Creates lists containing the neighbours of each
# atom
def find_neighbours(coordinates, atoms):
    neighbourList = [[] for i in range(len(coordinates))]

    atoms = np.array(atoms).flatten()

    for i in range(len(coordinates)-1):
        for j in range(i+1,len(coordinates)):
            li = coordinates[i][3] # Label of the atom i
            lj = coordinates[j][3] # Label of the atom j

            ii = np.where(atoms==li)[0]
            jj = np.where(atoms==lj)[0]

            ci = float(atoms[ii+1]) # Charge of the atom i
            cj = float(atoms[jj+1]) # Charge of the atom j

            if ci*cj < 0: # Checking if the charges have opposite signs
                d = distance(coordinates[i],coordinates[j])

                if d < float(atoms[ii+3]) and d < float(atoms[jj+3]):
                    neighbourList[i].append(j)
                    neighbourList[j].append(i)
    return neighbourList

# For each atom, finds if it has the correct number of neighbours,
# if not, modify its charge 
def evjen_charges(coordinates,atoms):
    neighbourList = find_neighbours(coordinates,atoms)

    atoms = np.array(atoms).flatten()

    charges = []

    for i in range(len(coordinates)):
        li = coordinates[i][3]
        ii = np.where(atoms==li)[0]

        nr = len(neighbourList[i])
        nt = int(atoms[ii+2])
        ci = float(atoms[ii+1])
        
        if nr > nt:
            print("Error : too much neighbours for atom n°%i, count %i neighbours where it should have a maximum of %i"%(i,nr,nt))
            sys.exit()
        charges.append(ci*nr/nt)

    return charges

# Computes the nuclear repulsion 
def nuclear_repulsion(coordinates,charges):

    rep = 0.0
    
    for i in range(len(coordinates)-1):
        for j in range(i+1,len(coordinates)):
            rij = distance(coordinates[i],coordinates[j])
            ci = charges[i]
            cj = charges[j]

            if(rij < 1):
                print(i,j,"\n",coordinates[i],"\n",coordinates[j],"\n",rij)

            rep += (ci*cj)/rij
    return rep

# Computes the symmetry in the whole system
def compute_symmetry(coordinates,charges,symmetry):
    symmetrizedCoordinates = []
    symmetrizedCharges = []
    uniqueIndexList = [] # The list containing the indexes of the unique atoms

    treated = [] # Will store the index of the atoms already treated

    symOp = []

    # Storing all the symmetry operations
    for s in symmetry:
        if s == 'C2x':
            symOp.append(np.array([1,-1,-1]))
        elif s == 'C2y':
            symOp.append(np.array([-1,1,-1]))
        elif s == 'C2z':
            symOp.append(np.array([-1,-1,1]))
        elif s == 'xOy':
            symOp.append(np.array([1,1,-1]))
        elif s == 'xOz':
            symOp.append(np.array([1,-1,1]))
        elif s == 'yOz':
            symOp.append(np.array([-1,1,1]))
        elif s == 'i':
            symOp.append(np.array([-1,-1,-1]))

    for i in range(len(coordinates)):
        if i in treated:
            continue
        
        treated.append(i)
        at1 = np.array(coordinates[i][:3])
        symmetrizedCoordinates.append(coordinates[i])
        symmetrizedCharges.append(charges[i])
        uniqueIndexList.append(len(symmetrizedCoordinates)-1)

        for j in range(len(coordinates)):
            if j in treated or coordinates[i][3] != coordinates[j][3]:
                continue

            at2 = np.array(coordinates[j][:3])

            for s in symOp:
                if distance(at1,at1*s) > 1e-4 and distance(at2,at1*s) < 1e-4: # Checking if op.at1 != at1 and that op.at2 = at1
                    p = at1*s
                    treated.append(j)
                    symmetrizedCoordinates.append([p[0],p[1],p[2],coordinates[i][3],coordinates[i][4]])
                    symmetrizedCharges.append(charges[i])
                    break

    return symmetrizedCoordinates,symmetrizedCharges,uniqueIndexList