Commit 287c707f authored by Carl Poelking's avatar Carl Poelking
Browse files

Beaker mods and abandoning non-nearest-img convention.

parent 4bf82771
#ifndef _SOAP_BOUNDARY_HPP
#define _SOAP_BOUNDARY_HPP
#include <cmath>
#include <boost/serialization/base_object.hpp>
#include <boost/serialization/export.hpp>
......@@ -41,6 +43,10 @@ public:
return r_j - r_i;
}
virtual std::vector<int> calculateRepetitions(double cutoff) {
std::vector<int> na_nb_nc = { 0, 0, 0 };
}
virtual eBoxType getBoxType() { return _type; }
template<class Archive>
......@@ -108,12 +114,24 @@ public:
BoundaryTriclinic(const matrix &box) {
_type = Boundary::typeTriclinic;
_box = box;
// Set-up inverse box
vec a = _box.getCol(0);
vec b = _box.getCol(1);
vec c = _box.getCol(2);
double V = this->BoxVolume();
vec a_inv = b ^ c / V;
vec b_inv = c ^ a / V;
vec c_inv = a ^ b / V;
_inv_box = matrix(a_inv, b_inv, c_inv);
}
BoundaryTriclinic() {
_type = Boundary::typeTriclinic;
_box.UnitMatrix();
}
vec connect(const vec &r_i, const vec &r_j) const {
/*
// This only works if a = (*,0,0), b = (*,*,0), c = (*,*,*) => e.g., GROMACS
vec r_tp, r_dp, r_sp, r_ij;
vec a = _box.getCol(0); vec b = _box.getCol(1); vec c = _box.getCol(2);
r_tp = r_j - r_i;
......@@ -121,12 +139,43 @@ public:
r_sp = r_dp - b*round(r_dp.getY()/b.getY());
r_ij = r_sp - a*round(r_sp.getX()/a.getX());
return r_ij;
*/
vec a = _box.getCol(0);
vec b = _box.getCol(1);
vec c = _box.getCol(2);
vec u = _inv_box.getCol(0);
vec v = _inv_box.getCol(1);
vec w = _inv_box.getCol(2);
vec dr = r_j - r_i;
dr = dr - std::floor(u*dr)*a;
dr = dr - std::floor(v*dr)*b;
dr = dr - std::floor(w*dr)*c;
vec dr_min = dr;
double d_min = soap::linalg::abs(dr);
for (int i=0; i < 2; ++i) {
for (int j=0; j < 2; ++j) {
for (int k=0; k < 2; ++k) {
vec dr_ijk = dr - i*a - j*b - k*c;
double d_ijk = soap::linalg::abs(dr_ijk);
if (d_ijk < d_min) {
d_min = d_ijk;
dr_min = dr_ijk;
}
}}}
return dr_min;
}
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
arch & boost::serialization::base_object<Boundary>(*this);
}
private:
matrix _inv_box;
};
}
......
......@@ -212,13 +212,31 @@ class KernelAdaptorSpecific(object):
return X, X_norm
class KernelAdaptorGlobalSpecific(object):
def __init__(self, options, *args, **kwargs):
if kwargs != {}: raise ValueError
if args != (): raise ValueError
def __init__(self, options, types_global=None):
self.types = types_global
self.S = len(types_global)
return
def getListAtomic(self, spectrum):
return [ spectrum.getGlobal() ]
def adapt(self, spectrum):
# EXTRACT A SET OF CENTER-BASED POWER EXPANSIONS
# Here: only global
IX = np.zeros((0,0), dtype='complex128')
dimX = -1
atomic_global = spectrum.getGlobal()
Xi_unnorm, Xi_norm = self.adaptScalar(atomic_global)
dimX = Xi_norm.shape[0]
IX = np.copy(Xi_norm)
IX.resize((1,dimX))
return IX
def adaptScalar(self, atomic):
xnklab_atomic = Xnklab(atomic, self.types)
X = xnklab_atomic.reduce()
x_norm = X/np.dot(X,X)**0.5
return X, X_norm
class KernelAdaptorGeneric(object):
def __init__(self, options):
def __init__(self, options, types_global=None):
return
def getListAtomic(self, spectrum):
return spectrum
......@@ -282,7 +300,7 @@ class KernelAdaptorGeneric(object):
return dX_dx, dX_dy, dX_dz
class KernelAdaptorGlobalGeneric(object):
def __init__(self, options):
def __init__(self, options, types_global=None):
return
def getListAtomic(self, spectrum):
return [ spectrum.getGlobal() ]
......@@ -482,7 +500,8 @@ KernelAdaptorFactory = {
'generic': KernelAdaptorGeneric,
'specific': KernelAdaptorSpecific,
'specific-unique': KernelAdaptorSpecificUnique,
'global-generic': KernelAdaptorGlobalGeneric
'global-generic': KernelAdaptorGlobalGeneric,
'global-specific': KernelAdaptorGlobalSpecific
}
KernelFunctionFactory = {
......
......@@ -3,12 +3,14 @@ import numpy.linalg as la
import kernel as kern
import soap
import multiprocessing as mp
import momo
import datasets.gdb
import datasets.soap
import libmatch.environments # Alchemy
import libmatch.structures # structure
try:
import datasets.gdb
import datasets.soap
import libmatch.environments # Alchemy
import libmatch.structures # structure
except ImportError:
pass
LAGRAPH_DEBUG = False
LAGRAPH_CHECK = True
......
......@@ -92,10 +92,24 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_ar
GLOG() << "Compute atomic spectrum for particle " << center->getId()
<< " (type " << center->getType() << ", targets " << targets.size() << ") ..." << std::endl;
vec box_a = _structure->getBoundary()->getBox().getCol(0);
vec box_b = _structure->getBoundary()->getBox().getCol(1);
vec box_c = _structure->getBoundary()->getBox().getCol(2);
double rc = _basis->getCutoff()->getCutoff();
int na_max = int(1 + rc/box_a.getX() - 0.5);
int nb_max = int(1 + rc/box_b.getY() - 0.5);
int nc_max = int(1 + rc/box_c.getZ() - 0.5);
GLOG() << box_a << " " << box_b << " " << box_c << std::endl;
GLOG() << rc << std::endl;
GLOG() << na_max << " " << nb_max << " " << nc_max << std::endl;
AtomicSpectrum *atomic_spectrum = new AtomicSpectrum(center, this->_basis);
Structure::particle_it_t pit;
for (pit = targets.begin(); pit != targets.end(); ++pit) {
for (pit = targets.begin(); pit != targets.end(); ++pit) { // TODO Consider images
// CHECK FOR EXCLUSIONS
if (_options->doExcludeTarget((*pit)->getType()) ||
......@@ -108,7 +122,7 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_ar
vec d = (r > 0.) ? dr/r : vec(0.,0.,1.);
// APPLY CUTOFF (= WEIGHT REDUCTION)
bool is_center = (*pit == center);
bool is_center = (*pit == center); // TODO Consider images
double weight0 = (*pit)->getWeight();
double weight_scale = _basis->getCutoff()->calculateWeight(r);
if (is_center) {
......@@ -119,7 +133,7 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_ar
bool gradients = (is_center) ? false : _options->get<bool>("spectrum.gradients");
BasisExpansion *nb_expansion = new BasisExpansion(this->_basis); // <- kept by AtomicSpectrum
nb_expansion->computeCoefficients(r, d, weight0, weight_scale, (*pit)->getSigma(), gradients);
atomic_spectrum->addQnlmNeighbour(*pit, nb_expansion);
atomic_spectrum->addQnlmNeighbour(*pit, nb_expansion); // TODO Consider images
}
return atomic_spectrum;
......
......@@ -168,6 +168,7 @@ public:
Particle &addParticle(Segment &seg);
// BOUNDARY CREATION & INTERFACE
Boundary *getBoundary() { return _box; }
void setBoundary(const matrix &box);
vec connect(const vec &r1, const vec &r2) { return _box->connect(r1, r2); /* 1->2 */ }
void setBoundaryNumeric(const boost::python::numeric::array &m);
......
......@@ -132,7 +132,7 @@ class AseConfig(object):
config_idx=None,
frame_idx=None,
config_file=None,
datastring=None):
datastring=''):
# COORDINATES
self.atoms = ase_config
# BOOK-KEEPING
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment