Develop

lammps_plugin/ML-UF3/pair_uf3.cpp

@@ -129,25 +129,45 @@ void PairUF3::coeff(int narg, char **arg)
 
 
     if (narg == 3){
-      for (int i = ilo; i <= ihi; i++) {
-        for (int j = MAX(jlo, i); j <= jhi; j++) {
-          if (comm->me == 0)
-            utils::logmesg(lmp, "\nUF3: Opening {} file\n", arg[2]);
-          uf3_read_pot_file(i,j,arg[2]);
+      if (utils::strmatch(arg[0],".*\\*.*") || utils::strmatch(arg[1],".*\\*.*")){
+        for (int i = ilo; i <= ihi; i++) {
+          for (int j = MAX(jlo, i); j <= jhi; j++) {
+            if (comm->me == 0)
+              utils::logmesg(lmp, "\nUF3: Opening {} file\n", arg[2]);
+            uf3_read_pot_file(i,j,arg[2]);
+          }
         }
       }
+
+      else{
+        int i = utils::inumeric(FLERR, arg[0], true, lmp);
+        int j = utils::inumeric(FLERR, arg[1], true, lmp);
+        if (comm->me == 0)
+          utils::logmesg(lmp, "\nUF3: Opening {} file\n", arg[2]);
+        uf3_read_pot_file(i,j,arg[2]);
+      }
     }
 
     if (narg == 4){
-      for (int i = ilo; i <= ihi; i++) {
-        for (int j = jlo; j <= jhi; j++) {
-          for (int k = MAX(klo, jlo); k <= khi; k++) {
-            if (comm->me == 0)
-              utils::logmesg(lmp, "\nUF3: Opening {} file\n", arg[3]);
-            uf3_read_pot_file(i,j,k,arg[3]);
+      if (utils::strmatch(arg[0],".*\\*.*") || utils::strmatch(arg[1],".*\\*.*") || utils::strmatch(arg[2],".*\\*.*")){
+        for (int i = ilo; i <= ihi; i++) {
+          for (int j = jlo; j <= jhi; j++) {
+            for (int k = MAX(klo, jlo); k <= khi; k++) {
+              if (comm->me == 0)
+                utils::logmesg(lmp, "\nUF3: Opening {} file\n", arg[3]);
+              uf3_read_pot_file(i,j,k,arg[3]);
+            }
           }
         }
       }
+      else{
+        if (comm->me == 0)
+          utils::logmesg(lmp, "\nUF3: Opening {} file\n", arg[3]);
+        int i = utils::inumeric(FLERR, arg[0], true, lmp);
+        int j = utils::inumeric(FLERR, arg[1], true, lmp);
+        int k = utils::inumeric(FLERR, arg[2], true, lmp);
+        uf3_read_pot_file(i,j,k,arg[3]);
+      }
     }
   }
 
@@ -281,6 +301,9 @@ void PairUF3::uf3_read_pot_file(int itype, int jtype, char *potf_name)
   utils::logmesg(lmp, "UF3: {} file should contain UF3 potential for {} {}\n", \
           potf_name, itype, jtype);
 
+  if (!platform::file_is_readable(potf_name))
+    error->all(FLERR, "UF3: {} file is not readable", potf_name);
+
   FILE *fp;
   fp = utils::open_potential(potf_name, lmp, nullptr);
 
@@ -400,6 +423,9 @@ void PairUF3::uf3_read_pot_file(int itype, int jtype, int ktype, char *potf_name
   utils::logmesg(lmp, "UF3: {} file should contain UF3 potential for {} {} {}\n",
           potf_name, itype, jtype, ktype);
 
+  if (!platform::file_is_readable(potf_name))
+    error->all(FLERR, "UF3: {} file is not readable", potf_name);
+
   FILE *fp;
   fp = utils::open_potential(potf_name, lmp, nullptr);
 

lammps_plugin/README.rst

@@ -58,9 +58,9 @@ To use UF3 potentials in lammps just add the following tags to the lammps input
 
 The 'uf3' keyword in :code:`pair_style` invokes the UF3 potentials in lammps. The number next to the :code:`uf3` keyword tells lammps whether the user wants to run the MD code with just 2-body or 2 and 3-body UF3 potentials. The last number of this line specifies the number of elemnts in the system. So in the above example, the user wants to run MD simulation with UF3 potentials containing both 2-body and 3-body interactions on a system containing only 1 element.
 
-The :code:`pair_coeff` tag is used to read in the user-provided UF3 lammps potential files. These files can be generated directly from the :code:`json` potential files of UF3. We recommend using the :code:`generate_uf3_lammps_pots.py` script (`found here <https://github.com/monk-04/uf3/tree/lammps_implementation/lammps_plugin/scripts>`_) for generating the UF3 lammps potential files. It will also additionally print lines that should be added to the lammps input file for using UF3 lammps potential files.
+The :code:`pair_coeff` tag is used to read in the user-provided UF3 lammps potential files. These files can be generated directly from the :code:`json` potential files of UF3. We recommend using the :code:`generate_uf3_lammps_pots.py` script (`found here </lammps_plugin/scripts>`_) for generating the UF3 lammps potential files. It will also additionally print lines that should be added to the lammps input file for using UF3 lammps potential files. **Note, nothing is inferred from the name of the UF3 lammps potential file. The name of the files can be completely arbitrary**
 
-After :code:`pair_coeff` specify the interactions (two numbers for 2-body, three numbers for 3-body) followed by the name of the potential file. The user can also use asterisks:code:`*` for wild-card characters. In this case the behaviour is similar to other LAMMPS :code:`pair_style` for example LJ. The user can also specify. Make sure these files are present in the current run directory or in directories where lammps can find them.
+After :code:`pair_coeff` specify the interactions (two numbers for 2-body, three numbers for 3-body) followed by the name of the potential file. The user can also use :code:`*` for wild-card characters. In this case the behaviour is similar to other LAMMPS :code:`pair_style` for example LJ. Make sure these files are present in the current run directory or in directories where lammps can find them.
 
 As an example for a multicomponet system containing elements 'A' and 'B' the above lines can be-
 
@@ -77,6 +77,26 @@ As an example for a multicomponet system containing elements 'A' and 'B' the abo
    pair_coeff 2 1 2 B_A_B
    pair_coeff 2 2 2 B_B_B
 
+
+If potential file is specified for 2-1 interaction (:code:`pair_coeff 2 1 FileName`), the potential for 1-2 interaction is automatically mapped. So, the following lines are also valid-
+
+.. code:: bash
+
+   pair_style uf3 3 2
+
+   pair_coeff 2 2 A_A
+   pair_coeff 1 1 B_B
+   pair_coeff 2 1 A_B
+
+   pair_coeff 2 2 2 A_A_A
+   pair_coeff 2 2 1 A_A_B
+   pair_coeff 2 1 1 A_B_B
+
+   pair_coeff 1 1 1 B_B_B
+   pair_coeff 1 2 1 B_A_B
+   pair_coeff 1 2 2 B_A_A
+
+
 Following format is also a valid for system containing elements 'A' and 'B'
 
 .. code:: bash

uf3/representation/angles.py

@@ -416,58 +416,44 @@ def generate_triplets(i_where: np.ndarray,
     i_groups = np.array_split(j_where, np.cumsum(group_sizes)[:-1])
     # generate j-k combinations
     for i in range(len(i_groups)):
-        tuples = np.array(np.meshgrid(i_groups[i],
-                                      i_groups[i])).T.reshape(-1, 2)
-        tuples = np.insert(tuples, 0, i_values[i], axis=1)
+        j_arr, k_arr = np.meshgrid(i_groups[i], i_groups[i])
+
+        # Pick out unique neighbor pairs of central atom i
+        # ex: With center atom 0 and its neighbors [2, 1, 3],
+        # j_arr = [[2, 1, 3],
+        #          [2, 1, 3],
+        #          [2, 1, 3]]
+        # k_arr = [[2, 2, 2],
+        #          [1, 1, 1],
+        #          [3, 3, 3]]
+        # j_indices = [1, 2, 1]
+        # k_indices = [2, 3, 3]
+        # => unique pairs: (1, 2), (2, 3), (1, 3)
+        # The unique_pair_mask has filtered out pairs like (1, 1), (2, 1), (3, 2), etc.
+        unique_pair_mask = (j_arr < k_arr)
+        j_indices = j_arr[unique_pair_mask]
+        k_indices = k_arr[unique_pair_mask]
+        tuples = np.vstack((i_values[i] * np.ones(len(j_indices), dtype=int),
+                            j_indices, k_indices)).T  # array of unique triplets
+
         comp_tuples = sup_composition[tuples]
-        
-        sort_indices = np.argsort(comp_tuples[:, 1:],axis=1)
+
+        sort_indices = np.argsort(comp_tuples[:, 1:],axis=1)  # to sort by atomic number
         comp_tuples_slice = np.take_along_axis(comp_tuples[:, 1:],sort_indices,axis=1)
         tuples_slice = np.take_along_axis(tuples[:, 1:],sort_indices,axis=1)
-
+
+        # sort comp_tuples and tuples the same way
         comp_tuples = np.hstack((comp_tuples[:, [0]], comp_tuples_slice))
         tuples = np.hstack((tuples[:, [0]], tuples_slice))
-        
+
         ijk_hash = composition.get_szudzik_hash(comp_tuples)
 
         grouped_triplets = [None] * n_hashes
         for j, hash_ in enumerate(hashes):
             ituples = tuples[ijk_hash == hash_]
-            icomp_tuples = comp_tuples[ijk_hash == hash_]
             if len(ituples) == 0:
                 grouped_triplets[j] = None
                 continue
-
-            # Check if same neighbouring elements
-            if icomp_tuples[0][1] == icomp_tuples[0][2]:
-                # element at j and k are same;
-                # remove redundant interactions
-                comparison_mask = (ituples[:, 1] < ituples[:, 2])
-                ituples = ituples[comparison_mask]
-                icomp_tuples = icomp_tuples[comparison_mask]
-
-            else:
-                # Elements at j and k not same
-                # cordinates of j and k are same; filter eg 011, 022, 033
-                # this comparison mask is redundant
-                comparison_mask = (ituples[:, 1] != ituples[:, 2])
-
-                ituples = ituples[comparison_mask]
-                icomp_tuples = icomp_tuples[comparison_mask]
-
-                if len(ituples) > 0:
-                    # Remove repetitive interactions
-                    ituples = np.unique(ituples,axis=0)
-                    icomp_tuples = np.unique(icomp_tuples,axis=0)
-
-                    # sort by electro-negativity as the interaction tuple is always
-                    # sorted by electro-negativity
-                    en_j = composition.reference_X[ase_symbols.chemical_symbols[icomp_tuples[0][1]]]
-                    en_k = composition.reference_X[ase_symbols.chemical_symbols[icomp_tuples[0][2]]]
-                    if  en_k < en_j:
-                        # Interchange columns 1 and 2
-                        ituples[:, [1, 2]] = ituples[:, [2, 1]]
-
             # extract distance tuples
             r_l = distance_matrix[ituples[:, 0], ituples[:, 1]]
             r_m = distance_matrix[ituples[:, 0], ituples[:, 2]]
@@ -484,11 +470,7 @@ def generate_triplets(i_where: np.ndarray,
             r_m = r_m[dist_mask]
             r_n = r_n[dist_mask]
             ituples = ituples[dist_mask]
-
-            if len(ituples) == 0:
-                grouped_triplets[j] = None
-            else:
-                grouped_triplets[j] = i, r_l, r_m, r_n, ituples
+            grouped_triplets[j] = i, r_l, r_m, r_n, ituples
         yield grouped_triplets