Fixing errors in the calculations
Contents
7.1. Fixing errors in the calculations#
import psi4
import py3Dmol
import pandas as pd
import sys
sys.path.append("..")
from helpers import *
psi4.set_num_threads(2)
psi4.set_memory('2 GB')
Threads set to 2 by Python driver.
Memory set to 1.863 GiB by Python driver.
2000000000
7.1.1. Allyl radical#
Exercise 1
Fix the error in the single-point wavefunction optimization for the allyl radical using the BLYP exchange-correlation functional.
allyl_radical = psi4.geometry("""
symmetry c1
1 2
C 1.307894 -0.292412 0.050945
C 0.000000 0.443183 0.200966
C -1.307894 -0.292412 0.050945
H 2.171777 0.390984 0.102354
H 1.440837 -1.062338 0.844942
H 1.356754 -0.838161 -0.914887
H 0.000000 1.488670 0.522326
H -2.171777 0.390984 0.102354
H -1.356754 -0.838161 -0.914887
H -1.440837 -1.062338 0.844942""")
c: [1.0]
fc: [1.0]
m: [2]
fm: [2]
---------------------------------------------------------------------------
ValidationError Traceback (most recent call last)
Cell In[3], line 1
----> 1 allyl_radical = psi4.geometry("""
2 symmetry c1
3 1 2
4 C 1.307894 -0.292412 0.050945
5 C 0.000000 0.443183 0.200966
6 C -1.307894 -0.292412 0.050945
7 H 2.171777 0.390984 0.102354
8 H 1.440837 -1.062338 0.844942
9 H 1.356754 -0.838161 -0.914887
10 H 0.000000 1.488670 0.522326
11 H -2.171777 0.390984 0.102354
12 H -1.356754 -0.838161 -0.914887
13 H -1.440837 -1.062338 0.844942""")
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/psi4/driver/molutil.py:259, in geometry(geom, name)
252 def geometry(geom: str, name: str = "default") -> core.Molecule:
253 """Function to create a molecule object of name *name* from the
254 geometry in string *geom*. Permitted for user use but deprecated
255 in driver in favor of explicit molecule-passing. Comments within
256 the string are filtered.
257
258 """
--> 259 molrec = qcel.molparse.from_string(
260 geom, enable_qm=True, missing_enabled_return_qm='minimal', enable_efp=True, missing_enabled_return_efp='none')
262 molecule = core.Molecule.from_dict(molrec['qm'])
263 if "geom" in molrec["qm"]:
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/qcelemental/molparse/from_string.py:276, in from_string(molstr, dtype, name, fix_com, fix_orientation, fix_symmetry, return_processed, enable_qm, enable_efp, missing_enabled_return_qm, missing_enabled_return_efp, verbose)
273 print(">>>\n")
275 # << 4 >> dict-->molspec
--> 276 molrec = from_input_arrays(
277 speclabel=True,
278 enable_qm=enable_qm,
279 enable_efp=enable_efp,
280 missing_enabled_return_qm=missing_enabled_return_qm,
281 missing_enabled_return_efp=missing_enabled_return_efp,
282 **molinit,
283 )
285 # replace from_arrays stamp with from_string stamp
286 if "qm" in molrec and molrec["qm"]:
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/qcelemental/molparse/from_arrays.py:97, in from_input_arrays(enable_qm, enable_efp, missing_enabled_return_qm, missing_enabled_return_efp, geom, elea, elez, elem, mass, real, elbl, name, units, input_units_to_au, fix_com, fix_orientation, fix_symmetry, fragment_separators, fragment_charges, fragment_multiplicities, molecular_charge, molecular_multiplicity, fragment_files, hint_types, geom_hints, geom_unsettled, variables, speclabel, tooclose, zero_ghost_fragments, nonphysical, mtol, copy, verbose)
95 if enable_qm:
96 dm = "qmvz" if geom_unsettled else "qm"
---> 97 processed = from_arrays(
98 domain=dm,
99 missing_enabled_return=missing_enabled_return_qm,
100 geom=geom,
101 elea=elea,
102 elez=elez,
103 elem=elem,
104 mass=mass,
105 real=real,
106 elbl=elbl,
107 name=name,
108 units=units,
109 input_units_to_au=input_units_to_au,
110 fix_com=fix_com,
111 fix_orientation=fix_orientation,
112 fix_symmetry=fix_symmetry,
113 fragment_separators=fragment_separators,
114 fragment_charges=fragment_charges,
115 fragment_multiplicities=fragment_multiplicities,
116 molecular_charge=molecular_charge,
117 molecular_multiplicity=molecular_multiplicity,
118 geom_unsettled=geom_unsettled,
119 variables=variables,
120 # processing details
121 speclabel=speclabel,
122 tooclose=tooclose,
123 zero_ghost_fragments=zero_ghost_fragments,
124 nonphysical=nonphysical,
125 mtol=mtol,
126 copy=copy,
127 verbose=1,
128 )
129 update_with_error(molinit, {"qm": processed})
130 if molinit["qm"] == {}:
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/qcelemental/molparse/from_arrays.py:380, in from_arrays(geom, elea, elez, elem, mass, real, elbl, name, units, input_units_to_au, fix_com, fix_orientation, fix_symmetry, fragment_separators, fragment_charges, fragment_multiplicities, molecular_charge, molecular_multiplicity, comment, provenance, connectivity, fragment_files, hint_types, geom_hints, geom_unsettled, variables, domain, missing_enabled_return, np_out, speclabel, tooclose, zero_ghost_fragments, nonphysical, mtol, copy, verbose)
377 update_with_error(molinit, processed)
379 Z_available = molinit["elez"] * molinit["real"] * 1.0
--> 380 processed = validate_and_fill_chgmult(
381 zeff=Z_available,
382 fragment_separators=molinit["fragment_separators"],
383 molecular_charge=molecular_charge,
384 fragment_charges=molinit["fragment_charges"],
385 molecular_multiplicity=molecular_multiplicity,
386 fragment_multiplicities=molinit["fragment_multiplicities"],
387 zero_ghost_fragments=zero_ghost_fragments,
388 verbose=verbose,
389 )
390 del molinit["fragment_charges"] # sometimes safe update is too picky about overwriting v_a_f_fragments values
391 del molinit["fragment_multiplicities"]
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/qcelemental/molparse/chgmult.py:501, in validate_and_fill_chgmult(zeff, fragment_separators, molecular_charge, fragment_charges, molecular_multiplicity, fragment_multiplicities, zero_ghost_fragments, verbose)
495 return fcgmp.format(final) if final == start else fcgmp.format("(" + str(int(final)) + ")")
497 # TODO could winnow down the exact_* lists a bit by ruling out
498 # independent values. do this if many-frag molecular systems take too
499 # long in the itertools.product
--> 501 c_final, fc_final, m_final, fm_final = reconcile(cgmp_exact_c, cgmp_exact_fc, cgmp_exact_m, cgmp_exact_fm)
503 c_text = stringify(molecular_charge, c_final)
504 fc_text = ", ".join((stringify(fs, ff) for fs, ff in zip(fragment_charges, fc_final)))
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/qcelemental/molparse/chgmult.py:491, in validate_and_fill_chgmult.<locals>.reconcile(exact_c, exact_fc, exact_m, exact_fm)
489 if verbose > -1:
490 print("\n\n" + "\n".join(text))
--> 491 raise ValidationError(err)
ValidationError: Inconsistent or unspecified chg/mult: sys chg: None, frag chg: [1.0], sys mult: None, frag mult: [2]
psi4.core.set_output_file('allyl_radical.log', False)
psi4.set_options({'reference':'uks'})
psi4.energy('blyp/3-21G', molecule=allyl_radical)
drawXYZ(allyl_radical)
7.1.2. Copper Water complex#
Exercise 2
Fix the error in the single-point wavefunction optimisation of the \(Cu[H_{2}O]^{2+}\) aqua complex (as observed in the gas phase) using the PBE exchange-correlation functional.
cu_h2o_2plus = psi4.geometry("""
2 2
Cu 0.009000 0.015000 0.004000
O 0.043000 1.451000 -0.430000
H 0.050000 1.450000 -1.397000
H -0.798000 1.844000 -0.158000
""")
drawXYZ(cu_h2o_2plus)
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
psi4.core.clean()
psi4.core.set_output_file('copper_water.log', False)
psi4.set_options({'reference':'uks'})
psi4.energy('pbe/3-21G', molecule=cu_h2o_2plus)
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
Cell In[15], line 3
1 psi4.core.set_output_file('copper_water.log', False)
2 psi4.set_options({'reference':'uks'})
----> 3 psi4.energy('pbe/3-21G', molecule=cu_h2o_2plus)
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/psi4/driver/driver.py:526, in energy(name, **kwargs)
524 logger.info(f"Compute energy(): method={lowername}, basis={core.get_global_option('BASIS').lower()}, molecule={molecule.name()}, nre={'w/EFP' if hasattr(molecule, 'EFP') else molecule.nuclear_repulsion_energy()}")
525 logger.debug("w/EFP" if hasattr(molecule, "EFP") else pp.pformat(molecule.to_dict()))
--> 526 wfn = procedures['energy'][lowername](lowername, molecule=molecule, **kwargs)
527 logger.info(f"Return energy(): {core.variable('CURRENT ENERGY')}")
529 for postcallback in hooks['energy']['post']:
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/psi4/driver/procrouting/proc.py:2571, in run_scf(name, **kwargs)
2568 if not core.has_global_option_changed('SCF_TYPE'):
2569 core.set_global_option('SCF_TYPE', 'DF')
-> 2571 scf_wfn = scf_helper(name, post_scf=False, **kwargs)
2572 returnvalue = scf_wfn.energy()
2574 ssuper = scf_wfn.functional()
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/psi4/driver/procrouting/proc.py:1767, in scf_helper(name, post_scf, **kwargs)
1764 core.print_out("\n ---------------------------------------------------------\n")
1765 core.print_out(" " + banner.center(58))
-> 1767 scf_wfn = scf_wavefunction_factory(name, base_wfn, core.get_option('SCF', 'REFERENCE'), **kwargs)
1769 # The wfn from_file routine adds the npy suffix if needed, but we add it here so that
1770 # we can use os.path.isfile to query whether the file exists before attempting to read
1771 read_filename = scf_wfn.get_scratch_filename(180) + '.npy'
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/psi4/driver/procrouting/proc.py:1435, in scf_wavefunction_factory(name, ref_wfn, reference, **kwargs)
1433 core.prepare_options_for_module("SCF")
1434 if reference in ["RHF", "RKS"]:
-> 1435 wfn = core.RHF(ref_wfn, superfunc)
1436 elif reference == "ROHF":
1437 wfn = core.ROHF(ref_wfn, superfunc)
RuntimeError:
Fatal Error: RHF: RHF reference is only for singlets.
Error occurred in file: /scratch/psilocaluser/conda-builds/psi4-multiout_1670993662927/work/psi4/src/psi4/libscf_solver/rhf.cc on line: 92
The most recent 5 function calls were:
7.1.3. N2 Geometry optimization#
Exercise 3
Fix the error in the geometry optimisation of \(N_2\) using the BP86 exchange-correlation functional.
n2 = psi4.geometry("""
0 1
N 0.000 0.000 0.000
N 0.000 0.000 1.098
""")
psi4.core.clean()
psi4.core.set_output_file('n2_opt.log', False)
psi4.energy('BP86/6-31+G*', molecule=n2)
7.1.4. \(^3 O\) Single point optimization#
Exercise 4
Fix the error in the single point calculation of \(^3 O\) using the B97-2 exchange-correlation functional.
psi4.core.clean()
psi4.core.clean_options()
o = psi4.geometry("""
0 3
O 0.0 0.0 0.0
""")
psi4.core.set_output_file('o3_singlept.log', False)
psi4.energy('b97-2/6-31+G*', molecule=o)
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
Cell In[11], line 2
1 psi4.core.set_output_file('o3_singlept.log', False)
----> 2 psi4.energy('b97-2/6-31+G*', molecule=o)
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/psi4/driver/driver.py:526, in energy(name, **kwargs)
524 logger.info(f"Compute energy(): method={lowername}, basis={core.get_global_option('BASIS').lower()}, molecule={molecule.name()}, nre={'w/EFP' if hasattr(molecule, 'EFP') else molecule.nuclear_repulsion_energy()}")
525 logger.debug("w/EFP" if hasattr(molecule, "EFP") else pp.pformat(molecule.to_dict()))
--> 526 wfn = procedures['energy'][lowername](lowername, molecule=molecule, **kwargs)
527 logger.info(f"Return energy(): {core.variable('CURRENT ENERGY')}")
529 for postcallback in hooks['energy']['post']:
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/psi4/driver/procrouting/proc.py:2571, in run_scf(name, **kwargs)
2568 if not core.has_global_option_changed('SCF_TYPE'):
2569 core.set_global_option('SCF_TYPE', 'DF')
-> 2571 scf_wfn = scf_helper(name, post_scf=False, **kwargs)
2572 returnvalue = scf_wfn.energy()
2574 ssuper = scf_wfn.functional()
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/psi4/driver/procrouting/proc.py:1767, in scf_helper(name, post_scf, **kwargs)
1764 core.print_out("\n ---------------------------------------------------------\n")
1765 core.print_out(" " + banner.center(58))
-> 1767 scf_wfn = scf_wavefunction_factory(name, base_wfn, core.get_option('SCF', 'REFERENCE'), **kwargs)
1769 # The wfn from_file routine adds the npy suffix if needed, but we add it here so that
1770 # we can use os.path.isfile to query whether the file exists before attempting to read
1771 read_filename = scf_wfn.get_scratch_filename(180) + '.npy'
File /opt/CompChem/psi4conda-1.7/lib/python3.10/site-packages/psi4/driver/procrouting/proc.py:1435, in scf_wavefunction_factory(name, ref_wfn, reference, **kwargs)
1433 core.prepare_options_for_module("SCF")
1434 if reference in ["RHF", "RKS"]:
-> 1435 wfn = core.RHF(ref_wfn, superfunc)
1436 elif reference == "ROHF":
1437 wfn = core.ROHF(ref_wfn, superfunc)
RuntimeError:
Fatal Error: RHF: RHF reference is only for singlets.
Error occurred in file: /scratch/psilocaluser/conda-builds/psi4-multiout_1670993662927/work/psi4/src/psi4/libscf_solver/rhf.cc on line: 92
The most recent 5 function calls were: