psi4 / psi4numpy

Combining Psi4 and Numpy for education and development.
BSD 3-Clause "New" or "Revised" License
330 stars 151 forks source link

Broken example: Fitting Lennard-Jones Parameters from Potential Energy Scan #139

Open KonstantinUshenin opened 5 months ago

KonstantinUshenin commented 5 months ago

I tried to repeat: https://github.com/psi4/psi4numpy/blob/master/Tutorials/01_Psi4NumPy-Basics/1b_molecule.ipynb

I got an error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[39], line 9
      5 mol = psi4.geometry(he_dimer.replace('**R**', 
      6                                      str(d)))
      8 # Compute the Counterpoise-Corrected interaction energy
----> 9 en = psi4.energy('MP2/aug-cc-pVDZ', 
     10                  molecule=mol, 
     11                  bsse_type='cp')
     13 # Place in a reasonable unit, Wavenumbers in this case
     14 en *= 219474.6

File ~/anaconda3/envs/p4env/lib/python3.10/site-packages/psi4/driver/driver.py:464, in energy(name, **kwargs)
    461 elif not isinstance(plan, AtomicComputer):
    462     # Advanced "Computer" active
    463     plan.compute()
--> 464     return plan.get_psi_results(return_wfn=return_wfn)
    466 else:
    467     # We have unpacked to an AtomicInput
    468     lowername = plan.method

File ~/anaconda3/envs/p4env/lib/python3.10/site-packages/psi4/driver/driver_nbody.py:1451, in ManyBodyComputer.get_psi_results(self, client, return_wfn)
   1424 def get_psi_results(
   1425     self,
   1426     client: Optional["qcportal.FractalClient"] = None,
   1427     *,
   1428     return_wfn: bool = False) -> EnergyGradientHessianWfnReturn:
   1429     """Called by driver to assemble results into ManyBody-flavored QCSchema,
   1430     then reshape and return them in the customary Psi4 driver interface: ``(e/g/h, wfn)``.
   1431 
   (...)
   1449 
   1450     """
-> 1451     nbody_model = self.get_results(client=client)
   1452     ret = nbody_model.return_result
   1454     wfn = core.Wavefunction.build(self.molecule, "def2-svp", quiet=True)

File ~/anaconda3/envs/p4env/lib/python3.10/site-packages/psi4/driver/driver_nbody.py:1368, in ManyBodyComputer.get_results(self, client)
   1365 core.print_out(info)
   1366 logger.info(info)
-> 1368 results = self.prepare_results(client=client)
   1369 ret_energy = results.pop("ret_energy")
   1370 ret_ptype = results.pop("ret_ptype")

File ~/anaconda3/envs/p4env/lib/python3.10/site-packages/psi4/driver/driver_nbody.py:1304, in ManyBodyComputer.prepare_results(self, results, client)
   1293 metadata = {
   1294     "quiet": self.quiet,
   1295     "nbodies_per_mc_level": self.nbodies_per_mc_level,
   (...)
   1301     "max_nbody": self.max_nbody,
   1302 }
   1303 if self.driver.name == "energy":
-> 1304     nbody_results = assemble_nbody_components("energy", trove["energy"], metadata.copy())
   1306 elif self.driver.name == "gradient":
   1307     nbody_results = assemble_nbody_components("energy", trove["energy"], metadata.copy())

File ~/anaconda3/envs/p4env/lib/python3.10/site-packages/psi4/driver/driver_nbody.py:648, in assemble_nbody_components(ptype, component_results, metadata)
    645         nocp_compute_list[len(w[0])].add(w)
    647 for nb in range(1, nbodies[-1] + 1):
--> 648     cp_by_level[nb] = _sum_cluster_ptype_data(
    649         ptype,
    650         component_results,
    651         cp_compute_list[nb],
    652         fragment_slice_dict,
    653         fragment_size_dict,
    654         mc_level_lbl=mc_level_lbl,
    655     )
    656     nocp_by_level[nb] = _sum_cluster_ptype_data(
    657         ptype,
    658         component_results,
   (...)
    662         mc_level_lbl=mc_level_lbl,
    663     )
    664     if nb in compute_dict["vmfc_levels"]:

File ~/anaconda3/envs/p4env/lib/python3.10/site-packages/psi4/driver/driver_nbody.py:307, in _sum_cluster_ptype_data(ptype, ptype_dict, compute_list, fragment_slice_dict, fragment_size_dict, mc_level_lbl, vmfc, nb)
    304         if vmfc:
    305             sign = ((-1)**(nb - len(frag)))
--> 307         ret += sign * ene
    309     return ret
    311 elif ptype == 'gradient':

TypeError: unsupported operand type(s) for *: 'int' and 'NoneType'
loriab commented 4 months ago

Thanks for the report. This looks like it should be an simple syntax fix, but I'm not able to reproduce the error on a recent v1.9 psi4. That is, the below runs to completion. Is this the same portion of the code that you hit the error with? Or had you added anything more to the script? Thanks!

import psi4

he_dimer = """
He
--
He 1 **R**
"""

distances = [2.875, 3.0, 3.125, 3.25, 3.375, 3.5, 3.75, 4.0, 4.5, 5.0, 6.0, 7.0]
energies = []
for d in distances:
    # Build a new molecule at each separation
    mol = psi4.geometry(he_dimer.replace('**R**', str(d)))

    # Compute the Counterpoise-Corrected interaction energy
    en = psi4.energy('MP2/aug-cc-pVDZ', molecule=mol, bsse_type='cp')

    # Place in a reasonable unit, Wavenumbers in this case
    en *= 219474.6

    # Append the value to our list
    energies.append(en)

print("Finished computing the potential!")