Develop

.travis.yml

@@ -1,12 +1,16 @@
 language: python
-
 sudo: required
 dist: trusty
-
 group: edge
-
 deploy:
-  provider: pages
+- provider: pypi
+  user: andersx
+  password:
+    secure: mXeDkazSavPo2Cx3Drhap9VkT42pK0XCDx2vK8IVsu/3i+fNLg/RtqT/Yg1hhMNWZQ76rg8TfzGWNEDvEa7JNNVpyTrA8CtWkS5KR3Oj/bkP+n7xcQKC2npjsmrReS97mBq4Sp3fEogcsYH9jCzhhMLnlfpoVlWEFDExcqeKwvXxTqvyv2ymdx6QLu1KkOQWCzsBoa6wIqH9kbzUSjqTVKIWgNZh3UHhDDv8X5rJZUUUulnwmbsRHZSTKFRVaB4gqv3s6gqgMZzmFjHdii3O059quEu+JrpisGNJ3EHSyaDZ1hk9DWqBqXK6LhNVjtwDPgHmABXLrZCKJVudm/7fuwI7GdTSghWjhtLDlaZYSVTBD9vD/aabuM397ZLt6ci8d7AD2AoCDWtc/8EjRVjdErVYL0dO0moJAu7mtxnvPKaT0ZWWPwfUQQlPt/vh17gFWOQhF/YT226J8K/nnnBHmhgHMIR7jBr3ve/YuQGLPUcURhTUnqfUS4O4UTkbXG23rOHPRPs45+F42dNlPb0d1d67C2NOtGhuqe14sQ47rip13FtXCq6HG09WJ6tOY7TvBEu+43DzVPam4qSlS5PoVbC2DwZQKEl8YVbXLCOAzaymk+1BAIofZ201wu4M7Y4HbRIdTLM4F2maKxo52c4Dn2flRA8iX5/OTR1osDIo4qg=
+  on:
+    branch: master
+    condition: ${TRAVIS_PYTHON_VERSION} = "2.7"
+- provider: pages
   local_dir: ${TRAVIS_BUILD_DIR}/docs/build/html
   skip_cleanup: true
   repo: qmlcode/qmlcode.github.io
@@ -15,47 +19,40 @@ deploy:
   on:
     branch: master
     condition: ${TRAVIS_PYTHON_VERSION} = "2.7"
-
 python:
-  - "2.7"
-  - "3.5"
-  - "3.6"
-
+- '2.7'
+- '3.5'
+- '3.6'
 before_install:
-  - sudo apt-get update -qq
-
+- sudo apt-get update -qq
 install:
- - sudo apt-get install -qq gcc gfortran libblas-dev liblapack-dev
- - sudo apt-get install -qq gcc-4.8 gfortran-4.8
- - |
-   if [ ${TRAVIS_PYTHON_VERSION:0:1} = 3 ]; then
-       sudo apt-get install python3-numpy
-       pip3 install scipy
-       python3 setup.py build
-       python3 setup.py install
-       pip3 install sphinx
-       pip3 install sphinx-rtd-theme
-       cd ${TRAVIS_BUILD_DIR}/docs
-       make html
-   elif [ ${TRAVIS_PYTHON_VERSION} = "2.7" ]; then
-       sudo apt-get install python-numpy
-       pip2 install scipy
-       python2 setup.py build
-       python2 setup.py install
-       pip2 install sphinx
-       pip2 install sphinx-rtd-theme
-       cd ${TRAVIS_BUILD_DIR}/docs
-       make html
-   else
-       echo "ERROR: Unknown Python version."
-   fi
-
-
+- sudo apt-get install -qq gcc gfortran libblas-dev liblapack-dev
+- sudo apt-get install -qq gcc-4.8 gfortran-4.8
+- |
+  if [ ${TRAVIS_PYTHON_VERSION:0:1} = 3 ]; then
+      sudo apt-get install python3-numpy
+      pip3 install scipy
+      python3 setup.py build
+      python3 setup.py install
+      pip3 install sphinx
+      pip3 install sphinx-rtd-theme
+      cd ${TRAVIS_BUILD_DIR}/docs
+      make html
+  elif [ ${TRAVIS_PYTHON_VERSION} = "2.7" ]; then
+      sudo apt-get install python-numpy
+      pip2 install scipy
+      python2 setup.py build
+      python2 setup.py install
+      pip2 install sphinx
+      pip2 install sphinx-rtd-theme
+      cd ${TRAVIS_BUILD_DIR}/docs
+      make html
+  else
+      echo "ERROR: Unknown Python version."
+  fi
 before_script:
- - cd ${TRAVIS_BUILD_DIR}/tests/
-
+- cd ${TRAVIS_BUILD_DIR}/tests/
 script:
- - nosetests -v
-
+- nosetests -v
 notifications:
   email: false

docs/source/examples.rst

@@ -1,4 +1,4 @@
-Examples
+xamples
 --------
 
 Generating representations using the ``Compound`` class
@@ -106,19 +106,25 @@ The easiest way to calculate the kernel matrix using an explicit, local represen
 .. code:: python
 
     import numpy as np
-    from qml.wrappers import get_atomic_kernels_gaussian
+    from qml.kernels import get_local_kernels_gaussian
 
     # Assume the QM7 dataset is loaded into a list of Compound()
     for compound in qm7:
 
         # Generate the desired representation for each compound
         compound.generate_atomic_coulomb_matrix(size=23, sort="row-norm")
 
+    # Make a big array with all the atomic representations
+    X = np.concatenate([mol.representation for mol in qm7])
+
+    # Make an array with the number of atoms in each compound
+    N = np.array([mol.natoms for mol in qm7])
+
     # List of kernel-widths
     sigmas = [50.0, 100.0, 200.0]
 
     # Calculate the kernel-matrix
-    K = get_atomic_kernels_gaussian(qm7, qm7, sigmas)
+    K = get_local_kernels_gaussian(X, X, N, N, sigmas)
 
     print(K.shape)
 

docs/source/installation.rst

@@ -70,8 +70,7 @@ using ``brew``:
     # Install GCC
     brew install gcc
 
-Note the Clang Fortran compiler from brew unfortunately does not support
-OpenMP.
+Note the Fortran compiler from brew (gfortran) unfortunately does not support OpenMP.
 Therefore parallelism via OpenMP is disabled as default for MacOS systems.
 
 Additionally, we found that some users have multiple version of the ``as`` assembler - this might happen if you have GCC from e.g. brew and macports at the same time. Look for the following error:

docs/source/qml.rst

@@ -52,6 +52,13 @@ qml\.arad module
     :members:
     :show-inheritance:
 
+qml\.fchl module
+----------------
+
+.. automodule:: qml.fchl
+    :members:
+    :show-inheritance:
+
 
 qml\.wrappers module
 --------------------

mkldiscover.py

@@ -0,0 +1,86 @@
+#!/usr/bin/env python
+# MIT License
+#
+# Copyright (c) 2017 Anders Steen Christensen
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+from __future__ import print_function
+
+import os
+
+def mkl_exists(verbose=False):
+
+    # Get environment variables
+    __MKLROOT__ = os.environ.get('MKLROOT')
+    __LD_LIBRARY_PATH__ = os.environ.get('LD_LIBRARY_PATH')
+
+    # Check if $MKLROOT is set
+
+    if __MKLROOT__ is None:
+
+        if verbose: 
+            print("MKL-discover: MKLROOT was not set")
+
+        return False
+
+    else:
+
+        if verbose: 
+            print("MKL-discover: MKLROOT was set to")
+            print(__MKLROOT__)
+
+    # Check if path exists
+    mklroot_exists = os.path.isdir(__MKLROOT__)
+
+    if not mklroot_exists:
+        if verbose: 
+            print("MKL-discover: MKLROOT path does not exist")
+
+        return False
+
+    found_libmkl_rt = False
+
+    # Check if libmkl_rt.so exists below $MKLROOT
+    for dirpath, dirnames, filenames in os.walk(__MKLROOT__, followlinks=True):
+
+        if "libmkl_rt.so" in filenames:
+
+            if verbose:
+                print("MKL-discover: Found libmkl_rt.so at ", dirpath)
+
+            # Check that the dirpath where libmkl_rt.so is in $LD_LIBRARY_PATH
+            if dirpath in __LD_LIBRARY_PATH__:
+
+                if verbose:
+                    print("MKL-discover: Found %s in $LD_LIBRARY_PATH" % dirpath)
+                    print("MKL-discover: Concluding that MKL can be used.")
+                found_libmkl_rt = True
+
+
+    return found_libmkl_rt
+
+if __name__ == "__main__":
+
+    mkl_present = mkl_exists(verbose=False)
+
+    if mkl_present:
+        print("MKL could be found")
+    else:
+        print("MKL could NOT be found")

qml/alchemy.py

@@ -0,0 +1,275 @@
+# MIT License
+#
+# Copyright (c) 2017 Anders Steen Christensen and Felix A. Faber
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+from __future__ import division
+from __future__ import print_function
+
+import numpy as np
+from copy import copy
+
+PTP = {\
+         1  :[1,1] ,2:  [1,8]#Row1
+
+        ,3  :[2,1] ,4:  [2,2]#Row2\
+        ,5  :[2,3] ,6:  [2,4] ,7  :[2,5] ,8  :[2,6] ,9  :[2,7] ,10 :[2,8]\
+
+        ,11 :[3,1] ,12: [3,2]#Row3\
+        ,13 :[3,3] ,14: [3,4] ,15 :[3,5] ,16 :[3,6] ,17 :[3,7] ,18 :[3,8]\
+
+        ,19 :[4,1] ,20: [4,2]#Row4\
+        ,31 :[4,3] ,32: [4,4] ,33 :[4,5] ,34 :[4,6] ,35 :[4,7] ,36 :[4,8]\
+        ,21 :[4,9] ,22: [4,10],23 :[4,11],24 :[4,12],25 :[4,13],26 :[4,14],27 :[4,15],28 :[4,16],29 :[4,17],30 :[4,18]\
+
+        ,37 :[5,1] ,38: [5,2]#Row5\
+        ,49 :[5,3] ,50: [5,4] ,51 :[5,5] ,52 :[5,6] ,53 :[5,7] ,54 :[5,8]\
+        ,39 :[5,9] ,40: [5,10],41 :[5,11],42 :[5,12],43 :[5,13],44 :[5,14],45 :[5,15],46 :[5,16],47 :[5,17],48 :[5,18]\
+
+        ,55 :[6,1] ,56: [6,2]#Row6\
+        ,81 :[6,3] ,82: [6,4] ,83 :[6,5] ,84 :[6,6] ,85 :[6,7] ,86 :[6,8]
+               ,72: [6,10],73 :[6,11],74 :[6,12],75 :[6,13],76 :[6,14],77 :[6,15],78 :[6,16],79 :[6,17],80 :[6,18]\
+        ,57 :[6,19],58: [6,20],59 :[6,21],60 :[6,22],61 :[6,23],62 :[6,24],63 :[6,25],64 :[6,26],65 :[6,27],66 :[6,28],67 :[6,29],68 :[6,30],69 :[6,31],70 :[6,32],71 :[6,33]\
+
+        ,87 :[7,1] ,88: [7,2]#Row7\
+        ,113:[7,3] ,114:[7,4] ,115:[7,5] ,116:[7,6] ,117:[7,7] ,118:[7,8]\
+               ,104:[7,10],105:[7,11],106:[7,12],107:[7,13],108:[7,14],109:[7,15],110:[7,16],111:[7,17],112:[7,18]\
+        ,89 :[7,19],90: [7,20],91 :[7,21],92 :[7,22],93 :[7,23],94 :[7,24],95 :[7,25],96 :[7,26],97 :[7,27],98 :[7,28],99 :[7,29],100:[7,30],101:[7,31],101:[7,32],102:[7,14],103:[7,33]}
+
+QtNm = {
+    #Row1
+    1  :[1,0,0,1./2.]
+    ,2:  [1,0,0,-1./2.]
+
+
+    #Row2
+    ,3  :[2,0,0,1./2.]
+    ,4: [2,0,0,-1./2.]
+
+    ,5  :[2,-1,1,1./2.],   6: [2,0,1,1./2.]  , 7  : [2,1,1,1./2.]
+    ,8  : [2,-1,1,-1./2.] ,9: [2,0,1,-1./2.] ,10 :[2,1,1,-1./2.]
+
+
+    #Row3
+    ,11 :[3,0,0,1./2.]
+    ,12: [3,0,0,-1./2.]
+
+    ,13 :[3,-1,1,1./2.] ,  14: [3,0,1,1./2.] ,   15 :[3,1,1,1./2.]
+    ,16 :[3,-1,1,-1./2.]  ,17 :[3,0,1,-1./2.]  ,18 :[3,1,1,-1./2.]
+
+
+    #Row3
+    ,19 :[4,0,0,1./2.]
+    ,20: [4,0,0,-1./2.]
+
+    ,31 :[4,-1,2,1./2.] , 32: [4,0,1,1./2.]  , 33 :[4,1,1,1./2.]
+    ,34 :[4,-1,1,-1./2.] ,35 :[4,0,1,-1./2.] ,36 :[4,1,1,-1./2.]
+
+    ,21 :[4,-2,2,1./2.],  22:[4,-1,2,1./2.],  23 :[4,0,2,1./2.], 24 :[4,1,2,1./2.], 25 :[4,2,2,1./2.]
+    ,26 :[4,-2,2,-1./2.], 27:[4,-1,2,-1./2.], 28 :[4,0,2,-1./2.],29 :[4,1,2,-1./2.],30 :[4,2,2,-1./2.]
+
+
+    #Row5
+    ,37 :[5,0,0,1./2.]
+    ,38: [5,0,0,-1./2.]
+
+    ,49 :[5,-1,1,1./2.] ,  50: [5,0,1,1./2.]  ,  51 :[5,1,1,1./2.]
+    ,52 :[5,-1,1,-1./2.]  ,53 :[5,0,1,-1./2.]  ,54 :[5,1,1,-1./2.]
+
+
+    ,39 :[5,-2,2,1./2.], 40:[5,-1,2,1./2.],   41 :[5,0,2,1./2.], 42 :[5,1,2,1./2.], 43 :[5,2,2,1./2.]
+    ,44 :[5,-2,2,-1./2.],45 :[5,-1,2,-1./2.],46 :[5,0,2,-1./2.],47 :[5,1,2,-1./2.],48 :[5,2,2,-1./2.]
+
+
+    #Row6
+    ,55 :[6,0,0,1./2.]
+    ,56: [6,0,0,-1./2.]
+
+    ,81 :[6,-1,1,1./2.] ,82: [6,0,1,1./2.] ,83: [6,1,1,1./2.]
+    ,84 :[6,-1,1,-1./2.] ,85 :[6,0,1,-1./2.] ,86 :[6,1,1,-1./2.]
+
+    ,71 :[6,-2,2,1./2.], 72: [6,-1,2,1./2.],  73 :[6,0,2,1./2.], 74 :[6,1,2,1./2.], 75 :[6,2,2,1./2.]
+    ,76 :[6,-2,2,-1./2.],77 :[6,-1,2,-1./2.],78 :[6,0,2,-1./2.],79 :[6,1,2,-1./2.],80 :[6,2,2,-1./2.]
+
+    ,57 :[6,-3,3,1./2.], 58: [6,-2,3,1./2.],  59 :[6,-1,3,1./2.], 60 :[6,0,3,1./2.],  61 :[6,1,3,1./2.], 62 :[6,2,3,1./2.], 63 :[6,3,3,1./2.]
+    ,64 :[6,-3,3,-1./2.],65 :[6,-2,3,-1./2.],66 :[6,-1,3,-1./2.],67 :[6,0,3,-1./2.],68 :[6,1,3,-1./2.],69 :[6,2,3,-1./2.],70 :[6,3,3,-1./2.]
+
+
+    #Row7
+    ,87 :[7,0,0,1./2.]
+    ,88: [7,0,0,-1./2.]
+
+    ,113:[7,-1,1,1./2.] , 114:[7,0,1,1./2.] ,  115:[7,1,1,1./2.]
+    ,116:[7,-1,1,-1./2.] ,117:[7,0,1,-1./2.] ,118:[7,1,1,-1./2.]
+
+    ,103:[7,-2,2,1./2.], 104:[7,-1,2,1./2.],  105:[7,0,2,1./2.], 106:[7,1,2,1./2.],  107:[7,2,2,1./2.]
+    ,108:[7,-2,2,-1./2.],109:[7,-1,2,-1./2.],110:[7,0,2,-1./2.],111:[7,1,2,-1./2.],112:[7,2,2,-1./2.]
+
+    ,89 :[7,-3,3,1./2.], 90: [7,-2,3,1./2.],  91 :[7,-1,3,1./2.], 92 :[7,0,3,1./2.],  93 :[7,1,3,1./2.],  94 :[7,2,3,1./2.],  95 :[7,3,3,1./2.]
+    ,96 :[7,-3,3,-1./2.],97 :[7,-2,3,-1./2.],98 :[7,-1,3,-1./2.],99 :[7,0,3,-1./2.],100:[7,1,3,-1./2.],101:[7,2,3,-1./2.],102:[7,3,3,-1./2.]}
+
+
+
+
+
+
+
+def QNum_distance(a,b, n_width, m_width, l_width, s_width):
+    """ Calculate stochiometric distance
+        a -- nuclear charge of element a
+        b -- nuclear charge of element b
+        r_width -- sigma in row-direction
+        c_width -- sigma in column direction
+    """
+
+    na = QtNm[int(a)][0]
+    nb = QtNm[int(b)][0]
+
+    ma = QtNm[int(a)][1]
+    mb = QtNm[int(b)][1]
+
+    la = QtNm[int(a)][2]
+    lb = QtNm[int(b)][2]
+
+    sa = QtNm[int(a)][3]
+    sb = QtNm[int(b)][3]
+
+    return  np.exp(-(na - nb)**2/(4 * n_width**2)
+                   -(ma - mb)**2/(4 * m_width**2)
+                   -(la - lb)**2/(4 * l_width**2)
+                   -(sa - sb)**2/(4 * s_width**2))
+
+def gen_QNum_distances(emax=100, n_width = 0.001, m_width = 0.001, l_width = 0.001, s_width = 0.001):
+    """ Generate stochiometric ditance matrix
+        emax -- Largest element
+        r_width -- sigma in row-direction
+        c_width -- sigma in column direction
+    """
+
+    pd = np.zeros((emax,emax))
+
+    for i in range(emax):
+        for j in range(emax):
+
+            pd[i,j] = QNum_distance(i+1, j+1, n_width, m_width, l_width, s_width)
+
+    return pd
+
+def periodic_distance(a, b, r_width, c_width):
+    """ Calculate stochiometric distance
+
+        a -- nuclear charge of element a
+        b -- nuclear charge of element b
+        r_width -- sigma in row-direction
+        c_width -- sigma in column direction
+    """
+
+    ra = PTP[int(a)][0]
+    rb = PTP[int(b)][0]
+    ca = PTP[int(a)][1]
+    cb = PTP[int(b)][1]
+
+    # return (r_width**2 * c_width**2) / ((r_width**2 + (ra - rb)**2) * (c_width**2 + (ca - cb)**2))
+
+    return np.exp(-(ra - rb)**2/(4 * r_width**2)-(ca - cb)**2/(4 * c_width**2))
+
+
+def gen_pd(emax=100, r_width=0.001, c_width=0.001):
+    """ Generate stochiometric ditance matrix
+
+        emax -- Largest element
+        r_width -- sigma in row-direction
+        c_width -- sigma in column direction
+    """
+
+    pd = np.zeros((emax,emax))
+
+    for i in range(emax):
+        for j in range(emax):
+
+            pd[i,j] = periodic_distance(i+1, j+1, r_width, c_width)
+
+    return pd
+
+
+def gen_custom(e_vec, emax=100):
+    """ Generate stochiometric ditance matrix
+        emax -- Largest element
+        r_width -- sigma in row-direction
+        c_width -- sigma in column direction
+    """
+
+    def check_if_unique(iterator):
+        return len(set(iterator)) == 1
+
+    num_dims = []
+
+    for k,v in e_vec.items():
+        assert isinstance(k,int), 'Error! Keys need to be int'
+        num_dims.append(len(v))
+
+    assert check_if_unique(num_dims), 'Error! Unequal number of dimensions'
+
+
+    tmp = np.zeros((emax,num_dims[0]))
+
+    for k,v in e_vec.items():
+        tmp[k,:] = copy(v)
+    pd = np.dot(tmp,tmp.T)
+
+    return pd
+
+def get_alchemy(alchemy, emax=100, r_width=0.001, c_width=0.001, elemental_vectors={}, \
+                n_width = 0.001, m_width = 0.001, l_width = 0.001, s_width = 0.001):
+
+    if (alchemy == "off"):
+
+        pd = np.eye(emax)
+        doalchemy = False
+
+        return doalchemy, pd
+
+    elif (alchemy == "periodic-table"):
+
+        pd = gen_pd(emax=emax, r_width=r_width, c_width=c_width)
+        doalchemy = True
+
+        return doalchemy, pd
+
+
+    elif (alchemy == "quantum-numbers"):
+        pd = gen_QNum_distances(emax=emax, n_width = n_width, m_width = m_width,
+                                          l_width = l_width, s_width = s_width)
+        doalchemy = True
+
+        return doalchemy, pd
+
+    elif (alchemy == "custom"):
+
+        pd = gen_custom(elemental_vectors,emax)
+        doalchemy = True
+
+
+        return doalchemy, pd
+
+    else:
+
+        print("QML ERROR: Unknown alchemical method specified:", alchemy)
+        exit(1)

qml/arad.py

@@ -24,6 +24,9 @@
 
 import numpy as np
 
+from .farad_kernels import fget_global_kernels_arad
+from .farad_kernels import fget_global_symmetric_kernels_arad
+
 from .farad_kernels import fget_local_kernels_arad
 from .farad_kernels import fget_local_symmetric_kernels_arad
 
@@ -131,6 +134,85 @@ def generate_arad_representation(coordinates, nuclear_charges, size=23, cut_dist
 
     return M
 
+def get_global_kernels_arad(X1, X2, sigmas, 
+        width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5):
+    """ Calculates the global Gaussian kernel matrix K for atomic ARAD
+        descriptors for a list of different sigmas. Each kernel element
+        is the sum of all kernel elements between pairs of atoms in two molecules.
+
+        K is calculated using an OpenMP parallel Fortran routine.
+
+        :param X1: ARAD descriptors for molecules in set 1.
+        :type X1: numpy array
+        :param X2: Array of ARAD descriptors for molecules in set 2.
+        :type X2: numpy array
+        :param sigmas: List of sigmas for which to calculate the Kernel matrices.
+        :type sigmas: list
+
+        :return: The kernel matrices for each sigma - shape (number_sigmas, number_molecules1, number_molecules2)
+        :rtype: numpy array
+    """
+
+    amax = X1.shape[1]
+
+    assert X1.shape[3] == amax, "ERROR: Check ARAD decriptor sizes! code = 1"
+    assert X2.shape[1] == amax, "ERROR: Check ARAD decriptor sizes! code = 2"
+    assert X2.shape[3] == amax, "ERROR: Check ARAD decriptor sizes! code = 3"
+
+    nm1 = X1.shape[0]
+    nm2 = X2.shape[0]
+
+    N1 = np.empty(nm1, dtype = np.int32)
+    Z1_arad = np.zeros((nm1, amax, 2))
+    for i in range(nm1):
+        N1[i] = len(np.where(X1[i,:,2,0] > 0)[0])
+        Z1_arad[i] = X1[i,:,1:3,0]
+
+    N2 = np.empty(nm2, dtype = np.int32)
+    Z2_arad = np.zeros((nm2, amax, 2))
+    for i in range(nm2):
+        N2[i] = len(np.where(X2[i,:,2,0] > 0)[0])
+        Z2_arad[i] = X2[i,:,1:3,0]
+
+    sigmas = np.array(sigmas)
+    nsigmas = sigmas.size
+
+    return fget_global_kernels_arad(X1, X2, Z1_arad, Z2_arad, N1, N2, sigmas, 
+                nm1, nm2, nsigmas, width, cut_distance, r_width, c_width)
+
+
+def get_global_symmetric_kernels_arad(X1, sigmas, 
+    width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5):
+    """ Calculates the global Gaussian kernel matrix K for atomic ARAD
+        descriptors for a list of different sigmas. Each kernel element
+        is the sum of all kernel elements between pairs of atoms in two molecules.
+
+        K is calculated using an OpenMP parallel Fortran routine.
+
+        :param X1: ARAD descriptors for molecules in set 1.
+        :type X1: numpy array
+        :param sigmas: List of sigmas for which to calculate the Kernel matrices.
+        :type sigmas: list
+
+        :return: The kernel matrices for each sigma - shape (number_sigmas, number_molecules1, number_molecules1)
+        :rtype: numpy array
+    """
+
+    nm1 = X1.shape[0]
+    amax = X1.shape[1]
+
+    N1 = np.empty(nm1, dtype = np.int32)
+    Z1_arad = np.zeros((nm1, amax, 2))
+    for i in range(nm1):
+        N1[i] = len(np.where(X1[i,:,2,0] > 0)[0])
+        Z1_arad[i] = X1[i,:,1:3,0]
+
+    sigmas = np.array(sigmas)
+    nsigmas = sigmas.size
+
+    return fget_global_symmetric_kernels_arad(X1, Z1_arad, N1, sigmas, 
+                nm1, nsigmas, width, cut_distance, r_width, c_width)
+
 
 def get_local_kernels_arad(X1, X2, sigmas, 
         width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5):

qml/compound.py

@@ -1,6 +1,6 @@
 # MIT License
 #
-# Copyright (c) 2016 Anders Steen Christensen, Felix Faber
+# Copyright (c) 2016-2017 Anders Steen Christensen, Felix Faber, Lars Andersen Bratholm
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
@@ -70,7 +70,7 @@ def __init__(self, xyz = None):
         if xyz is not None:
             self.read_xyz(xyz)
 
-    def generate_coulomb_matrix(self, size = 23, sorting = "row-norm"):
+    def generate_coulomb_matrix(self, size = 23, sorting = "row-norm", indices = None):
         """ Creates a Coulomb Matrix representation of a molecule.
             A matrix :math:`M` is constructed with elements
 
@@ -132,7 +132,8 @@ def generate_eigenvalue_coulomb_matrix(self, size = 23):
                 self.nuclear_charges, self.coordinates, size = size)
 
     def generate_atomic_coulomb_matrix(self, size = 23, sorting = "row-norm", 
-            central_cutoff = 1e6, central_decay = -1, interaction_cutoff = 1e6, interaction_decay = -1):
+            central_cutoff = 1e6, central_decay = -1, interaction_cutoff = 1e6, interaction_decay = -1,
+            indices = None):
         """ Creates a Coulomb Matrix representation of the local environment of a central atom.
             For each central atom :math:`k`, a matrix :math:`M` is constructed with elements
 
@@ -178,6 +179,11 @@ def generate_atomic_coulomb_matrix(self, size = 23, sorting = "row-norm",
 
             The upper triangular of M, including the diagonal, is concatenated to a 1D
             vector representation.
+
+            The representation can be calculated for a subset by either specifying
+            ``indices = [0,1,...]``, where :math:`[0,1,...]` are the requested atom indices,
+            or by specifying ``indices = 'C'`` to only calculate central carbon atoms.
+
             The representation is calculated using an OpenMP parallel Fortran routine.
 
             :param size: The size of the largest molecule supported by the representation
@@ -194,6 +200,8 @@ def generate_atomic_coulomb_matrix(self, size = 23, sorting = "row-norm",
             :type interaction_cutoff: float
             :param interaction_decay: The distance over which the the coulomb interaction decays from full to none
             :type interaction_decay: float
+            :param indices: Subset indices or atomtype
+            :type indices: Nonetype/array/string
 
 
             :return: nD representation - shape (:math:`N_{atoms}`, size(size+1)/2)
@@ -202,11 +210,11 @@ def generate_atomic_coulomb_matrix(self, size = 23, sorting = "row-norm",
 
 
         self.representation = generate_atomic_coulomb_matrix(
-            self.nuclear_charges, self.coordinates, size = size,
+            self.nuclear_charges, self.coordinates, size = size, indices = indices,
             sorting = sorting, central_cutoff = central_cutoff, central_decay = central_decay,
             interaction_cutoff = interaction_cutoff, interaction_decay = interaction_decay)
 
-    def generate_bob(self, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}):
+    def generate_bob(self, size=23, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}):
         """ Creates a Bag of Bonds (BOB) representation of a molecule.
             The representation expands on the coulomb matrix representation.
             For each element a bag (vector) is constructed for self interactions
@@ -256,6 +264,7 @@ def generate_slatm(self, mbtypes,
 
     A version that works for periodic boundary conditions will be released soon.
 
+
     NOTE: You will need to run the ``get_slatm_mbtypes()`` function to get the ``mbtypes`` input (or generate it manually).
 
     :param mbtypes: Many-body types for the whole dataset, including 1-, 2- and 3-body types. Could be obtained by calling ``get_slatm_mbtypes()``.

qml/fcho_solve.f90

@@ -34,35 +34,117 @@ subroutine fcho_solve(A,y,x)
 
     call dpotrf("U", na, A, na, info)
     if (info > 0) then
-        write (*,*) "WARNING: Cholesky decomposition DPOTRF() exited with error code:", info
+        write (*,*) "WARNING: Error in LAPACK Cholesky decomposition DPOTRF()."
+        write (*,*) "WARNING: The", info, "-th leading order is not positive definite."
+    else if (info < 0) then
+        write (*,*) "WARNING: Error in LAPACK Cholesky decomposition DPOTRF()."
+        write (*,*) "WARNING: The", -info, "-th argument had an illegal value."
     endif
 
     x(:na) = y(:na)
 
     call dpotrs("U", na, 1, A, na, x, na, info)
-    if (info > 0) then
-        write (*,*) "WARNING: Cholesky solve DPOTRS() exited with error code:", info
+    if (info < 0) then
+        write (*,*) "WARNING: Error in LAPACK Cholesky solver DPOTRS()."
+        write (*,*) "WARNING: The", -info, "-th argument had an illegal value."
     endif
 
 end subroutine fcho_solve
 
-! subroutine fcho_invert(A)
-! 
-!     implicit none
-! 
-!     double precision, dimension(:,:), intent(inout) :: A
-!     integer :: info, na
-! 
-!     na = size(A, dim=1)
-! 
-!     call dpotrf("L", na, A , na, info)
-!     if (info > 0) then
-!         write (*,*) "WARNING: Cholesky decomposition DPOTRF() exited with error code:", info
-!     endif
-! 
-!     call dpotri("L", na, A , na, info )
-!     if (info > 0) then
-!         write (*,*) "WARNING: Cholesky inversion DPOTRI() exited with error code:", info
-!     endif
-! 
-! end subroutine fcho_invert
+subroutine fcho_invert(A)
+
+    implicit none
+
+    double precision, dimension(:,:), intent(inout) :: A
+    integer :: info, na
+
+    na = size(A, dim=1)
+
+    call dpotrf("L", na, A , na, info)
+    if (info > 0) then
+        write (*,*) "WARNING: Cholesky decomposition DPOTRF() exited with error code:", info
+    endif
+
+    call dpotri("L", na, A , na, info )
+    if (info > 0) then
+        write (*,*) "WARNING: Cholesky inversion DPOTRI() exited with error code:", info
+    endif
+
+end subroutine fcho_invert
+
+
+subroutine fbkf_invert(A)
+
+    implicit none
+
+    double precision, dimension(:,:), intent(inout) :: A
+    integer :: info, na, nb
+
+    integer, dimension(size(A,1)) :: ipiv   ! pivot indices
+    integer :: ilaenv
+
+    integer :: lwork
+
+    double precision, allocatable, dimension(:) :: work
+
+    na = size(A, dim=1)
+
+    nb = ilaenv( 1, 'DSYTRF', "L", na, -1, -1, -1 )
+
+    lwork = na*nb
+
+    allocate(work(lwork))
+
+    ! call dpotrf("L", na, A , na, info)
+    call dsytrf("L", na, A, na, ipiv, work, lwork, info)
+    if (info > 0) then
+        write (*,*) "WARNING: Bunch-Kaufman factorization DSYTRI() exited with error code:", info
+    endif
+
+    ! call dpotri("L", na, A , na, info )
+    call dsytri( "L", na, a, na, ipiv, work, info )
+    if (info > 0) then
+        write (*,*) "WARNING: BKF inversion DPOTRI() exited with error code:", info
+    endif
+
+    deallocate(work)
+
+end subroutine fbkf_invert
+
+
+subroutine fbkf_solve(A,y,x)
+
+    implicit none
+
+    double precision, dimension(:,:), intent(in) :: A
+    double precision, dimension(:), intent(in) :: y
+    double precision, dimension(:), intent(inout) :: x
+
+    double precision, allocatable, dimension(:) :: work
+    integer :: ilaenv
+
+    integer, dimension(size(A,1)) :: ipiv   ! pivot indices
+    integer :: info, na, nb, lwork
+
+    na = size(A, dim=1)
+
+    nb = ilaenv( 1, 'DSYTRF', "L", na, -1, -1, -1 )
+
+    lwork = na*nb
+    allocate(work(lwork))
+
+    call dsytrf("L", na, A, na, ipiv, work, lwork, info)
+    if (info > 0) then
+        write (*,*) "WARNING: Bunch-Kaufman factorization DSYTRI() exited with error code:", info
+    endif
+
+    x(:na) = y(:na)
+
+    call dsytrs("L", na, 1, A, na, ipiv, x, na, info )
+
+    if (info > 0) then
+        write (*,*) "WARNING: Bunch-Kaufman solver DSYTRS() exited with error code:", info
+    endif
+
+    deallocate(work)
+end subroutine fbkf_solve

qml/fkernels.f90

@@ -1,6 +1,6 @@
 ! MIT License
 !
-! Copyright (c) 2016 Anders Steen Christensen
+! Copyright (c) 2016 Anders Steen Christensen, Lars A. Bratholm, Felix A. Faber
 !
 ! Permission is hereby granted, free of charge, to any person obtaining a copy
 ! of this software and associated documentation files (the "Software"), to deal
@@ -20,6 +20,181 @@
 ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 ! SOFTWARE.
 
+subroutine fget_local_kernels_gaussian(q1, q2, n1, n2, sigmas, &
+        & nm1, nm2, nsigmas, kernels)
+
+    implicit none
+
+    double precision, dimension(:,:), intent(in) :: q1
+    double precision, dimension(:,:), intent(in) :: q2
+
+    ! List of numbers of atoms in each molecule
+    integer, dimension(:), intent(in) :: n1
+    integer, dimension(:), intent(in) :: n2
+
+    ! Sigma in the Gaussian kernel
+    double precision, dimension(:), intent(in) :: sigmas
+
+    ! Number of molecules
+    integer, intent(in) :: nm1
+    integer, intent(in) :: nm2
+
+    ! Number of sigmas
+    integer, intent(in) :: nsigmas
+
+    ! -1.0 / sigma^2 for use in the kernel
+    double precision, dimension(nsigmas) :: inv_sigma2
+
+    ! Resulting alpha vector
+    double precision, dimension(nsigmas,nm1,nm2), intent(out) :: kernels
+
+    ! Internal counters
+    integer :: a, b, i, j, k, ni, nj
+
+    ! Temporary variables necessary for parallelization
+    double precision, allocatable, dimension(:,:) :: atomic_distance
+
+    integer, allocatable, dimension(:) :: i_starts
+    integer, allocatable, dimension(:) :: j_starts
+
+    allocate(i_starts(nm1))
+    allocate(j_starts(nm2))
+
+    !$OMP PARALLEL DO
+    do i = 1, nm1
+        i_starts(i) = sum(n1(:i)) - n1(i)
+    enddo
+    !$OMP END PARALLEL DO
+
+    !$OMP PARALLEL DO
+    do j = 1, nm2
+        j_starts(j) = sum(n2(:j)) - n2(j)
+    enddo
+    !$OMP END PARALLEL DO
+
+    inv_sigma2(:) = -0.5d0 / (sigmas(:))**2
+    kernels(:,:,:) = 0.0d0
+
+    allocate(atomic_distance(maxval(n1), maxval(n2)))
+    atomic_distance(:,:) = 0.0d0
+
+    !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj)
+     do a = 1, nm1
+        ni = n1(a)
+        do b = 1, nm2
+            nj = n2(b)
+
+            atomic_distance(:,:) = 0.0d0
+            do i = 1, ni
+                do j = 1, nj
+
+                    atomic_distance(i, j) = sum((q1(:,i + i_starts(a)) - q2(:,j + j_starts(b)))**2)
+
+                enddo
+            enddo
+
+            do k = 1, nsigmas
+                kernels(k, a, b) =  sum(exp(atomic_distance(:ni,:nj) * inv_sigma2(k)))
+            enddo
+
+        enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    deallocate(atomic_distance)
+    deallocate(i_starts)
+    deallocate(j_starts)
+
+end subroutine fget_local_kernels_gaussian
+
+subroutine fget_local_kernels_laplacian(q1, q2, n1, n2, sigmas, &
+        & nm1, nm2, nsigmas, kernels)
+
+    implicit none
+
+    double precision, dimension(:,:), intent(in) :: q1
+    double precision, dimension(:,:), intent(in) :: q2
+
+    ! List of numbers of atoms in each molecule
+    integer, dimension(:), intent(in) :: n1
+    integer, dimension(:), intent(in) :: n2
+
+    ! Sigma in the Gaussian kernel
+    double precision, dimension(:), intent(in) :: sigmas
+
+    ! Number of molecules
+    integer, intent(in) :: nm1
+    integer, intent(in) :: nm2
+
+    ! Number of sigmas
+    integer, intent(in) :: nsigmas
+
+    ! -1.0 / sigma^2 for use in the kernel
+    double precision, dimension(nsigmas) :: inv_sigma2
+
+    ! Resulting alpha vector
+    double precision, dimension(nsigmas,nm1,nm2), intent(out) :: kernels
+
+    ! Internal counters
+    integer :: a, b, i, j, k, ni, nj
+
+    ! Temporary variables necessary for parallelization
+    double precision, allocatable, dimension(:,:) :: atomic_distance
+
+    integer, allocatable, dimension(:) :: i_starts
+    integer, allocatable, dimension(:) :: j_starts
+
+    allocate(i_starts(nm1))
+    allocate(j_starts(nm2))
+
+    !$OMP PARALLEL DO
+    do i = 1, nm1
+        i_starts(i) = sum(n1(:i)) - n1(i)
+    enddo
+    !$OMP END PARALLEL DO
+
+    !$OMP PARALLEL DO
+    do j = 1, nm2
+        j_starts(j) = sum(n2(:j)) - n2(j)
+    enddo
+    !$OMP END PARALLEL DO
+
+    inv_sigma2(:) = -1.0d0 / sigmas(:)
+    kernels(:,:,:) = 0.0d0
+
+    allocate(atomic_distance(maxval(n1), maxval(n2)))
+    atomic_distance(:,:) = 0.0d0
+
+    !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj)
+     do a = 1, nm1
+        ni = n1(a)
+        do b = 1, nm2
+            nj = n2(b)
+
+            atomic_distance(:,:) = 0.0d0
+            do i = 1, ni
+                do j = 1, nj
+
+                    atomic_distance(i, j) = sum(abs(q1(:,i + i_starts(a)) - q2(:,j + j_starts(b))))
+
+                enddo
+            enddo
+
+
+            do k = 1, nsigmas
+                kernels(k, a, b) =  sum(exp(atomic_distance(:ni,:nj) * inv_sigma2(k)))
+            enddo
+
+        enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    deallocate(atomic_distance)
+    deallocate(i_starts)
+    deallocate(j_starts)
+
+end subroutine fget_local_kernels_laplacian
+
 subroutine fget_vector_kernels_laplacian(q1, q2, n1, n2, sigmas, &
         & nm1, nm2, nsigmas, kernels)
 
@@ -225,6 +400,30 @@ subroutine flaplacian_kernel(a, na, b, nb, k, sigma)
 end subroutine flaplacian_kernel
 
 
+subroutine flinear_kernel(a, na, b, nb, k)
+
+    implicit none
+
+    double precision, dimension(:,:), intent(in) :: a
+    double precision, dimension(:,:), intent(in) :: b
+
+    integer, intent(in) :: na, nb
+
+    double precision, dimension(:,:), intent(inout) :: k
+
+    integer :: i, j
+
+!$OMP PARALLEL DO
+    do i = 1, nb
+        do j = 1, na
+            k(j,i) = dot_product(a(:,j), b(:,i))
+        enddo
+    enddo
+!$OMP END PARALLEL DO
+
+end subroutine flinear_kernel
+
+
 subroutine fmatern_kernel_l2(a, na, b, nb, k, sigma, order)
 
     implicit none

qml/frepresentations.f90

@@ -1,6 +1,6 @@
 ! MIT License
 !
-! Copyright (c) 2016 Anders Steen Christensen
+! Copyright (c) 2016-2017 Anders Steen Christensen, Lars Andersen Bratholm
 !
 ! Permission is hereby granted, free of charge, to any person obtaining a copy
 ! of this software and associated documentation files (the "Software"), to deal
@@ -204,18 +204,21 @@ subroutine fgenerate_unsorted_coulomb_matrix(atomic_charges, coordinates, nmax,
 
 end subroutine fgenerate_unsorted_coulomb_matrix
 
-subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, nmax, &
-        & cent_cutoff, cent_decay, int_cutoff, int_decay, cm)
+subroutine fgenerate_local_coulomb_matrix(central_atom_indices, central_natoms, &
+        & atomic_charges, coordinates, natoms, nmax, cent_cutoff, cent_decay, &
+        & int_cutoff, int_decay, cm)
 
     implicit none
 
+    integer, intent(in) :: central_natoms
+    integer, dimension(:), intent(in) :: central_atom_indices
     double precision, dimension(:), intent(in) :: atomic_charges
     double precision, dimension(:,:), intent(in) :: coordinates
     integer,intent(in) :: natoms
     integer, intent(in) :: nmax
     double precision, intent(inout) :: cent_cutoff, cent_decay, int_cutoff, int_decay
 
-    double precision, dimension(natoms, ((nmax + 1) * nmax) / 2), intent(out):: cm
+    double precision, dimension(central_natoms, ((nmax + 1) * nmax) / 2), intent(out):: cm
 
     integer :: idx
 
@@ -230,12 +233,12 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
 
     double precision, allocatable, dimension(:, :, :) :: pair_distance_matrix
     double precision, allocatable, dimension(:, :) :: distance_matrix
-    double precision, allocatable, dimension(:, :) :: distance_matrix_tmp
 
-    integer i, j, m, n, k
+    integer i, j, m, n, k, l
 
     double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
 
+
     if (size(coordinates, dim=1) /= size(atomic_charges, dim=1)) then
         write(*,*) "ERROR: Coulomb matrix generation"
         write(*,*) size(coordinates, dim=1), "coordinates, but", &
@@ -287,25 +290,28 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
     enddo
     !$OMP END PARALLEL DO
 
-    do i = 1, natoms
-        if (cutoff_count(i) > nmax) then
+    do i = 1, central_natoms
+        j = central_atom_indices(i)
+        if (cutoff_count(j) > nmax) then
             write(*,*) "ERROR: Coulomb matrix generation"
             write(*,*) nmax, "size set, but", &
-                & cutoff_count(i), "size needed!"
+                & cutoff_count(j), "size needed!"
             stop
         endif
     enddo
 
     ! Allocate temporary
-    allocate(pair_distance_matrix(natoms, natoms, natoms))
-    allocate(row_norms(natoms, natoms))
+    allocate(pair_distance_matrix(natoms, natoms, central_natoms))
+    allocate(row_norms(natoms, central_natoms))
 
     pair_distance_matrix = 0.0d0
     row_norms = 0.0d0
 
-    !$OMP PARALLEL DO PRIVATE(pair_norm, prefactor) REDUCTION(+:row_norms) COLLAPSE(2)
+
+    !$OMP PARALLEL DO PRIVATE(pair_norm, prefactor, k) REDUCTION(+:row_norms) COLLAPSE(2)
     do i = 1, natoms
-        do k = 1, natoms
+        do l = 1, central_natoms
+            k = central_atom_indices(l)
             ! self interaction
             if (distance_matrix(i,k) > cent_cutoff) then
                 cycle
@@ -318,8 +324,8 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
             endif
 
             pair_norm = prefactor * prefactor * 0.5d0 * atomic_charges(i) ** 2.4d0
-            pair_distance_matrix(i,i,k) = pair_norm
-            row_norms(i,k) = row_norms(i,k) + pair_norm * pair_norm
+            pair_distance_matrix(i,i,l) = pair_norm
+            row_norms(i,l) = row_norms(i,l) + pair_norm * pair_norm
 
             do j = i+1, natoms
                 if (distance_matrix(j,k) > cent_cutoff) then
@@ -344,32 +350,34 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
                         & * (distance_matrix(j,k) - cent_cutoff + cent_decay) / cent_decay) + 1)
                 endif
 
-                pair_distance_matrix(i, j, k) = pair_norm
-                pair_distance_matrix(j, i, k) = pair_norm
+                pair_distance_matrix(i, j, l) = pair_norm
+                pair_distance_matrix(j, i, l) = pair_norm
                 pair_norm = pair_norm * pair_norm
-                row_norms(i,k) = row_norms(i,k) + pair_norm
-                row_norms(j,k) = row_norms(j,k) + pair_norm
+                row_norms(i,l) = row_norms(i,l) + pair_norm
+                row_norms(j,l) = row_norms(j,l) + pair_norm
             enddo
         enddo
     enddo
     !$OMP END PARALLEL DO
 
     ! Allocate temporary
-    allocate(sorted_atoms_all(natoms, natoms))
+    allocate(sorted_atoms_all(natoms, central_natoms))
 
-    !$OMP PARALLEL DO
-        do k = 1, natoms
-            row_norms(k,k) = huge_double
+    !$OMP PARALLEL DO PRIVATE(k)
+        do l = 1, central_natoms
+            k = central_atom_indices(l)
+            row_norms(k,l) = huge_double
         enddo
     !$OMP END PARALLEL DO
 
-    !$OMP PARALLEL DO PRIVATE(j)
-        do k = 1, natoms
+    !$OMP PARALLEL DO PRIVATE(j,k)
+        do l = 1, central_natoms
+            k = central_atom_indices(l)
             !$OMP CRITICAL
                 do i = 1, cutoff_count(k)
-                    j = maxloc(row_norms(:,k), dim=1)
-                    sorted_atoms_all(i, k) = j
-                    row_norms(j,k) = 0.0d0
+                    j = maxloc(row_norms(:,l), dim=1)
+                    sorted_atoms_all(i, l) = j
+                    row_norms(j,l) = 0.0d0
                 enddo
             !$OMP END CRITICAL
         enddo
@@ -383,14 +391,15 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
     ! Fill coulomb matrix according to sorted row-2-norms
     cm = 0.0d0
 
-    !$OMP PARALLEL DO PRIVATE(i, j, idx)
-        do k = 1, natoms
+    !$OMP PARALLEL DO PRIVATE(i, j, k, idx)
+        do l = 1, central_natoms
+            k = central_atom_indices(l)
             do m = 1, cutoff_count(k)
-                i = sorted_atoms_all(m, k)
+                i = sorted_atoms_all(m, l)
                 idx = (m*m+m)/2 - m
                 do n = 1, m
-                    j = sorted_atoms_all(n, k)
-                    cm(k, idx+n) = pair_distance_matrix(i,j,k)
+                    j = sorted_atoms_all(n, l)
+                    cm(l, idx+n) = pair_distance_matrix(i,j,l)
                 enddo
             enddo
         enddo
@@ -403,18 +412,20 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
 
 end subroutine fgenerate_local_coulomb_matrix
 
-subroutine fgenerate_atomic_coulomb_matrix(atomic_charges, coordinates, natoms, nmax, &
-        & cent_cutoff, cent_decay, int_cutoff, int_decay, cm)
+subroutine fgenerate_atomic_coulomb_matrix(central_atom_indices, central_natoms, atomic_charges, &
+        & coordinates, natoms, nmax, cent_cutoff, cent_decay, int_cutoff, int_decay, cm)
 
     implicit none
 
+    integer, dimension(:), intent(in) :: central_atom_indices
+    integer, intent(in) :: central_natoms
     double precision, dimension(:), intent(in) :: atomic_charges
     double precision, dimension(:,:), intent(in) :: coordinates
     integer,intent(in) :: natoms
     integer, intent(in) :: nmax
     double precision, intent(inout) :: cent_cutoff, cent_decay, int_cutoff, int_decay
 
-    double precision, dimension(natoms, ((nmax + 1) * nmax) / 2), intent(out):: cm
+    double precision, dimension(central_natoms, ((nmax + 1) * nmax) / 2), intent(out):: cm
 
     integer :: idx
 
@@ -430,7 +441,7 @@ subroutine fgenerate_atomic_coulomb_matrix(atomic_charges, coordinates, natoms,
     double precision, allocatable, dimension(:, :) :: distance_matrix
     double precision, allocatable, dimension(:, :) :: distance_matrix_tmp
 
-    integer i, j, m, n, k
+    integer i, j, m, n, k, l
 
     double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
 
@@ -485,11 +496,12 @@ subroutine fgenerate_atomic_coulomb_matrix(atomic_charges, coordinates, natoms,
     enddo
     !$OMP END PARALLEL DO
 
-    do i = 1, natoms
-        if (cutoff_count(i) > nmax) then
+    do i = 1, central_natoms
+        k = central_atom_indices(i)
+        if (cutoff_count(k) > nmax) then
             write(*,*) "ERROR: Coulomb matrix generation"
             write(*,*) nmax, "size set, but", &
-                & cutoff_count(i), "size needed!"
+                & cutoff_count(k), "size needed!"
             stop
         endif
     enddo
@@ -520,20 +532,22 @@ subroutine fgenerate_atomic_coulomb_matrix(atomic_charges, coordinates, natoms,
     !$OMP END PARALLEL DO
 
     ! Allocate temporary
-    allocate(sorted_atoms_all(natoms, natoms))
+    allocate(distance_matrix_tmp(natoms, natoms))
+    allocate(sorted_atoms_all(natoms, central_natoms))
 
     distance_matrix_tmp = distance_matrix
     !Generate sorted list of atom ids by distance matrix
-    !$OMP PARALLEL DO PRIVATE(j)
-        do k = 1, natoms
-            !$OMP CRITICAL
-                do i = 1, cutoff_count(k)
-                    j = minloc(distance_matrix_tmp(:,k), dim=1)
-                    sorted_atoms_all(i, k) = j
-                    distance_matrix_tmp(j, k) = huge_double
-                enddo
-            !$OMP END CRITICAL
-        enddo
+    !$OMP PARALLEL DO PRIVATE(j, k)
+    do l = 1, central_natoms
+        k = central_atom_indices(l)
+        !$OMP CRITICAL
+            do i = 1, cutoff_count(k)
+                j = minloc(distance_matrix_tmp(:,k), dim=1)
+                sorted_atoms_all(i, l) = j
+                distance_matrix_tmp(j, k) = huge_double
+            enddo
+        !$OMP END CRITICAL
+    enddo
     !$OMP END PARALLEL DO
 
     ! Clean up
@@ -544,34 +558,36 @@ subroutine fgenerate_atomic_coulomb_matrix(atomic_charges, coordinates, natoms,
 
     pair_norm = 0.0d0
 
-        do k = 1, natoms
-            do m = 1, cutoff_count(k)
-                i = sorted_atoms_all(m, k)
+    !$OMP PARALLEL DO PRIVATE(i, prefactor, idx, j, pair_norm, k)
+    do l = 1, central_natoms
+        k = central_atom_indices(l)
+        do m = 1, cutoff_count(k)
+            i = sorted_atoms_all(m, l)
 
-                if (distance_matrix(i,k) > cent_cutoff) then
-                    cycle
-                endif
-                prefactor = 1.0d0
-                if (distance_matrix(i,k) > cent_cutoff - cent_decay) then
-                    prefactor = 0.5d0 * (cos(pi &
-                        & * (distance_matrix(i,k) - cent_cutoff + cent_decay) &
-                        & / cent_decay) + 1.0d0)
-                endif
+            if (distance_matrix(i,k) > cent_cutoff) then
+                cycle
+            endif
+            prefactor = 1.0d0
+            if (distance_matrix(i,k) > cent_cutoff - cent_decay) then
+                prefactor = 0.5d0 * (cos(pi &
+                    & * (distance_matrix(i,k) - cent_cutoff + cent_decay) &
+                    & / cent_decay) + 1.0d0)
+            endif
 
-                idx = (m*m+m)/2 - m
-                do n = 1, m
-                    j = sorted_atoms_all(n, k)
-
-                    pair_norm = prefactor * pair_distance_matrix(i, j)
-                    if (distance_matrix(j,k) > cent_cutoff - cent_decay) then
-                        pair_norm = pair_norm * 0.5d0 * (cos(pi &
-                            & * (distance_matrix(j,k) - cent_cutoff + cent_decay) &
-                            & / cent_decay) + 1)
-                    endif
-                    cm(k, idx+n) = pair_norm
-                enddo
+            idx = (m*m+m)/2 - m
+            do n = 1, m
+                j = sorted_atoms_all(n, l)
+
+                pair_norm = prefactor * pair_distance_matrix(i, j)
+                if (distance_matrix(j,k) > cent_cutoff - cent_decay) then
+                    pair_norm = pair_norm * 0.5d0 * (cos(pi &
+                        & * (distance_matrix(j,k) - cent_cutoff + cent_decay) &
+                        & / cent_decay) + 1)
+                endif
+                cm(l, idx+n) = pair_norm
             enddo
         enddo
+    enddo
 
     ! Clean up
     deallocate(distance_matrix)

qml/kernels.py

@@ -1,6 +1,6 @@
 # MIT License
 #
-# Copyright (c) 2016 Anders Steen Christensen, Felix Faber
+# Copyright (c) 2016 Anders Steen Christensen, Felix A. Faber, Lars A. Bratholm
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
@@ -26,9 +26,13 @@
 
 from .fkernels import fgaussian_kernel
 from .fkernels import flaplacian_kernel
+from .fkernels import flinear_kernel
 from .fkernels import fsargan_kernel
 from .fkernels import fmatern_kernel_l2
 
+from .fkernels import fget_local_kernels_gaussian
+from .fkernels import fget_local_kernels_laplacian
+
 def laplacian_kernel(A, B, sigma):
     """ Calculates the Laplacian kernel matrix K, where :math:`K_{ij}`:
 
@@ -37,9 +41,9 @@ def laplacian_kernel(A, B, sigma):
         Where :math:`A_{i}` and :math:`B_{j}` are representation vectors.
         K is calculated using an OpenMP parallel Fortran routine.
 
-        :param A: 2D array of descriptors - shape (N, representation size).
+        :param A: 2D array of representations - shape (N, representation size).
         :type A: numpy array
-        :param B: 2D array of descriptors - shape (M, representation size).
+        :param B: 2D array of representations - shape (M, representation size).
         :type B: numpy array
         :param sigma: The value of sigma in the kernel matrix.
         :type sigma: float
@@ -66,9 +70,9 @@ def gaussian_kernel(A, B, sigma):
         Where :math:`A_{i}` and :math:`B_{j}` are representation vectors.
         K is calculated using an OpenMP parallel Fortran routine.
 
-        :param A: 2D array of descriptors - shape (N, representation size).
+        :param A: 2D array of representations - shape (N, representation size).
         :type A: numpy array
-        :param B: 2D array of descriptors - shape (M, representation size).
+        :param B: 2D array of representations - shape (M, representation size).
         :type B: numpy array
         :param sigma: The value of sigma in the kernel matrix.
         :type sigma: float
@@ -87,6 +91,34 @@ def gaussian_kernel(A, B, sigma):
 
     return K
 
+def linear_kernel(A, B):
+    """ Calculates the linear kernel matrix K, where :math:`K_{ij}`:
+
+            :math:`K_{ij} = A_i \cdot B_j`
+
+        VWhere :math:`A_{i}` and :math:`B_{j}` are  representation vectors. 
+
+        K is calculated using an OpenMP parallel Fortran routine.
+
+        :param A: 2D array of representations - shape (N, representation size).
+        :type A: numpy array
+        :param B: 2D array of representations - shape (M, representation size).
+        :type B: numpy array
+
+        :return: The Gaussian kernel matrix - shape (N, M)
+        :rtype: numpy array
+    """
+
+    na = A.shape[0]
+    nb = B.shape[0]
+
+    K = np.empty((na, nb), order='F')
+
+    # Note: Transposed for Fortran
+    flinear_kernel(A.T, na, B.T, nb, K)
+
+    return K
+
 def sargan_kernel(A, B, sigma, gammas):
     """ Calculates the Sargan kernel matrix K, where :math:`K_{ij}`:
 
@@ -95,9 +127,9 @@ def sargan_kernel(A, B, sigma, gammas):
         Where :math:`A_{i}` and :math:`B_{j}` are representation vectors.
         K is calculated using an OpenMP parallel Fortran routine.
 
-        :param A: 2D array of descriptors - shape (N, representation size).
+        :param A: 2D array of representations - shape (N, representation size).
         :type A: numpy array
-        :param B: 2D array of descriptors - shape (M, representation size).
+        :param B: 2D array of representations - shape (M, representation size).
         :type B: numpy array
         :param sigma: The value of sigma in the kernel matrix.
         :type sigma: float
@@ -137,9 +169,9 @@ def matern_kernel(A, B, sigma, order = 0, metric = "l1"):
 
         K is calculated using an OpenMP parallel Fortran routine.
 
-        :param A: 2D array of descriptors - shape (N, representation size).
+        :param A: 2D array of representations - shape (N, representation size).
         :type A: numpy array
-        :param B: 2D array of descriptors - shape (M, representation size).
+        :param B: 2D array of representations - shape (M, representation size).
         :type B: numpy array
         :param sigma: The value of sigma in the kernel matrix.
         :type sigma: float
@@ -183,3 +215,84 @@ def matern_kernel(A, B, sigma, order = 0, metric = "l1"):
 
     return K
 
+def get_local_kernels_gaussian(A, B, na, nb, sigmas):
+    """ Calculates the Gaussian kernel matrix K, for a local representation where :math:`K_{ij}`:
+
+            :math:`K_{ij} = \sum_{a \in i} \sum_{b \in j} \\exp \\big( -\\frac{\\|A_a - B_b\\|_2^2}{2\sigma^2} \\big)`
+
+        Where :math:`A_{a}` and :math:`B_{b}` are representation vectors.
+
+        Note that the input array is one big 2D array with all atoms concatenated along the same axis.
+        Further more a series of kernels is produced (since calculating the distance matrix is expensive
+        but getting the resulting kernels elements for several sigmas is not.)
+
+        K is calculated using an OpenMP parallel Fortran routine.
+
+        :param A: 2D array of descriptors - shape (total atoms A, representation size).
+        :type A: numpy array
+        :param B: 2D array of descriptors - shape (total atoms B, representation size).
+        :type B: numpy array
+        :param na: 1D array containing numbers of atoms in each compound.
+        :type na: numpy array
+        :param nb: 1D array containing numbers of atoms in each compound.
+        :type nb: numpy array
+        :param sigma: The value of sigma in the kernel matrix.
+        :type sigma: float
+
+        :return: The Gaussian kernel matrix - shape (nsigmas, N, M)
+        :rtype: numpy array
+    """
+
+    assert np.sum(na) == A.shape[0], "Error in A input"
+    assert np.sum(nb) == B.shape[0], "Error in B input"
+
+    assert A.shape[1] == B.shape[1], "Error in representation sizes"
+
+    nma = len(na)
+    nmb = len(nb)
+
+    sigmas = np.asarray(sigmas)
+    nsigmas = len(sigmas)
+
+    return fget_local_kernels_gaussian(A.T, B.T, na, nb, sigmas, nma, nmb, nsigmas)
+
+def get_local_kernels_laplacian(A, B, na, nb, sigmas):
+    """ Calculates the Local Laplacian kernel matrix K, for a local representation where :math:`K_{ij}`:
+
+            :math:`K_{ij} = \sum_{a \in i} \sum_{b \in j} \\exp \\big( -\\frac{\\|A_a - B_b\\|_1}{\sigma} \\big)`
+
+        Where :math:`A_{a}` and :math:`B_{b}` are representation vectors.
+
+        Note that the input array is one big 2D array with all atoms concatenated along the same axis.
+        Further more a series of kernels is produced (since calculating the distance matrix is expensive
+        but getting the resulting kernels elements for several sigmas is not.)
+        
+        K is calculated using an OpenMP parallel Fortran routine.
+
+        :param A: 2D array of descriptors - shape (N, representation size).
+        :type A: numpy array
+        :param B: 2D array of descriptors - shape (M, representation size).
+        :type B: numpy array
+        :param na: 1D array containing numbers of atoms in each compound.
+        :type na: numpy array
+        :param nb: 1D array containing numbers of atoms in each compound.
+        :type nb: numpy array
+        :param sigmas: List of the sigmas.
+        :type sigmas: list
+
+        :return: The Laplacian kernel matrix - shape (nsigmas, N, M)
+        :rtype: numpy array
+    """
+
+    assert np.sum(na) == A.shape[0], "Error in A input"
+    assert np.sum(nb) == B.shape[0], "Error in B input"
+
+    assert A.shape[1] == B.shape[1], "Error in representation sizes"
+
+    nma = len(na)
+    nmb = len(nb)
+
+    sigmas = np.asarray(sigmas)
+    nsigmas = len(sigmas)
+
+    return fget_local_kernels_laplacian(A.T, B.T, na, nb, sigmas, nma, nmb, nsigmas)

qml/math.py

@@ -22,45 +22,56 @@
 
 import numpy as np
 
+from copy import deepcopy
+
 from .fcho_solve import fcho_solve
-# from .fcho_solve import fcho_invert
-
-
-# Disabled due to bug.
-# def cho_invert(A):
-#     """ Solves [A x = y] for x using a Cholesky decomposition
-#         via calls to LAPACK dpotrf and dpotri in the F2PY module.
-# 
-#         Arguments:
-#         ==============
-#         A -- the A-matrix (symmetric and positive definite).
-# 
-#         Returns:
-#         ==============
-#         A -- the inverted A-matrix
-#     """
-# 
-#     B = np.asfortranarray(A)
-#     fcho_invert(B)
-# 
-#     return B
+from .fcho_solve import fcho_invert
+from .fcho_solve import fbkf_solve
+from .fcho_solve import fbkf_invert
+
+
+def cho_invert(A):
+    """ Returns the inverse of a positive definite matrix, using a Cholesky decomposition
+        via calls to LAPACK dpotrf and dpotri in the F2PY module.
+
+        :param A: Matrix (symmetric and positive definite, left-hand side).
+        :type A: numpy array
+
+        :return: The inverse matrix
+        :rtype: numpy array
+    """
+
+    if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
+        raise ValueError('expected square matrix')
+
+    I = np.asfortranarray(A)
+
+    fcho_invert(I)
+
+    # Matrix to store the inverse
+    i_lower = np.tril_indices_from(A)
+
+    # Copy lower triangle to upper
+    I.T[i_lower] = I[i_lower]
+
+    return I
 
 
 def cho_solve(A, y):
     """ Solves the equation
 
-            :math:`A x = y` 
+            :math:`A x = y`
 
-        for x using a Cholesky decomposition  via calls to LAPACK dpotrf and dpotrs in the F2PY module.
+        for x using a Cholesky decomposition  via calls to LAPACK dpotrf and dpotrs in the F2PY module. Preserves the input matrix A.
 
         :param A: Matrix (symmetric and positive definite, left-hand side).
-        :type A: numpy array 
+        :type A: numpy array
         :param y: Vector (right-hand side of the equation).
-        :type y: numpy array 
+        :type y: numpy array
 
         :return: The solution vector.
-        :rtype: numpy array 
-    """
+        :rtype: numpy array
+        """
 
     if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
         raise ValueError('expected square matrix')
@@ -70,7 +81,84 @@ def cho_solve(A, y):
 
     n = A.shape[0]
 
+    # Backup diagonal before Cholesky-decomposition
+    A_diag = A[np.diag_indices_from(A)]
+
     x = np.zeros(n)
     fcho_solve(A, y, x)
 
+    # Reset diagonal after Cholesky-decomposition
+    A[np.diag_indices_from(A)] = A_diag
+
+    # Copy lower triangle to upper
+    i_lower = np.tril_indices_from(A)
+    A.T[i_lower] = A[i_lower]
+
+    return x
+
+
+def bkf_invert(A):
+    """ Returns the inverse of a positive definite matrix, using a Cholesky decomposition
+        via calls to LAPACK dpotrf and dpotri in the F2PY module.
+
+        :param A: Matrix (symmetric and positive definite, left-hand side).
+        :type A: numpy array
+
+        :return: The inverse matrix
+        :rtype: numpy array
+    """
+
+    if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
+        raise ValueError('expected square matrix')
+
+    I = np.asfortranarray(A)
+
+    fbkf_invert(I)
+
+    # Matrix to store the inverse
+    i_lower = np.tril_indices_from(A)
+
+    # Copy lower triangle to upper
+    I.T[i_lower] = I[i_lower]
+
+    return I
+
+
+def bkf_solve(A, y):
+    """ Solves the equation
+
+            :math:`A x = y`
+
+        for x using a Cholesky decomposition  via calls to LAPACK dpotrf and dpotrs in the F2PY module. Preserves the input matrix A.
+
+        :param A: Matrix (symmetric and positive definite, left-hand side).
+        :type A: numpy array
+        :param y: Vector (right-hand side of the equation).
+        :type y: numpy array
+
+        :return: The solution vector.
+        :rtype: numpy array
+        """
+
+    if len(A.shape) != 2 or A.shape[0] != A.shape[1]:
+        raise ValueError('expected square matrix')
+
+    if len(y.shape) != 1 or y.shape[0] != A.shape[1]:
+        raise ValueError('expected matrix and vector of same stride size')
+
+    n = A.shape[0]
+
+    # Backup diagonal before Cholesky-decomposition
+    A_diag = A[np.diag_indices_from(A)]
+
+    x = np.zeros(n)
+    fbkf_solve(A, y, x)
+
+    # Reset diagonal after Cholesky-decomposition
+    A[np.diag_indices_from(A)] = A_diag
+
+    # Copy lower triangle to upper
+    i_lower = np.tril_indices_from(A)
+    A.T[i_lower] = A[i_lower]
+
     return x

qml/representations.py

@@ -1,6 +1,6 @@
 # MIT License
 #
-# Copyright (c) 2017 Anders Steen Christensen, Lars A. Bratholm and Bing Huang
+# Copyright (c) 2017 Anders Steen Christensen, Lars Andersen Bratholm and Bing Huang
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
@@ -30,14 +30,41 @@
 from .frepresentations import fgenerate_local_coulomb_matrix
 from .frepresentations import fgenerate_atomic_coulomb_matrix
 from .frepresentations import fgenerate_eigenvalue_coulomb_matrix
-from .frepresentations import fgenerate_bob
 
 from .data import NUCLEAR_CHARGE
 
 from .slatm import get_boa
 from .slatm import get_sbop
 from .slatm import get_sbot
 
+def vector_to_matrix(v):
+    """ Converts a representation from 1D vector to 2D square matrix.
+    :param v: 1D input representation.
+    :type v: numpy array 
+    :return: Square matrix representation.
+    :rtype: numpy array 
+    """
+
+    if not (np.sqrt(8*v.shape[0]+1) == int(np.sqrt(8*v.shape[0]+1))):
+        print("ERROR: Can not make a square matrix.")
+        exit(1)
+
+    n = v.shape[0]
+    l = (-1 + int(np.sqrt(8*n+1)))//2
+    M = np.empty((l,l))
+
+    index = 0
+    for i in range(l):
+        for j in range(l):
+            if j > i:
+                continue
+
+            M[i,j] = v[index]
+            M[j,i] = M[i,j]
+
+            index += 1
+    return M
+
 def generate_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "row-norm"):
     """ Creates a Coulomb Matrix representation of a molecule.
         Sorting of the elements can either be done by ``sorting="row-norm"`` or ``sorting="unsorted"``.
@@ -90,9 +117,9 @@ def generate_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "
         print("ERROR: Unknown sorting scheme requested")
         raise SystemExit
 
-
 def generate_atomic_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "distance",
-            central_cutoff = 1e6, central_decay = -1, interaction_cutoff = 1e6, interaction_decay = -1):
+            central_cutoff = 1e6, central_decay = -1, interaction_cutoff = 1e6, interaction_decay = -1,
+            indices = None):
     """ Creates a Coulomb Matrix representation of the local environment of a central atom.
         For each central atom :math:`k`, a matrix :math:`M` is constructed with elements
 
@@ -138,6 +165,11 @@ def generate_atomic_coulomb_matrix(nuclear_charges, coordinates, size = 23, sort
 
         The upper triangular of M, including the diagonal, is concatenated to a 1D
         vector representation.
+
+        The representation can be calculated for a subset by either specifying
+        ``indices = [0,1,...]``, where :math:`[0,1,...]` are the requested atom indices,
+        or by specifying ``indices = 'C'`` to only calculate central carbon atoms.
+
         The representation is calculated using an OpenMP parallel Fortran routine.
 
         :param nuclear_charges: Nuclear charges of the atoms in the molecule
@@ -158,20 +190,40 @@ def generate_atomic_coulomb_matrix(nuclear_charges, coordinates, size = 23, sort
         :type interaction_cutoff: float
         :param interaction_decay: The distance over which the the coulomb interaction decays from full to none
         :type interaction_decay: float
+        :param indices: Subset indices or atomtype
+        :type indices: Nonetype/array/string
 
 
         :return: nD representation - shape (:math:`N_{atoms}`, size(size+1)/2)
         :rtype: numpy array
     """
 
 
+    if indices == None:
+        nindices = len(nuclear_charges)
+        indices = np.arange(1,1+nindices, 1, dtype = int)
+    elif type("") == type(indices):
+        if indices in NUCLEAR_CHARGE:
+            indices = np.where(nuclear_charges == NUCLEAR_CHARGE[indices])[0] + 1
+            nindices = indices.size
+            if nindices == 0:
+                return np.zeros((0,0))
+
+        else:
+            print("ERROR: Unknown value %s given for 'indices' variable" % indices)
+            raise SystemExit
+    else:
+        indices = np.asarray(indices, dtype = int) + 1
+        nindices = indices.size
+
+
     if (sorting == "row-norm"):
-        return fgenerate_local_coulomb_matrix(nuclear_charges,
+        return fgenerate_local_coulomb_matrix(indices, nindices, nuclear_charges,
             coordinates, nuclear_charges.size, size,
             central_cutoff, central_decay, interaction_cutoff, interaction_decay)
 
     elif (sorting == "distance"):
-        return fgenerate_atomic_coulomb_matrix(nuclear_charges,
+        return fgenerate_atomic_coulomb_matrix(indices, nindices, nuclear_charges,
             coordinates, nuclear_charges.size, size, 
             central_cutoff, central_decay, interaction_cutoff, interaction_decay)
 
@@ -209,7 +261,7 @@ def generate_eigenvalue_coulomb_matrix(nuclear_charges, coordinates, size = 23):
     return fgenerate_eigenvalue_coulomb_matrix(nuclear_charges,
         coordinates, size)
 
-def generate_bob(nuclear_charges, coordinates, atomtypes, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}):
+def generate_bob(nuclear_charges, coordinates, atomtypes, size=23, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}):
     """ Creates a Bag of Bonds (BOB) representation of a molecule.
         The representation expands on the coulomb matrix representation.
         For each element a bag (vector) is constructed for self interactions
@@ -234,12 +286,45 @@ def generate_bob(nuclear_charges, coordinates, atomtypes, asize = {"O":3, "C":7,
         :type nuclear_charges: numpy array
         :param coordinates: 3D Coordinates of the atoms in the molecule
         :type coordinates: numpy array
+        :param size: The maximum number of atoms in the representation
+        :type size: integer
         :param asize: The maximum number of atoms of each element type supported by the representation
-        :type size: dictionary
+        :type asize: dictionary
 
         :return: 1D representation
         :rtype: numpy array
     """
+    natoms = len(nuclear_charges)
+
+    coulomb_matrix = fgenerate_unsorted_coulomb_matrix(nuclear_charges, coordinates, size)
+
+    coulomb_matrix = vector_to_matrix(coulomb_matrix)
+    descriptor = []
+    atomtypes = np.asarray(atomtypes)
+    for atom1, size1 in sorted(asize.items()):
+        pos1 = np.where(atomtypes == atom1)[0]
+        feature_vector = np.zeros(size1)
+        feature_vector[:pos1.size] = np.diag(coulomb_matrix)[pos1]
+        feature_vector.sort()
+        descriptor.append(feature_vector[:])
+        for atom2, size2 in sorted(asize.items()):
+            if atom1 > atom2:
+                continue
+            if atom1 == atom2:
+                size = size1*(size1-1)//2
+                feature_vector = np.zeros(size)
+                sub_matrix = coulomb_matrix[np.ix_(pos1,pos1)]
+                feature_vector[:pos1.size*(pos1.size-1)//2] = sub_matrix[np.triu_indices(pos1.size, 1)]
+                feature_vector.sort()
+                descriptor.append(feature_vector[:])
+            else:
+                pos2 = np.where(atomtypes == atom2)[0]
+                feature_vector = np.zeros(size1*size2)
+                feature_vector[:pos1.size*pos2.size] = coulomb_matrix[np.ix_(pos1,pos2)].ravel()
+                feature_vector.sort()
+                descriptor.append(feature_vector[:])
+
+    return np.concatenate(descriptor)
 
     n = 0
     atoms = sorted(asize, key=asize.get)
@@ -399,7 +484,7 @@ def generate_slatm(coordinates, nuclear_charges, mbtypes,
                     #print ' 001, pbc = ', pbc
                     mbsi = get_sbop(mbtype, obj, iloc=iloc, ia=ia, \
                                     sigma=sigmas[0], dgrid=dgrids[0], rcut=rcut, \
-                                    pbc=pbc, rpower=rpower)[1]
+                                    pbc=pbc, rpower=rpower)
                     mbsi *= 0.5 # only for the two-body parts, local rpst
                     #print ' 002'
                     if alchemy:
@@ -417,7 +502,8 @@ def generate_slatm(coordinates, nuclear_charges, mbtypes,
                         mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 )
                 else: # len(mbtype) == 3:
                     mbsi = get_sbot(mbtype, obj, iloc=iloc, ia=ia, \
-                                    sigma=sigmas[1], dgrid=dgrids[1], rcut=rcut, pbc=pbc)[1]
+                                    sigma=sigmas[1], dgrid=dgrids[1], rcut=rcut, pbc=pbc)
+
                     if alchemy:
                         n3 = len(mbsi)
                         n3_0 = mbs_ia.shape[0]
@@ -457,7 +543,8 @@ def generate_slatm(coordinates, nuclear_charges, mbtypes,
                     mbs = np.concatenate( (mbs, mbsi), axis=0 )
             elif len(mbtype) == 2:
                 mbsi = get_sbop(mbtype, obj, sigma=sigmas[0], \
-                                dgrid=dgrids[0], rcut=rcut, rpower=rpower)[1]
+                                dgrid=dgrids[0], rcut=rcut, rpower=rpower)
+
 
                 if alchemy:
                     n2 = len(mbsi)
@@ -474,7 +561,8 @@ def generate_slatm(coordinates, nuclear_charges, mbtypes,
                     mbs = np.concatenate( (mbs, mbsi), axis=0 )
             else: # len(mbtype) == 3:
                 mbsi = get_sbot(mbtype, obj, sigma=sigmas[1], \
-                        dgrid=dgrids[1], rcut=rcut)[1]
+                        dgrid=dgrids[1], rcut=rcut)
+
                 if alchemy:
                     n3 = len(mbsi)
                     n3_0 = mbs.shape[0]

qml/slatm.py

@@ -20,45 +20,28 @@
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 # SOFTWARE.
 
+
+from __future__ import print_function
+
 import scipy.spatial.distance as ssd
 import itertools as itl
 import numpy as np
 
-def get_pbc(obj, d0 = 3.6):
-    """
-    automatically tell if an compound object is periodic or not
-
-    :param d0: the threshhold value to tell if two cells are adjacent
-    :type d0: float
-    """
-
-    pbc = []
-
-    zs, ps, c = obj
-    na = len(zs); idxs = np.arange(na)
-
-    for i in range(3):
-        psx = ps[:,i]; xmin = min(psx)
-        idxs_i = idxs[ psx == xmin ]
-        ps1 = ps[idxs_i[0]] + c[i]
-        if np.min( ssd.cdist([ps1,], ps)[0] ) < d0:
-            pbc.append( '1' )
-        else:
-            pbc.append( '0' )
-
-    return ''.join(pbc)
+from .fslatm import fget_sbot
+from .fslatm import fget_sbot_local
+from .fslatm import fget_sbop
+from .fslatm import fget_sbop_local
 
 def update_m(obj, ia, rcut=9.0, pbc=None):
     """
     retrieve local structure around atom `ia
     for periodic systems (or very large system)
     """
-
     zs, coords, c = obj
     v1, v2, v3 = c
     vs = ssd.norm(c, axis=0)
 
-    nns = []; ns = []
+    nns = []
     for i,vi in enumerate(vs):
         n1_doulbe = rcut/li
         n1 = int(n1_doulbe)
@@ -83,7 +66,8 @@ def update_m(obj, ia, rcut=9.0, pbc=None):
     n123s = np.array(n123s, np.float)
 
     na = len(zs)
-    ai = m[ia]; cia = coords[ia]
+    cia = coords[ia]
+
     if na == 1:
         ds = np.array([[0.]])
     else:
@@ -116,8 +100,9 @@ def get_boa(z1, zs_):
     return z1*np.array( [(zs_ == z1).sum(), ])
     #return -0.5*z1**2.4*np.array( [(zs_ == z1).sum(), ])
 
+
 def get_sbop(mbtype, obj, iloc=False, ia=None, normalize=True, sigma=0.05, \
-             rcut=4.8, dgrid=0.03, ipot=True, pbc='000', rpower=6):
+             rcut=4.8, dgrid=0.03, pbc='000', rpower=6):
     """
     two-body terms
 
@@ -139,77 +124,22 @@ def get_sbop(mbtype, obj, iloc=False, ia=None, normalize=True, sigma=0.05, \
         # after update of `m, the query atom `ia will become the first atom
         ia = 0
 
-    na = len(zs)
-    ds_triu = ssd.pdist(coords) # upper triangle part of distance matrix
-    ds = ssd.squareform(ds_triu)
-    dmin = 0.25 # Angstrom
-    assert np.all(ds_triu.ravel() > dmin), '#ERROR: two small distance detected!'
-
-    ias = np.arange(na)
-    ias1 = ias[zs == z1]
-    ias2 = ias[zs == z2]
-
-    if z1 == z2:
-        ias12 = list( itl.combinations(ias1,2) )
-    else:
-        ias12 = itl.product(ias1,ias2)
-
-    if iloc:
-        dsu = []; icnt = 0
-        for j1,j2 in ias12:
-            if ia == j1 or ia == j2:
-                dsu.append( ds[j1,j2] )
-            icnt += 1
-    else:
-        dsu = [ ds[i,j] for (i,j) in ias12 ]
-
-    dsu = np.array(dsu)
-
     # bop potential distribution
     r0 = 0.1
-    nx = (rcut - r0)/dgrid + 1
-    xs = np.linspace(r0, rcut, nx)
-    ys0 = np.zeros(xs.shape)
+    nx = int((rcut - r0)/dgrid) + 1
 
-    # update dsu by exluding d > 6.0
-    nr = dsu.shape[0]
-    if nr > 0:
-        dsu = dsu[ dsu <= rcut ]
-        nr = len(dsu)
-
-    #print ' -- dsu = ', dsu
 
     coeff = 1/np.sqrt(2*sigma**2*np.pi) if normalize else 1.0
-    #print ' -- now calculating 2-body terms...'
-    if ipot:
-        # get distribution of 2-body potentials
-        # unit of x: Angstrom
-        c0 = (z1%1000)*(z2%1000)*coeff
-        ys = ys0
-        for i in range(nr):
-            ys += ( c0/(xs**rpower) )*np.exp( -0.5*((xs-dsu[i])/sigma)**2 )
-        ys *= dgrid
-    else:
-        # print distribution of distances
-        c0 = coeff
-        ys = ys0
-        for i in range(nr):
-            ys += c0*np.exp( -0.5*((xs-dsu[i])/sigma)**2 )
 
-    return xs, ys
-
-def vang(u,v):
-    cost = np.dot(u,v)/(np.linalg.norm(u) * np.linalg.norm(v))
-    # sometimes, cost might be 1.00000000002, then np.arccos(cost)
-    # does not exist!
-    u = cost if abs(cost) <= 1 else 1.0
-    return np.arccos( u )
+    if iloc:
+        ys = fget_sbop_local(coords, zs, ia, z1, z2, rcut, nx, dgrid, sigma, coeff, rpower)
+    else:
+        ys = fget_sbop(coords, zs, z1, z2, rcut, nx, dgrid, sigma, coeff, rpower)
 
-def cvang(u,v):
-    return np.dot(u,v)/np.sqrt(np.dot(u,u)*np.dot(v,v))
+    return ys
 
-def get_sbot(mbtype, obj, iloc=False, ia=None, normalize=True, sigma=0.05, label=None, \
-             rcut=4.8, dgrid=0.0262, ipot=True, pbc='000'):
+def get_sbot(mbtype, obj, iloc=False, ia=None, normalize=True, sigma=0.05, \
+             rcut=4.8, dgrid=0.0262, pbc='000'):
 
     """
     sigma -- standard deviation of gaussian distribution centered on a specific angle
@@ -231,78 +161,19 @@ def get_sbot(mbtype, obj, iloc=False, ia=None, normalize=True, sigma=0.05, label
         # after update of `m, the query atom `ia will become the first atom
         ia = 0
 
-    na = len(zs)
-    ds = ssd.squareform( ssd.pdist(coords) )
-    #get_date(' ds matrix calc done ')
-
-    ias = np.arange(na)
-    ias1 = ias[zs == z1]; n1 = len(ias1)
-    ias2 = ias[zs == z2]; n2 = len(ias2)
-    ias3 = ias[zs == z3]; n3 = len(ias3)
-    tas = []
-
-    for ia1 in ias1:
-        ias2u = ias2[ np.logical_and( ds[ia1,ias2] > 0, ds[ia1,ias2] <= rcut ) ]
-        for ia2 in ias2u:
-            filt1 = np.logical_and( ds[ia1,ias3] > 0, ds[ia1,ias3] <= rcut )
-            filt2 = np.logical_and( ds[ia2,ias3] > 0, ds[ia2,ias3] <= rcut )
-            ias3u = ias3[ np.logical_and(filt1, filt2) ]
-            for ia3 in ias3u:
-                tasi = [ia1,ia2,ia3]
-                iok1 = (tasi not in tas)
-                iok2 = (tasi[::-1] not in tas)
-                if iok1 and iok2:
-                    tas.append( tasi )
-
-    if iloc:
-        tas_u = []
-        for tas_i in tas:
-            if ia == tas_i[1]:
-                tas_u.append( tas_i )
-        tas = tas_u
+    # for a normalized gaussian distribution, u should multiply this coeff
+    coeff = 1/np.sqrt(2*sigma**2*np.pi) if normalize else 1.0
 
+    # Setup grid in Python
     d2r = np.pi/180 # degree to rad
-    a0 = -20.0*d2r; a1 = np.pi + 20.0*d2r
+    a0 = -20.0*d2r
+    a1 = np.pi + 20.0*d2r
     nx = int((a1-a0)/dgrid) + 1
-    xs = np.linspace(a0, a1, nx)
-    ys0 = np.zeros(nx, np.float)
-    nt = len(tas)
 
-    # u actually have considered the same 3-body term for
-    # three times, so rescale it
-    prefactor = 1.0/3
+    if iloc:
+        ys = fget_sbot_local(coords, zs, ia, z1, z2, z3, rcut, nx, dgrid, sigma, coeff)
+    else:
+        ys = fget_sbot(coords, zs, z1, z2, z3, rcut, nx, dgrid, sigma, coeff)
 
-    # for a normalized gaussian distribution, u should multiply this coeff
-    coeff = 1/np.sqrt(2*sigma**2*np.pi) if normalize else 1.0
+    return ys
 
-    tidxs = np.array(tas, np.int)
-    if ipot:
-        # get distribution of 3-body potentials
-        # unit of x: Angstrom
-        c0 = prefactor*(z1%1000)*(z2%1000)*(z3%1000)*coeff
-        ys = ys0
-        for it in range(nt):
-            i,j,k = tas[it]
-            # angle spanned by i <-- j --> k, i.e., vector ji and jk
-            u = coords[i]-coords[j]; v = coords[k] - coords[j]
-            ang = vang( u, v ) # ang_j
-            #print ' -- (i,j,k) = (%d,%d,%d),  ang = %.2f'%(i,j,k, ang)
-            cak = cvang( coords[j]-coords[k], coords[i]-coords[k] ) # cos(ang_k)
-            cai = cvang( coords[k]-coords[i], coords[j]-coords[i] ) # cos(ang_i)
-            ys += c0*( (1.0 + 1.0*np.cos(xs)*cak*cai)/(ds[i,j]*ds[i,k]*ds[j,k])**3 )*\
-                                ( np.exp(-(xs-ang)**2/(2*sigma**2)) )
-        ys *= dgrid
-    else:
-        # print distribution of angles (unit: degree)
-        sigma = sigma/d2r
-        xs = xs/d2r
-        c0 = 1
-
-        ys = ys0
-        for it in range(nt):
-            i,j,k = tas[it]
-            # angle spanned by i <-- j --> k, i.e., vector ji and jk
-            ang = vang( coords[i]-coords[j], coords[k]-coords[j] )/d2r
-            ys += c0*np.exp( -(xs-ang)**2/(2*sigma**2) )
-
-    return xs, ys

qml/wrappers.py

@@ -24,6 +24,7 @@
 
 from .fkernels import fget_vector_kernels_gaussian
 from .fkernels import fget_vector_kernels_laplacian
+from .fkernels import fget_local_kernels_gaussian
 
 from .arad import get_local_kernels_arad
 from .arad import get_local_symmetric_kernels_arad
@@ -55,7 +56,6 @@ def get_atomic_kernels_laplacian(mols1, mols2, sigmas):
     x1 = np.swapaxes(x1, 0, 2)
     x2 = np.swapaxes(x2, 0, 2)
 
-
     sigmas = np.asarray(sigmas, dtype=np.float64)
     nsigmas = sigmas.size
 

setup.py

@@ -1,6 +1,8 @@
 import sys
 from numpy.distutils.core import Extension, setup
 
+from mkldiscover import mkl_exists
+
 __author__ = "Anders S. Christensen"
 __copyright__ = "Copyright 2016"
 __credits__ = ["Anders S. Christensen (2016) https://github.com/qmlcode/qml"]
@@ -16,11 +18,16 @@
 FORTRAN = "f90"
 
 # GNU (default)
-COMPILER_FLAGS = ["-O3", "-fopenmp", "-m64", "-march=native", "-fPIC", 
+COMPILER_FLAGS = ["-O3", "-fopenmp", "-m64", "-march=native", "-fPIC",
                     "-Wno-maybe-uninitialized", "-Wno-unused-function", "-Wno-cpp"]
 LINKER_FLAGS = ["-lgomp"]
 MATH_LINKER_FLAGS = ["-lblas", "-llapack"]
 
+# UNCOMMENT TO FORCE LINKING TO MKL with GNU compilers:
+if mkl_exists(verbose=True):
+    LINKER_FLAGS = ["-lgomp", " -lpthread", "-lm", "-ldl"]
+    MATH_LINKER_FLAGS = ["-L${MKLROOT}/lib/intel64", "-lmkl_rt"]
+
 # For clang without OpenMP: (i.e. most Apple/mac system)
 if sys.platform == "darwin" and all(["gnu" not in arg for arg in sys.argv]):
     COMPILER_FLAGS = ["-O3", "-m64", "-march=native", "-fPIC"]
@@ -35,11 +42,38 @@
     MATH_LINKER_FLAGS = ["-L${MKLROOT}/lib/intel64", "-lmkl_rt"]
 
 
-# UNCOMMENT TO FORCE LINKING TO MKL with GNU compilers:
-# LINKER_FLAGS = ["-lgomp", " -lpthread", "-lm", "-ldl"]
-# MATH_LINKER_FLAGS = ["-L${MKLROOT}/lib/intel64", "-lmkl_rt"]
 
 
+ext_ffchl_module = Extension(name = 'ffchl_module',
+                          sources = [
+                                'qml/ffchl_module.f90',
+                                'qml/ffchl_scalar_kernels.f90',
+                            ],
+                          extra_f90_compile_args = COMPILER_FLAGS,
+                          extra_f77_compile_args = COMPILER_FLAGS,
+                          extra_compile_args = COMPILER_FLAGS ,
+                          extra_link_args = LINKER_FLAGS + MATH_LINKER_FLAGS,
+                          language = FORTRAN,
+                          f2py_options=['--quiet'])
+
+ext_ffchl_scalar_kernels = Extension(name = 'ffchl_scalar_kernels',
+                          sources = ['qml/ffchl_scalar_kernels.f90'],
+                          extra_f90_compile_args = COMPILER_FLAGS,
+                          extra_f77_compile_args = COMPILER_FLAGS,
+                          extra_compile_args = COMPILER_FLAGS,
+                          extra_link_args = LINKER_FLAGS + MATH_LINKER_FLAGS,
+                          language = FORTRAN,
+                          f2py_options=['--quiet'])
+
+ext_ffchl_vector_kernels = Extension(name = 'ffchl_vector_kernels',
+                          sources = ['qml/ffchl_vector_kernels.f90'],
+                          extra_f90_compile_args = COMPILER_FLAGS,
+                          extra_f77_compile_args = COMPILER_FLAGS,
+                          extra_compile_args = COMPILER_FLAGS,
+                          extra_link_args = LINKER_FLAGS + MATH_LINKER_FLAGS,
+                          language = FORTRAN,
+                          f2py_options=['--quiet'])
+
 ext_farad_kernels = Extension(name = 'farad_kernels',
                           sources = ['qml/farad_kernels.f90'],
                           extra_f90_compile_args = COMPILER_FLAGS,
@@ -85,6 +119,15 @@
                           language = FORTRAN,
                           f2py_options=['--quiet'])
 
+ext_fslatm = Extension(name = 'fslatm',
+                          sources = ['qml/fslatm.f90'],
+                          extra_f90_compile_args = COMPILER_FLAGS,
+                          extra_f77_compile_args = COMPILER_FLAGS,
+                          extra_compile_args = COMPILER_FLAGS,
+                          extra_link_args = LINKER_FLAGS,
+                          language = FORTRAN,
+                          f2py_options=['--quiet'])
+
 # use README.md as long description
 def readme():
     with open('README.md') as f:
@@ -112,10 +155,12 @@ def setup_pepytools():
 
         ext_package = 'qml',
         ext_modules = [
+              ext_ffchl_module,
               ext_farad_kernels,
               ext_fcho_solve,
               ext_fdistance,
               ext_fkernels,
+              ext_fslatm,
               ext_frepresentations,
         ],
 )

tests/data/eigenvalue_coulomb_matrix_representation.txt

@@ -0,0 +1,10 @@
+1.250726879e+02  4.540938889e+01  1.698020351e+01  1.431822726e+01  3.906495595e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  -3.867537429e-02  -2.820199687e-01  -4.573365589e-01  -7.079149287e-01  -9.099281091e-01  
+1.151615868e+02  4.881885871e+01  2.395834183e+01  1.287122167e+01  9.394697151e+00  1.293529945e-01  0.000000000e+00  -7.025128458e-02  -7.125064768e-02  -1.246546160e-01  -1.535039349e-01  -1.696710371e-01  -3.258540252e-01  -5.140916554e-01  -6.765120088e-01  -8.950223396e-01  -1.042119461e+00  
+1.189056060e+02  5.103243144e+01  2.301273537e+01  8.614411823e+00  3.808768061e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  -1.157426401e-03  -1.092926673e-01  -2.639914590e-01  -6.810228753e-01  -1.027360084e+00  
+1.192849858e+02  5.186714311e+01  2.372163699e+01  5.841982021e+00  3.176851825e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  -6.895415313e-02  -5.484217301e-01  -9.840956374e-01  
+1.167538495e+02  4.839740821e+01  2.418859595e+01  1.304863135e+01  6.400942465e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  -3.308768787e-02  -6.862271937e-02  -7.419078992e-02  -1.408210145e-01  -2.242849025e-01  -4.522285057e-01  -5.781734135e-01  -8.726026226e-01  -1.054287648e+00  
+1.188911764e+02  5.071123148e+01  2.326905421e+01  9.448295930e+00  3.069163128e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  -6.826059037e-02  -6.829293604e-02  -2.652259416e-01  -6.507763346e-01  -1.045237094e+00  
+1.206708186e+02  4.251229257e+01  2.395892026e+01  1.585796038e+01  5.759085611e+00  5.265555652e-03  0.000000000e+00  0.000000000e+00  0.000000000e+00  -1.121612652e-02  -5.757858887e-02  -1.537427408e-01  -1.705489260e-01  -4.850240279e-01  -7.287242773e-01  -8.282990379e-01  -1.038081066e+00  
+1.222328141e+02  4.556024708e+01  2.404558706e+01  1.090681002e+01  2.506181904e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  0.000000000e+00  -7.460984648e-03  -6.869045466e-02  -1.141822672e-01  -7.397383525e-01  -1.030439883e+00  
+1.195308461e+02  4.221302919e+01  2.410927785e+01  1.623264179e+01  8.120913601e+00  7.113442198e-02  0.000000000e+00  -2.667959114e-02  -8.606617808e-02  -1.094679848e-01  -1.350172283e-01  -1.672342174e-01  -2.583822389e-01  -5.128070226e-01  -8.160522987e-01  -8.478784498e-01  -1.027129517e+00  
+1.239972946e+02  3.074763292e+01  2.399739643e+01  2.399499946e+01  7.471295657e+00  3.472090476e-02  1.681756968e-02  0.000000000e+00  -2.301165158e-02  -1.182229236e-01  -1.326115900e-01  -1.581530629e-01  -1.830420736e-01  -7.075857805e-01  -8.165534186e-01  -8.204997319e-01  -1.009349091e+00  

tests/data/y_cho_solve.txt

@@ -0,0 +1,200 @@
+-1.620680000000000064e+03
+-1.623029999999999973e+03
+-1.314670000000000073e+03
+-1.757589999999999918e+03
+-1.423589999999999918e+03
+-1.503670000000000073e+03
+-1.606660000000000082e+03
+-1.374029999999999973e+03
+-1.551160000000000082e+03
+-1.565059999999999945e+03
+-1.362329999999999927e+03
+-1.751089999999999918e+03
+-1.274170000000000073e+03
+-1.484069999999999936e+03
+-1.619380000000000109e+03
+-1.354339999999999918e+03
+-1.122490000000000009e+03
+-1.618779999999999973e+03
+-1.015860000000000014e+03
+-1.597410000000000082e+03
+-1.631990000000000009e+03
+-1.215059999999999945e+03
+-1.475160000000000082e+03
+-1.223019999999999982e+03
+-1.481029999999999973e+03
+-1.468859999999999900e+03
+-1.419259999999999991e+03
+-1.427509999999999991e+03
+-1.552740000000000009e+03
+-8.690960000000000036e+02
+-1.319079999999999927e+03
+-1.550259999999999991e+03
+-1.312400000000000091e+03
+-1.201980000000000018e+03
+-1.749150000000000091e+03
+-1.274930000000000064e+03
+-1.246619999999999891e+03
+-1.018419999999999959e+03
+-1.408779999999999973e+03
+-1.892740000000000009e+03
+-1.339890000000000100e+03
+-1.618609999999999900e+03
+-1.615450000000000045e+03
+-8.299379999999999882e+02
+-1.427359999999999900e+03
+-1.218000000000000000e+03
+-1.387460000000000036e+03
+-1.617269999999999982e+03
+-1.405730000000000018e+03
+-1.485869999999999891e+03
+-1.412180000000000064e+03
+-1.542200000000000045e+03
+-1.292559999999999945e+03
+-1.160720000000000027e+03
+-1.214240000000000009e+03
+-1.604859999999999900e+03
+-1.156960000000000036e+03
+-1.427509999999999991e+03
+-1.190119999999999891e+03
+-1.593500000000000000e+03
+-1.193900000000000091e+03
+-1.483720000000000027e+03
+-1.567079999999999927e+03
+-1.232490000000000009e+03
+-9.582459999999999809e+02
+-1.313950000000000045e+03
+-1.324789999999999964e+03
+-1.236809999999999945e+03
+-1.304039999999999964e+03
+-1.564019999999999982e+03
+-1.366839999999999918e+03
+-1.246589999999999918e+03
+-1.358829999999999927e+03
+-1.089440000000000055e+03
+-1.094210000000000036e+03
+-1.262619999999999891e+03
+-1.254670000000000073e+03
+-1.073349999999999909e+03
+-1.259480000000000018e+03
+-1.460960000000000036e+03
+-1.548589999999999918e+03
+-1.065210000000000036e+03
+-1.489920000000000073e+03
+-1.502990000000000009e+03
+-1.613589999999999918e+03
+-1.247680000000000064e+03
+-1.063880000000000109e+03
+-7.240489999999999782e+02
+-1.487309999999999945e+03
+-1.354900000000000091e+03
+-1.223349999999999909e+03
+-1.314420000000000073e+03
+-1.208299999999999955e+03
+-1.316880000000000109e+03
+-1.075410000000000082e+03
+-1.357829999999999927e+03
+-1.242680000000000064e+03
+-1.692319999999999936e+03
+-1.621019999999999982e+03
+-9.469809999999999945e+02
+-1.614200000000000045e+03
+-1.294980000000000018e+03
+-1.239869999999999891e+03
+-1.491369999999999891e+03
+-1.221779999999999973e+03
+-1.274450000000000045e+03
+-1.192150000000000091e+03
+-1.102049999999999955e+03
+-1.626779999999999973e+03
+-1.472829999999999927e+03
+-1.746630000000000109e+03
+-1.621529999999999973e+03
+-1.622410000000000082e+03
+-1.417619999999999891e+03
+-8.701760000000000446e+02
+-1.268380000000000109e+03
+-1.603259999999999991e+03
+-1.540599999999999909e+03
+-1.225049999999999955e+03
+-1.766819999999999936e+03
+-1.617650000000000091e+03
+-1.472019999999999982e+03
+-1.174049999999999955e+03
+-1.455650000000000091e+03
+-1.230190000000000055e+03
+-1.419059999999999945e+03
+-1.307369999999999891e+03
+-1.316980000000000018e+03
+-1.213740000000000009e+03
+-1.226049999999999955e+03
+-1.341470000000000027e+03
+-1.332970000000000027e+03
+-1.228190000000000055e+03
+-1.384039999999999964e+03
+-1.228829999999999927e+03
+-1.628349999999999909e+03
+-1.489329999999999927e+03
+-1.451309999999999945e+03
+-1.506309999999999945e+03
+-1.422140000000000100e+03
+-1.007370000000000005e+03
+-1.319789999999999964e+03
+-1.024740000000000009e+03
+-1.429380000000000109e+03
+-1.153009999999999991e+03
+-1.232880000000000109e+03
+-1.269279999999999973e+03
+-1.057710000000000036e+03
+-1.182099999999999909e+03
+-1.194319999999999936e+03
+-1.613680000000000064e+03
+-1.167599999999999909e+03
+-1.348190000000000055e+03
+-1.374549999999999955e+03
+-1.224140000000000100e+03
+-9.713200000000000500e+02
+-1.158740000000000009e+03
+-1.629089999999999918e+03
+-1.758170000000000073e+03
+-1.104319999999999936e+03
+-1.121480000000000018e+03
+-1.697569999999999936e+03
+-1.265269999999999982e+03
+-1.364519999999999982e+03
+-1.689690000000000055e+03
+-9.284790000000000418e+02
+-1.612400000000000091e+03
+-1.124779999999999973e+03
+-1.098869999999999891e+03
+-1.462680000000000064e+03
+-9.756960000000000264e+02
+-1.348420000000000073e+03
+-1.541660000000000082e+03
+-1.242730000000000018e+03
+-1.255970000000000027e+03
+-1.408680000000000064e+03
+-1.487529999999999973e+03
+-1.480559999999999945e+03
+-8.064950000000000045e+02
+-1.106400000000000091e+03
+-1.549269999999999982e+03
+-1.369900000000000091e+03
+-1.510900000000000091e+03
+-1.603650000000000091e+03
+-1.378970000000000027e+03
+-1.325579999999999927e+03
+-1.355849999999999909e+03
+-1.684599999999999909e+03
+-1.243000000000000000e+03
+-1.546710000000000036e+03
+-1.185190000000000055e+03
+-1.757579999999999927e+03
+-1.611410000000000082e+03
+-1.471180000000000064e+03
+-1.616240000000000009e+03
+-1.295740000000000009e+03
+-1.748529999999999973e+03
+-9.805030000000000427e+02
+-1.454150000000000091e+03
+-1.121940000000000055e+03

tests/test_arad.py

@@ -13,6 +13,9 @@
 from qml.arad import get_local_kernels_arad
 from qml.arad import get_local_symmetric_kernels_arad
 
+from qml.arad import get_global_kernels_arad
+from qml.arad import get_global_symmetric_kernels_arad
+
 from qml.arad import get_atomic_kernels_arad
 from qml.arad import get_atomic_symmetric_kernels_arad
 
@@ -70,7 +73,11 @@ def test_arad():
     assert np.allclose(K_local_symm, K_local_asymm), "Symmetry error in local kernels"
     assert np.invert(np.all(np.isnan(K_local_asymm))), "ERROR: ARAD local symmetric kernel contains NaN"
 
-    K_local_asymm = get_local_kernels_arad(X1[-4:], X1[:6], sigmas)
+    K_global_asymm = get_global_kernels_arad(X1, X1, sigmas)
+    K_global_symm = get_global_symmetric_kernels_arad(X1, sigmas)
+
+    assert np.allclose(K_global_symm, K_global_asymm), "Symmetry error in global kernels"
+    assert np.invert(np.all(np.isnan(K_global_asymm))), "ERROR: ARAD global symmetric kernel contains NaN"
 
     molid = 5
     X1 = generate_arad_representation(mols[molid].coordinates,

tests/test_energy_krr_atomic_cmat.py

@@ -0,0 +1,222 @@
+# MIT License
+#
+# Copyright (c) 2017 Anders Steen Christensen
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+from __future__ import print_function
+
+import os
+import numpy as np
+import qml
+from qml.kernels import laplacian_kernel
+from qml.math import cho_solve
+
+from qml.representations import get_slatm_mbtypes
+
+from qml.kernels import get_local_kernels_gaussian
+from qml.kernels import get_local_kernels_laplacian
+
+from qml.wrappers import get_atomic_kernels_gaussian
+from qml.wrappers import get_atomic_kernels_laplacian
+
+def get_energies(filename):
+    """ Returns a dictionary with heats of formation for each xyz-file.
+    """
+
+    f = open(filename, "r")
+    lines = f.readlines()
+    f.close()
+
+    energies = dict()
+
+    for line in lines:
+        tokens = line.split()
+
+        xyz_name = tokens[0]
+        hof = float(tokens[1])
+
+        energies[xyz_name] = hof
+
+    return energies
+
+def test_krr_gaussian_local_cmat():
+
+    test_dir = os.path.dirname(os.path.realpath(__file__))
+
+    # Parse file containing PBE0/def2-TZVP heats of formation and xyz filenames
+    data = get_energies(test_dir + "/data/hof_qm7.txt")
+
+    # Generate a list of qml.Compound() objects"
+    mols = []
+
+
+    for xyz_file in sorted(data.keys())[:1000]:
+
+        # Initialize the qml.Compound() objects
+        mol = qml.Compound(xyz=test_dir + "/qm7/" + xyz_file)
+
+        # Associate a property (heat of formation) with the object
+        mol.properties = data[xyz_file]
+
+        # This is a Molecular Coulomb matrix sorted by row norm
+        mol.generate_atomic_coulomb_matrix(size=23, sorting="row-norm")
+
+        mols.append(mol)
+
+    # Shuffle molecules
+    np.random.seed(666)
+    np.random.shuffle(mols)
+
+    # Make training and test sets
+    n_test  = 100
+    n_train = 200
+
+    training = mols[:n_train]
+    test  = mols[-n_test:]
+
+    X = np.concatenate([mol.representation for mol in training])
+    Xs = np.concatenate([mol.representation for mol in test])
+
+    N = np.array([mol.natoms for mol in training])
+    Ns = np.array([mol.natoms for mol in test])
+
+    # List of properties
+    Y = np.array([mol.properties for mol in training])
+    Ys = np.array([mol.properties for mol in test])
+
+    # Set hyper-parameters
+    sigma = 724.0
+    llambda = 10**(-6.5)
+
+    K = get_local_kernels_gaussian(X, X, N, N, [sigma])[0]
+    assert np.allclose(K, K.T), "Error in local Gaussian kernel symmetry"
+
+    K_test = np.loadtxt(test_dir + "/data/K_local_gaussian.txt")
+    assert np.allclose(K, K_test), "Error in local Gaussian kernel (vs. reference)"
+
+    K_test = get_atomic_kernels_gaussian(training, training, [sigma])[0]
+    assert np.allclose(K, K_test), "Error in local Gaussian kernel (vs. wrapper)"
+
+    # Solve alpha
+    K[np.diag_indices_from(K)] += llambda
+    alpha = cho_solve(K,Y)
+
+    # Calculate prediction kernel
+    Ks = get_local_kernels_gaussian(Xs, X, Ns, N, [sigma])[0]
+
+    Ks_test = np.loadtxt(test_dir + "/data/Ks_local_gaussian.txt")
+    # Somtimes a few coulomb matrices differ because of parallel sorting and numerical error
+    # Allow up to 5 molecules to differ from the supplied reference.
+    differences_count = len(set(np.where(Ks - Ks_test > 1e-7)[0]))
+    assert differences_count < 5, "Error in local Laplacian kernel (vs. reference)"
+    # assert np.allclose(Ks, Ks_test), "Error in local Gaussian kernel (vs. reference)"
+
+    Ks_test = get_atomic_kernels_gaussian(test, training, [sigma])[0]
+    assert np.allclose(Ks, Ks_test), "Error in local Gaussian kernel (vs. wrapper)"
+
+    Yss = np.dot(Ks, alpha)
+
+    mae = np.mean(np.abs(Ys - Yss))
+    assert abs(19.0 - mae) < 1.0, "Error in local Gaussian kernel-ridge regression"
+
+def test_krr_laplacian_local_cmat():
+
+    test_dir = os.path.dirname(os.path.realpath(__file__))
+
+    # Parse file containing PBE0/def2-TZVP heats of formation and xyz filenames
+    data = get_energies(test_dir + "/data/hof_qm7.txt")
+
+    # Generate a list of qml.Compound() objects"
+    mols = []
+
+
+    for xyz_file in sorted(data.keys())[:1000]:
+
+        # Initialize the qml.Compound() objects
+        mol = qml.Compound(xyz=test_dir + "/qm7/" + xyz_file)
+
+        # Associate a property (heat of formation) with the object
+        mol.properties = data[xyz_file]
+
+        # This is a Molecular Coulomb matrix sorted by row norm
+        mol.generate_atomic_coulomb_matrix(size=23, sorting="row-norm")
+
+        mols.append(mol)
+
+    # Shuffle molecules
+    np.random.seed(666)
+    np.random.shuffle(mols)
+
+    # Make training and test sets
+    n_test  = 100
+    n_train = 200
+
+    training = mols[:n_train]
+    test  = mols[-n_test:]
+
+    X = np.concatenate([mol.representation for mol in training])
+    Xs = np.concatenate([mol.representation for mol in test])
+
+    N = np.array([mol.natoms for mol in training])
+    Ns = np.array([mol.natoms for mol in test])
+
+    # List of properties
+    Y = np.array([mol.properties for mol in training])
+    Ys = np.array([mol.properties for mol in test])
+
+    # Set hyper-parameters
+    sigma = 10**(3.6)
+    llambda = 10**(-12.0)
+
+    K = get_local_kernels_laplacian(X, X, N, N, [sigma])[0]
+    assert np.allclose(K, K.T), "Error in local Laplacian kernel symmetry"
+
+    K_test = np.loadtxt(test_dir + "/data/K_local_laplacian.txt")
+    assert np.allclose(K, K_test), "Error in local Laplacian kernel (vs. reference)"
+
+    K_test = get_atomic_kernels_laplacian(training, training, [sigma])[0]
+    assert np.allclose(K, K_test), "Error in local Laplacian kernel (vs. wrapper)"
+
+    # Solve alpha
+    K[np.diag_indices_from(K)] += llambda
+    alpha = cho_solve(K,Y)
+
+    # Calculate prediction kernel
+    Ks = get_local_kernels_laplacian(Xs, X, Ns, N, [sigma])[0]
+
+    Ks_test = np.loadtxt(test_dir + "/data/Ks_local_laplacian.txt")
+
+    # Somtimes a few coulomb matrices differ because of parallel sorting and numerical error
+    # Allow up to 5 molecules to differ from the supplied reference.
+    differences_count = len(set(np.where(Ks - Ks_test > 1e-7)[0]))
+    assert differences_count < 5, "Error in local Laplacian kernel (vs. reference)"
+
+    Ks_test = get_atomic_kernels_laplacian(test, training, [sigma])[0]
+    assert np.allclose(Ks, Ks_test), "Error in local Laplacian kernel (vs. wrapper)"
+
+    Yss = np.dot(Ks, alpha)
+
+    mae = np.mean(np.abs(Ys - Yss))
+    assert abs(8.7 - mae) < 1.0, "Error in local Laplacian kernel-ridge regression"
+
+if __name__ == "__main__":
+
+    test_krr_gaussian_local_cmat()
+    test_krr_laplacian_local_cmat()

tests/test_energy_krr_bob.py

@@ -0,0 +1,117 @@
+# MIT License
+#
+# Copyright (c) 2017 Anders Steen Christensen
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+from __future__ import print_function
+
+import os
+import numpy as np
+import qml
+from qml.kernels import laplacian_kernel
+from qml.math import cho_solve
+
+from qml.representations import get_slatm_mbtypes
+
+
+def get_energies(filename):
+    """ Returns a dictionary with heats of formation for each xyz-file.
+    """
+
+    f = open(filename, "r")
+    lines = f.readlines()
+    f.close()
+
+    energies = dict()
+
+    for line in lines:
+        tokens = line.split()
+
+        xyz_name = tokens[0]
+        hof = float(tokens[1])
+
+        energies[xyz_name] = hof
+
+    return energies
+
+def test_krr_bob():
+
+    test_dir = os.path.dirname(os.path.realpath(__file__))
+
+    # Parse file containing PBE0/def2-TZVP heats of formation and xyz filenames
+    data = get_energies(test_dir + "/data/hof_qm7.txt")
+
+    # Generate a list of qml.Compound() objects
+    mols = []
+
+    for xyz_file in sorted(data.keys())[:1000]:
+
+        # Initialize the qml.Compound() objects
+        mol = qml.Compound(xyz=test_dir + "/qm7/" + xyz_file)
+
+        # Associate a property (heat of formation) with the object
+        mol.properties = data[xyz_file]
+
+        # This is a Molecular Coulomb matrix sorted by row norm
+        mol.generate_bob()
+
+        mols.append(mol)
+
+    # Shuffle molecules
+    np.random.seed(666)
+    np.random.shuffle(mols)
+
+    # Make training and test sets
+    n_test  = 300
+    n_train = 700
+
+    training = mols[:n_train]
+    test  = mols[-n_test:]
+
+    # List of representations
+    X  = np.array([mol.representation for mol in training])
+    Xs = np.array([mol.representation for mol in test])
+
+    # List of properties
+    Y = np.array([mol.properties for mol in training])
+    Ys = np.array([mol.properties for mol in test])
+
+    # Set hyper-parameters
+    sigma = 26214.40
+    llambda = 1e-10
+
+    # Generate training Kernel
+    K = laplacian_kernel(X, X, sigma)
+
+    # Solve alpha
+    K[np.diag_indices_from(K)] += llambda
+    alpha = cho_solve(K,Y)
+
+    # Calculate prediction kernel
+    Ks = laplacian_kernel(X, Xs, sigma)
+    Yss = np.dot(Ks.transpose(), alpha)
+
+    mae = np.mean(np.abs(Ys - Yss))
+    print(mae)
+    assert mae < 2.6, "ERROR: Too high MAE!"
+
+if __name__ == "__main__":
+
+    test_krr_bob()

tests/test_energy_krr_cmat.py

@@ -109,4 +109,9 @@ def test_krr_cmat():
     Yss = np.dot(Ks.transpose(), alpha)
 
     mae = np.mean(np.abs(Ys - Yss))
-    print(mae)
+
+    assert mae < 6.0, "ERROR: Too high MAE!"
+
+if __name__ == "__main__":
+
+    test_krr_cmat()