From ac5c3a437a47efd32464f701416a38cb9c5eede9 Mon Sep 17 00:00:00 2001 From: hspark1212 Date: Fri, 24 May 2024 12:08:31 +0100 Subject: [PATCH] v0.1.14 --- fast_grid/libs/distance_matrix.pyx | 59 +++++++++--------- fast_grid/libs/potential.pyx | 97 ++++++++++++++++++------------ setup.py | 4 +- tutorial.ipynb | 47 +++------------ 4 files changed, 98 insertions(+), 109 deletions(-) diff --git a/fast_grid/libs/distance_matrix.pyx b/fast_grid/libs/distance_matrix.pyx index c9b1192..10d0947 100644 --- a/fast_grid/libs/distance_matrix.pyx +++ b/fast_grid/libs/distance_matrix.pyx @@ -8,31 +8,6 @@ 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): @@ -40,19 +15,43 @@ def distance_matrix_triclinic_cython(np.ndarray[np.float64_t, ndim=2] pos1, 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 \ No newline at end of file diff --git a/fast_grid/libs/potential.pyx b/fast_grid/libs/potential.pyx index 13e1aa7..89cc3be 100644 --- a/fast_grid/libs/potential.pyx +++ b/fast_grid/libs/potential.pyx @@ -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, @@ -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) @@ -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] @@ -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, @@ -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) @@ -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) diff --git a/setup.py b/setup.py index 59b610e..e4a94e0 100644 --- a/setup.py +++ b/setup.py @@ -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", diff --git a/tutorial.ipynb b/tutorial.ipynb index ef01b49..6cb1611 100644 --- a/tutorial.ipynb +++ b/tutorial.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 3, "metadata": {}, "outputs": [ { @@ -12,35 +12,6 @@ "Visualizing energy grid | supercell [1 1 1]...\n" ] }, - { - "data": { - "text/html": [ - " \n", - " " - ] - }, - "metadata": {}, - "output_type": "display_data" - }, { "data": { "application/vnd.plotly.v1+json": { @@ -110268,9 +110239,9 @@ } }, "text/html": [ - "