7.1. Fixing errors in the calculations#
When you hand in the homework for this exercise (at least for all problems in this page), you should give the modified code for each question and mark the position where you modified the code. Only the modified part is necessary. Thanks!
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: