SIMS-LJ tree systems.

test/test_forcefield/0_explore_lj.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, SimSpaceTree, SimSpaceJoint
+from simspace import LJRepulsive
+from simspace import optimize_node, multi_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
+
+# 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)
+"""
+
+
+
+
+
+
+
+osio.cd('out.files')
+tree = SimSpaceTree(options)
+tree.seed(struct)
+
+N_generations = 2
+N_spawn = [3,3]
+
+N_generations = 4
+N_spawn = [2,2,2,2]
+pot_options_dotlj = soap.Options()
+pot_options_dotlj.set('kernel.adaptor', 'global-generic')
+pot_options_dotlj.set('kernel.type', 'dot-lj')
+pot_options_dotlj.set('kernel.alpha', +1.)
+pot_options_dotlj.set('kernel.delta', 1.)
+pot_options_dotlj.set('kernel.xi', 2.)
+pot_options_dotlj.set('kernel.lj_sigma', 0.1) # 0.2
+pot_options_dotlj.set('kernel.lj_eps_cap', 0.001)
+
+for n in range(N_generations):
+    gen_id = n+1
+    gen_prefix = 'out.gen_%d' % gen_id
+    print "Spawn generation %d" % gen_id    
+
+    new_joints = tree.spawn(N_spawn[n])
+    tree.summarize()    
+    tree.writeParentChildPairs('%s.tree.txt' % gen_prefix)
+    
+    # POTENTIALS
+    tree.clearPotentials()
+    tree.addPairPotential(pot_options_dotlj)
+    
+    # OPTIMIZE 1st GENERATION
+    nodes_opt = [ joint.node for joint in tree.generations[-1] ]
+    x0_opt = []
+    for node in nodes_opt:
+        x0 = node.perturbPositions(scale=0.1, zero_pids=[], compute=True)
+        x0_opt.append(x0)
+    x0_opt = np.array(x0_opt).flatten()
+    multi_optimize_node(nodes_opt, x0_opt)
+    
+    # OUTPUT
+    tree.top.writeData(prefix=gen_prefix)
+    
+    raw_input('...')
+
+osio.root()
+
+
+
+
+
+
+

test/test_forcefield/simspace.py

@@ -97,6 +97,18 @@ class SimSpaceNode(object):
     def createPotential(self, nb_node, options): # TODO Simplify. Only one potential needed per pair.
         self.potentials.append(SimSpacePotential(self, nb_node, options))
         return
+    def clearPotentials(self):
+        del self.potentials
+        del self.potentials_self
+        self.potentials = []
+        self.potentials_self = []
+    def evaluatePotentialEnergy(self):
+        self.energy = 0.
+        for pot in self.potentials:
+            self.energy += pot.computeEnergy()
+        for pot in self.potentials_self:
+            self.energy += pot.computeEnergy()
+        return self.energy
     def getListAtomic(self):
         return self.adaptor.getListAtomic(self.spectrum)
     def getPidX(self, pid):
@@ -123,6 +135,11 @@ class SimSpaceTopology(object):
         if not self.IX.any(): self.compileIX()
         K = self.kernelfct.computeBlock(self.IX, return_distance)
         return K
+    def computePotentialEnergy(self):
+        self.energy = 0.
+        for node in self.nodes:
+            self.energy += node.computeEnergy()
+        return self.energy
     def createNode(self, structure):
         node_id = len(self.nodes)+1
         node = SimSpaceNode(node_id, structure, self.basis, self.options)
@@ -194,6 +211,84 @@ class SimSpacePotential(object):
     def computeGradients(self, verbose=False):
         return -1 * self.computeForces(verbose)
 
+class SimSpaceJoint(object):
+    def __init__(self, top, par_joint=None, structure=None):
+        self.top = top
+        self.par_joint = par_joint # TODO Extend to several parents => structural cross-over
+        if self.par_joint:
+            assert structure == None
+            self.node = top.createNode(par_joint.node.structure)
+        else:
+            self.node = top.createNode(structure)
+        self.child_joints = []
+    def hasChildren(self):
+        return len(self.child_joints) > 0
+    def spawn(self, n):
+        new_joints = []
+        if self.hasChildren():
+            for child_joint in self.child_joints:
+                new_joints = new_joints + child_joint.spawn(n)
+        else:
+            for i in range(n):
+                child_joint = SimSpaceJoint(self.top, self)
+                self.child_joints.append(child_joint)
+                new_joints.append(child_joint)
+        return new_joints
+    def summarize(self, indent=''):
+        print '%sID=%d: %d children' % (indent, self.node.id, len(self.child_joints))
+        for child in self.child_joints:
+            child.summarize(indent=indent+'  ')
+        return
+    def writeParentChildPairs(self, ofs):
+        for child in self.child_joints:
+            ofs.write('%d %d\n' % (self.node.id, child.node.id))
+            child.writeParentChildPairs(ofs)
+        return
+
+class SimSpaceTree(object):
+    def __init__(self, options):
+        self.top = SimSpaceTopology(options)
+        self.joints = []
+        self.generations = []
+        self.root_joint = None
+    def seed(self, structure):
+        self.root_joint = SimSpaceJoint(self.top, par_joint=None, structure=structure)
+        self.joints.append(self.root_joint)
+        self.generations.append([self.root_joint])
+        return self.root_joint
+    def spawn(self, n):
+        new_joints = self.root_joint.spawn(n)
+        self.joints = self.joints + new_joints
+        self.generations.append(new_joints)
+        return new_joints
+    def summarize(self):
+        n_generations = len(self.generations)
+        n_joints = len(self.joints)
+        print "Generations # = %d   Joints # = %d" % (n_generations, n_joints)
+        for idx, gen in enumerate(self.generations):
+            print "  Gen %d : Joints %d" % (idx, len(gen))
+        self.root_joint.summarize()
+    def writeParentChildPairs(self, outfile='out.tree.txt'):
+        ofs = open(outfile, 'w')
+        self.root_joint.writeParentChildPairs(ofs)
+        ofs.close()
+        return
+    def addPairPotential(self, options):
+        n_joints = len(self.joints)
+        for i in range(n_joints):
+            node_i = self.joints[i].node
+            for j in range(i+1, n_joints):
+                node_j = self.joints[j].node
+                print "%d:%d" % (node_i.id, node_j.id),
+                node_i.createPotential(node_j, options)
+                node_j.createPotential(node_i, options)
+            print ""
+        return
+    def clearPotentials(self):
+        for joint in self.joints:
+            joint.node.clearPotentials()
+        return
+
 class LJRepulsive(object):
     def __init__(self, node, options):
         self.sigma = float(options.get('potentials.lj.sigma'))

test/test_forcefield/symmetrize.py

@@ -18,7 +18,7 @@ from dimred import dimred_matrix
 
 method = sys.argv[1]
 
-prj_dimension = 2
+prj_dimension = 3
 
 if method == 'mds':
 	dimred_matrix('mds', kmat=K, distmat=Dsym, outfile='tmp.txt', prj_dimension=prj_dimension)

test/test_kernelpot/kernel.py

@@ -185,10 +185,25 @@ class KernelFunctionDotHarmonic(object):
     def computeBlockDot(self, IX, return_distance=False):
         return self.kfctdot.computeBlock(IX, return_distance)
         
-        
+class KernelFunctionDotLj(object):
+    def __init__(self, options):
+        self.sigma = float(options.get('kernel.lj_sigma'))
+        self.eps_cap = float(options.get('kernel.lj_eps_cap'))
+        self.kfctdot = KernelFunctionDot(options)
+    def compute(self, IX, X):
+        C = self.kfctdot.compute(IX, X)
+        D = (1.-C+self.eps_cap)**0.5
+        return (self.sigma/D)**12 - (self.sigma/D)**6
+    def computeDerivativeOuter(self, IX, X):
+        C = self.kfctdot.compute(IX, X)
+        D = (1.-C+self.eps_cap)**0.5
+        IC = self.kfctdot.computeDerivativeOuter(IX, X)
+        return (-6.*self.sigma**12*(1./D)**7 + 3.*self.sigma**6*(1./D)**4)*IC
+    def computeBlockDot(self, IX, return_distance=False):
+        return self.kfctdot.computeBlock(IX, return_distance)        
 
 KernelAdaptorFactory = { 'generic': KernelAdaptorGeneric, 'global-generic': KernelAdaptorGlobalGeneric }     
-KernelFunctionFactory = { 'dot':KernelFunctionDot, 'dot-harmonic': KernelFunctionDotHarmonic }
+KernelFunctionFactory = { 'dot':KernelFunctionDot, 'dot-harmonic': KernelFunctionDotHarmonic, 'dot-lj': KernelFunctionDotLj }
 
 class KernelPotential(object):
     def __init__(self, options):