Develop

qml/aglaia/aglaia.py

@@ -594,7 +594,7 @@ def _set_acsf_parameters(self, params):
         """
 
         self.acsf_parameters = {'rcut': 5.0, 'acut': 5.0, 'nRs2': 5, 'nRs3': 5, 'nTs': 5,
-                                      'zeta': 220.127, 'eta': 30.8065}
+                                      'zeta': 220.127, 'eta': 30.8065, 'bin_min': 0.8}
 
         if params is not None:
             for key, value in params.items():
@@ -842,6 +842,10 @@ def _check_acsf_values(self):
         if is_numeric_array(self.acsf_parameters['zeta']):
             raise InputError("Expecting a scalar value for zeta. Got %s." % (self.acsf_parameters['zeta']))
 
+        if not is_positive_or_zero(self.acsf_parameters['bin_min']):
+            raise InputError(
+                "Expected positive or zero float for variable 'bin_min'. Got %s." % str(self.acsf_parameters['bin_min']))
+
     def _get_msize(self, pad = 0):
         """
         Gets the maximum number of atoms in a single molecule. To support larger molecules
@@ -1695,7 +1699,8 @@ def _generate_acsf_tf(self, xyz, classes):
                                           nRs2=self.acsf_parameters['nRs2'],
                                           nRs3=self.acsf_parameters['nRs3'],
                                           nTs=self.acsf_parameters['nTs'], eta=self.acsf_parameters['eta'],
-                                          zeta=self.acsf_parameters['zeta'])
+                                          zeta=self.acsf_parameters['zeta'],
+                                          bin_min=self.acsf_parameters['bin_min'])
 
         sess = tf.Session()
         sess.run(tf.global_variables_initializer())
@@ -1765,7 +1770,8 @@ def _generate_acsf_fortran(self, xyz, classes):
                                   nTs=self.acsf_parameters['nTs'],
                                   eta2=self.acsf_parameters['eta'],
                                   eta3=self.acsf_parameters['eta'],
-                                  zeta=self.acsf_parameters['zeta'])
+                                  zeta=self.acsf_parameters['zeta'],
+                                  bin_min=self.acsf_parameters['bin_min'])
 
                 padded_g = np.zeros((initial_natoms, g.shape[-1]))
                 padded_g[:g.shape[0], :] = g
@@ -2121,7 +2127,7 @@ def _check_representation_parameters(self, parameters):
         elif self.representation_name == "acsf":
 
             acsf_parameters =  {'rcut': 5.0, 'acut': 5.0, 'nRs2': 5, 'nRs3': 5, 'nTs': 5,
-                                      'zeta': 220.127, 'eta': 30.8065}
+                                      'zeta': 220.127, 'eta': 30.8065, 'bin_min':0.8}
 
             for key, value in parameters.items():
                 try:
@@ -2423,7 +2429,8 @@ def _build_model_from_xyz(self, n_atoms, element_weights, element_biases):
                                                     nRs3=self.acsf_parameters['nRs3'],
                                                     nTs=self.acsf_parameters['nTs'],
                                                     eta=self.acsf_parameters['eta'],
-                                                    zeta=self.acsf_parameters['zeta'])
+                                                    zeta=self.acsf_parameters['zeta'],
+                                                    bin_min=self.acsf_parameters['bin_min'])
 
         with tf.name_scope("Model_pred"):
             batch_energies_nn = self._model(batch_representation, batch_zs, element_weights, element_biases)

qml/aglaia/np_symm_funct.py

@@ -188,7 +188,7 @@ def acsf_ang(xyzs, Zs, element_pairs, angular_cutoff, angular_rs, theta_s, zeta,
     return np.asarray(total_descriptor)
 
 def generate_acsf_np(xyzs, Zs, elements, element_pairs, rcut, acut, nRs2,
-                     nRs3, nTs, zeta, eta):
+                     nRs3, nTs, zeta, eta, bin_min):
     """
     This function calculates the symmetry functions used in the tensormol paper.
 
@@ -203,11 +203,12 @@ def generate_acsf_np(xyzs, Zs, elements, element_pairs, rcut, acut, nRs2,
     :param theta_s: list of all the thetas parameters. Numpy array of shape (n_thetas,)
     :param zeta: parameter. scalar.
     :param eta: parameter. scalar.
+    :param bin_min: value at which to start the binning of the distances
     :return: numpy array of shape (n_samples, n_atoms, n_rad_rs*n_elements + n_ang_rs*n_thetas*n_element_pairs)
     """
 
-    radial_rs = np.linspace(0, rcut, nRs2)
-    angular_rs = np.linspace(0, acut, nRs3)
+    radial_rs = np.linspace(bin_min, rcut, nRs2)
+    angular_rs = np.linspace(bin_min, acut, nRs3)
     theta_s = np.linspace(0, np.pi, nTs)
 
     rad_term = acsf_rad(xyzs, Zs, elements, rcut, radial_rs, eta)

qml/aglaia/symm_funct.py

@@ -381,7 +381,7 @@ def sum_ang(pre_sumterm, Zs, element_pairs_list, angular_rs, theta_s):
     return clean_final_term
 
 def generate_acsf_tf(xyzs, Zs, elements, element_pairs, rcut, acut,
-                     nRs2, nRs3, nTs, zeta, eta):
+                     nRs2, nRs3, nTs, zeta, eta, bin_min):
     """
     This function generates the atom centred symmetry function as used in the Tensormol paper. Currently only tested for
     single systems with many conformations. It requires the coordinates of all the atoms in each data sample, the atomic
@@ -410,13 +410,15 @@ def generate_acsf_tf(xyzs, Zs, elements, element_pairs, rcut, acut,
     :type zeta: scalar float
     :param eta: parameter in the exponential terms
     :type eta: scalar float
+    :param bin_min: the value at which to start binning the distances
+    :type bin_min: positive float
 
     :return: the atom centred symmetry functions
     :rtype: a tf tensor of shape a tf tensor of shape (n_samples, n_atoms, nRs2 * n_elements + nRs3 * nTs * n_elementpairs)
     """
 
-    radial_rs = np.linspace(0, rcut, nRs2)
-    angular_rs = np.linspace(0, acut, nRs3)
+    radial_rs = np.linspace(bin_min, rcut, nRs2)
+    angular_rs = np.linspace(bin_min, acut, nRs3)
     theta_s = np.linspace(0, np.pi, nTs)
 
     with tf.name_scope("acsf_params"):

qml/data/compound.py

@@ -295,7 +295,7 @@ def generate_slatm(self, mbtypes,
         self.representation = slatm
 
     def generate_acsf(self, elements = [1,6,7,8,16], nRs2 = 3, nRs3 = 3, nTs = 3, eta2 = 1,
-            eta3 = 1, zeta = 1, rcut = 5, acut = 5, gradients = False):
+                      eta3 = 1, zeta = 1, rcut = 5, acut = 5, bin_min=0.8, gradients = False):
         """
         Generate the variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J
 
@@ -317,12 +317,14 @@ def generate_acsf(self, elements = [1,6,7,8,16], nRs2 = 3, nRs3 = 3, nTs = 3, et
         :type rcut: float
         :param acut: Cut-off radius of the three-body terms
         :type acut: float
+        :param bin_min: the value at which to start binning the distances
+        :type bin_min: positive float
         :param gradients: To return gradients or not
         :type gradients: boolean
         """
 
-        Rs2 = np.linspace(0, rcut, nRs2)
-        Rs3 = np.linspace(0, acut, nRs3)
+        Rs2 = np.linspace(bin_min, rcut, nRs2)
+        Rs3 = np.linspace(bin_min, acut, nRs3)
         Ts = np.linspace(0, np.pi, nTs)
         n_elements = len(elements)
         natoms = len(self.coordinates)

qml/representations/representations.py

@@ -549,8 +549,8 @@ def generate_slatm(coordinates, nuclear_charges, mbtypes,
 
     return mbs
 
-def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 = 3, nRs3 = 3, nTs = 3, eta2 = 1, 
-        eta3 = 1, zeta = 1, rcut = 5, acut = 5, gradients = False):
+def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 = 3, nRs3 = 3, nTs = 3, eta2 = 1,
+                  eta3 = 1, zeta = 1, rcut = 5, acut = 5, bin_min=0.8, gradients = False):
     """
     Generate the variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J
 
@@ -576,14 +576,16 @@ def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 =
     :type rcut: float
     :param acut: Cut-off radius of the three-body terms
     :type acut: float
+    :param bin_min: the value at which to start binning the distances
+    :type bin_min: positive float
     :param gradients: To return gradients or not
     :type gradients: boolean
     :return: Atom-centered symmetry functions representation
     :rtype: numpy array
     """
 
-    Rs2 = np.linspace(0, rcut, nRs2)
-    Rs3 = np.linspace(0, acut, nRs3)
+    Rs2 = np.linspace(bin_min, rcut, nRs2)
+    Rs3 = np.linspace(bin_min, acut, nRs3)
     Ts = np.linspace(0, np.pi, nTs)
     n_elements = len(elements)
     natoms = len(coordinates)

test/test_acsf.py

@@ -49,6 +49,7 @@ def test_acsf_1():
     acut = 5
     zeta = 220.127
     eta = 30.8065
+    bin_min = 0.0
 
     input_data = test_dir + "/data/data_test_acsf.npz"
     data = np.load(input_data)
@@ -65,13 +66,13 @@ def test_acsf_1():
         zs_tf = tf.placeholder(shape=[n_samples, n_atoms], dtype=tf.int32, name="zs")
         xyz_tf = tf.placeholder(shape=[n_samples, n_atoms, 3], dtype=tf.float32, name="xyz")
 
-    acsf_tf_t = symm_funct.generate_acsf_tf(xyz_tf, zs_tf, elements, element_pairs, rcut, acut, nRs2, nRs3, nTs, zeta, eta)
+    acsf_tf_t = symm_funct.generate_acsf_tf(xyz_tf, zs_tf, elements, element_pairs, rcut, acut, nRs2, nRs3, nTs, zeta, eta, bin_min)
 
     sess = tf.Session()
     sess.run(tf.global_variables_initializer())
     acsf_tf = sess.run(acsf_tf_t, feed_dict={xyz_tf: xyzs, zs_tf: zs})
 
-    acsf_np = np_symm_funct.generate_acsf_np(xyzs, zs, elements, element_pairs, rcut, acut, nRs2, nRs3, nTs, zeta, eta)
+    acsf_np = np_symm_funct.generate_acsf_np(xyzs, zs, elements, element_pairs, rcut, acut, nRs2, nRs3, nTs, zeta, eta, bin_min)
 
     n_samples = xyzs.shape[0]
     n_atoms = xyzs.shape[1]
@@ -97,6 +98,7 @@ def test_acsf_2():
     acut = 5
     zeta = 220.127
     eta = 30.8065
+    bin_min = 0.0
 
     input_data = test_dir + "/data/qm7_testdata.npz"
     data = np.load(input_data)
@@ -113,13 +115,13 @@ def test_acsf_2():
         zs_tf = tf.placeholder(shape=[n_samples, max_n_atoms], dtype=tf.int32, name="zs")
         xyz_tf = tf.placeholder(shape=[n_samples, max_n_atoms, 3], dtype=tf.float32, name="xyz")
 
-    acsf_tf_t = symm_funct.generate_acsf_tf(xyz_tf, zs_tf, elements, element_pairs, rcut, acut, nRs2, nRs3, nTs, zeta, eta)
+    acsf_tf_t = symm_funct.generate_acsf_tf(xyz_tf, zs_tf, elements, element_pairs, rcut, acut, nRs2, nRs3, nTs, zeta, eta, bin_min)
 
     sess = tf.Session()
     sess.run(tf.global_variables_initializer())
     acsf_tf = sess.run(acsf_tf_t, feed_dict={xyz_tf: xyzs, zs_tf: zs})
 
-    acsf_np = np_symm_funct.generate_acsf_np(xyzs, zs, elements, element_pairs, rcut, acut, nRs2, nRs3, nTs, zeta, eta)
+    acsf_np = np_symm_funct.generate_acsf_np(xyzs, zs, elements, element_pairs, rcut, acut, nRs2, nRs3, nTs, zeta, eta, bin_min)
 
     for i in range(n_samples):
         for j in range(max_n_atoms):
@@ -130,7 +132,6 @@ def test_acsf_2():
                 acsf_tf_sort = np.sort(acsf_tf[i][j])
                 np.testing.assert_array_almost_equal(acsf_np_sort, acsf_tf_sort, decimal=4)
 
-
 if __name__ == "__main__":
     test_acsf_1()
     test_acsf_2()

test/test_symm_funct.py

@@ -63,7 +63,7 @@ def test_acsf():
              "qm7/0110.xyz"]
 
 
-    path = test_dir = os.path.dirname(os.path.realpath(__file__))
+    path = os.path.dirname(os.path.realpath(__file__))
 
     mols = []
     for xyz_file in files:
@@ -85,7 +85,7 @@ def fort_acsf(mols, path, elements):
     # Generate atom centered symmetry functions representation
     # using the Compound class
     for i, mol in enumerate(mols):
-        mol.generate_acsf(elements = elements)
+        mol.generate_acsf(elements = elements, bin_min=0.0)
 
     X_test = np.concatenate([mol.representation for mol in mols])
     X_ref = np.loadtxt(path + "/data/acsf_representation.txt")
@@ -96,8 +96,8 @@ def fort_acsf(mols, path, elements):
     rep = []
     for i, mol in enumerate(mols):
         rep.append(generate_acsf(coordinates = mol.coordinates,
-                nuclear_charges = mol.nuclear_charges, 
-                elements = elements))
+                                 nuclear_charges = mol.nuclear_charges,
+                                 elements = elements, bin_min=0.0))
 
     X_test = np.concatenate(rep)
     X_ref = np.loadtxt(path + "/data/acsf_representation.txt")
@@ -111,6 +111,7 @@ def tf_acsf(mols, path, elements):
     n_theta_s = 3
     zeta = 1.0
     eta = 1.0
+    bin_min=0.0
 
     element_pairs = []
     for i, ei in enumerate(elements):
@@ -128,7 +129,7 @@ def tf_acsf(mols, path, elements):
         zs_tf = tf.placeholder(shape=[n_samples, max_n_atoms], dtype=tf.int32, name="zs")
         xyz_tf = tf.placeholder(shape=[n_samples, max_n_atoms, 3], dtype=tf.float32, name="xyz")
 
-    acsf_tf_t = symm_funct.generate_acsf_tf(xyz_tf, zs_tf, elements, element_pairs, radial_cutoff, angular_cutoff, n_radial_rs, n_angular_rs, n_theta_s, zeta, eta)
+    acsf_tf_t = symm_funct.generate_acsf_tf(xyz_tf, zs_tf, elements, element_pairs, radial_cutoff, angular_cutoff, n_radial_rs, n_angular_rs, n_theta_s, zeta, eta, bin_min)
 
     sess = tf.Session()
     sess.run(tf.global_variables_initializer())
@@ -142,7 +143,7 @@ def fort_acsf_gradients(mols, path, elements):
     # Generate atom centered symmetry functions representation
     # and gradients using the Compound class
     for i, mol in enumerate(mols):
-        mol.generate_acsf(elements = elements, gradients = True)
+        mol.generate_acsf(elements = elements, gradients = True, bin_min=0.0)
 
     X_test = np.concatenate([mol.representation for mol in mols])
     X_ref = np.loadtxt(path + "/data/acsf_representation.txt")
@@ -159,7 +160,7 @@ def fort_acsf_gradients(mols, path, elements):
     grad = []
     for i, mol in enumerate(mols):
         r, g = generate_acsf(coordinates = mol.coordinates, nuclear_charges = mol.nuclear_charges,
-                elements = elements, gradients = True)
+                             elements = elements, gradients = True, bin_min=0.0)
         rep.append(r)
         grad.append(g)