Grid.cpp 2.41 KB
Newer Older
1
/*
2
# Copyright 2016-2018 Ruben Jesus Garcia Hernandez
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
 # You may obtain a copy of the License at
 #
 #     http://www.apache.org/licenses/LICENSE-2.0
 #
 # Unless required by applicable law or agreed to in writing, software
 # distributed under the License is distributed on an "AS IS" BASIS,
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
*/


18 19
#include "Grid.h"
#include "atoms.hpp" //for radius
20
#include <math.h>
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

grid::grid (float *m, float *M, int dims, float s):scale(s) {
	content = new std::vector<float*> [dims*dims*dims];
	for (int i=0;i<3;i++) {
		this->m[i]=m[i];
		this->M[i]=M[i];
		dif[i]=M[i]-m[i];
	}
	this->dims=dims;
	maxradius=0;
}

grid::~grid()
{
	delete [] content;
}

void grid::coordinates(const float pos[3], int c[3])
{
	for (int i=0;i<3;i++) {
41
		c[i]=static_cast<int>(floor((pos[i]-m[i])/dif[i]*dims));
42 43 44 45 46 47 48 49 50 51 52 53 54
		if (c[i]>=dims)
			c[i]=dims-1;
		else if (c[i]<0)
			c[i]=0;
	}
}

void grid::add (float *p) //compatible with the atoms xyzr
{
	int pos[3];
	coordinates (p, pos);
	
	content[pos[0]*dims*dims + pos[1]*dims+pos[2]].push_back(p);
55
	float ar=atomRadius(static_cast<int>(p[3]));
56 57 58 59 60 61 62 63
	if (ar>maxradius)
		maxradius=ar;
}

bool grid::compare (float *a, float *b)
{
	if (a<=b) //already returned when searching a beforehand
		return false;
64
	float sqd=atomRadius(static_cast<int>(a[3]))+atomRadius(static_cast<int>(b[3]));
65
	sqd*=sqd;
66 67 68
	float d=0;
	for (int i=0;i<3;i++)
		d+=(a[i]-b[i])*(a[i]-b[i]);
69
	if (d*scale < sqd)
70
		return true;
71
	return false;
72 73 74 75 76 77 78 79 80 81 82 83 84
}

std::vector<float*> grid::find (float *p) 
{
	std::vector<float*> result;
	//search a sphere centered in p, of radius (p[3]+maxradius) *scale

	//start by searching a cube
	int mc[3];
	int Mc[3];
	float mp[3];
	float Mp[3];
	for (int i=0;i<3;i++) {
85 86
		mp[i]=p[i]-(atomRadius(static_cast<int>(p[3]))+maxradius)/scale;
		Mp[i]=p[i]+(atomRadius(static_cast<int>(p[3]))+maxradius)/scale;
87 88 89
	}
	coordinates(mp, mc);
	coordinates(Mp, Mc);
90 91 92 93

	for (int x=mc[0];x<=Mc[0];x++)
		for (int y=mc[1];y<=Mc[1];y++)
			for (int z=mc[2];z<=Mc[2];z++) {
94 95 96 97 98 99 100 101 102
				const int c=x*dims*dims + y*dims+z;
				for (int i=0;i<content[c].size();i++) {
					if (compare (content[c][i], p))
						result.push_back(content[c][i]);
				}
			}
	return result;			
}