Repository: quantumlib/OpenFermion
Pull request #544 (2019-11-29 06:43:55)

Add functionality to retrieve specific types of structure info from PubChem
Currently, the function **geometry_from_pubchem** in **utils/_pubchem.py/** constructs the geometry of a molecule based on 3D structure info if available, otherwise the 2D structure info is used. If both 3D and 2D info are not available in PubChem, None is returned.

In this pull request I have implemented a new argument for this functiion: structure. With this argument users of the function will be able to request _only 3D_ or _only 2D_ information. The default behavior of the function is unchanged, i.e. if the 'structure' argument is not specified, the function will behave as before this pull request, for backwards compatibility.

The reason for the implementation of this argument is, that for a bunch of molecules, the 2D structure information has nothing to do with the actual geometric sizes of the molecule. And, when 3D is not available (e.g for H2), the default behavior will return incorrect geometry based on the non-accurate 2D structure info. It is sometimes better to return nothing, than return incorrect data. With this new option users will be able to force the function to use e.g. only 3D structure info. 

PubChem support confirms, that the 2D structure info is not guaranteed to reflect actual information about molecules. Here is a snipped from their reply:

> PubChem 2-D structures have coordinates generated by a layout algorithm.  They do not reflect anything other than that to help provide a depiction that is a suitable distance from each other.  Yes, they could be normalized in some way (e.g., 0 and 1 or -0.5 and 0.5) but they are arbitrary, in a sense.
[Re_ [Pubchem-help] 2D coordinates of compounds_a.pdf](https://github.com/quantumlib/OpenFermion/files/3904024/Re_.Pubchem-help.2D.coordinates.of.compounds_a.pdf)


My full email conversation with PubChem support is attached!

googlebot: 
Thanks for your pull request. It looks like this may be your first contribution to a Google open source project (if not, look below for help). Before we can look at your pull request, you'll need to sign a Contributor License Agreement (CLA).

:memo: **Please visit <https://cla.developers.google.com/> to sign.**

Once you've signed (or fixed any issues), please reply here with `@googlebot I signed it!` and we'll verify it.

----

#### What to do if you already signed the CLA

##### Individual signers

*   It's possible we don't have your GitHub username or you're using a different email address on your commit. Check [your existing CLA data](https://cla.developers.google.com/clas) and verify that your [email is set on your git commits](https://help.github.com/articles/setting-your-email-in-git/).

##### Corporate signers

*   Your company has a Point of Contact who decides which employees are authorized to participate. Ask your POC to be added to the group of authorized contributors. If you don't know who your Point of Contact is, direct the Google project maintainer to [go/cla#troubleshoot](http://go/cla#troubleshoot) ([Public version](https://opensource.google/docs/cla/#troubleshoot)).
*   The email used to register you as an authorized contributor must be the email used for the Git commit. Check [your existing CLA data](https://cla.developers.google.com/clas) and verify that your [email is set on your git commits](https://help.github.com/articles/setting-your-email-in-git/).
*   The email used to register you as an authorized contributor must also be [attached to your GitHub account](https://github.com/settings/emails).
		

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F544) for more info**.

<!-- need_sender_cla -->
hay-k: @googlebot I signed it!
googlebot: CLAs look good, thanks!

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F544) for more info**.

<!-- ok -->
ncrubin: @hay-k Thanks for the submission and clarifying this!  This is really helpful to get cleared up and make sure people are using pubchem in a way that is reasonable.  I've only requested minor changes to make.  
hay-k: Hi @ncrubin, I have commited the changes.

Patch:
diff --git a/src/openfermion/utils/_pubchem.py b/src/openfermion/utils/_pubchem.py
index f1270e0..6ede28e 100644
--- a/src/openfermion/utils/_pubchem.py
+++ b/src/openfermion/utils/_pubchem.py
@@ -11,13 +11,20 @@
 #   limitations under the License.
 
 
-def geometry_from_pubchem(name):
+def geometry_from_pubchem(name, structure=None):
     """Function to extract geometry using the molecule's name from the PubChem
-    database.
+    database. The 'structure' argument can be used to specify which structure
+    info to use to extract the geometry. If structure=None, the geometry will
+    be constructed based on 3D info, if available, otherwise on 2D (to keep
+    backwards compatibility with the times when the argument 'structure'
+    was not implemented).
 
     Args:
         name: a string giving the molecule's name as required by the PubChem
             database.
+        structure: a string '2d' or '3d', to specify a specific structure
+            information to be retrieved from pubchem. The default is None.
+            Recommended value is '3d'.
 
     Returns:
         geometry: a list of tuples giving the coordinates of each atom with
@@ -25,29 +32,29 @@ def geometry_from_pubchem(name):
     """
     import pubchempy
 
-    pubchempy_2d_molecule = pubchempy.get_compounds(name, 'name',
-                                                    record_type='2d')
-
-    # Check if 2-D geometry is available. If not then no geometry is.
-    if not pubchempy_2d_molecule:
-        print('Unable to find molecule in the PubChem database.')
+    if structure in ['2d', '3d']:
+        pubchempy_molecule = pubchempy.get_compounds(name, 'name',
+                                                     record_type=structure)
+    elif structure is None:
+        # Ideally get the 3-D geometry if available.
+        pubchempy_molecule = pubchempy.get_compounds(name, 'name',
+                                                        record_type='3d')
+
+        # If the 3-D geometry isn't available, get the 2-D geometry instead.
+        if not pubchempy_molecule:
+            pubchempy_molecule = pubchempy.get_compounds(name, 'name',
+                                                        record_type='2d')
+    else:
+        raise ValueError('Incorrect value for the argument structure=%s' % structure)
+
+    # Check if pubchempy_molecule is an empty list or None
+    if not pubchempy_molecule:
+        print('Unable to find structure info in the PubChem database for the specified molecule "%s".' % name)
         return None
 
-    # Ideally get the 3-D geometry if available.
-    pubchempy_3d_molecule = pubchempy.get_compounds(name, 'name',
-                                                    record_type='3d')
-
-    if pubchempy_3d_molecule:
-        pubchempy_geometry = \
-            pubchempy_3d_molecule[0].to_dict(properties=['atoms'])['atoms']
-        geometry = [(atom['element'], (atom['x'], atom['y'], atom['z']))
-                    for atom in pubchempy_geometry]
-        return geometry
-
-    # If the 3-D geometry isn't available, get the 2-D geometry instead.
     pubchempy_geometry = \
-        pubchempy_2d_molecule[0].to_dict(properties=['atoms'])['atoms']
-    geometry = [(atom['element'], (atom['x'], atom['y'], 0))
+        pubchempy_molecule[0].to_dict(properties=['atoms'])['atoms']
+    geometry = [(atom['element'], (atom['x'], atom['y'], atom.get('z', 0)))
                 for atom in pubchempy_geometry]
 
     return geometry


Test patch:
diff --git a/src/openfermion/utils/_pubchem_test.py b/src/openfermion/utils/_pubchem_test.py
index 5e42173..76b3d96 100644
--- a/src/openfermion/utils/_pubchem_test.py
+++ b/src/openfermion/utils/_pubchem_test.py
@@ -76,3 +76,16 @@ class OpenFermionPubChemTest(unittest.TestCase):
         none_geometry = geometry_from_pubchem('none')
 
         self.assertIsNone(none_geometry)
+
+    def test_water_2d(self):
+        water_geometry = geometry_from_pubchem('water', structure='2d')
+        self.water_natoms = len(water_geometry)
+
+        water_natoms = 3
+        self.assertEqual(water_natoms, self.water_natoms)
+
+        self.oxygen_z_1 = water_geometry[0][1][2]
+        self.oxygen_z_2 = water_geometry[1][1][2]
+        z = 0
+        self.assertEqual(z, self.oxygen_z_1)
+        self.assertEqual(z, self.oxygen_z_2)