Hi,
(Please assume that I have 0 idea of what I’m doing and go easy on me, I’ve never made a post on a forum like this so I’m sorry if it’s all over the place. Happy to follow up with additional info).
I’m trying to optimize/calc SPEs for a variety of different palladium based organometallic complexes using psi4 (in jupyter lab). I eventually want to use the following method for those calculations: MN15/6-311++G(2d,p)/SDD(Pd) // MN15L/6-31+G(d)/LANL2DZ(Pd). I noticed SDD is not an option but seeing as I can’t get past the initial set-up and geometry optimization, I can pin that for later!
(I’ve been able to do basic calculations on organic species using psi4)
My background is strictly in g16 so this is pretty new to me overall. It’s very rare where I have to consider the details of the output of a g16.log aside from the values I’m most interested in. I have a basic handle of python (mostly for data viz/statistical learning).
Computer specifications:
Apple M3 Max Silicon
16 cores (4 efficiency) /16 threads
64 GB RAM
I’ve created a conda environment called “psi4dev” and installed as follows:
conda install psi4 python=3.10 -c conda-forge
Followed by installation of jupyter lab
conda install jupyter
Some previous posts suggested checking the version of libint (and that its different than libint2). When I check, I get:
(psi4dev) orz@MacBook-Pro-de-Brandon ~ % conda list libint
# packages in environment at /Users/orz/anaconda3/envs/psi4dev:
#
# Name Version Build Channel
libint 2.8.2 h078ce10_0 conda-forge
and
(psi4dev) orz@MacBook-Pro-de-Brandon ~ % conda list libint2
# packages in environment at /Users/orz/anaconda3/envs/psi4dev:
#
# Name Version Build Channel
So libint2 is not installed (I learn about this after the fact so lets jump to the code with the error):
first I import psi4 in jupyter lab followed by setting up my molecule:
psi4.set_memory('8 GB') #probably wrong
psi4.set_num_threads(8) #probably wrong
molecule_name = "dppp"
dppp = psi4.geometry("""
0 1
C 1.90867300 0.20404100 -2.22940800
C 0.90779200 1.32772800 -2.49907600
C -0.57874600 0.96855400 -2.48175900
H 2.88587000 0.46782500 -2.63824100
H 1.08685200 2.16412900 -1.82417400
H 1.12752700 1.70920500 -3.49874100
H -1.13068100 1.70097800 -3.07330300
H -0.72910900 -0.00578800 -2.94585500
H 1.60652400 -0.70655400 -2.74709800
Pd 0.04300400 -0.48287200 0.28465700
P -1.34236200 0.91077000 -0.78055300
P 2.12719800 -0.30603600 -0.45297400
C -1.12006700 2.66623900 -0.27752200
C -0.14490300 2.97423900 0.66816200
C -1.83073800 3.70379300 -0.88027900
C 0.12770700 4.29464600 0.99677100
H 0.40722500 2.16765000 1.13216000
C -1.56819100 5.02239300 -0.54177800
H -2.59949700 3.47782400 -1.60765200
C -0.58601700 5.32115300 0.39524600
H 0.89698100 4.51597800 1.72464600
H -2.12969000 5.81991100 -1.01110500
H -0.37910900 6.35188800 0.65312500
C -3.13904400 0.85363000 -1.15241600
C -3.65464300 0.37047500 -2.35431000
C -4.03264500 1.21193900 -0.13922500
C -5.02656900 0.26077700 -2.54274500
H -2.99098700 0.07610000 -3.15499500
C -5.39907300 1.11434500 -0.33256300
H -3.64770000 1.57877700 0.80426400
C -5.90370000 0.63547900 -1.53693600
H -5.40854900 -0.11453900 -3.48349000
H -6.07548800 1.40760100 0.45989000
H -6.97219900 0.55197800 -1.68556600
C 3.01141300 1.10875600 0.31145800
C 3.57192700 2.16468900 -0.40305300
C 3.06428900 1.14248500 1.70614600
C 4.16536600 3.23104800 0.25967800
H 3.54702700 2.16868900 -1.48332300
C 3.67060900 2.19730100 2.36780200
H 2.60544700 0.33818000 2.26795000
C 4.21817600 3.25047900 1.64504200
H 4.58657800 4.04874700 -0.31050900
H 3.70467500 2.20532400 3.44938100
H 4.67997900 4.08270900 2.16012300
C 3.45626700 -1.56247700 -0.66031500
C 4.81284800 -1.24595600 -0.69646000
C 3.06719600 -2.89355500 -0.80264400
C 5.76088100 -2.24350600 -0.87548400
H 5.13031400 -0.21907300 -0.57413400
C 4.01400300 -3.88734200 -0.99902900
H 2.01444400 -3.14496900 -0.74437400
C 5.36452800 -3.56486900 -1.03084300
H 6.81257900 -1.98789000 -0.89278900
H 3.69919400 -4.91688100 -1.10985000
H 6.10530400 -4.34141800 -1.17020700
""")
psi4.basis_helper("""
assign 6-31G_D_P_ ## I know this is the wrong basis set for the method I mentioned earlier
assign pd LANL2DZ
""",
name='MIXED1',
key='BASIS',
set_option = True)
This structure was previously optimized using g16 so I was curious to see the SPE which I performed as follows:
energy = psi4.energy('scf')
This calculation converges but I have a few questions about the output. First, the loading of the basis set as specified seems correct:
Scratch directory: /tmp/
gradient() will perform analytic gradient computation.
=> Libint2 <=
Primary basis highest AM E, G, H: 5, 4, 3
Auxiliary basis highest AM E, G, H: 6, 5, 4
Onebody basis highest AM E, G, H: 6, 5, 4
Solid Harmonics ordering: gaussian
*** tstart() called on MacBook-Pro-de-Brandon.local
*** at Sat Feb 24 18:50:07 2024
=> Loading Basis Set <=
Name: MIXED1
Role: ORBITAL
Keyword: BASIS
atoms 1-3, 13-16, 18, 20, 24-27, 29, 31, 35-38, 40, 42, 46-49, 51, 53 entry C line 115 file /Users/orz/anaconda3/envs/psi4-2/share/psi4/basis/6-31g_d_p_.gbs
atoms 4-9, 17, 19, 21-23, 28, 30, 32-34, 39, 41, 43-45, 50, 52, 54-56 entry H line 44 file /Users/orz/anaconda3/envs/psi4-2/share/psi4/basis/6-31g_d_p_.gbs
atoms 10 entry **PD** line 922 (ECP: line 2407) file /Users/orz/anaconda3/envs/psi4-2/share/psi4/basis/**lanl2dz**.gbs
atoms 11-12 entry P line 296 file /Users/orz/anaconda3/envs/psi4-2/share/psi4/basis/6-31g_d_p_.gbs
!!! WARNING: ECP capability is in beta. Please check occupations closely. !!!
But then later in the output, it says this:
Basis Set: MIXED1
Blend: 6-31G_D_P_ + LANL2DZ
Number of shells: 264
Number of basis functions: 597
Number of Cartesian functions: 597
Spherical Harmonics?: false
Max angular momentum: 2
Core potential: MIXED1
Number of shells: 4
Number of ECP primitives: 20
Number of ECP core electrons: 28
Max angular momentum: 3
=> Loading Basis Set <=
**Name: (MIXED1 AUX)**
Role: JKFIT
Keyword: DF_BASIS_SCF
atoms 1-3, 13-16, 18, 20, 24-27, 29, 31, 35-38, 40, 42, 46-49, 51, 53 entry C line 121 file /Users/orz/anaconda3/envs/psi4dev/share/psi4/basis/cc-pvdz-jkfit.gbs
atoms 4-9, 17, 19, 21-23, 28, 30, 32-34, 39, 41, 43-45, 50, 52, 54-56 entry H line 51 file /Users/orz/anaconda3/envs/psi4dev/share/psi4/basis/cc-pvdz-jkfit.gbs
atoms 10 entry PD line 4258 file /Users/orz/anaconda3/envs/psi4dev/share/psi4/basis/def2-universal-jkfit.gbs
atoms 11-12 entry P line 519 file /Users/orz/anaconda3/envs/psi4dev/share/psi4/basis/cc-pvdz-jkfit.gbs
What is MIXED1 AUX? I see the basis are different than the ones I specified and I could just be misunderstanding the theory here…
I’m going to ignore my likely misunderstanding of how to specify the mixed basis set correctly in python and move on to trying to geometry optimize the start complex by calling:
psi4.optimize('MN15L/MIXED1', molecule=dppp)
Which again, I couldn’t find any case studies in my search for how to use the mixed basis set for an optimization. When I try to run this, I get an error:
RuntimeError Traceback (most recent call last)
Cell In[4], line 1
----> 1 psi4.optimize('MN15L/MIXED1', molecule=dppp)
File ~/anaconda3/envs/psi4dev/lib/python3.10/site-packages/psi4/driver/driver.py:1245, in optimize(name, **kwargs)
1242 opt_calcs = opt_object.calculations_needed() # tuple of strings ('energy', 'gradient', etc)
1244 # Compute the gradient - no longer need to worry about opt_data being wiped
-> 1245 G, wfn = gradient(lowername, return_wfn=True, molecule=molecule, **kwargs)
1246 thisenergy = core.variable('CURRENT ENERGY')
1247 opt_object.E = thisenergy
File ~/anaconda3/envs/psi4dev/lib/python3.10/site-packages/psi4/driver/driver.py:638, in gradient(name, **kwargs)
636 logger.info(f"Compute gradient(): method={lowername}, basis={core.get_global_option('BASIS').lower()}, molecule={molecule.name()}, nre={'w/EFP' if hasattr(molecule, 'EFP') else molecule.nuclear_repulsion_energy()}")
637 logger.debug("w/EFP" if hasattr(molecule, "EFP") else pp.pformat(molecule.to_dict()))
--> 638 wfn = procedures['gradient'][lowername](lowername, molecule=molecule, **kwargs)
639 logger.info(f"Return gradient(): {core.variable('CURRENT ENERGY')}")
640 logger.info(nppp(wfn.gradient().np))
File ~/anaconda3/envs/psi4dev/lib/python3.10/site-packages/psi4/driver/procrouting/proc.py:96, in select_scf_gradient(name, **kwargs)
94 return
95 else:
---> 96 return func(name, **kwargs)
File ~/anaconda3/envs/psi4dev/lib/python3.10/site-packages/psi4/driver/procrouting/proc.py:2697, in run_scf_gradient(name, **kwargs)
2694 disp_grad = ref_wfn._disp_functor.compute_gradient(ref_wfn.molecule(), ref_wfn)
2695 ref_wfn.set_variable("-D Gradient", disp_grad)
-> 2697 grad = core.scfgrad(ref_wfn)
2699 ref_wfn.set_gradient(grad)
2701 ref_wfn.set_variable("SCF TOTAL GRADIENT", grad) # P::e SCF
RuntimeError: Engine::lmax_exceeded -- angular momentum limit exceeded
I’m not sure what the next steps here and would appreciate any pointers in all parts of my code/handling as I’m really motivated to learn this software well for future projects. Please let me know if there are any other details I can include to clarify!
If we can get this working, I’m more than happy to have it included as a worked example. Many of my colleagues have expressed interest in using psi4 and work heavily in the area of organometallic catalysis. We rely almost exclusively on g16 for DFT studies with some people using ORCA. The integration of psi4 with python is extremely attractive for workflows in parallel chemistry and chemoinformatics.
Thank you so much!
Brandon