Skip to content

Commit

Permalink
Merge pull request #16 from hspark1212/develop
Browse files Browse the repository at this point in the history 
v0.1.14
  • Loading branch information
@hspark1212
hspark1212 committed May 24, 2024
2 parents a4390e4 + ac5c3a4 commit 81c8821
Show file tree
Hide file tree
Showing 4 changed files with 98 additions and 109 deletions.
59 changes: 29 additions & 30 deletions fast_grid/libs/distance_matrix.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,51 +8,50 @@ cimport numpy as np
from libc.math cimport sqrt


def minimum_image_triclinic(np.ndarray[np.float64_t, ndim=1] dx, np.ndarray[np.float64_t, ndim=2] box):
cdef int ix, iy, iz
cdef float rx, ry0, ry1, rz0, rz1, rz2, dsq
cdef float dsq_min = np.finfo(np.float64).max
cdef np.ndarray[np.float64_t, ndim=1] dx_min = np.zeros(3, dtype=np.float64)

for ix in range(-1, 2):
rx = dx[0] + box[0, 0] * ix
for iy in range(-1, 2):
ry0 = rx + box[1, 0] * iy
ry1 = dx[1] + box[1, 1] * iy
for iz in range(-1, 2):
rz0 = ry0 + box[2, 0] * iz
rz1 = ry1 + box[2, 1] * iz
rz2 = dx[2] + box[2, 2] * iz
dsq = rz0 * rz0 + rz1 * rz1 + rz2 * rz2
if dsq < dsq_min:
dsq_min = dsq
dx_min[0] = rz0
dx_min[1] = rz1
dx_min[2] = rz2

dx[:] = dx_min


def distance_matrix_triclinic_cython(np.ndarray[np.float64_t, ndim=2] pos1,
np.ndarray[np.float64_t, ndim=2] pos2,
np.ndarray[np.float64_t, ndim=2] box):
cdef int i, j
cdef int n = pos1.shape[0]
cdef int m = pos2.shape[0]
cdef double r2

cdef int ix, iy, iz
cdef float rx, ry0, ry1, rz0, rz1, rz2, dsq
cdef float dsq_min = np.finfo(np.float64).max

cdef np.ndarray[np.float64_t, ndim=2] distances = np.zeros((n, m), dtype=np.float64)
cdef np.ndarray[np.float64_t, ndim=1] diff = np.zeros(3, dtype=np.float64)
cdef np.ndarray[np.float64_t, ndim=1] dx_min = np.zeros(3, dtype=np.float64)

for i in prange(n, nogil=True):
for j in range(m):
diff[0] = pos2[j, 0] - pos1[i, 0]
diff[1] = pos2[j, 1] - pos1[i, 1]
diff[2] = pos2[j, 2] - pos1[i, 2]

with gil:
minimum_image_triclinic(diff, box)
r2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]

dsq_min = 10000.
dx_min[0] = 0
dx_min[1] = 0
dx_min[2] = 0

for ix in range(-1, 2):
rx = diff[0] + box[0, 0] * ix
for iy in range(-1, 2):
ry0 = rx + box[1, 0] * iy
ry1 = diff[1] + box[1, 1] * iy
for iz in range(-1, 2):
rz0 = ry0 + box[2, 0] * iz
rz1 = ry1 + box[2, 1] * iz
rz2 = diff[2] + box[2, 2] * iz
dsq = rz0 * rz0 + rz1 * rz1 + rz2 * rz2
if dsq < dsq_min:
dsq_min = dsq
dx_min[0] = rz0
dx_min[1] = rz1
dx_min[2] = rz2

r2 = dx_min[0] * dx_min[0] + dx_min[1] * dx_min[1] + dx_min[2] * dx_min[2]
distances[i, j] = sqrt(r2)

return distances
return distances
97 changes: 58 additions & 39 deletions fast_grid/libs/potential.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,34 +8,6 @@ cimport numpy as np
from libc.math cimport exp


def minimum_image_triclinic(np.ndarray[np.float64_t, ndim=1] dx, np.ndarray[np.float64_t, ndim=2] box):
cdef int ix, iy, iz
cdef float rx, ry0, ry1, rz0, rz1, rz2, dsq
cdef float dsq_min = np.finfo(np.float64).max
cdef np.ndarray[np.float64_t, ndim=1] dx_min = np.zeros(3, dtype=np.float64)

for ix in range(-1, 2):
rx = dx[0] + box[0, 0] * ix
for iy in range(-1, 2):
ry0 = rx + box[1, 0] * iy
ry1 = dx[1] + box[1, 1] * iy
for iz in range(-1, 2):
rz0 = ry0 + box[2, 0] * iz
rz1 = ry1 + box[2, 1] * iz
rz2 = dx[2] + box[2, 2] * iz
dsq = rz0 * rz0 + rz1 * rz1 + rz2 * rz2
if dsq < dsq_min:
dsq_min = dsq
dx_min[0] = rz0
dx_min[1] = rz1
dx_min[2] = rz2

dx[:] = dx_min


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def lj_potential_cython(np.ndarray[np.float64_t, ndim=2] pos1,
np.ndarray[np.float64_t, ndim=2] pos2,
np.ndarray[np.float64_t, ndim=2] box,
Expand All @@ -48,6 +20,11 @@ def lj_potential_cython(np.ndarray[np.float64_t, ndim=2] pos1,
cdef double r2, lj6, lj12, inv_r2, inv_r6, inv_r12, e, s, s6, s12, energy
cdef double threshold = 1e-10
cdef double cutoff_squared = cutoff * cutoff

cdef int ix, iy, iz
cdef float rx, ry0, ry1, rz0, rz1, rz2, dsq
cdef float dsq_min = np.finfo(np.float64).max
cdef np.ndarray[np.float64_t, ndim=1] dx_min = np.zeros(3, dtype=np.float64)

cdef np.ndarray[np.float64_t, ndim=1] energy_grid = np.zeros((n), dtype=np.float64)
cdef np.ndarray[np.float64_t, ndim=1] diff = np.zeros(3, dtype=np.float64)
Expand All @@ -58,11 +35,31 @@ def lj_potential_cython(np.ndarray[np.float64_t, ndim=2] pos1,
diff[0] = pos2[j, 0] - pos1[i, 0]
diff[1] = pos2[j, 1] - pos1[i, 1]
diff[2] = pos2[j, 2] - pos1[i, 2]

with gil:
minimum_image_triclinic(diff, box)
r2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]

dsq_min = 10000.
dx_min[0] = 0
dx_min[1] = 0
dx_min[2] = 0

# minimum_image_triclinic
for ix in range(-1, 2):
rx = diff[0] + box[0, 0] * ix
for iy in range(-1, 2):
ry0 = rx + box[1, 0] * iy
ry1 = diff[1] + box[1, 1] * iy
for iz in range(-1, 2):
rz0 = ry0 + box[2, 0] * iz
rz1 = ry1 + box[2, 1] * iz
rz2 = diff[2] + box[2, 2] * iz
dsq = rz0 * rz0 + rz1 * rz1 + rz2 * rz2
if dsq < dsq_min:
dsq_min = dsq
dx_min[0] = rz0
dx_min[1] = rz1
dx_min[2] = rz2

r2 = dx_min[0] * dx_min[0] + dx_min[1] * dx_min[1] + dx_min[2] * dx_min[2]

if r2 < cutoff_squared and r2 > threshold:
e = epsilon[j]
s = sigma[j]
Expand All @@ -81,9 +78,6 @@ def lj_potential_cython(np.ndarray[np.float64_t, ndim=2] pos1,
return energy_grid


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def gaussian_cython(np.ndarray[np.float64_t, ndim=2] pos1,
np.ndarray[np.float64_t, ndim=2] pos2,
np.ndarray[np.float64_t, ndim=2] box,
Expand All @@ -97,6 +91,11 @@ def gaussian_cython(np.ndarray[np.float64_t, ndim=2] pos1,
cdef double threshold = 1e-10
cdef double width_squared = width * width
cdef double cutoff_squared = cutoff * cutoff

cdef int ix, iy, iz
cdef float rx, ry0, ry1, rz0, rz1, rz2, dsq
cdef float dsq_min = np.finfo(np.float64).max
cdef np.ndarray[np.float64_t, ndim=1] dx_min = np.zeros(3, dtype=np.float64)

cdef np.ndarray[np.float64_t, ndim=1] energy_grid = np.zeros((n), dtype=np.float64)
cdef np.ndarray[np.float64_t, ndim=1] diff = np.zeros(3, dtype=np.float64)
Expand All @@ -107,10 +106,30 @@ def gaussian_cython(np.ndarray[np.float64_t, ndim=2] pos1,
diff[0] = pos2[j, 0] - pos1[i, 0]
diff[1] = pos2[j, 1] - pos1[i, 1]
diff[2] = pos2[j, 2] - pos1[i, 2]

with gil:
minimum_image_triclinic(diff, box)
r2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]

dsq_min = 10000.
dx_min[0] = 0
dx_min[1] = 0
dx_min[2] = 0

# minimum_image_triclinic
for ix in range(-1, 2):
rx = diff[0] + box[0, 0] * ix
for iy in range(-1, 2):
ry0 = rx + box[1, 0] * iy
ry1 = diff[1] + box[1, 1] * iy
for iz in range(-1, 2):
rz0 = ry0 + box[2, 0] * iz
rz1 = ry1 + box[2, 1] * iz
rz2 = diff[2] + box[2, 2] * iz
dsq = rz0 * rz0 + rz1 * rz1 + rz2 * rz2
if dsq < dsq_min:
dsq_min = dsq
dx_min[0] = rz0
dx_min[1] = rz1
dx_min[2] = rz2

r2 = dx_min[0] * dx_min[0] + dx_min[1] * dx_min[1] + dx_min[2] * dx_min[2]

if r2 < cutoff_squared and r2 > threshold:
energy += height * exp(-1 * r2 / width_squared)
Expand Down
4 changes: 2 additions & 2 deletions setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,10 @@

setup(
name="fast-grid",
version="0.1.13",
version="0.1.14",
description="Fast grid calculation",
author="Hyunsoo Park",
author_email="hpark@ic.ac.uk",
author_email="phs68660888@gmail.com",
url="https://github.com/hspark1212/fast-grid.git",
install_requires=[
"ase",
Expand Down
47 changes: 9 additions & 38 deletions tutorial.ipynb
View file

Large diffs are not rendered by default.

0 comments on commit 81c8821

Please sign in to comment.