Update data format (#4)

Author: jarrodmcc

examples/generate_diatomic.py

@@ -0,0 +1,86 @@
+# FermiLib plugin to interface with Psi4
+# Copyright (C) 2017 FermiLib Developers
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""This is a simple script for generating data."""
+import os
+
+from fermilib.utils import MolecularData
+
+from fermilibpluginpsi4 import run_psi4
+
+if __name__ == '__main__':
+
+    # Set chemical parameters.
+    element_names = ['H', 'H']
+    basis = 'sto-3g'
+    charge = 0
+    multiplicity = 1
+
+    # Single point at equilibrium for testing
+    spacings = [0.7414]
+
+    # Add points for a full dissociation curve from 0.1 to 3.0 angstroms
+    spacings += [0.1 * r for r in range(1, 31)]
+
+    # Set run options
+    run_scf = 1
+    run_mp2 = 1
+    run_cisd = 1
+    run_ccsd = 1
+    run_fci = 1
+    verbose = 1
+    tolerate_error = 1
+
+    # Run Diatomic Curve
+    for spacing in spacings:
+        description = "{}".format(spacing)
+        geometry = [[element_names[0], [0, 0, 0]],
+                    [element_names[1], [0, 0, spacing]]]
+        molecule = MolecularData(geometry,
+                                 basis,
+                                 multiplicity,
+                                 charge,
+                                 description)
+
+        molecule = run_psi4(molecule,
+                            run_scf=run_scf,
+                            run_mp2=run_mp2,
+                            run_cisd=run_cisd,
+                            run_ccsd=run_ccsd,
+                            run_fci=run_fci,
+                            verbose=verbose,
+                            tolerate_error=tolerate_error)
+        molecule.save()
+
+    # Run Li H single point
+    description = "1.45"
+    geometry = [['Li', [0, 0, 0]],
+                ['H', [0, 0, 1.45]]]
+    molecule = MolecularData(geometry,
+                             basis,
+                             multiplicity,
+                             charge,
+                             description)
+
+    molecule = run_psi4(molecule,
+                        run_scf=run_scf,
+                        run_mp2=run_mp2,
+                        run_cisd=run_cisd,
+                        run_ccsd=run_ccsd,
+                        run_fci=run_fci,
+                        verbose=verbose,
+                        tolerate_error=tolerate_error)
+    molecule.save()

examples/psi4_demo.ipynb

@@ -191,7 +191,7 @@
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 2.0
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
@@ -202,5 +202,5 @@
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
-}
+ "nbformat_minor": 0
+}

fermilibpluginpsi4/_psi4_conversion_functions.py

@@ -209,49 +209,49 @@ def beta_unoccupied(i):
     # Store singles
     for entry in T1IA_Amps:
         i, a, value = entry
-        single_amplitudes[alpha_occupied(i),
-                          alpha_unoccupied(a)] = value
+        single_amplitudes[alpha_unoccupied(a),
+                          alpha_occupied(i)] = value
         if (restricted):
-            single_amplitudes[beta_occupied(i),
-                              beta_unoccupied(a)] = value
+            single_amplitudes[beta_unoccupied(a),
+                              beta_occupied(i)] = value
 
     for entry in T1ia_Amps:
         i, a, value = entry
-        single_amplitudes[beta_occupied(i),
-                          beta_unoccupied(a)] = value
+        single_amplitudes[beta_unoccupied(a),
+                          beta_occupied(i)] = value
 
     # Store doubles, include factor of 1/2 for convention
     for entry in T2IJAB_Amps:
         i, j, a, b, value = entry
-        double_amplitudes[alpha_occupied(i),
-                          alpha_unoccupied(a),
-                          alpha_occupied(j),
-                          alpha_unoccupied(b)] = -value / 2.
+        double_amplitudes[alpha_unoccupied(a),
+                          alpha_occupied(i),
+                          alpha_unoccupied(b),
+                          alpha_occupied(j)] = value / 2.
         if (restricted):
-            double_amplitudes[beta_occupied(i),
-                              beta_unoccupied(a),
-                              beta_occupied(j),
-                              beta_unoccupied(b)] = -value / 2.
+            double_amplitudes[beta_unoccupied(a),
+                              beta_occupied(i),
+                              beta_unoccupied(b),
+                              beta_occupied(j)] = value / 2.
 
     for entry in T2ijab_Amps:
         i, j, a, b, value = entry
-        double_amplitudes[beta_occupied(i),
-                          beta_unoccupied(a),
-                          beta_occupied(j),
-                          beta_unoccupied(b)] = -value / 2.
+        double_amplitudes[beta_unoccupied(a),
+                          beta_occupied(i),
+                          beta_unoccupied(b),
+                          beta_occupied(j)] = value / 2.
 
     for entry in T2IjAb_Amps:
         i, j, a, b, value = entry
-        double_amplitudes[alpha_occupied(i),
-                          alpha_unoccupied(a),
-                          beta_occupied(j),
-                          beta_unoccupied(b)] = -value / 2.
+        double_amplitudes[alpha_unoccupied(a),
+                          alpha_occupied(i),
+                          beta_unoccupied(b),
+                          beta_occupied(j)] = value / 2.
 
         if (restricted):
-            double_amplitudes[beta_occupied(i),
-                              beta_unoccupied(a),
-                              alpha_occupied(j),
-                              alpha_unoccupied(b)] = -value / 2.
+            double_amplitudes[beta_unoccupied(a),
+                              beta_occupied(i),
+                              alpha_unoccupied(b),
+                              alpha_occupied(j)] = value / 2.
 
     # Package into InteractionOperator.
     molecule = InteractionOperator(0.0,

fermilibpluginpsi4/_psi4_template

@@ -6,6 +6,7 @@ from fermilib.config import *
 from fermilib.ops._interaction_tensor import (one_body_basis_change,
                                               two_body_basis_change)
 from fermilib.utils import MolecularData
+from numpy import array
 
 sys.path.append('&THIS_DIRECTORY')
 from _psi4_conversion_functions import *
@@ -22,7 +23,8 @@ molecule = MolecularData(&geometry,
                          '&basis',
                          &multiplicity,
                          &charge,
-                         _description)
+                         _description,
+                          filename='&mol_filename')
 
 # Set molecular geometry and symmetry.
 molecule mol {
@@ -45,11 +47,11 @@ set globals {
     basis &basis
     freeze_core false
     fail_on_maxiter true
-    df_scf_guess true
+    df_scf_guess false
     opdm true
     tpdm true
     soscf false
-    scf_type df
+    scf_type pk
     maxiter 1e6
     num_amps_print 1e6
     r_convergence 1e-6
@@ -93,8 +95,7 @@ if &run_scf:
         two_body_integrals = numpy.einsum('psqr', two_body_integrals)
         two_body_integrals = two_body_basis_change(
             two_body_integrals, molecule.canonical_orbitals)
-        integrals_name = molecule.filename + '_eri'
-        numpy.save(integrals_name, two_body_integrals)
+        molecule.two_body_integrals = two_body_integrals
         molecule.save()
 
 
@@ -157,8 +158,7 @@ if &run_cisd:
 
         # Store 1-RDM in molecule file, 2-RDM separately in other file.
         molecule.cisd_one_rdm = cisd_one_rdm
-        cisd_rdm_name = molecule.filename + '_cisd_rdm'
-        numpy.save(cisd_rdm_name, cisd_two_rdm)
+        molecule.cisd_two_rdm = cisd_two_rdm
         molecule.save()
 
 
@@ -203,8 +203,7 @@ if &run_fci:
 
         # Store 1-RDM in molecule file, 2-RDM separately in other file.
         molecule.fci_one_rdm = fci_one_rdm
-        fci_rdm_name = molecule.filename + '_fci_rdm'
-        numpy.save(fci_rdm_name, fci_two_rdm)
+        molecule.fci_two_rdm = fci_two_rdm
         molecule.save()
 
 

fermilibpluginpsi4/_run_psi4.py

@@ -57,7 +57,8 @@ def generate_psi4_input(molecule,
                         run_fci,
                         verbose,
                         tolerate_error,
-                        memory):
+                        memory,
+                        template_file):
     """This function creates and saves a psi4 input file.
 
     Args:
@@ -70,6 +71,7 @@ def generate_psi4_input(molecule,
         verbose: Boolean whether to print calculation results to screen.
         tolerate_error: Whether to fail or merely warn when Psi4 fails.
         memory: Int giving amount of memory to allocate in MB.
+        template_file(str): Specify the filename of a Psi4 template
 
     Returns:
         input_file: A string giving the name of the saved input file.
@@ -81,7 +83,8 @@ def generate_psi4_input(molecule,
     psi4_directory = os.path.dirname(os.path.realpath(__file__))
 
     # Parse input template.
-    template_file = psi4_directory + '/_psi4_template'
+    if template_file is None:
+        template_file = psi4_directory + '/_psi4_template'
     input_template = []
     with open(template_file, 'r') as stream:
         for line in stream:
@@ -103,6 +106,8 @@ def generate_psi4_input(molecule,
                      for line in input_content]
     input_content = [re.sub('&description', str(molecule.description), line)
                      for line in input_content]
+    input_content = [re.sub('&mol_filename', str(molecule.filename), line)
+                     for line in input_content]
     input_content = [re.sub('&geo_string', geo_string, line)
                      for line in input_content]
 
@@ -158,7 +163,8 @@ def run_psi4(molecule,
              tolerate_error=False,
              delete_input=True,
              delete_output=False,
-             memory=8000):
+             memory=8000,
+             template_file=None):
     """This function runs a Psi4 calculation.
 
     Args:
@@ -173,6 +179,7 @@ def run_psi4(molecule,
         delete_input: Optional boolean to delete psi4 input file.
         delete_output: Optional boolean to delete psi4 output file.
         memory: Optional int giving amount of memory to allocate in MB.
+        template_file(str): Path to Psi4 template file
 
     Returns:
         molecule: The updated MolecularData object.
@@ -189,7 +196,8 @@ def run_psi4(molecule,
                                      run_fci,
                                      verbose,
                                      tolerate_error,
-                                     memory)
+                                     memory,
+                                     template_file)
 
     # Run psi4.
     output_file = molecule.filename + '.out'