Commit 167d9294 authored by Carl Poelking's avatar Carl Poelking
Browse files

Option for 2l+1 norm: on/off.

parent 67817423
......@@ -28,6 +28,7 @@ public:
RadialBasis *getRadBasis() { return _radbasis; }
AngularBasis *getAngBasis() { return _angbasis; }
CutoffFunction *getCutoff() { return _cutoff; }
Options *getOptions() { return _options; }
const int &N() { return _radbasis->N(); }
const int &L() { return _angbasis->L(); }
......
......@@ -8,6 +8,7 @@ Options::Options() :
// Set defaults
this->set("spectrum.gradients", false);
this->set("spectrum.2l1_norm", true);
this->set("radialbasis.type", "gaussian");
this->set("radialbasis.mode", "equispaced");
this->set("radialbasis.N", 9);
......
......@@ -54,11 +54,16 @@ public:
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
arch & _key_value_map;
arch & _exclude_center;
arch & _exclude_target;
arch & _exclude_center_list;
arch & _exclude_target_list;
arch & _exclude_center_id;
arch & _exclude_target_id;
arch & _exclude_center_id_list;
arch & _exclude_target_id_list;
return;
}
......
......@@ -13,6 +13,9 @@ PowerExpansion::PowerExpansion(Basis *basis) :
_has_scalars = true;
_coeff = coeff_zero_t(_N*_N, _L+1);
_with_sqrt_2l_1_norm = _basis->getOptions()->get<bool>("spectrum.2l1_norm"); // With/without normalization sqrt(8\pi^2/(2l+1))
}
void PowerExpansion::computeCoefficientsGradients(BasisExpansion *basex1, BasisExpansion *basex2, bool same_types) {
......@@ -83,9 +86,10 @@ void PowerExpansion::computeCoefficientsGradients(BasisExpansion *basex1, BasisE
c_nkl_dz += dqnlm_dz(n,lm)*std::conj(qnlm(k,lm))
+ qnlm(n,lm)*std::conj(dqnlm_dz(k,lm));
}
_coeff_grad_x(n*_N+k, l) = 2.*sqrt(2.)*M_PI/sqrt(2.*l+1)*c_nkl_dx; // Normalization = sqrt(8\pi^2/(2l+1))
_coeff_grad_y(n*_N+k, l) = 2.*sqrt(2.)*M_PI/sqrt(2.*l+1)*c_nkl_dy;
_coeff_grad_z(n*_N+k, l) = 2.*sqrt(2.)*M_PI/sqrt(2.*l+1)*c_nkl_dz;
double prefac_2l_1 = (_with_sqrt_2l_1_norm) ? 2.*sqrt(2.)*M_PI/sqrt(2.*l+1) : 1.; // Normalization = sqrt(8\pi^2/(2l+1))
_coeff_grad_x(n*_N+k, l) = prefac_2l_1*c_nkl_dx;
_coeff_grad_y(n*_N+k, l) = prefac_2l_1*c_nkl_dy;
_coeff_grad_z(n*_N+k, l) = prefac_2l_1*c_nkl_dz;
}
}
}
......@@ -112,9 +116,10 @@ void PowerExpansion::computeCoefficientsGradients(BasisExpansion *basex1, BasisE
c_nkl_dz += qnlm(n,lm)*std::conj(dqnlm_dz(k,lm));
}
}
_coeff_grad_x(n*_N+k, l) = 2.*sqrt(2.)*M_PI/sqrt(2.*l+1)*c_nkl_dx; // Normalization = sqrt(8\pi^2/(2l+1))
_coeff_grad_y(n*_N+k, l) = 2.*sqrt(2.)*M_PI/sqrt(2.*l+1)*c_nkl_dy;
_coeff_grad_z(n*_N+k, l) = 2.*sqrt(2.)*M_PI/sqrt(2.*l+1)*c_nkl_dz;
double prefac_2l_1 = (_with_sqrt_2l_1_norm) ? 2.*sqrt(2.)*M_PI/sqrt(2.*l+1) : 1.; // Normalization = sqrt(8\pi^2/(2l+1))
_coeff_grad_x(n*_N+k, l) = prefac_2l_1*c_nkl_dx;
_coeff_grad_y(n*_N+k, l) = prefac_2l_1*c_nkl_dy;
_coeff_grad_z(n*_N+k, l) = prefac_2l_1*c_nkl_dz;
}
}
}
......@@ -138,7 +143,8 @@ void PowerExpansion::computeCoefficients(BasisExpansion *basex1, BasisExpansion
//std::cout << m << " " << std::flush;
c_nkl += coeff1(n, l*l+l+m)*std::conj(coeff2(k, l*l+l+m));
}
_coeff(n*_N+k, l) = 2.*sqrt(2.)*M_PI/sqrt(2.*l+1)*c_nkl; // Normalization = sqrt(8\pi^2/(2l+1))
double prefac_2l_1 = (_with_sqrt_2l_1_norm) ? 2.*sqrt(2.)*M_PI/sqrt(2.*l+1) : 1.; // Normalization = sqrt(8\pi^2/(2l+1))
_coeff(n*_N+k, l) = prefac_2l_1*c_nkl;
//std::cout << std::endl;
}
}
......
......@@ -47,15 +47,25 @@ public:
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
arch & _basis;
arch & _N;
arch & _L;
arch & _with_sqrt_2l_1_norm;
arch & _has_scalars;
arch & _coeff;
arch & _has_gradients;
arch & _coeff_grad_x;
arch & _coeff_grad_y;
arch & _coeff_grad_z;
}
private:
Basis *_basis;
int _N;
int _L;
bool _with_sqrt_2l_1_norm;
bool _has_scalars;
coeff_t _coeff; // access via (N*n+k, l) with shape (N*N, L+1)
......
......@@ -134,7 +134,7 @@ class KernelPotential:
#print self.IX
logging.info("Acquired %d environments." % n_acqu)
return
def computeEnergy(self, structure):
def computeEnergy(self, structure, return_prj_mat=False):
logging.info("Start energy ...")
energy = 0.0
spectrum = soap.Spectrum(structure, self.options, self.basis)
......@@ -145,13 +145,18 @@ class KernelPotential:
n_acqu = IX_acqu.shape[0]
dim_acqu = IX_acqu.shape[1]
logging.info("Compute energy from %d atomic environments ..." % n_acqu)
projection_matrix = []
for n in range(n_acqu):
X = IX_acqu[n]
#print "X_MAG", np.dot(X,X)
ic = self.kernelfct.compute(self.IX, X)
energy += self.alpha.dot(ic)
#print "Projection", ic
return energy
#print "Projection", ic
projection_matrix.append(ic)
if return_prj_mat:
return energy, projection_matrix
else:
return energy
def computeForces(self, structure, verbose=False):
logging.info("Start forces ...")
if verbose:
......
......@@ -15,9 +15,12 @@ logging.basicConfig(
level=logging.ERROR)
verbose = False
exclude_center_pids = []
options = soap.Options()
options.excludeCenters([])
options.excludeTargets([])
options.excludeCenterIds(exclude_center_pids)
options.excludeTargetIds([])
options.set('radialbasis.type', 'gaussian')
options.set('radialbasis.mode', 'adaptive')
options.set('radialbasis.N', 9) # 9
......@@ -26,16 +29,17 @@ options.set('radialbasis.integration_steps', 15)
options.set('radialcutoff.Rc', 6.8)
options.set('radialcutoff.Rc_width', 0.5)
options.set('radialcutoff.type', 'heaviside')
options.set('radialcutoff.center_weight', 1.)
options.set('radialcutoff.center_weight', 0.5)
options.set('angularbasis.type', 'spherical-harmonic')
options.set('angularbasis.L', 6) # 6
options.set('densitygrid.N', 20)
options.set('densitygrid.dx', 0.15)
options.set('spectrum.gradients', True)
options.set('kernel.adaptor', 'generic')
options.set('spectrum.2l1_norm', False) # <- pull here
options.set('kernel.adaptor', 'generic') # <- pull here
options.set('kernel.type', 'dot')
options.set('kernel.delta', 1.)
options.set('kernel.xi', 4.)
options.set('kernel.xi', 4.) # <- pull here
# STRUCTURE
xyzfile = 'config.xyz'
......@@ -48,9 +52,10 @@ for part in structure.particles:
positions_0 = [ part.pos for part in structure ]
# KERNEL
soap.silence()
kernelpot = KernelPotential(options)
kernelpot.acquire(structure, 1.)
kernelpot.acquire(structure, 1.)
soap.silence()
positions = [
np.array([0., 0., 0.]),
......@@ -88,15 +93,27 @@ positions = [
restore_positions(structure, positions)
print "D", kernelpot.computeEnergy(structure)
positions = [
np.array([0., 0., 0.]),
np.array([1.5, 0., 0.]),
np.array([100.,100.,100.]) ]
restore_positions(structure, positions)
#kernelpot.acquire(structure, 1.)
print "E", kernelpot.computeEnergy(structure)
# MAP
Nx = 20
Ny = 20
Nx = 40
Ny = 40
dx = 0.25
dy = 0.25
xs = []
ys = []
zs = []
ofs = open('out', 'w')
for i in range(Nx+1):
for j in range(Ny+1):
......@@ -108,9 +125,27 @@ for i in range(Nx+1):
np.array([1.5, 0., 0.]),
np.array([x, y, 0.]) ]
restore_positions(structure, positions)
e_in = kernelpot.computeEnergy(structure)
e_in, prj_mat = kernelpot.computeEnergy(structure, True)
prj_mat = np.array(prj_mat)
prj_mat = prj_mat.reshape(prj_mat.shape[0]*prj_mat.shape[1])
prj_mat_str = ('%s' % prj_mat).replace('[','').replace(']','').replace('\n','')
xs.append(x)
ys.append(y)
zs.append(prj_mat)
print prj_mat_str
ofs.write('%+1.7e %+1.7e %+1.7e\n' % (x, y, e_in))
#f = kernelpot.computeForces(structure)[2]
#ofs.write('%+1.7e %+1.7e %+1.7e %+1.7e %+1.7e %+1.7e\n' % (x, y, e_in, f[0], f[1], f[2]))
ofs.write('\n')
ofs.close()
xs = np.array(xs)
ys = np.array(ys)
zs = np.array(zs)
np.savetxt('plot/np.x.txt', xs)
np.savetxt('plot/np.y.txt', ys)
np.savetxt('plot/np.z.txt', zs)
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