Metafield tests.

src/soap/atomicspectrum.cpp

@@ -384,6 +384,7 @@ void AtomicSpectrum::registerPython() {
         .def("getPower", &AtomicSpectrum::getPower, return_value_policy<reference_existing_object>())
         .def("getPowerGradGeneric", &AtomicSpectrum::getPowerGradGeneric, return_value_policy<reference_existing_object>())
         .def("getCenter", &AtomicSpectrum::getCenter, return_value_policy<reference_existing_object>())
+        .def("getCenterId", &AtomicSpectrum::getCenterId)
         .def("getCenterType", &AtomicSpectrum::getCenterType, return_value_policy<reference_existing_object>())
         .def("getCenterPos", &AtomicSpectrum::getCenterPos, return_value_policy<reference_existing_object>());
 }

src/soap/atomicspectrum.hpp

@@ -48,6 +48,7 @@ public:
 
     // CENTER & BASIS METHODS
     Particle *getCenter() { return _center; }
+    int getCenterId() { return _center_id; }
 	std::string &getCenterType() { return _center_type; }
 	vec &getCenterPos() { return _center_pos; }
 	Basis *getBasis() { return _basis; }

src/soap/globals.cpp

@@ -5,7 +5,7 @@ namespace soap {
 Logger GLOG;
 
 void GLOG_SILENCE() {
-    GLOG() << "Silencing logger ..." << std::endl;
+    //GLOG() << "Silencing logger ..." << std::endl;
     GLOG.silence();
     return;
 }

src/soap/structure.cpp

@@ -19,6 +19,17 @@ void Particle::null() {
 	_segment = NULL;
 }
 
+void Particle::model(Particle &model) {
+    // ID assigned internally, so is Segment pointer
+    _pos = model._pos;
+    _name = model._name;
+    _type_id = model._type_id;
+    _type = model._type;
+    _mass = model._mass;
+    _weight = model._weight;
+    _sigma = model._sigma;
+}
+
 boost::python::numeric::array Particle::getPosNumeric() {
 	boost::python::numeric::array pos(boost::python::make_tuple(_pos.x(), _pos.y(), _pos.z())); return pos;
 }
@@ -70,6 +81,27 @@ Structure::Structure() {
 	this->null();
 }
 
+Structure::Structure(const Structure &structure) {
+    this->null();
+    assert(false && "COPY CONSTRUCTOR NOT AVAILABLE. USE ::model(Structure &) instead.");
+}
+
+void Structure::model(Structure &structure) {
+    this->null();
+    // ID, label, box
+    _id = structure._id;
+    _label = structure._label;
+    this->setBoundary(structure._box->getBox());
+    // Segments, particles 
+    for (segment_it_t sit = structure.beginSegments(); sit != structure.endSegments(); ++sit) {
+        Segment &new_seg = this->addSegment();
+        for (particle_it_t pit = (*sit)->beginParticles(); pit != (*sit)->endParticles(); ++pit) {
+            Particle &new_part = this->addParticle(new_seg);
+            new_part.model(*(*pit));
+        }
+    }
+}
+
 void Structure::null() {
 	_id = -1;
 	_label = "?";
@@ -157,6 +189,7 @@ boost::python::numeric::array Structure::getBoundaryNumeric() {
 void Structure::registerPython() {
 	using namespace boost::python;
 	class_<Structure>("Structure", init<std::string>())
+       .def("model", &Structure::model)
 	   .def("addSegment", &Structure::addSegment, return_value_policy<reference_existing_object>())
 	   .def("getSegment", &Structure::getSegment, return_value_policy<reference_existing_object>())
 	   .def("addParticle", &Structure::addParticle, return_value_policy<reference_existing_object>())

src/soap/structure.hpp

@@ -20,10 +20,11 @@ class Structure;
 class Particle
 {
 public:
-    Particle(int id) { this->null(); _id = id; }
+    explicit Particle(int id) { this->null(); _id = id; }
     Particle() { this->null(); }
    ~Particle() {;}
     void null();
+    void model(Particle &model);
     // Position
     void setPos(vec &pos) { _pos = pos; }
     void setPos(double x, double y, double z) { _pos = vec(x,y,z); }
@@ -91,7 +92,7 @@ public:
 	typedef std::vector<Particle*> particle_array_t;
 	typedef particle_array_t::iterator particle_it_t;
 
-    Segment(int id) : _id(id), _name("?"), _type("?") { ; }
+    explicit Segment(int id) : _id(id), _name("?"), _type("?") { ; }
     Segment() : _id(-1), _name("?"), _type("?") { ; }
    ~Segment() { _particles.clear(); }
 
@@ -138,10 +139,12 @@ public:
 	typedef std::vector<Segment*> segment_array_t;
 	typedef std::vector<Segment*>::iterator segment_it_t;
 
-    Structure(std::string label);
+    explicit Structure(std::string label);
+    Structure(const Structure &structure);
     Structure();
    ~Structure();
     void null();
+    void model(Structure &structure);
 
     // PARTICLE CONTAINER
     particle_array_t &particles() { return _particles; }

test/test_forcefield/0_interact.py

+#! /usr/bin/env python
+import soap
+import soap.tools
+soap.silence()
+
+import os
+import sys
+import numpy as np
+import logging
+import io
+
+from momo import osio, endl, flush
+from kernel import KernelAdaptorFactory, KernelFunctionFactory, TrajectoryLogger
+from simspace import SimSpaceTopology, SimSpaceNode, SimSpacePotential
+from simspace import LJRepulsive
+from simspace import optimize_node
+from dimred import dimred
+
+logging.basicConfig(
+    format='[%(asctime)s] %(message)s',
+    datefmt='%I:%M:%S',
+    level=logging.ERROR)
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+select_struct_label = 'config.xyz'
+exclude_centers = []
+adaptor_type = 'global-generic'
+kernel_type = 'dot' # 'dot-harmonic'
+alpha = +1.
+mu = 0.9
+sc = 0.1 # 1.5 # 0.1
+average_energy = True
+lj_sigma = 0.02 # 0.00001 # 0.02
+# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+# NOTE adaptor_type = 'generic' => exclude center if not constrained in optimization
+
+def sample_random_x0(struct, exclude_pids_rnd, opt_pids, b0, b1):
+    # Sample
+    x0 = [ part.pos for part in struct ]
+    x0 = random_positions(struct, exclude_pid=exclude_pids_rnd, b0=b0, b1=b1)
+    x0 = np.array(x0)    
+    # Short-list
+    opt_pids = np.array(opt_pids)
+    opt_pidcs = opt_pids-1
+    x0 = x0[opt_pidcs]
+    x0 = x0.flatten()    
+    return x0
+
+def evaluate_energy_local(positions, structure, kernelpot, opt_pids, verbose=False, ofs=None, average=False):
+    if verbose: print "Energy"
+    # Impose positions
+    pid_pos = positions.reshape((opt_pids.shape[0],3))
+    for pidx, pid in enumerate(opt_pids):
+        pos = pid_pos[pidx,:]
+        particle = structure.getParticle(pid)        
+        particle.pos = pos        
+    for part in structure:
+        if verbose: print part.id, part.type, part.pos
+    # Evaluate energy function
+    energy = kernelpot.computeEnergy(structure)
+    if average: energy /= kernelpot.nConfigs()
+    # Evaluate additional potentials
+    lj = LJRepulsive()
+    energy_add = lj.computeEnergy(struct)
+    energy += energy_add
+    # Log
+    if ofs: ofs.logFrame(structure)
+    if verbose: print energy
+    print energy
+    return energy
+
+def evaluate_energy_gradient_local(positions, structure, kernelpot, opt_pids, verbose=False, ofs=None, average=False):
+    if verbose: print "Forces"
+    # Adjust positions
+    pid_pos = positions.reshape((opt_pids.shape[0],3))
+    for pidx, pid in enumerate(opt_pids):
+        pos = pid_pos[pidx,:]
+        particle = structure.getParticle(pid)        
+        particle.pos = pos        
+    for part in structure:
+        if verbose: print part.id, part.type, part.pos
+    # Evaluate forces
+    forces = kernelpot.computeForces(structure)
+    gradients = -1.*np.array(forces)
+    if average:
+        gradients = gradients/kernelpot.nConfigs()
+    # Evaluate additional potentials
+    lj = LJRepulsive()
+    gradients_add = lj.computeGradient(struct)
+    gradients = gradients + gradients_add
+    # Short-list
+    opt_pidcs = opt_pids-1
+    gradients_short = gradients[opt_pidcs]
+    gradients_short = gradients_short.flatten()
+    if verbose: print gradients_short
+    #forces[2] = 0. # CONSTRAIN TO Z=0 PLANE
+    return gradients_short
+
+# OPTIONS
+options = soap.Options()
+options.excludeCenters(['H'])
+options.excludeTargets(['H'])
+options.excludeCenterIds(exclude_centers)
+options.excludeTargetIds([])
+# Spectrum
+options.set('spectrum.gradients', True)
+options.set('spectrum.2l1_norm', False)                     # <- pull here (True/False)
+# Basis
+options.set('radialbasis.type', 'gaussian')
+options.set('radialbasis.mode', 'adaptive')
+options.set('radialbasis.N', 9)
+options.set('radialbasis.sigma', 0.5)
+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)
+# Kernel
+options.set('kernel.adaptor', adaptor_type)                 # <- pull here (generic/global-generic)
+options.set('kernel.type', kernel_type)
+options.set('kernel.delta', 1.)
+options.set('kernel.xi', 2.)                                # <- pull here (1., ..., 4., ...)
+options.set('kernel.mu', mu)
+# Cube files
+options.set('densitygrid.N', 20)
+options.set('densitygrid.dx', 0.15)
+
+
+
+
+
+
+# LOAD STRUCTURES
+structures = [ 
+    soap.tools.setup_structure_ase(c.config_file, c.atoms)
+    for c in soap.tools.ase_load_all('in.configs')
+]
+struct_dict = {}
+for struct in structures:
+    print struct.label
+    struct_dict[struct.label] = struct
+
+# SELECT STRUCTURE
+struct = struct_dict[select_struct_label]
+
+# DEFINE INTERACTION TEMPLATES
+pot_options_dot = soap.Options()
+pot_options_dot.set('kernel.adaptor', 'global-generic')
+pot_options_dot.set('kernel.type', 'dot')
+pot_options_dot.set('kernel.alpha', -1.)
+pot_options_dot.set('kernel.delta', 1.)
+pot_options_dot.set('kernel.xi', 2.)
+pot_options_dot.set('kernel.mu', 0.5)
+pot_options_harm = soap.Options()
+pot_options_harm.set('kernel.adaptor', 'global-generic')
+pot_options_harm.set('kernel.type', 'dot-harmonic')
+pot_options_harm.set('kernel.alpha', +1.)
+pot_options_harm.set('kernel.delta', 1.)
+pot_options_harm.set('kernel.xi', 2.)
+pot_options_harm.set('kernel.mu', 0.5)
+pot_options_lj = soap.Options()
+pot_options_lj.set('potentials.lj.sigma', 0.02)
+
+# TOPOLOGY
+top = SimSpaceTopology(options)
+
+# DEFINE NODES
+node_root = top.createNode(struct)
+
+"""
+node_leaf = top.createNode(struct)
+top.summarize()
+
+# META-INTERACTIONS
+potentials = []
+pot = SimSpacePotential(node_leaf, node_root, pot_options_harm)
+potentials.append(pot)
+
+# SELF-INTERACTIONS
+potentials_self = []
+pot = LJRepulsive(node_leaf, pot_options_lj)
+potentials_self.append(pot)
+
+# OPTIMIZE
+x0 = node_leaf.randomizePositions(scale=sc, zero_pids=[0], compute=True)
+optimize_node(node_leaf, potentials, potentials_self, np.array([1,2]), x0)
+top.writeData()
+"""
+
+for n in range(6):
+    node_leaf = top.createNode(struct)    
+    node_root = top.nodes[0]
+    
+    # POTENTIALS
+    potentials = []
+    potentials_self = []
+    
+    # DOT-HARMONIC
+    pot_options_harm.set('kernel.adaptor', 'global-generic')
+    pot_options_harm.set('kernel.type', 'dot-harmonic')
+    pot_options_harm.set('kernel.alpha', +10.)
+    pot_options_harm.set('kernel.delta', 1.)
+    pot_options_harm.set('kernel.xi', 2.)
+    pot_options_harm.set('kernel.mu', 0.95)
+    pot = SimSpacePotential(node_leaf, node_root, pot_options_harm)
+    potentials.append(pot)
+    print "<dot-harm, ij>", node_leaf.id, node_root.id
+    
+    # DOT
+    for j in range(1, len(top.nodes)-1):
+        pot_options_dot.set('kernel.adaptor', 'global-generic')
+        pot_options_dot.set('kernel.type', 'dot-harmonic')
+        pot_options_dot.set('kernel.alpha', +1.)
+        pot_options_dot.set('kernel.delta', 1.)
+        pot_options_dot.set('kernel.xi', 2.)
+        pot_options_dot.set('kernel.mu', 0.95)
+        pot = SimSpacePotential(node_leaf, top.nodes[j], pot_options_dot)
+        potentials.append(pot)
+        print "<dot, ij>", node_leaf.id, top.nodes[j].id
+    
+    # LJ
+    pot_options_lj.set('potentials.lj.sigma', 0.02)
+    pot = LJRepulsive(node_leaf, pot_options_lj)
+    potentials_self.append(pot)
+    
+    if n == 0:
+        x0 = node_leaf.randomizePositions(scale=sc, zero_pids=[0], compute=True)
+    else:
+        x0 = np.array([ part.pos for part in node_leaf.structure ])
+    optimize_node(node_leaf, potentials, potentials_self, np.array([1,2]), x0)
+
+
+top.writeData()
+    
+
+
+
+
+
+
+
+

test/test_forcefield/in.configs/config.xyz

+3
+
+C 0.0 0.0 0.0
+C 1.5 0.0 0.0
+C 0.0 1.5 0.0

test/test_forcefield/simspace.py

+#! /usr/bin/env python
+import soap
+import soap.tools
+
+import os
+import sys
+import numpy as np
+import logging
+import io
+
+import scipy.optimize
+
+from kernel import KernelAdaptorFactory, KernelFunctionFactory, TrajectoryLogger
+
+class SimSpaceNode(object):
+    def __init__(self, node_id, structure, basis, options, compute=True):
+        self.id = node_id
+        # Copy structure
+        self.structure = soap.Structure("?")
+        self.structure.model(structure)
+        # Basis, options
+        self.basis = basis
+        self.options = options
+        # Stored spectra
+        self.reset()
+        # Spectrum adaptor
+        self.adaptor = KernelAdaptorFactory[options.get('kernel.adaptor')](options)
+        if compute: self.acquire()
+        return
+    def size(self):
+        return self.IX.shape[0]
+    def reset(self):
+        self.IX = None
+        self.dimX = None
+        self.pid_X_unnorm = {}
+        self.pid_X_norm = {}
+        self.pid_nbpid_dX = {}
+    def assignPositions(self, x0, compute=True):
+        for idx, part in enumerate(self.structure):
+            part.pos = x0[idx]
+        if compute: self.acquire()
+        return
+    def randomizePositions(self, scale=0.1, zero_pids=[0], compute=True):
+        x0 = np.random.uniform(-1.*scale, +1.*scale, (self.structure.n_particles, 3))
+        for pid in zero_pids:
+            x0[pid-1,:] = 0.
+        self.assignPositions(x0, compute)
+        return x0
+    def acquire(self, reset=True):
+        if reset: self.reset()
+        # Compute spectrum
+        self.spectrum = soap.Spectrum(self.structure, self.options, self.basis)
+        self.spectrum.compute()
+        self.spectrum.computePower()
+        self.spectrum.computePowerGradients()
+        self.spectrum.computeGlobal()
+        # PID-resolved storage for X, dX
+        for atomic in self.adaptor.getListAtomic(self.spectrum):
+            pid = atomic.getCenterId()
+            X_unnorm, X_norm = self.adaptor.adaptScalar(atomic) # TODO Done again below in ::adapt => simplify
+            self.pid_X_unnorm[pid] = X_unnorm
+            self.pid_X_norm[pid] = X_norm
+            self.pid_nbpid_dX[pid] = {}
+            # NB-PID-resolved derivatives
+            nb_pids = atomic.getNeighbourPids()
+            for nb_pid in nb_pids:
+                dX_dx, dX_dy, dX_dz = self.adaptor.adaptGradients(atomic, nb_pid, X_unnorm)
+                self.pid_nbpid_dX[pid][nb_pid] = (dX_dx, dX_dy, dX_dz)
+        # New X's
+        IX_acqu = self.adaptor.adapt(self.spectrum)
+        n_acqu = IX_acqu.shape[0]
+        dim_acqu = IX_acqu.shape[1]
+        if not self.dimX:
+            # First time ...
+            self.dimX = dim_acqu
+            self.IX = IX_acqu
+        else:
+            # Check and extend ...
+            assert self.dimX == dim_acqu # Acquired descr. should match linear dim. of previous descr.'s
+            I = self.IX.shape[0]
+            self.IX.resize((I+n_acqu, self.dimX))
+            self.IX[I:I+n_acqu,:] = IX_acqu
+        return
+    def getListAtomic(self):
+        return self.adaptor.getListAtomic(self.spectrum)
+    def getPidX(self, pid):
+        return self.pid_X_unnorm[pid], self.pid_X_norm[pid]
+    def getPidGradX(self, pid, nb_pid):
+        return self.pid_nbpid_dX[pid][nb_pid]
+        
+class SimSpaceTopology(object):
+    def __init__(self, options):
+        self.basis = soap.Basis(options)
+        self.options = options
+        self.nodes = []
+        self.kernelfct = KernelFunctionFactory[options.get('kernel.type')](options)
+        self.IX = np.zeros((0))
+        return
+    def compileIX(self):
+        n_nodes = len(self.nodes)
+        dim = self.nodes[0].IX.shape[1]
+        self.IX = np.zeros((n_nodes,dim))
+        for idx, node in enumerate(self.nodes):
+            self.IX[idx] = node.IX
+        return
+    def computeKernelMatrix(self, return_distance=False):
+        if not self.IX.any(): self.compileIX()
+        K = self.kernelfct.computeBlock(self.IX, return_distance)
+        return K
+    def createNode(self, structure):
+        node_id = len(self.nodes)+1
+        node = SimSpaceNode(node_id, structure, self.basis, self.options)
+        self.nodes.append(node)
+        return node
+    def summarize(self):
+        for idx, node in enumerate(self.nodes):
+            print "Node %d" % (idx+1)            
+            for part in node.structure:
+                print part.id, part.type, part.pos
+        return
+    def writeData(self, prefix='out.top'):
+        # Structures
+        ofs = TrajectoryLogger('%s.xyz' % prefix)
+        for node in self.nodes:
+            ofs.logFrame(node.structure)
+        ofs.close()
+        # Descriptors
+        self.compileIX()
+        np.savetxt('%s.ix.txt' % prefix, self.IX)
+        # Kernel
+        K = self.computeKernelMatrix()
+        np.savetxt('%s.kernelmatrix.txt' % prefix, K)
+        return
+
+class SimSpacePotential(object):
+    def __init__(self, target, source, options):
+        self.target = target
+        self.source = source
+        self.alpha = np.array([float(options.get('kernel.alpha'))])
+        # INTERACTION FUNCTION
+        self.kernelfct = KernelFunctionFactory[options.get('kernel.type')](options)
+        # Spectrum adaptor
+        self.adaptor = KernelAdaptorFactory[options.get('kernel.adaptor')](options)   
+        return
+    def computeEnergy(self, return_prj_mat=False):
+        energy = 0.0
+        projection_matrix = []
+        for n in range(self.target.size()):
+            X = self.target.IX[n]
+            ic = self.kernelfct.compute(self.source.IX, X)
+            energy += self.alpha.dot(ic)
+            projection_matrix.append(ic)
+        if return_prj_mat:
+            return energy, projection_matrix
+        else:
+            return energy
+    def computeForces(self, verbose=False):
+        forces = [ np.zeros((3)) for i in range(self.target.structure.n_particles) ]        
+        for atomic in self.target.getListAtomic():
+            pid = atomic.getCenterId()
+            nb_pids = atomic.getNeighbourPids()
+            if verbose: print "  Center %d" % (pid)
+            # neighbour-pid-independent kernel "prevector" (outer derivative)
+            X_unnorm, X_norm = self.target.getPidX(pid)
+            dIC = self.kernelfct.computeDerivativeOuter(self.source.IX, X_norm)
+            alpha_dIC = self.alpha.dot(dIC)
+            for nb_pid in nb_pids:
+                # Force on neighbour
+                if verbose: print "    -> Nb %d" % (nb_pid)
+                dX_dx, dX_dy, dX_dz = self.target.getPidGradX(pid, nb_pid)
+                force_x = -alpha_dIC.dot(dX_dx)
+                force_y = -alpha_dIC.dot(dX_dy)
+                force_z = -alpha_dIC.dot(dX_dz)                
+                forces[nb_pid-1][0] += force_x
+                forces[nb_pid-1][1] += force_y
+                forces[nb_pid-1][2] += force_z
+        return np.array(forces)
+    def computeGradients(self, verbose=False):
+        return -1 * self.computeForces(verbose)
+
+class LJRepulsive(object):
+    def __init__(self, node, options):
+        self.sigma = float(options.get('potentials.lj.sigma'))
+        self.struct = node.structure
+        return
+    def computeEnergy(self):
+        coords = np.zeros((self.struct.n_particles, 3))
+        for idx, part in enumerate(self.struct):
+            coords[idx,:] = part.pos
+        E = 0.
+        for i in range(self.struct.n_particles):
+            for j in range(i):
+                r = np.sqrt(np.sum((coords[i, :] - coords[j, :]) ** 2))
+                E += (self.sigma/r)**12
+        return E
+    def computeGradient(self):
+        coords = np.zeros((self.struct.n_particles, 3))
+        for idx, part in enumerate(self.struct):
+            coords[idx,:] = part.pos            
+        grad = np.zeros(coords.shape)
+        for i in range(self.struct.n_particles):
+            for j in range(i):
+                dr = coords[i, :] - coords[j, :]
+                r = np.sqrt(np.sum(dr**2))
+                g = 12./self.sigma*(self.sigma/r)**13
+                grad[i, :] += -g * dr / r
+                grad[j, :] += g * dr / r
+        return grad
+
+def evaluate_potential_energy(x0, node, potentials, potentials_self, opt_pidcs, trjlog):
+    #print "Energy at", x0
+    # Adjust positions
+    pid_pos = x0.reshape((opt_pidcs.shape[0],3))
+    for rel_idx, abs_idx in enumerate(opt_pidcs):
+        pos = pid_pos[rel_idx,:]
+        particle = node.structure.getParticle(abs_idx+1)        
+        particle.pos = pos
+    node.acquire()
+    # Compute energy
+    energy = 0.
+    for pot in potentials:
+        energy += pot.computeEnergy()
+    # Structure self-energy
+    energy_self = 0.