Commit 758d7722 authored by Carl Poelking's avatar Carl Poelking
Browse files

Global spectrum toy analysis.

parent 7f3aa2d8
......@@ -107,6 +107,9 @@ void RadialBasisGaussian::configure(Options &options) {
GLOG() << "Adjusted radial cutoff to " << _Rc
<< " based on sigma_0 = " << sigma_0 << ", L = " << L << ", stride = " << sigma_stride_factor << std::endl;
}
else if (_mode == "incremental") {
throw std::runtime_error("Not implemented.");
}
else {
throw std::runtime_error("Not implemented.");
}
......
......@@ -7,7 +7,7 @@ import numpy as np
import logging
from momo import osio, endl, flush
class KernelAdaptorGeneric:
class KernelAdaptorGeneric(object):
def __init__(self, options):
return
def adapt(self, spectrum):
......@@ -69,6 +69,46 @@ class KernelAdaptorGeneric:
dX_dz = dX_dz/mag_X - np.dot(X, dX_dz)/mag_X**3 * X
return dX_dx, dX_dy, dX_dz
class KernelAdaptorGlobalGeneric(object):
def __init__(self, options):
return
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):
# EXTRACT POWER EXPANSION FROM ATOMIC SPECTRUM
# Here: type "":"" (= generic)
X = np.real(atomic.getPower("","").array)
dimX = X.shape[0]*X.shape[1]
X = np.real(X.reshape((dimX)))
# Normalize
X_norm = X/np.dot(X,X)**0.5
return X, X_norm
def adaptGradients(self, atomic, nb_pid, X):
# NOTE X is not normalized (=> X = X')
dxnkl_pid = atomic.getPowerGradGeneric(nb_pid)
dX_dx = dxnkl_pid.getArrayGradX()
dX_dy = dxnkl_pid.getArrayGradY()
dX_dz = dxnkl_pid.getArrayGradZ()
dimX = dX_dx.shape[0]*dX_dx.shape[1]
dX_dx = np.real(dX_dx.reshape((dimX)))
dX_dy = np.real(dX_dy.reshape((dimX)))
dX_dz = np.real(dX_dz.reshape((dimX)))
# Normalize
mag_X = np.dot(X,X)**0.5
dX_dx = dX_dx/mag_X - np.dot(X, dX_dx)/mag_X**3 * X
dX_dy = dX_dy/mag_X - np.dot(X, dX_dy)/mag_X**3 * X
dX_dz = dX_dz/mag_X - np.dot(X, dX_dz)/mag_X**3 * X
return dX_dx, dX_dy, dX_dz
class KernelFunctionDot(object):
def __init__(self, options):
self.delta = float(options.get('kernel.delta'))
......@@ -82,10 +122,10 @@ class KernelFunctionDot(object):
c = self.computeDot(IX, X, self.xi-1, self.delta)
return self.xi*np.diag(c).dot(IX)
KernelAdaptorFactory = { 'generic': KernelAdaptorGeneric }
KernelAdaptorFactory = { 'generic': KernelAdaptorGeneric, 'global-generic': KernelAdaptorGlobalGeneric }
KernelFunctionFactory = { 'dot':KernelFunctionDot }
class KernelPotential:
class KernelPotential(object):
def __init__(self, options):
logging.info("Construct kernel potential ...")
self.basis = soap.Basis(options)
......@@ -100,6 +140,7 @@ class KernelPotential:
# ADAPTOR
logging.info("Choose adaptor ...")
self.adaptor = KernelAdaptorFactory[options.get('kernel.adaptor')](options)
self.use_global_spectrum = True if options.get('kernel.adaptor') == 'global-generic' else False
# INCLUSIONS / EXCLUSIONS
# -> Already be enforced when computing spectra
return
......@@ -109,6 +150,8 @@ class KernelPotential:
spectrum.compute()
spectrum.computePower()
spectrum.computePowerGradients()
if self.use_global_spectrum:
spectrum.computeGlobal()
# New X's
logging.info("Adapt spectrum ...")
IX_acqu = self.adaptor.adapt(spectrum)
......@@ -141,6 +184,9 @@ class KernelPotential:
spectrum.compute()
spectrum.computePower()
spectrum.computePowerGradients()
# TODO The kernel policy should take care of this:
if self.use_global_spectrum:
spectrum.computeGlobal()
IX_acqu = self.adaptor.adapt(spectrum)
n_acqu = IX_acqu.shape[0]
dim_acqu = IX_acqu.shape[1]
......@@ -168,9 +214,14 @@ class KernelPotential:
spectrum.compute()
spectrum.computePower()
spectrum.computePowerGradients()
if self.use_global_spectrum:
atomic_global = spectrum.computeGlobal()
spectrum_iter = [ atomic_global ]
else:
spectrum_iter = spectrum
# Extract & compute force
for atomic in spectrum:
pid = atomic.getCenter().id
for atomic in spectrum_iter:
pid = atomic.getCenter().id if not use_global_spectrum else -1
#if not pid in self.pid_list_force:
# logging.debug("Skip forces derived from environment with pid = %d" % pid)
# continue
......
#! /usr/bin/env python
import soap
import soap.tools
import os
import numpy as np
import logging
from momo import osio, endl, flush
from kernel import KernelPotential, restore_positions, apply_force_step
logging.basicConfig(
format='[%(asctime)s] %(message)s',
datefmt='%I:%M:%S',
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
options.set('radialbasis.sigma', 0.5) # 0.9
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', 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('spectrum.2l1_norm', False) # <- pull here
options.set('kernel.adaptor', 'global-generic') # <- pull here
options.set('kernel.type', 'dot')
options.set('kernel.delta', 1.)
options.set('kernel.xi', 4.) # <- pull here
# STRUCTURE
xyzfile = 'config.xyz'
config = soap.tools.ase_load_single(xyzfile)
structure = soap.tools.setup_structure_ase(config.config_file, config.atoms)
for part in structure.particles:
print part.id, part.type, part.pos
positions_0 = [ part.pos for part in structure ]
# KERNEL
kernelpot = KernelPotential(options)
kernelpot.acquire(structure, 1.)
soap.silence()
positions = [
np.array([0., 0., 0.]),
np.array([1.75, 0., 0.]),
np.array([0., 1.75, 0.]) ]
restore_positions(structure, positions)
#kernelpot.acquire(structure, 1.)
print "O'", kernelpot.computeEnergy(structure)
positions = [
np.array([0., 0., 0.]),
np.array([1.5, 0., 0.]),
np.array([1.5, 0., 0.]) ]
restore_positions(structure, positions)
print "A", kernelpot.computeEnergy(structure)
positions = [
np.array([0., 0., 0.]),
np.array([1.5, 0., 0.]),
np.array([0., 1.5, 0.]) ]
restore_positions(structure, positions)
print "B", kernelpot.computeEnergy(structure)
positions = [
np.array([0., 0., 0.]),
np.array([0., 1.5, 0.]),
np.array([0., 0., 1.5]) ]
restore_positions(structure, positions)
print "C", kernelpot.computeEnergy(structure)
positions = [
np.array([0., 0., 0.]),
np.array([0., 0., 1.5]),
np.array([1.5, 0., 0.]) ]
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 = 40
Ny = 40
dx = 0.25
dy = 0.25
xs = []
ys = []
zs = []
ofs = open('out_global', 'w')
for i in range(Nx+1):
for j in range(Ny+1):
x = i*dx - 0.5*Nx*dx
y = j*dy - 0.5*Ny*dy
print i, j
positions = [
np.array([0., 0., 0.]),
np.array([1.5, 0., 0.]),
np.array([x, y, 0.]) ]
restore_positions(structure, positions)
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_global/np.x.txt', xs)
np.savetxt('plot_global/np.y.txt', ys)
np.savetxt('plot_global/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