Skip to content

Instantly share code, notes, and snippets.

@hlscalon
Created August 25, 2018 11:41
Show Gist options
  • Save hlscalon/ccd0ce084532973d25c8835baffa19a7 to your computer and use it in GitHub Desktop.
Save hlscalon/ccd0ce084532973d25c8835baffa19a7 to your computer and use it in GitHub Desktop.
Nearest neighbor search not cubic unit cell - ase, aboria
48
Lattice="5.105310960166873 0.0 0.0 0.0 8.842657971447274 0.0 0.0 0.0 12.505406830647296" Properties=species:S:1:pos:R:3:Z:I:1 pbc="T T F"
Cu 0.00000000 0.00000000 0.00000000 29
Cu 0.00000000 1.47377633 2.08423447 29
Cu 0.00000000 2.94755266 4.16846894 29
Cu 1.27632774 0.73688816 4.16846894 29
Cu 1.27632774 2.21066449 -0.00000000 29
Cu 1.27632774 3.68444082 2.08423447 29
Cu 2.55265548 -0.00000000 -0.00000000 29
Cu 2.55265548 1.47377633 2.08423447 29
Cu 2.55265548 2.94755266 4.16846894 29
Cu 3.82898322 0.73688816 4.16846894 29
Cu 3.82898322 2.21066449 -0.00000000 29
Cu 3.82898322 3.68444082 2.08423447 29
Cu 0.00000000 4.42132899 0.00000000 29
Cu 0.00000000 5.89510531 2.08423447 29
Cu 0.00000000 7.36888164 4.16846894 29
Cu 1.27632774 5.15821715 4.16846894 29
Cu 1.27632774 6.63199348 -0.00000000 29
Cu 1.27632774 8.10576981 2.08423447 29
Cu 2.55265548 4.42132899 -0.00000000 29
Cu 2.55265548 5.89510531 2.08423447 29
Cu 2.55265548 7.36888164 4.16846894 29
Cu 3.82898322 5.15821715 4.16846894 29
Cu 3.82898322 6.63199348 -0.00000000 29
Cu 3.82898322 8.10576981 2.08423447 29
Cu 0.00000000 0.00000000 6.25270342 29
Cu 0.00000000 1.47377633 8.33693789 29
Cu 0.00000000 2.94755266 10.42117236 29
Cu 1.27632774 0.73688816 10.42117236 29
Cu 1.27632774 2.21066449 6.25270342 29
Cu 1.27632774 3.68444082 8.33693789 29
Cu 2.55265548 0.00000000 6.25270342 29
Cu 2.55265548 1.47377633 8.33693789 29
Cu 2.55265548 2.94755266 10.42117236 29
Cu 3.82898322 0.73688816 10.42117236 29
Cu 3.82898322 2.21066449 6.25270342 29
Cu 3.82898322 3.68444082 8.33693789 29
Cu 0.00000000 4.42132899 6.25270342 29
Cu 0.00000000 5.89510531 8.33693789 29
Cu 0.00000000 7.36888164 10.42117236 29
Cu 1.27632774 5.15821715 10.42117236 29
Cu 1.27632774 6.63199348 6.25270342 29
Cu 1.27632774 8.10576981 8.33693789 29
Cu 2.55265548 4.42132899 6.25270342 29
Cu 2.55265548 5.89510531 8.33693789 29
Cu 2.55265548 7.36888164 10.42117236 29
Cu 3.82898322 5.15821715 10.42117236 29
Cu 3.82898322 6.63199348 6.25270342 29
Cu 3.82898322 8.10576981 8.33693789 29
36
Lattice="8.850000000000001 0.0 0.0 -4.424999999999998 7.664324823492283 0.0 0.0 0.0 9.369200000000001" Properties=species:S:1:pos:R:3:Z:I:1 pbc="T T F"
Ti 0.00000000 0.00000000 0.00000000 22
Ti 0.00000000 1.70318329 2.34230000 22
Ti 2.95000000 0.00000000 0.00000000 22
Ti 2.95000000 1.70318329 2.34230000 22
Ti 5.90000000 0.00000000 0.00000000 22
Ti 5.90000000 1.70318329 2.34230000 22
Ti -1.47500000 2.55477494 0.00000000 22
Ti -1.47500000 4.25795824 2.34230000 22
Ti 1.47500000 2.55477494 0.00000000 22
Ti 1.47500000 4.25795824 2.34230000 22
Ti 4.42500000 2.55477494 0.00000000 22
Ti 4.42500000 4.25795824 2.34230000 22
Ti -2.95000000 5.10954988 0.00000000 22
Ti -2.95000000 6.81273318 2.34230000 22
Ti 0.00000000 5.10954988 0.00000000 22
Ti 0.00000000 6.81273318 2.34230000 22
Ti 2.95000000 5.10954988 0.00000000 22
Ti 2.95000000 6.81273318 2.34230000 22
Ti 0.00000000 0.00000000 4.68460000 22
Ti 0.00000000 1.70318329 7.02690000 22
Ti 2.95000000 0.00000000 4.68460000 22
Ti 2.95000000 1.70318329 7.02690000 22
Ti 5.90000000 0.00000000 4.68460000 22
Ti 5.90000000 1.70318329 7.02690000 22
Ti -1.47500000 2.55477494 4.68460000 22
Ti -1.47500000 4.25795824 7.02690000 22
Ti 1.47500000 2.55477494 4.68460000 22
Ti 1.47500000 4.25795824 7.02690000 22
Ti 4.42500000 2.55477494 4.68460000 22
Ti 4.42500000 4.25795824 7.02690000 22
Ti -2.95000000 5.10954988 4.68460000 22
Ti -2.95000000 6.81273318 7.02690000 22
Ti 0.00000000 5.10954988 4.68460000 22
Ti 0.00000000 6.81273318 7.02690000 22
Ti 2.95000000 5.10954988 4.68460000 22
Ti 2.95000000 6.81273318 7.02690000 22
#include "Aboria.h"
#include <vector>
#include <iostream>
#include <set>
#include <fstream>
using namespace Aboria;
using Pair = std::pair<double, int>;
using particle_t = Particles<std::tuple<>, 3, std::vector, Kdtree>;
using position = particle_t::position;
int main(int argc, char *argv[]) {
if (argc != 5) {
std::cerr << "1: Filename; 2: random point x; 3: random point y; 4: random point z.\n";
return -1;
}
std::ifstream file(argv[1], std::ios::binary);
std::string line;
int totalAtoms = 0;
file >> totalAtoms;
std::getline(file, line); // skip second line
std::getline(file, line); // skip second line
particle_t particles(totalAtoms);
int i = 0;
double x, y, z; std::string nop;
while (file >> nop >> x >> y >> z >> nop) {
get<position>(particles)[i] = vdouble3(x, y, z);
i++;
}
vdouble3 min = vdouble3(0.0, 0.0, 0.0);
vdouble3 max = vdouble3(5.105310960166873, 8.842657971447274, 12.505406830647296);
vbool3 periodic = vbool3(true, true, false);
particles.init_neighbour_search(min, max, periodic);
vdouble3 rp = vdouble3(std::stod(argv[2]), std::stod(argv[3]), std::stod(argv[4]));
std::set<Pair> orderedPoints;
double radius = 3;
int count = 0;
for (auto i = euclidean_search(particles.get_query(), rp, radius); i != false; ++i) {
count++;
orderedPoints.emplace(Pair{(i.dx()).norm(), get<id>(*i)});
std::cout << "Found a particle with dx = " << i.dx() << " and id = " << get<id>(*i) << "\n";
}
for (const auto & p : orderedPoints) {
std::cout << p.second << "\t" << p.first << "\n";
}
std::cout << "There are " << count << " particles.\n";
return 0;
}
import sys
from ase import Atom
from ase.io import read
from operator import itemgetter
def getAllDistancesFromPoint(slab, atomic_number, x, y, z):
slab.append(Atom(atomic_number, (x, y, z))) # append phantom atom
idxAtom = len(slab) - 1
indices = range(0, idxAtom)
allDistances = slab.get_distances(idxAtom, indices, mic=True)
slab.pop() # remove phantom
return allDistances
def printNClosestNeighborsFromPoint(slab, n, x, y, z, atomic_number):
allDistances = getAllDistancesFromPoint(slab, atomic_number, x, y, z)
distances = dict((idx, d) for idx, d in enumerate(allDistances))
nClosest = sorted(distances.items(), key=itemgetter(1))[:n]
for i in nClosest:
print i[0], "\t:\t", i[1]
def main():
filename = sys.argv[1]
x, y, z = float(sys.argv[2]), float(sys.argv[3]), float(sys.argv[4]) # random point
slab = read(filename)
atomic_number = slab.get_atomic_numbers()[0]
n = 10 # number of neighbors
printNClosestNeighborsFromPoint(slab, n, x, y, z, atomic_number)
if __name__ == "__main__":
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment