[WIP] Development branch with many body feature and energy mapping

.gitignore

@@ -11,7 +11,9 @@ tags
 pwscf.save/
 /docs/source/_build/html
 /docs/build/
+/docs/doctrees/
 /datasets/
+/tests/otf_data/
 /tests/test_ase/out/
 /tests/test_ase/Al/
 /tests/test_ase/__pycache__/
@@ -30,6 +32,17 @@ test_files/qe_output_*.out
 /tests/pwscf.xml
 /tests/*-RESTART.kp*
 /tests/*.wfn*
+/tests/*.pwi*
+/tests/*.pwo*
 /tests/*-bak
 /tests/cp2k.in*
 /tests/*.in_run
+/tests/tmp*
+/tests/grid3*npy
+/tests/AgI_Molten_15.txt
+/tests/log.lammps
+/tests/*.txt
+/tests/*.mgp
+/tests/*.log
+/tests/POSCAR*
+.idea/

.travis.yml

@@ -9,10 +9,9 @@ cache:
 
 before_install:
     - sudo apt update
-    - sudo apt install -y quantum-espresso tree
+    - sudo apt install -y quantum-espresso 
     - pip install codecov pytest pytest-cov
     - pwd
-    - tree
     - wget http://folk.uio.no/anjohan/lmp
     - "wget https://github.com/cp2k/cp2k/releases/download/\
       v6.1.0/cp2k-6.1-Linux-x86_64.sopt"
@@ -22,10 +21,10 @@ before_install:
 
 script:
     - pwd
-    - tree
     - cd tests
-    - PWSCF_COMMAND=pw.x CP2K_COMMAND=../cp2k-6.1-Linux-x86_64.sopt lmp=../lmp
-      pytest --cov=../flare/
+    - PWSCF_COMMAND=pw.x CP2K_COMMAND=../cp2k-6.1-Linux-x86_64.sopt 
+      lmp=$(pwd)/../lmp
+      pytest --show-capture=all -vv --durations=0 --cov=../flare/
     - coverage xml
 
 after_success:

docs/source/FAQs/all.rst

@@ -35,7 +35,7 @@ Gaussian Processes and OTF
 
         * `maxiter`: maximal number of iterations, usually set to ~10 to prevent training for too long.
 
-        * `par`: if `True`, then parallelization is used in optimization. 
+        * `parallel`: if `True`, then parallelization is used in optimization. 
           The serial version could be very slow.
 
         * `output` (in `train` function): set up an output file for monitoring optimization iterations.

docs/source/flare/gp.rst

@@ -6,7 +6,7 @@ Gaussian Process Force Fields
 
    struc
    env
-   kernels/kernels
+   kernels/sc
    kernels/mc_simple
    kernels/mc_sephyps
    cutoffs

docs/source/flare/kernels/kernels.rst → docs/source/flare/kernels/sc.rst

@@ -1,5 +1,5 @@
 Single-element Kernels
 ======================
 
-.. automodule:: flare.kernels.kernels
+.. automodule:: flare.kernels.sc
     :members:

docs/source/flare/kernels/util.rst

@@ -0,0 +1,5 @@
+Kernel interface (for developer)
+======================
+
+.. automodule:: flare.kernels.utils
+    :members:

docs/source/install.rst

@@ -85,8 +85,8 @@ Second, in the initialization of the GP class and OTF class, turn on the GP para
 
 .. code-block:: python
 
-    gp_model = GaussianProcess(..., n_cpus=2)
-    otf_instance = OTF(..., par=False, n_cpus=2)
+    gp_model = GaussianProcess(..., parallel=True, per_atom_par=False, n_cpus=2)
+    otf_instance = OTF(..., par, n_cpus=2)
 
 Third, set the number of threads for MKL before running your python script.
 
@@ -100,7 +100,7 @@ Third, set the number of threads for MKL before running your python script.
    If these numbers are larger than the actual CPUs number, it can lead to an overload of the machine.
 
 .. note::
-   If OTF.par=True and GaussianProcess.n_cpus>1, it is equivalent to run with n_cpu**2 threads
+   If gp_model.per_atom_par=True and NUM_OMB_THREAD>1, it is equivalent to run with NUM_OMB_THREAD*otf.n_cpus threads
    because the MKL calls are nested in the multiprocessing code. 
 
 The current version of FLARE can only support parallel calculations within one compute node.

docs/source/tutorials/after_training.py

@@ -8,13 +8,13 @@
 # -------------------------------------------------------------------------
 #                  reconstruct gp model from otf snippet
 # -------------------------------------------------------------------------
-import flare.kernels.mc_simple as mc_simple
 
-kernel = mc_simple.two_plus_three_body_mc
-kernel_grad = mc_simple.two_plus_three_body_mc_grad
-gp_model = otf_object.make_gp(kernel=kernel, kernel_grad=kernel_grad,
+
+
+
+gp_model = otf_object.make_gp(kernel_name="2+3_mc",
                               hyp_no=hyp_no)
-gp_model.par = True
+gp_model.parallel = True
 gp_model.hyp_labels = ['sig2', 'ls2', 'sig3', 'ls3', 'noise']
 
 # write model to a binary file

docs/source/tutorials/gpfa.rst

@@ -18,7 +18,7 @@ Step 1: Setting up a Gaussian Process Object
 Our goal is to train a GP, which first must be instantiated with a set of parameters.
 
 For the sake of this example, which is a molecule, we will use a two-plus-three body kernel. 
-We must provide the kernel and the kernel gradient as ``callable`` s to the GP. 
+We must provide the kernel name to the GP, "two-plus-three-mc" or "2+3mc'.
 Our initial guesses for the hyperparameters are not important. 
 The hyperparameter labels are included below for later output.
 The system contains a small number of atoms, so we choose a relatively 
@@ -30,7 +30,6 @@ At the header of a file, include the following imports:
 .. code-block:: python
 
 	from flare.gp import GaussianProcess
-	from flare.kernels.mc_simple import two_plus_three_body_mc, two_plus_three_body_mc_grad
 
 We will then set up the ``GaussianProcess`` object.
 
@@ -47,17 +46,15 @@ We will then set up the ``GaussianProcess`` object.
 * | The two cutoff values correspond to the functions which set up 
   | the two- and three-body Atomic Environments. Since Methanol is a small 
   | molecule, 7 Angstrom each will be sufficent.
-* | The kernels which facilitate these comparisons must be imported as Python  ``callable`` s. 
+* | The kernel name must contain the terms you want to use. 
 * | Here, we will use the ``two_plus_three_body_mc`` kernel, which 
   | uses two-body and three-body comparisons. ``mc`` means multi-component, 
   | indicating that it can handle multiple atomic species being present.
-* | We must also import the gradient of the kernel, which is
-  | ``two_plus_three_body_mc_grad``.
 
 
 .. code-block:: python
 
-	gp = GaussianProcess(kernel=two_plus_three_body_mc, kernel_grad=two_plus_three_body_mc_grad,
+	gp = GaussianProcess(kernel_name="2+3mc", 
 	hyps=[0.01, 0.01, 0.01, 0.01, 0.01],
 	cutoffs = (7,7),
 	hyp_labels=['Two-Body Signal Variance','Two-Body Length Scale','Three-Body Signal Variance',

docs/source/tutorials/otf_al.rst

@@ -8,42 +8,35 @@ Step 1: Set up a GP Model
 -------------------------
 
 Let's start up with the GP model with three-body kernel function. 
-(See :doc:`kernels.py <../flare/kernels/kernels>` (single component)
+(See :doc:`kernels.py <../flare/kernels/sc>` (single component)
 or :doc:`mc_simple.py <../flare/kernels/mc_simple>` (multi-component) for more options.)
 
 .. code-block:: python
     :linenos:
 
     # make gp model
-    kernel = en.three_body
-    kernel_grad = en.three_body_grad
     hyps = np.array([0.1, 1, 0.01])
     hyp_labels = ['Signal Std', 'Length Scale', 'Noise Std']
     cutoffs = np.array([3.9, 3.9])
-    energy_force_kernel = en.three_body_force_en
 
     gp = \
-        GaussianProcess(kernel=kernel,
-                        kernel_grad=kernel_grad,
+        GaussianProcess(kernel_name='3b',
                         hyps=hyps,
                         cutoffs=cutoffs,
                         hyp_labels=hyp_labels,
-                        energy_force_kernel=energy_force_kernel,
                         maxiter=50)
 
 
 **Some Explanation about the parameters:**
 
-* ``kernel``: set to be the kernel function  
+* ``kernel_name``: set to be the name of kernel functions
 
-    * import from :doc:`kernels.py <../flare/kernels/kernels>` (single-component system) 
+    * import from :doc:`sc.py <../flare/kernels/sc>` (single-component system) 
       or :doc:`mc_simple.py <../flare/kernels/mc_simple>` (multi-component system). 
     * Currently we have the choices of two-body, three-body and two-plus-three-body kernel functions.
     * Two-plus-three-body kernel function is simply the summation of two-body and three-body kernels,
       and is tested to have best performance.
 
-* ``kernel_grad``: set to be the gradient of kernel function used for hyperparameter training. 
-
 * ``hyps``: the array of hyperparameters, whose names are shown in ``hyp_labels``.
 
     * For two-body kernel function, an array of length 3 is needed, ``hyps=[sigma_2, ls_2, sigma_n]``;
@@ -53,21 +46,13 @@ or :doc:`mc_simple.py <../flare/kernels/mc_simple>` (multi-component) for more o
 * ``cutoffs``: consists of two values. The 1st is the cutoff of two-body and the 2nd is for three-body kernel. 
   Usually we will set a larger one for two-body.
 
-* ``energy_force_kernel``: set to calculate local energy for each atom based on the integral of forces. 
-
 * ``maxiter``: set to constrain the number of steps in training hyperparameters. 
 
 
 **Note:**
 
 1. See :doc:`GaussianProcess <../flare/gp>` for complete description of arguments of ``GaussianProcess`` class.
 
-2. In any case, you need to set up ``kernel``. 
-
-    a. If you want to train the hyperparameters, you need to set up ``kernel_grad``;
-    b. If you want the output of local energy, you need to set up ``energy_force_kernel`` and ``energy_kernel``, 
-       otherwise you don't need them.
-
 
 Step 2: Set up DFT Calculator
 -----------------------------

flare/ase/calculator.py

@@ -37,8 +37,6 @@ def get_property(self, name, atoms=None, allow_calculation=True):
 
 
     def get_potential_energy(self, atoms=None, force_consistent=False):
-        if self.use_mapping:
-            print('MGP energy mapping not implemented, temporarily set to 0')
         return self.get_property('energy', atoms)
 
 
@@ -99,32 +97,30 @@ def calculate_mgp_serial(self, atoms):
                                atoms.get_atomic_numbers(),
                                atoms.positions)
 
-        forces = np.zeros((nat, 3))
-        stress = np.zeros((nat, 6))
-        stds = np.zeros((nat, 3))
+        self.results['forces'] = np.zeros((nat, 3))
+        self.results['stresses'] = np.zeros((nat, 6))
+        self.results['stds'] = np.zeros((nat, 3))
+        self.results['local_energies'] = np.zeros(nat)
         for n in range(nat):
             chemenv = AtomicEnvironment(struc_curr, n,
                                         self.mgp_model.cutoffs)
-            f, v, vir = self.mgp_model.predict(chemenv, mean_only=False)
-            forces[n] = f
-            stress[n] = vir
-            stds[n] = np.sqrt(np.absolute(v))
+            f, v, vir, e = self.mgp_model.predict(chemenv, mean_only=False)
+            self.results['forces'][n] = f
+            self.results['stresses'][n] = vir
+            self.results['stds'][n] = np.sqrt(np.absolute(v))
+            self.results['local_energies'][n] = e
 
-        self.results['forces'] = forces
-        self.results['stds'] = stds
-        self.results['stresses'] = stress
-        self.results['stress'] = np.sum(stress, axis=0)
-
-        # TODO: implement energy mapping
-        self.results['local_energies'] = np.zeros(forces.shape)
-        self.results['energy'] = 0
+        volume = atoms.get_volume()
+        total_stress = np.sum(self.results['stresses'], axis=0) 
+        self.results['stress'] = total_stress / volume
+        self.results['energy'] = np.sum(self.results['local_energies'])
 
         atoms.get_uncertainties = self.get_uncertainties
-        return forces
 
 
     def calculate_mgp_par(self, atoms):
-        return self.calculate_mgp_serial(atoms)
+        # TODO: to be done
+        self.calculate_mgp_serial(atoms)
 
 
     def calculation_required(self, atoms, quantities):