Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Floating point division by zero when running SAPT with H+ as a compound #3156

Open
dspoel opened this issue Apr 19, 2024 · 3 comments
Open

Comments

@dspoel
Copy link

dspoel commented Apr 19, 2024

When using H+ as the compound in a SAPT calculation it crashes with a divide by zero error:

  File "/var/spool/slurm/d/job333055/slurm_script", line 21, in <module>
    myener = psi4.energy("sapt2+(ccd)dmp2")
  File "/home/spoel/miniconda3/lib/python3.8/site-packages/psi4/driver/driver.py", line 525, in energy
    wfn = procedures['energy'][lowername](lowername, molecule=molecule, **kwargs)
  File "/home/spoel/miniconda3/lib/python3.8/site-packages/psi4/driver/procrouting/proc.py", line 4581, in run_sapt
    monomerA_wfn = scf_helper('RHF', molecule=monomerA, **kwargs)
  File "/home/spoel/miniconda3/lib/python3.8/site-packages/psi4/driver/procrouting/proc.py", line 1887, in scf_helper
    e_scf = scf_wfn.compute_energy()
  File "/home/spoel/miniconda3/lib/python3.8/site-packages/psi4/driver/procrouting/scf_proc/scf_iterator.py", line 99, in scf_compute_energy
    scf_energy = self.finalize_energy()
  File "/home/spoel/miniconda3/lib/python3.8/site-packages/psi4/driver/procrouting/scf_proc/scf_iterator.py", line 692, in scf_finalize_energy
    self.print_energies()
  File "/home/spoel/miniconda3/lib/python3.8/site-packages/psi4/driver/procrouting/scf_proc/scf_iterator.py", line 804, in scf_print_energies
    self.set_variable("HF VIRIAL RATIO", - potential / ke)  # P::e SCF
ZeroDivisionError: float division by zero

The script to reproduce this is below.

#!/usr/bin/env python3
import os
import numpy as np
import psi4 as psi4
psi4.core.set_num_threads(4)
psi4.set_options({"cachelevel": 1, "print": 1})
psi4.set_memory(12000000000)
psi4_io = psi4.core.IOManager.shared_object()
psi4.core.set_output_file('hcl-sapt.log', False)
geometry= """
 1 1
 H -0.0 0.0 -1.36
 --
 -1 1
 Cl 0.0 0.0 1
"""
geom = psi4.geometry(geometry)
psi4.basis_helper("""
assign aug-cc-pvtz
""")
myener = psi4.energy("sapt2+(ccd)dmp2")
for ener in [ 'SAPT ELST ENERGY', 'SAPT EXCH ENERGY', 'SAPT IND ENERGY', 'SAPT DISP ENERGY', 'SAPT TOTAL ENERGY' ]:
    print("%s %g" % ( ener, psi4.variable(ener) ))
@JonathonMisiewicz
Copy link
Contributor

Ah yes, you can't have an electron kinetic energy if you don't have electrons...

@dspoel
Copy link
Author

dspoel commented Apr 19, 2024

That makes sense, but does that mean this cannot be fixed with an if statement? This is only printing after all.

@JonathonMisiewicz
Copy link
Contributor

No, this isn't only printing. This is a Psi variable that can't be defined due to the exceptional circumstance of having zero electrons. I defer to @loriab for deciding what to do about this one.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants