Three-body terms.

test/test_forcefield/0_explore_3body.py

@@ -26,7 +26,7 @@ select_struct_label = 'config.xyz'
 exclude_centers = []
 adaptor_type = 'global-generic'
 kernel_type = 'dot'
-sc = 0.1
+sc = 0.05
 lj_sigma = 0.02
 N_network = 3
 # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -69,28 +69,48 @@ 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.alpha', +1.) # TODO
 pot_options_harm.set('kernel.delta', 1.)
 pot_options_harm.set('kernel.xi', 2.)
-pot_options_harm.set('kernel.mu', 0.95)
+pot_options_harm.set('kernel.mu', 0.95) # TODO
 
 pot_options_3harm = soap.Options()
 pot_options_3harm.set('kernel.adaptor', 'global-generic')
 pot_options_3harm.set('kernel.type', 'dot-3-harmonic')
-pot_options_3harm.set('kernel.alpha', +1.)
+pot_options_3harm.set('kernel.alpha', +5.) # TODO
 pot_options_3harm.set('kernel.delta', 1.)
 pot_options_3harm.set('kernel.xi', 2.)
 
-pot_options_lj = soap.Options()
-pot_options_lj.set('potentials.lj.sigma', 0.02)
+
+hex_c = 0.1
+pot_options_harm_dist = soap.Options()
+pot_options_harm_dist.set('kernel.adaptor', 'global-generic')
+pot_options_harm_dist.set('kernel.type', 'dot-harmonic-dist')
+pot_options_harm_dist.set('kernel.alpha', +1.) # TODO
+pot_options_harm_dist.set('kernel.delta', 1.)
+pot_options_harm_dist.set('kernel.xi', 2.)
+pot_options_harm_dist.set('kernel.dotharmdist_d0', hex_c) # TODO
+pot_options_harm_dist.set('kernel.dotharmdist_eps', 0.001)
+
+pot_options_3harm_dist = soap.Options()
+pot_options_3harm_dist.set('kernel.adaptor', 'global-generic')
+pot_options_3harm_dist.set('kernel.type', 'dot-3-harmonic-dist')
+pot_options_3harm_dist.set('kernel.alpha', +1.) # TODO
+pot_options_3harm_dist.set('kernel.delta', 1.)
+pot_options_3harm_dist.set('kernel.xi', 2.)
+pot_options_3harm_dist.set('kernel.dot3harmdist_d0', 3.**0.5*hex_c)
+pot_options_3harm_dist.set('kernel.dot3harmdist_eps', 0.001)
+
+#pot_options_lj = soap.Options()
+#pot_options_lj.set('potentials.lj.sigma', 0.02)
 
 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.alpha', +0.1) # TODO
 pot_options_dotlj.set('kernel.delta', 1.)
 pot_options_dotlj.set('kernel.xi', 2.)
-pot_options_dotlj.set('kernel.lj_sigma', 0.125) # 0.1 # 0.2 too large => TODO Better force capping required
+pot_options_dotlj.set('kernel.lj_sigma', 0.1) # TODO
 pot_options_dotlj.set('kernel.lj_eps_cap', 0.001)
 
 
@@ -160,7 +180,32 @@ class SimSpaceNetwork(object):
         self.spawned = self.spawned + spawned
         self.gen_spawned.append(spawned)
         return spawned
+    def addNonbondedPairPotential(self, options, exclude_nbs=True):
+        print "Add non-bonded pair potential ..."
+        n_spawns = len(self.spawned)
+        for i in range(n_spawns):
+            for j in range(i+1, n_spawns):
+                vert_i = self.spawned[i]
+                vert_j = self.spawned[j]
+                if exclude_nbs:
+                    excluded = False
+                    if vert_i in vert_j.nbs: excluded = True
+                    if vert_j in vert_i.nbs: assert excluded
+                    if excluded: continue
+                    for nb in vert_i.nbs:
+                        if vert_j in nb.nbs: excluded = True
+                    for nb in vert_j.nbs:
+                        if vert_i in nb.nbs: assert excluded
+                    if excluded: continue
+                node_i = self.spawned[i].node
+                node_j = self.spawned[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 addPairPotential(self, options):
+        print "Add bonded pair potential ..."
         for spawn in self.spawned:
             for nb in spawn.nbs:
                 if nb.isInitialized():
@@ -169,6 +214,7 @@ class SimSpaceNetwork(object):
             print ""
         return
     def addThreePotential(self, options):
+        print "Add three-body potential ..."
         for spawn in self.spawned:
             if len(spawn.nbs) != 3: continue
             nb1 = spawn.nbs[0]
@@ -276,7 +322,7 @@ for struct in structures:
     struct_dict[struct.label] = struct
 
 # CHANGE TO WORKING DIRECTORY
-osio.cd('out.files')
+osio.cd('out.files.hexcomb')
 
 # SELECT STRUCTURE
 struct = struct_dict[select_struct_label]
@@ -287,10 +333,10 @@ struct = struct_dict[select_struct_label]
 hexcomb = SimSpaceNetwork(options)
 hexcomb.createHexComb(c=1.4, N=N_network)
 hexcomb.createNeighbours(cutoff=1.41)
+hexcomb.writeXyz()
 
 root_idx = 2*((2*N_network+1)*2*N_network + N_network)
 """
-hexcomb.writeXyz()
 hexcomb.printInfo()
 hexcomb.printConnectivity(root_idx=root_idx, levels=3)
 """
@@ -298,7 +344,7 @@ hexcomb.printConnectivity(root_idx=root_idx, levels=3)
 # SEED
 hexcomb.seed(root_idx, struct)
 
-N_generations = 2
+N_generations = 3
 for i in range(N_generations):
     gen_id = i+1
     gen_prefix = "out.gen_%d" % gen_id
@@ -311,8 +357,11 @@ for i in range(N_generations):
     
     # POTENTIALS
     hexcomb.clearPotentials()
-    hexcomb.addPairPotential(pot_options_harm)
-    hexcomb.addThreePotential(pot_options_3harm)
+    #hexcomb.addPairPotential(pot_options_harm)
+    #hexcomb.addThreePotential(pot_options_3harm)
+    hexcomb.addPairPotential(pot_options_harm_dist)
+    hexcomb.addThreePotential(pot_options_3harm_dist)
+    #hexcomb.addNonbondedPairPotential(pot_options_dotlj)
     
     # OPTIMIZE 1st GENERATION
     nodes_opt = [ vertex.node for vertex in hexcomb.gen_spawned[-1] ]
@@ -322,11 +371,11 @@ for i in range(N_generations):
         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)
+    multi_optimize_node(nodes_opt, x0_opt)
     
     # OUTPUT
     hexcomb.top.writeData(prefix=gen_prefix)
-    raw_input('...')
+    #raw_input('...')
 
 # RETURN
 osio.root()

test/test_kernelpot/kernel.py

@@ -185,6 +185,22 @@ class KernelFunctionDotHarmonic(object):
     def computeBlockDot(self, IX, return_distance=False):
         return self.kfctdot.computeBlock(IX, return_distance)
 
+class KernelFunctionDotHarmonicDist(object):
+    def __init__(self, options):
+        self.d0 = float(options.get('kernel.dotharmdist_d0'))
+        self.eps = float(options.get('kernel.dotharmdist_eps'))
+        self.kfctdot = KernelFunctionDot(options)
+        return
+    def compute(self, IX, X):
+        C = self.kfctdot.compute(IX, X)
+        D = (1.-C+self.eps)**0.5
+        return (D-self.d0)**2
+    def computeDerivativeOuter(self, IX, X):
+        C = self.kfctdot.compute(IX, X)
+        D = (1.-C+self.eps)**0.5
+        IC = self.kfctdot.computeDerivativeOuter(IX, X)
+        return 2*(D-self.d0)*0.5/D*(-1.)*IC
+        
 class KernelFunctionDotLj(object):
     def __init__(self, options):
         self.sigma = float(options.get('kernel.lj_sigma'))
@@ -224,6 +240,36 @@ class KernelFunctionDot3Harmonic(object):
         IC_B = self.kfctdot.computeDerivativeOuter(IX_B, X)
         return 2*(C_A - C_B)*(IC_A - IC_B) + 2*(C_A - C_AB)*IC_A + 2*(C_B - C_AB)*IC_B
 
+class KernelFunctionDot3HarmonicDist(object):
+    def __init__(self, options):
+        self.kfctdot = KernelFunctionDot(options)
+        self.d0 = float(options.get('kernel.dot3harmdist_d0'))
+        self.eps = float(options.get('kernel.dot3harmdist_eps'))
+    def compute(self, IX_A, IX_B, X):
+        C_A = self.kfctdot.compute(IX_A, X)
+        C_B = self.kfctdot.compute(IX_B, X)
+        # TODO Generalize this to multi-environment representation:
+        assert IX_A.shape[0] == 1
+        assert IX_B.shape[0] == 1
+        C_AB = self.kfctdot.compute(IX_A, IX_B[0])
+        D_A = (1. - C_A + self.eps)**0.5
+        D_B = (1. - C_B + self.eps)**0.5
+        D_AB = (1. - C_AB + self.eps)**0.5
+        return (D_A - D_B)**2 + (D_A - D_AB)**2 + (D_B - D_AB)**2 + (D_A - self.d0)**2 + (D_B - self.d0)**2
+    def computeDerivativeOuter(self, IX_A, IX_B, X):
+        C_A = self.kfctdot.compute(IX_A, X)
+        C_B = self.kfctdot.compute(IX_B, X)
+        # TODO Generalize this to multi-environment representation:
+        assert IX_A.shape[0] == 1
+        assert IX_B.shape[0] == 1
+        C_AB = self.kfctdot.compute(IX_A, IX_B[0])
+        D_A = (1. - C_A + self.eps)**0.5
+        D_B = (1. - C_B + self.eps)**0.5
+        D_AB = (1. - C_AB + self.eps)**0.5
+        IC_A = self.kfctdot.computeDerivativeOuter(IX_A, X)
+        IC_B = self.kfctdot.computeDerivativeOuter(IX_B, X)
+        return 2*(D_A - D_B)*(0.5/D_A*(-1.)*IC_A - 0.5/D_B*(-1.)*IC_B) + 2*(D_A - D_AB)*0.5/D_A*(-1.)*IC_A + 2*(D_B - D_AB)*0.5/D_B*(-1.)*IC_B + 2*(D_A-self.d0)*1./D_A*(-1.)*IC_A + 2*(D_B-self.d0)*1./D_B*(-1.)*IC_B
+
 KernelAdaptorFactory = {
 'generic': KernelAdaptorGeneric, 
 'global-generic': KernelAdaptorGlobalGeneric 
@@ -232,8 +278,10 @@ KernelAdaptorFactory = {
 KernelFunctionFactory = { 
 'dot': KernelFunctionDot, 
 'dot-harmonic': KernelFunctionDotHarmonic, 
+'dot-harmonic-dist': KernelFunctionDotHarmonicDist,
 'dot-lj': KernelFunctionDotLj,
-'dot-3-harmonic': KernelFunctionDot3Harmonic
+'dot-3-harmonic': KernelFunctionDot3Harmonic, 
+'dot-3-harmonic-dist': KernelFunctionDot3HarmonicDist
 }
 
 class KernelPotential(object):