Commit dd3edadb authored by Malthe Kjær Bisbo's avatar Malthe Kjær Bisbo
Browse files

added compiled code for descriptor and repulsive prior

parent 4eb061b1
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.
import numpy as np
from itertools import product
from ase.ga.utilities import closest_distances_generator
cimport numpy as np
from libc.math cimport *
from cymem.cymem cimport Pool
cimport cython
# Custom functions
ctypedef struct Point:
double coord[3]
cdef Point subtract(Point p1, Point p2):
cdef Point p
p.coord[0] = p1.coord[0] - p2.coord[0]
p.coord[1] = p1.coord[1] - p2.coord[1]
p.coord[2] = p1.coord[2] - p2.coord[2]
return p
cdef Point add(Point p1, Point p2):
cdef Point p
p.coord[0] = p1.coord[0] + p2.coord[0]
p.coord[1] = p1.coord[1] + p2.coord[1]
p.coord[2] = p1.coord[2] + p2.coord[2]
return p
cdef double norm(Point p):
return sqrt(p.coord[0]*p.coord[0] + p.coord[1]*p.coord[1] + p.coord[2]*p.coord[2])
cdef double euclidean(Point p1, Point p2):
return norm(subtract(p1,p2))
def convert_atom_types(num):
cdef int Natoms = len(num)
cdef list atomic_types = sorted(list(set(num)))
cdef int Ntypes = len(atomic_types)
cdef list num_converted = [0]*Natoms
cdef int i, j
for i in range(Natoms):
for j in range(Ntypes):
if num[i] == atomic_types[j]:
num_converted[i] = j
return num_converted
class delta():
def __init__(self, atoms, rcut=5, ratio_of_covalent_radii=0.7):
self.rcut = rcut
self.pbc = atoms.get_pbc()
self.cell = atoms.get_cell()
self.cell_displacements = self.__get_neighbour_cells_displacement(self.pbc, self.cell)
self.Ncells = len(self.cell_displacements)
num = atoms.get_atomic_numbers()
atomic_types = sorted(list(set(num)))
self.Ntypes = len(atomic_types)
blmin_dict = closest_distances_generator(num, ratio_of_covalent_radii=ratio_of_covalent_radii)
self.blmin = -np.ones((self.Ntypes, self.Ntypes))
for i1, type1 in enumerate(atomic_types):
for i2, type2 in enumerate(atomic_types):
self.blmin[i1,i2] = blmin_dict[(type1, type2)]
self.blmin = self.blmin.reshape(-1).tolist()
def energy(self, a):
# Memory allocation pool
cdef Pool mem
mem = Pool()
cdef int Natoms = a.get_number_of_atoms()
cdef double rcut = self.rcut
cdef list x_np = a.get_positions().tolist()
cdef Point *x
x = <Point*>mem.alloc(Natoms, sizeof(Point))
cdef int m
for m in range(Natoms):
x[m].coord[0] = x_np[m][0]
x[m].coord[1] = x_np[m][1]
x[m].coord[2] = x_np[m][2]
# Get neighbourcells and convert to Point-struct
cdef int Ncells = self.Ncells
cdef list cell_displacements_old = self.cell_displacements
cdef Point *cell_displacements
cell_displacements = <Point*>mem.alloc(Ncells, sizeof(Point))
for m in range(Ncells):
cell_displacements[m].coord[0] = cell_displacements_old[m][0]
cell_displacements[m].coord[1] = cell_displacements_old[m][1]
cell_displacements[m].coord[2] = cell_displacements_old[m][2]
# Convert 2body bondtype list into c-array
cdef int Ntypes = self.Ntypes
cdef list blmin_old = self.blmin
cdef double *blmin
blmin = <double*>mem.alloc(Ntypes*Ntypes, sizeof(double))
for m in range(Ntypes*Ntypes):
blmin[m] = blmin_old[m]
# Get converted atom Ntypes
cdef list num_converted_old = convert_atom_types(a.get_atomic_numbers())
cdef int *num_converted
num_converted = <int*>mem.alloc(Natoms, sizeof(int))
for m in range(Natoms):
num_converted[m] = num_converted_old[m]
cdef double E=0, rmin, radd
cdef int i, j
for i in range(Natoms):
xi = x[i]
for cell_index in range(Ncells):
displacement = cell_displacements[cell_index]
for j in range(Natoms):
rmin = blmin[num_converted[i]*Ntypes + num_converted[j]]
radd = 1 - rmin
xj = add(x[j], displacement)
r = euclidean(xi, xj)
if r > 1e-5 and r < rcut:
r_scaled = r + radd
E += 1/pow(r_scaled,12)
return E/2
def forces(self, a):
# Memory allocation pool
cdef Pool mem
mem = Pool()
cdef int Natoms = a.get_number_of_atoms()
cdef int dim = 3
cdef double rcut = self.rcut
cdef list x_np = a.get_positions().tolist()
cdef Point *x
x = <Point*>mem.alloc(Natoms, sizeof(Point))
cdef int m
for m in range(Natoms):
x[m].coord[0] = x_np[m][0]
x[m].coord[1] = x_np[m][1]
x[m].coord[2] = x_np[m][2]
# Get neighbourcells and convert to Point-struct
cdef int Ncells = self.Ncells
cdef list cell_displacements_old = self.cell_displacements
cdef Point *cell_displacements
cell_displacements = <Point*>mem.alloc(Ncells, sizeof(Point))
for m in range(Ncells):
cell_displacements[m].coord[0] = cell_displacements_old[m][0]
cell_displacements[m].coord[1] = cell_displacements_old[m][1]
cell_displacements[m].coord[2] = cell_displacements_old[m][2]
# Convert 2body bondtype list into c-array
cdef int Ntypes = self.Ntypes
cdef list blmin_old = self.blmin
cdef double *blmin
blmin = <double*>mem.alloc(Ntypes*Ntypes, sizeof(double))
for m in range(Ntypes*Ntypes):
blmin[m] = blmin_old[m]
# Get converted atom Ntypes
cdef list num_converted_old = convert_atom_types(a.get_atomic_numbers())
cdef int *num_converted
num_converted = <int*>mem.alloc(Natoms, sizeof(int))
for m in range(Natoms):
num_converted[m] = num_converted_old[m]
#cdef double rmin = 0.7 * self.cov_dist
#cdef double radd = 1 - rmin
# Initialize Force object
cdef double *dE
dE = <double*>mem.alloc(Natoms * dim, sizeof(double))
cdef int i, j, k
for i in range(Natoms):
xi = x[i]
for cell_index in range(Ncells):
displacement = cell_displacements[cell_index]
for j in range(Natoms):
rmin = blmin[num_converted[i]*Ntypes + num_converted[j]]
radd = 1 - rmin
xj = add(x[j], displacement)
r = euclidean(xi,xj)
if r > 1e-5 and r < rcut:
r_scaled = r + radd
rijVec = subtract(xi,xj)
for k in range(dim):
dE[dim*i + k] += -12*rijVec.coord[k] / (pow(r_scaled,13)*r)
dE[dim*j + k] += 12*rijVec.coord[k] / (pow(r_scaled,13)*r)
dE_np = np.zeros((Natoms, dim))
for i in range(Natoms):
for k in range(dim):
dE_np[i,k] = dE[dim*i + k]
return - dE_np.reshape(-1)/2
def __get_neighbour_cells_displacement(self, pbc, cell):
# Calculate neighbour cells
cell_vec_norms = np.linalg.norm(cell, axis=0)
neighbours = []
for i in range(3):
if pbc[i]:
ncellmax = int(np.ceil(abs(self.rcut/cell_vec_norms[i]))) + 1 # +1 because atoms can be outside unitcell.
neighbours.append(range(-ncellmax,ncellmax+1))
else:
neighbours.append([0])
neighbourcells_disp = []
for x,y,z in product(*neighbours):
xyz = (x,y,z)
displacement = np.dot(xyz, cell)
neighbourcells_disp.append(list(displacement))
return neighbourcells_disp
import numpy as np
from scipy.spatial.distance import euclidean
from itertools import product
from ase.ga.utilities import closest_distances_generator
class delta():
def __init__(self, atoms, rcut=5, ratio_of_covalent_radii=0.7):
self.rcut = rcut
self.pbc = atoms.get_pbc()
self.cell = atoms.get_cell()
self.cell_displacements = self.__get_neighbour_cells_displacement(self.pbc, self.cell)
self.Ncells = len(self.cell_displacements)
num = atoms.get_atomic_numbers()
self.blmin = closest_distances_generator(num, ratio_of_covalent_radii=ratio_of_covalent_radii)
def energy(self, a):
num = a.get_atomic_numbers()
x = a.get_positions()
E = 0
for i, xi in enumerate(x):
for cell_index in range(self.Ncells):
displacement = self.cell_displacements[cell_index]
for j, xj in enumerate(x):
key = (num[i], num[j])
rmin = self.blmin[key]
radd = 1 - rmin
r = euclidean(xi, xj + displacement)
if r > 1e-5 and r < self.rcut:
r_scaled = r + radd
E += 1/r_scaled**12
# Devide by two decause every pair is counted twice
return E/2
def forces(self, a):
num = a.get_atomic_numbers()
x = a.get_positions()
Natoms, dim = x.shape
dE = np.zeros((Natoms, dim))
for i, xi in enumerate(x):
for cell_index in range(self.Ncells):
displacement = self.cell_displacements[cell_index]
for j, xj in enumerate(x):
key = (num[i], num[j])
rmin = self.blmin[key]
radd = 1 - rmin
xj_pbc = xj+displacement
r = euclidean(xi,xj_pbc)
if r > 1e-5 and r < self.rcut:
r_scaled = r + radd
rijVec = xi-xj_pbc
dE[i] += 12*rijVec*(-1 / (r_scaled**13*r))
dE[j] += -12*rijVec*(-1 / (r_scaled**13*r))
# Devide by two decause every pair is counted twice
return - dE.reshape(-1) / 2
def __get_neighbour_cells_displacement(self, pbc, cell):
# Calculate neighbour cells
cell_vec_norms = np.linalg.norm(cell, axis=0)
neighbours = []
for i in range(3):
if pbc[i]:
ncellmax = int(np.ceil(abs(self.rcut/cell_vec_norms[i]))) + 1 # +1 because atoms can be outside unitcell.
neighbours.append(range(-ncellmax,ncellmax+1))
else:
neighbours.append([0])
neighbourcells_disp = []
for x,y,z in product(*neighbours):
xyz = (x,y,z)
displacement = np.dot(xyz, cell)
neighbourcells_disp.append(list(displacement))
return neighbourcells_disp
from distutils.core import setup
from Cython.Build import cythonize
import numpy
setup(ext_modules=cythonize('delta.pyx'), include_dirs=[numpy.get_include()])
import numpy as np
import matplotlib.pyplot as plt
from angular_fingerprintFeature import Angular_Fingerprint
from featureCalculators_multi.angular_fingerprintFeature_cy import Angular_Fingerprint as Angular_Fingerprint_cy
from custom_calculators import doubleLJ_calculator
from ase import Atoms
from ase.visualize import view
from ase.io import read, write
from time import time
dim = 3
L = 1
d = 1
pbc = [0,0,0]
"""
N = 2
x = np.array([0.7*L, 0.9*L, d/2,
0.7*L, 0.5*L, d/2])
positions = x.reshape((-1, dim))
atomtypes = ['He', 'He']
a = Atoms(atomtypes,
positions=positions,
cell=[L,L,d],
pbc=pbc)
N = 3
x = np.array([1.5*L, 0.2*L, d/2,
0.5*L, 0.9*L, d/2,
-0.5*L, 0.5*L, d/2,])
positions = x.reshape((-1,dim))
atomtypes = ['He', 'H', 'H']
a = Atoms(atomtypes,
positions=positions,
cell=[L,L,d],
pbc=pbc)
x = np.array([1, 0, 0, 2, 0, 0, 3, 0, 0, 1.5, 1, 0])
positions = x.reshape((-1,dim))
a = Atoms('H4',
positions=positions,
cell=[4,2,1],
pbc=[0, 0, 0])
N = 4
x = np.array([0.2*L, 0.7*L, d/2,
0.3*L, 0.2*L, d/2,
0.7*L, 0.9*L, d/2,
0.7*L, 0.5*L, d/2])
positions = x.reshape((-1,dim))
atomtypes = ['H', 'H', 'He', 'He']
a = Atoms(atomtypes,
positions=positions,
cell=[L,L,d],
pbc=pbc)
N = 5
x = np.array([0.2*L, 0.7*L, d/2,
0.3*L, 0.2*L, d/2,
0.7*L, 0.9*L, d/2,
0.7*L, 0.5*L, d/2,
0.9*L, 0.1*L, d/2])
positions = x.reshape((-1,dim))
atomtypes = ['H', 'He', 'O', 'H', 'H']
a = Atoms(atomtypes,
positions=positions,
cell=[L,L,d],
pbc=pbc)
"""
atoms = read('graphene_all2.traj', index=':')
a = atoms[100]
atomtypes = a.get_atomic_numbers()
N = len(a.get_atomic_numbers())
x = a.get_positions().reshape(-1)
"""
calc = doubleLJ_calculator()
a.set_calculator(calc)
print('1:', a.get_potential_energy())
print(a.get_positions())
print('2:', a.get_potential_energy())
print(a.get_positions())
#a.set_positions(a.get_scaled_positions())
a.wrap()
print('3:', a.get_potential_energy())
print(a.get_positions())
"""
#view(a)
from delta import delta as delta_cy
from delta_py import delta as delta_py
from ase.data import covalent_radii
from ase.ga.utilities import closest_distances_generator
dcy = delta_cy(atoms=a)
dpy = delta_py(atoms=a)
print(dcy.energy(a))
print(dpy.energy(a))
print(dcy.forces(a))
print(dpy.forces(a))
"""
print('pbc check:')
print('before wrap')
print(dpy.energy(a))
a.wrap()
print('after wrap')
print(dpy.energy(a))
print(dcy.energy(a))
print('pbc check:')
print('before wrap')
print(dpy.forces(a))
a.wrap()
print('after wrap')
print(dpy.forces(a))
print(dcy.forces(a))
"""
import numpy as np
import matplotlib.pyplot as plt
from angular_fingerprintFeature import Angular_Fingerprint
from featureCalculators_multi.angular_fingerprintFeature_cy import Angular_Fingerprint as Angular_Fingerprint_cy
from custom_calculators import doubleLJ_calculator
from ase import Atoms
from ase.visualize import view
from ase.io import read, write
from time import time
dim = 3
L = 2
d = 1
pbc = [1,1,0]
"""
N = 2
x = np.array([0.7*L, 0.9*L, d/2,
0.7*L, 0.5*L, d/2])
positions = x.reshape((-1, dim))
atomtypes = ['He', 'He']
a = Atoms(atomtypes,
positions=positions,
cell=[L,L,d],
pbc=pbc)
N = 3
x = np.array([1.5*L, 0.2*L, d/2,
0.5*L, 0.9*L, d/2,
-0.5*L, 0.5*L, d/2,])
positions = x.reshape((-1,dim))
atomtypes = ['He', 'H', 'H']
a = Atoms(atomtypes,
positions=positions,
cell=[L,L,d],
pbc=pbc)
x = np.array([1, 0, 0, 2, 0, 0, 3, 0, 0, 1.5, 1, 0])
positions = x.reshape((-1,dim))
a = Atoms('H4',
positions=positions,
cell=[4,2,1],
pbc=[0, 0, 0])
N = 4
x = np.array([0.2*L, 0.7*L, d/2,
0.3*L, 0.2*L, d/2,
0.7*L, 0.9*L, d/2,
0.7*L, 0.5*L, d/2])
positions = x.reshape((-1,dim))
atomtypes = ['H', 'H', 'He', 'He']
a = Atoms(atomtypes,
positions=positions,
cell=[L,L,d],
pbc=pbc)
"""
N = 5
x = np.array([0.2*L, 0.7*L, d/2,
0.3*L, 0.2*L, d/2,
0.7*L, 0.9*L, d/2,
0.7*L, 0.5*L, d/2,
0.9*L, 0.1*L, d/2])
positions = x.reshape((-1,dim))
atomtypes = ['H', 'He', 'O', 'H', 'H']
a = Atoms(atomtypes,
positions=positions,
cell=[L,L,d],
pbc=pbc)
"""
atoms = read('graphene_data/graphene_all2.traj', index=':')
a = atoms[100]
atomtypes = a.get_atomic_numbers()
N = len(a.get_atomic_numbers())
x = a.get_positions().reshape(-1)
"""
#view(a)
Rc1 = 4
binwidth1 = 0.1
sigma1 = 0.2
Rc2 = 3
Nbins2 = 50
sigma2 = 0.2
eta = 30
gamma = 2
use_angular = True
featureCalculator = Angular_Fingerprint(a, Rc1=Rc1, Rc2=Rc2, binwidth1=binwidth1, Nbins2=Nbins2, sigma1=sigma1, sigma2=sigma2,
eta=eta, gamma=gamma, use_angular=use_angular)
t0 = time()
fingerprint = featureCalculator.get_feature(a)
fingerprint_grad = featureCalculator.get_featureGradient(a)
runtime = time() - t0
featureCalculator_cy = Angular_Fingerprint_cy(a, Rc1=Rc1, Rc2=Rc2, binwidth1=binwidth1, Nbins2=Nbins2, sigma1=sigma1, sigma2=sigma2,
eta=eta, gamma=gamma, use_angular=use_angular)
t0_cy = time()
fingerprint_cy = featureCalculator_cy.get_feature(a)
fingerprint_grad_cy = featureCalculator_cy.get_featureGradient(a)
runtime_cy = time() - t0_cy
print(featureCalculator_cy.Nelements)
print('runtime python:', runtime)
print('runtime cython:', runtime_cy)
print('times faster:', runtime/runtime_cy)
Nelements = featureCalculator.Nelements
Nbins1 = int(np.ceil(Rc1/binwidth1))
if use_angular:
Nbins = Nbins1 + Nbins2
else:
Nbins = Nbins1
Nbondtypes = int(Nelements/Nbins)
r = np.linspace(0,Rc1*Nbondtypes, Nelements)
print('py:', fingerprint.shape)
print('cy:', fingerprint_cy.shape)
plt.figure(1)
plt.plot(r, fingerprint)
plt.plot(r, fingerprint_cy, linestyle=':', color='k')
plt.figure(2)
plt.plot(r, fingerprint_grad.T)
plt.plot(r, fingerprint_grad_cy.T, linestyle=':', color='k')
plt.show()
import numpy as np
from ase import Atoms