Add extrapolation level calcualtion

.github/workflows/cuda-build.yml

@@ -9,7 +9,8 @@ jobs:
   build:
     runs-on: ubuntu-latest
     container: 
-      image: nvidia/cuda:11.2.2-devel-ubuntu20.04  # 使用适合你的CUDA版本的官方NVIDIA容器
+      image: nvidia/cuda:12.2.2-devel-ubuntu22.04  # 使用适合你的CUDA版本的官方NVIDIA容器
+
 
     steps:
       - name: Checkout code

.github/workflows/cuda_build_windows.yml

@@ -21,12 +21,12 @@ jobs:
 
     - name: Install CUDA      
       run: |
-        choco install cuda --version=11.2.2.46133
+        choco install cuda --version=12.2.0.53625
 
     - name: Build project
       shell: cmd
       run: |
-        set PATH=%PATH%;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v11.2\bin
+        set PATH=%PATH%;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.2\bin
         nvcc --version
         cd src && make
 

doc/bibliography.rst

@@ -273,3 +273,14 @@ Bibliography
    | Computational Materials Science **207**, 111275 (2022)
    | DOI: `10.1016/j.commatsci.2022.111275 <https://doi.org/10.1016/j.commatsci.2022.111275>`_
 
+.. [Podryabinkin2023]
+   | Evgeny Podryabinkin, Kamil Garifullin, Alexander Shapeev, and Ivan Novikov
+   | *MLIP-3: Active learning on atomic environments with moment tensor potentials*
+   | J. Chem. Phys. **159**, 084112 (2023)
+   | DOI: `10.1063/5.0155887 <https://doi.org/10.1063/5.0155887>`_
+
+.. [Lysogorskiy2023]
+   | Yury Lysogorskiy, Anton Bochkarev, Matous Mrovec, and Ralf Drautz
+   | *Active learning strategies for atomic cluster expansion models*
+   | Phys. Rev. M **7**, 043801 (2023)
+   | DOI: `10.1103/PhysRevMaterials.7.043801 <https://doi.org/10.1103/PhysRevMaterials.7.043801>`_

doc/gpumd/input_parameters/compute_extrapolation.rst

@@ -0,0 +1,44 @@
+.. _kw_compute_extrapolation:
+.. index::
+   single: compute_extrapolation (keyword in run.in)
+
+:attr:`compute_extrapolation`
+=============================
+
+This keyword is used to compute the extrapolation grade of structures in an NEP potential.
+
+The extrapolation grade `gamma` can be considered as the uncertainty of a structure relative to the training set.
+
+A structure with large `gamma` tends to have higher energy and force errors.
+
+Similiar methods have been applied to MTP ([Podryabinkin2023]_) and ACE ([Lysogorskiy2023]_). You can refer to their papers for more details.
+
+Before computing `gamma`, you need to obtain an `active set` from your training set. There are some Python scripts to do it <https://github.com/psn417/nep_active>.
+
+There are also some Python scripts to perform active learning automatically <https://github.com/psn417/nep_maker>.
+
+Syntax
+------
+
+This keyword is used as follows::
+
+  compute_extrapolation asi_file <asi_file> gamma_low <gamma_low> gamma_high <gamma_high> check_interval <check_interval> dump_interval <dump_interval>
+
+:attr:`asi_file` is the name of the Active Set Inversion (ASI) file. This file is generate by the Python script in <https://github.com/psn417/nep_active>.
+
+:attr:`gamma_low`: Only if the max gamma value of a structure exceeds `gamma_low`, then the structure will be dumped into `extrapolation_dump.xyz` file. The default value is `0`.
+
+:attr:`gamma_high`: If the max gamma value of a structure exceeds `gamma_high`, then the simulation will stop. The default value is very large so it will never stop.
+
+:attr:`check_interval`: Since calculating gamma value is slow, you can check the gamma value every `check_interval` steps. The default value is `1` (check every step).
+
+:attr:`dump_interval`: You can set the minimum interval between dumps to `dump_interval` steps. The default value is `1`.
+
+Example
+-------
+
+.. code::
+
+    compute_extrapolation asi_file active_set.asi gamma_low 5 gamma_high 10 check_interval 10 dump_interval 10
+
+This means that the structures with max gamma between 5-10 will be dumped. The gamma value will be checked every 10 steps.

doc/gpumd/input_parameters/index.rst

@@ -15,6 +15,7 @@ Below you can find a listing of keywords for the ``run.in`` input file.
    velocity
    correct_velocity
    potential
+   compute_extrapolation
    dftd3
    change_box
    deform

src/force/nep.cu

@@ -27,11 +27,11 @@ heat transport, Phys. Rev. B. 104, 104309 (2021).
 #include "utilities/error.cuh"
 #include "utilities/gpu_macro.cuh"
 #include "utilities/nep_utilities.cuh"
+#include <cstring>
 #include <fstream>
 #include <iostream>
 #include <string>
 #include <vector>
-#include <cstring>
 
 const std::string ELEMENTS[NUM_ELEMENTS] = {
   "H",  "He", "Li", "Be", "B",  "C",  "N",  "O",  "F",  "Ne", "Na", "Mg", "Al", "Si", "P",  "S",
@@ -357,6 +357,7 @@ NEP::NEP(const char* file_potential, const int num_atoms)
 #endif
 
   initialize_dftd3();
+  B_projection_size = annmb.num_neurons1 * (annmb.dim + 2);
 }
 
 NEP::~NEP(void)
@@ -587,7 +588,10 @@ static __global__ void find_descriptor(
   double* g_pe,
   float* g_Fp,
   double* g_virial,
-  float* g_sum_fxyz)
+  float* g_sum_fxyz,
+  bool need_B_projection,
+  double* B_projection,
+  int B_projection_size)
 {
   int n1 = blockIdx.x * blockDim.x + threadIdx.x + N1;
   if (n1 < N2) {
@@ -695,7 +699,8 @@ static __global__ void find_descriptor(
         accumulate_s(paramb.L_max, d12, x12, y12, z12, gn12, s);
 #endif
       }
-      find_q(paramb.L_max, paramb.num_L, paramb.n_max_angular + 1, n, s, q + (paramb.n_max_radial + 1));
+      find_q(
+        paramb.L_max, paramb.num_L, paramb.n_max_angular + 1, n, s, q + (paramb.n_max_radial + 1));
       for (int abc = 0; abc < NUM_OF_ABC; ++abc) {
         g_sum_fxyz[(n * NUM_OF_ABC + abc) * N + n1] = s[abc];
       }
@@ -747,16 +752,29 @@ static __global__ void find_descriptor(
         F,
         Fp);
     } else {
-      apply_ann_one_layer(
-        annmb.dim,
-        annmb.num_neurons1,
-        annmb.w0[t1],
-        annmb.b0[t1],
-        annmb.w1[t1],
-        annmb.b1,
-        q,
-        F,
-        Fp);
+      if (!need_B_projection)
+        apply_ann_one_layer(
+          annmb.dim,
+          annmb.num_neurons1,
+          annmb.w0[t1],
+          annmb.b0[t1],
+          annmb.w1[t1],
+          annmb.b1,
+          q,
+          F,
+          Fp);
+      else
+        apply_ann_one_layer(
+          annmb.dim,
+          annmb.num_neurons1,
+          annmb.w0[t1],
+          annmb.b0[t1],
+          annmb.w1[t1],
+          annmb.b1,
+          q,
+          F,
+          Fp,
+          B_projection + n1 * B_projection_size);
     }
     g_pe[n1] += F;
 
@@ -979,7 +997,18 @@ static __global__ void find_partial_force_angular(
           g_gn_angular[index_left_all] * weight_left + g_gn_angular[index_right_all] * weight_right;
         float gnp12 = g_gnp_angular[index_left_all] * weight_left +
                       g_gnp_angular[index_right_all] * weight_right;
-        accumulate_f12(paramb.L_max, paramb.num_L, n, paramb.n_max_angular + 1, d12, r12, gn12, gnp12, Fp, sum_fxyz, f12);
+        accumulate_f12(
+          paramb.L_max,
+          paramb.num_L,
+          n,
+          paramb.n_max_angular + 1,
+          d12,
+          r12,
+          gn12,
+          gnp12,
+          Fp,
+          sum_fxyz,
+          f12);
       }
 #else
       float fc12, fcp12;
@@ -1007,7 +1036,18 @@ static __global__ void find_partial_force_angular(
           gn12 += fn12[k] * annmb.c[c_index];
           gnp12 += fnp12[k] * annmb.c[c_index];
         }
-        accumulate_f12(paramb.L_max, paramb.num_L, n, paramb.n_max_angular + 1, d12, r12, gn12, gnp12, Fp, sum_fxyz, f12);
+        accumulate_f12(
+          paramb.L_max,
+          paramb.num_L,
+          n,
+          paramb.n_max_angular + 1,
+          d12,
+          r12,
+          gn12,
+          gnp12,
+          Fp,
+          sum_fxyz,
+          f12);
       }
 #endif
       g_f12x[index] = f12[0];
@@ -1234,7 +1274,10 @@ void NEP::compute_large_box(
     potential_per_atom.data(),
     nep_data.Fp.data(),
     virial_per_atom.data(),
-    nep_data.sum_fxyz.data());
+    nep_data.sum_fxyz.data(),
+    need_B_projection,
+    B_projection,
+    B_projection_size);
   GPU_CHECK_KERNEL
 
   bool is_dipole = paramb.model_type == 1;
@@ -1416,7 +1459,10 @@ void NEP::compute_small_box(
     potential_per_atom.data(),
     nep_data.Fp.data(),
     virial_per_atom.data(),
-    nep_data.sum_fxyz.data());
+    nep_data.sum_fxyz.data(),
+    need_B_projection,
+    B_projection,
+    B_projection_size);
   GPU_CHECK_KERNEL
 
   bool is_dipole = paramb.model_type == 1;
@@ -1703,7 +1749,8 @@ static __global__ void find_descriptor(
         accumulate_s(paramb.L_max, d12, x12, y12, z12, gn12, s);
 #endif
       }
-      find_q(paramb.L_max, paramb.num_L, paramb.n_max_angular + 1, n, s, q + (paramb.n_max_radial + 1));
+      find_q(
+        paramb.L_max, paramb.num_L, paramb.n_max_angular + 1, n, s, q + (paramb.n_max_radial + 1));
       for (int abc = 0; abc < NUM_OF_ABC; ++abc) {
         g_sum_fxyz[(n * NUM_OF_ABC + abc) * N + n1] = s[abc];
       }

src/force/nep_small_box.cuh

@@ -16,8 +16,8 @@
 #include "model/box.cuh"
 #include "nep.cuh"
 #include "utilities/common.cuh"
-#include "utilities/nep_utilities.cuh"
 #include "utilities/gpu_macro.cuh"
+#include "utilities/nep_utilities.cuh"
 
 #ifdef USE_KEPLER
 #if defined(__CUDA_ARCH__) && (__CUDA_ARCH__ < 600)
@@ -166,7 +166,10 @@ static __global__ void find_descriptor_small_box(
   double* g_pe,
   float* g_Fp,
   double* g_virial,
-  float* g_sum_fxyz)
+  float* g_sum_fxyz,
+  bool need_B_projection,
+  double* B_projection,
+  int B_projection_size)
 {
   int n1 = blockIdx.x * blockDim.x + threadIdx.x + N1;
   if (n1 < N2) {
@@ -263,7 +266,8 @@ static __global__ void find_descriptor_small_box(
         accumulate_s(paramb.L_max, d12, r12[0], r12[1], r12[2], gn12, s);
 #endif
       }
-      find_q(paramb.L_max, paramb.num_L, paramb.n_max_angular + 1, n, s, q + (paramb.n_max_radial + 1));
+      find_q(
+        paramb.L_max, paramb.num_L, paramb.n_max_angular + 1, n, s, q + (paramb.n_max_radial + 1));
       for (int abc = 0; abc < NUM_OF_ABC; ++abc) {
         g_sum_fxyz[(n * NUM_OF_ABC + abc) * N + n1] = s[abc];
       }
@@ -311,16 +315,29 @@ static __global__ void find_descriptor_small_box(
         F,
         Fp);
     } else {
-      apply_ann_one_layer(
-        annmb.dim,
-        annmb.num_neurons1,
-        annmb.w0[t1],
-        annmb.b0[t1],
-        annmb.w1[t1],
-        annmb.b1,
-        q,
-        F,
-        Fp);
+      if (!need_B_projection)
+        apply_ann_one_layer(
+          annmb.dim,
+          annmb.num_neurons1,
+          annmb.w0[t1],
+          annmb.b0[t1],
+          annmb.w1[t1],
+          annmb.b1,
+          q,
+          F,
+          Fp);
+      else
+        apply_ann_one_layer(
+          annmb.dim,
+          annmb.num_neurons1,
+          annmb.w0[t1],
+          annmb.b0[t1],
+          annmb.w1[t1],
+          annmb.b1,
+          q,
+          F,
+          Fp,
+          B_projection + n1 * B_projection_size);
     }
     g_pe[n1] += F;
 
@@ -452,7 +469,8 @@ static __global__ void find_descriptor_small_box(
         accumulate_s(paramb.L_max, d12, r12[0], r12[1], r12[2], gn12, s);
 #endif
       }
-      find_q(paramb.L_max, paramb.num_L, paramb.n_max_angular + 1, n, s, q + (paramb.n_max_radial + 1));
+      find_q(
+        paramb.L_max, paramb.num_L, paramb.n_max_angular + 1, n, s, q + (paramb.n_max_radial + 1));
       for (int abc = 0; abc < NUM_OF_ABC; ++abc) {
         g_sum_fxyz[(n * NUM_OF_ABC + abc) * N + n1] = s[abc];
       }
@@ -665,7 +683,18 @@ static __global__ void find_force_angular_small_box(
           g_gn_angular[index_left_all] * weight_left + g_gn_angular[index_right_all] * weight_right;
         float gnp12 = g_gnp_angular[index_left_all] * weight_left +
                       g_gnp_angular[index_right_all] * weight_right;
-        accumulate_f12(paramb.L_max, paramb.num_L, n, paramb.n_max_angular + 1, d12, r12, gn12, gnp12, Fp, sum_fxyz, f12);
+        accumulate_f12(
+          paramb.L_max,
+          paramb.num_L,
+          n,
+          paramb.n_max_angular + 1,
+          d12,
+          r12,
+          gn12,
+          gnp12,
+          Fp,
+          sum_fxyz,
+          f12);
       }
 #else
       float fc12, fcp12;
@@ -692,7 +721,18 @@ static __global__ void find_force_angular_small_box(
           gn12 += fn12[k] * annmb.c[c_index];
           gnp12 += fnp12[k] * annmb.c[c_index];
         }
-        accumulate_f12(paramb.L_max, paramb.num_L, n, paramb.n_max_angular + 1, d12, r12, gn12, gnp12, Fp, sum_fxyz, f12);
+        accumulate_f12(
+          paramb.L_max,
+          paramb.num_L,
+          n,
+          paramb.n_max_angular + 1,
+          d12,
+          r12,
+          gn12,
+          gnp12,
+          Fp,
+          sum_fxyz,
+          f12);
       }
 #endif
       double s_sxx = 0.0;