Commit 579137a5 authored by Carl Poelking's avatar Carl Poelking
Browse files

Smooth expansion, numerical integration.

parent fdfe189d
......@@ -76,6 +76,7 @@ public:
}
}
}
_coeff *= weight;
return;
}
void add(BasisExpansion &other) {
......@@ -178,8 +179,9 @@ public:
vec dr_center_particle = structure->connect(center->getPos(), (*pit)->getPos());
vec dr_particle_target = dr - dr_center_particle;
double r_particle_target = soap::linalg::abs(dr_particle_target);
double sigma = 0.5;
density_dr += pow(1./(2.*M_PI*sigma*sigma), 1.5)
double sigma = (*pit)->getSigma();
double weight = (*pit)->getWeight();
density_dr += weight*pow(1./(2.*M_PI*sigma*sigma), 1.5)
* exp(-r_particle_target*r_particle_target/(2*sigma*sigma));
}
int_density_dr += density_dr*dx*dy*dz;
......
......@@ -14,6 +14,8 @@ public:
_N_recip(8), _L_recip(6), _Rc_recip(5.),
_b1(vec(1,0,0)), _b2(vec(0,1,0)), _b3(vec(0,0,1)),
_center_w0(1.), _center_excludes(boost::python::list()) {
// SET DEFAULTS
this->set("integration_steps", 15);
}
~Options() {}
......
......@@ -19,6 +19,7 @@ namespace ub = boost::numeric::ublas;
void RadialBasis::configure(Options &options) {
_N = options.get<int>("radialbasis.N");
_Rc = options.get<double>("radialbasis.Rc");
_integration_steps = options.get<int>("radialbasis.integration_steps");
}
RadialCoefficients RadialBasis::computeCoefficients(double r) {
......@@ -172,6 +173,34 @@ struct SphericalGaussian
double _norm_g_dV;
};
struct ModifiedSphericalBessel1stKind
{
static std::vector<double> eval(int degree, double r) {
std::vector<double> il;
if (r < RADZERO) {
il.push_back(1.);
il.push_back(0.);
}
else {
il.push_back(sinh(r)/r);
il.push_back(cosh(r)/r - sinh(r)/(r*r));
}
for (int l = 2; l <= degree; ++l) {
if (r < RADZERO) {
il.push_back(0.);
}
else {
if (il[l-1] < SPHZERO) il.push_back(0.);
else il.push_back( il[l-2] - (2*(l-1)+1)/r*il[l-1] );
}
}
return il;
}
static const double RADZERO = 1e-10;
static const double SPHZERO = 1e-4;
};
void RadialBasisGaussian::computeCoefficients(double r, double particle_sigma, radcoeff_t &save_here) {
// Delta-type expansion =>
......@@ -185,29 +214,151 @@ void RadialBasisGaussian::computeCoefficients(double r, double particle_sigma, r
save_here(n, l) = gn_at_r;
}
}
save_here = ub::prod(_Tij, save_here);
}
else {
throw soap::base::NotImplemented("...");
// std::cout << "sigma > 0." << std::endl;
// int degree_sph_in = save_here.size2();
// std::vector<double> sph_il = ModifiedSphericalBessel1stKind::eval(degree_sph_in, 1e-4);
// for (int l = 0; l != save_here.size2(); ++l) {
// std::cout << "l=" << l << " il=" << sph_il[l] << std::endl;
// }
//
// throw soap::base::NotImplemented("...");
// Particle properties
double ai = 1./(2*particle_sigma*particle_sigma);
double ri = r;
SphericalGaussian gi_sph(vec(0,0,0), particle_sigma); // <- position should not matter, as only normalization used here
double norm_g_dV_sph_i = gi_sph._norm_g_dV;
basis_it_t it;
int k = 0;
basis_it_t it;
for (it = _basis.begin(), k = 0; it != _basis.end(); ++it, ++k) {
// Radial Gaussian properties
double ak = (*it)->_alpha;
double rk = (*it)->_r0;
double norm_r2_g2_dr_rad_k = (*it)->_norm_r2_g2_dr;
// Combined properties
double beta_ik = ai+ak;
double rho_ik = ak*rk/beta_ik;
double bessel_arg_ik = 2*ai*ri*rho_ik;
double sigma_ik = sqrt(0.5/beta_ik);
RadialGaussian gik_rad(rho_ik, sigma_ik);
double norm_r2_g_dr_rad_ik = gik_rad._norm_r2_g_dr;
double prefac =
4*M_PI *
norm_r2_g2_dr_rad_k*norm_g_dV_sph_i/norm_r2_g_dr_rad_ik *
exp(-ai*ri*ri) *
exp(-ak*rk*rk*(1-ak/beta_ik));
// // DELTA APPROXIMATION
// std::cout << "r " << r;
// std::cout << " beta_ik " << beta_ik;
// std::cout << " rho_ik " << rho_ik;
// std::cout << std::endl;
// double bessel_arg_ik = 2*ai*ri*rho_ik;
// int degree_sph = save_here.size2(); // <- L+1
// std::vector<double> sph_il = ModifiedSphericalBessel1stKind::eval(degree_sph, bessel_arg_ik);
// for (int l = 0; l != degree_sph; ++l) {
// save_here(k, l) = prefac*sph_il[l];
// }
// std::cout << std::endl;
// std::cout << "DELTA" << std::endl;
// std::cout << "k = " << k;
// for (int l = 0; l != save_here.size2(); ++l) {
// std::cout << " " << save_here(k, l);
// }
// std::cout << std::endl;
// NUMERICAL INTEGRATION
// ZERO COEFFS
for (int l = 0; l != save_here.size2(); ++l) {
save_here(k, l) = 0.0;
}
// // INTEGRATE
// double delta_r = 0.01;
// double r_min = 0.0;
// double r_max = 10.;
// int steps = int((r_max-r_min)/delta_r)+1;
// std::cout << "r " << r;
// std::cout << " steps " << steps;
// std::cout << " beta_ik " << beta_ik;
// std::cout << " rho_ik " << rho_ik;
// std::cout << std::endl;
//
// for (int s = 0; s <= steps; ++s) {
// // Current radius
// double r_step = r_min + s*delta_r;
// // Calculate ModSphBessels
// double arg = 2*ai*ri*r_step;
// std::vector<double> sph_il = ModifiedSphericalBessel1stKind::eval(save_here.size2(), arg);
// // Add increments
// for (int l = 0; l != save_here.size2(); ++l) {
// save_here(k,l) +=
// prefac*
// delta_r*r_step*r_step*
// sph_il[l]*
// exp(-beta_ik*(r_step-rho_ik)*(r_step-rho_ik))*norm_r2_g_dr_rad_ik;
// }
// }
// SIMPSON'S RULE
double r_min = rho_ik - 4*sigma_ik;
double r_max = rho_ik + 4*sigma_ik;
if (r_min < 0.) {
r_max -= r_min;
r_min = 0.;
}
int n_steps = this->_integration_steps;
double delta_r_step = (r_max-r_min)/(n_steps-1);
int n_sample = 2*n_steps+1;
double delta_r_sample = 0.5*delta_r_step;
// Compute samples for all l's
// For each l, store integrand at r_sample
ub::matrix<double> integrand_l_at_r = ub::zero_matrix<double>(save_here.size2(), n_sample);
for (int s = 0; s < n_sample; ++s) {
// f0 f1 f2 f3 .... f-3 f-2 f-1
// |-----||----||----||-------|
double r_sample = r_min - delta_r_sample + s*delta_r_sample;
// ... Generate Bessels
std::vector<double> sph_il = ModifiedSphericalBessel1stKind::eval(save_here.size2(), 2*ai*ri*r_sample);
// ... Compute & store integrands
for (int l = 0; l != save_here.size2(); ++l) {
integrand_l_at_r(l, s) =
prefac*
r_sample*r_sample*
sph_il[l]*
exp(-beta_ik*(r_sample-rho_ik)*(r_sample-rho_ik))*norm_r2_g_dr_rad_ik;
}
}
// Apply Simpson's rule
for (int s = 0; s < n_steps; ++s) {
for (int l = 0; l != save_here.size2(); ++l) {
save_here(k,l) +=
delta_r_step/6.*(
integrand_l_at_r(l, 2*s)+
4*integrand_l_at_r(l, 2*s+1)+
integrand_l_at_r(l, 2*s+2)
);
}
}
// std::cout << "NUMERICAL" << std::endl;
// std::cout << "k = " << k;
// for (int l = 0; l != save_here.size2(); ++l) {
// std::cout << " " << save_here(k, l);
// }
// std::cout << std::endl;
}
save_here = ub::prod(_Tij, save_here);
}
return;
}
......
......@@ -50,6 +50,7 @@ protected:
std::string _type;
int _N;
double _Rc;
int _integration_steps;
static const double RADZERO = 1e-10;
};
......
......@@ -52,13 +52,12 @@ void Spectrum::computeAtomic(Particle *center) {
Structure::particle_it_t pit;
for (pit = _structure->beginParticles(); pit != _structure->endParticles(); ++pit) {
// TODO Cut-off check
// TODO Cut-off check && cut-off smoothing (= weight reduction)
vec dr = _structure->connect(center->getPos(), (*pit)->getPos());
double r = soap::linalg::abs(dr);
vec d = dr/r;
std::string type_other = (*pit)->getType();
// TODO Account for weight
BasisExpansion nb_expansion(this->_basis);
nb_expansion.computeCoefficients(r, d, (*pit)->getWeight(), (*pit)->getSigma());
nbhood_expansion->add(nb_expansion);
......
......@@ -6,6 +6,19 @@ import os
import numpy as np
from momo import osio, endl, flush
element_vdw_radius = {
'C':1.70,
'N':1.55,
'O':1.52,
'H':1.20
}
element_mass = {
'C':12,
'N':14,
'O':16,
'H':1
}
ase_config_list = soap.tools.ase_load_all('configs')
config = ase_config_list[0]
osio << config.atoms << endl
......@@ -16,18 +29,21 @@ options.set('radialbasis.type', 'gaussian')
options.set('radialbasis.N', 9)
options.set('radialbasis.Rc', 4.0)
options.set('radialbasis.sigma', 0.5)
options.set('radialbasis.integration_steps', 15)
options.set('angularbasis.type', 'spherical-harmonic')
options.set('angularbasis.L', 6)
options.set('densitygrid.N', 20)
options.set('densitygrid.dx', 0.1)
options.set('densitygrid.dx', 0.10)
# LOAD STRUCTURE
structure = soap.tools.setup_structure_ase(config.config_file, config.atoms)
osio << osio.mg << structure.label << endl
for atom in structure:
atom.sigma = 0.4*element_vdw_radius[atom.type]
atom.weight = 1. #element_mass[atom.type]
osio << atom.name << atom.type << atom.weight << atom.sigma << atom.pos << endl
atom.sigma = 0.0
# COMPUTE SPECTRUM
spectrum = soap.Spectrum(structure, options)
......
Supports Markdown
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