Hi All, I am an old chemist trying to learn new tricks. I am new to DFT and Psi4. I am playing with acetone, which is a symmetrical molecule but regardless of the basis set/theory I choose (6-31G, 6-311G, 6-311+G(d,p)/B3LYP, M06-2X, wB97M-V) the charges are not symmetrical. I have a set of data generated on Gaussian that shows the charges being symmetrical. What am I doing wrong? Any help is greatly appreciated.
Can you supply an example output file? There’s lots of information (what exactly is your input file, how big is the asymmetry, what Psi4 version are you using, how are you getting the charges) missing right now.
Hi and thank you for your reply. Unfortunately I cannot upload a file yet. Below is the end of the file. I did a stepwise geometry optimization starting with B3LYP, then M06-2X, then to wB97M-V to save time. I then do an energy calculation and generate a .fchk file for Multiwfn input to get Fukui and Dual Descriptor data.
Properties will be evaluated at 0.000000, 0.000000, 0.000000 [a0]
OEProp: No title given, name of density matrix used for the following properties is ‘SCF alpha density’
Mulliken Charges: (a.u.)
Center Symbol Alpha Beta Spin Total
1 C 3.30975 3.30975 0.00000 -0.61950
2 C 2.84006 2.84006 -0.00000 0.31989
3 O 4.16817 4.16817 -0.00000 -0.33634
4 C 3.31982 3.31982 0.00000 -0.63963
5 H 0.39081 0.39081 -0.00000 0.21838
6 H 0.39572 0.39572 -0.00000 0.20855
7 H 0.39572 0.39572 -0.00000 0.20855
8 H 0.40479 0.40479 0.00000 0.19041
9 H 0.38758 0.38758 0.00000 0.22485
10 H 0.38758 0.38758 0.00000 0.22485
Total alpha = 16.00000, Total beta = 16.00000, Total charge = -0.00000
Molecular point group: cs
Full point group: Cs
Geometry (in Angstrom), charge = 0, multiplicity = 1:
Center X Y Z Mass
------------ ----------------- ----------------- ----------------- -----------------
C -1.395701764705 -0.447897468896 0.000000000000 12.000000000000
C 0.020306659865 0.080054290412 0.000000000000 12.000000000000
O 0.239185037104 1.297231183844 0.000000000000 15.994914619570
C 1.143154899748 -0.922904695516 0.000000000000 12.000000000000
H 2.099146397741 -0.406643599480 0.000000000000 1.007825032230
H 1.075840526310 -1.568550784773 0.878882249439 1.007825032230
H 1.075840526310 -1.568550784773 -0.878882249439 1.007825032230
H -1.436526767216 -1.535462590696 0.000000000000 1.007825032230
H -1.922548229315 -0.070039347489 0.877170446261 1.007825032230
H -1.922548229315 -0.070039347489 -0.877170446261 1.007825032230
Variable Map:
“-D ENERGY” => 0.000000000000
“CURRENT DIPOLE X” => -0.682205556212
“CURRENT DIPOLE Y” => -3.074052686072
“CURRENT DIPOLE Z” => 0.000000000000
“CURRENT ENERGY” => -193.052125652775
“CURRENT REFERENCE ENERGY” => -193.052125652775
“DFT FUNCTIONAL TOTAL ENERGY” => -193.052125652775
“DFT TOTAL ENERGY” => -193.052125652775
“DFT VV10 ENERGY” => 0.126792917804
“DFT XC ENERGY” => -18.989540218941
“NUCLEAR REPULSION ENERGY” => 119.000086808206
“ONE-ELECTRON ENERGY” => -496.185502041880
“OPTIMIZATION ITERATIONS” => 3.000000000000
“PCM POLARIZATION ENERGY” => 0.000000000000
“SCF DIPOLE X” => -0.682205556212
“SCF DIPOLE Y” => -3.074052686072
“SCF DIPOLE Z” => 0.000000000000
“SCF ITERATION ENERGY” => -193.052125652775
“SCF ITERATIONS” => 7.000000000000
“SCF TOTAL ENERGY” => -193.052125652775
“TWO-ELECTRON ENERGY” => 202.996036882037
Psi4 stopped on: Tuesday, 23 February 2021 08:35PM
Psi4 wall time for execution: 1:41:42.52
*** Psi4 exiting successfully. Buy a developer a beer!
Below is the input portion:
-----------------------------------------------------------------------
Psi4: An Open-Source Ab Initio Electronic Structure Package
Psi4 1.3.2 release
Git: Rev {HEAD} ecbda83
R. M. Parrish, L. A. Burns, D. G. A. Smith, A. C. Simmonett,
A. E. DePrince III, E. G. Hohenstein, U. Bozkaya, A. Yu. Sokolov,
R. Di Remigio, R. M. Richard, J. F. Gonthier, A. M. James,
H. R. McAlexander, A. Kumar, M. Saitow, X. Wang, B. P. Pritchard,
P. Verma, H. F. Schaefer III, K. Patkowski, R. A. King, E. F. Valeev,
F. A. Evangelista, J. M. Turney, T. D. Crawford, and C. D. Sherrill,
J. Chem. Theory Comput. 13(7) pp 3185--3197 (2017).
(doi: 10.1021/acs.jctc.7b00174)
Additional Contributions by
P. Kraus, H. Kruse, M. H. Lechner, M. C. Schieber, R. A. Shaw,
A. Alenaizan, R. Galvelis, Z. L. Glick, S. Lehtola, and J. P. Misiewicz
-----------------------------------------------------------------------
Psi4 started on: Tuesday, 23 February 2021 06:53PM
Process ID: 27677
Host: euawsvdcpcapp01
PSIDATADIR: /usr/local/miniconda/share/psi4
Memory: 500.0 MiB
Threads: 1
==> Input File <==
memory 15000 mb
set_num_threads (4)
molecule mol {
0 1
C -1.393852143065 -0.448182806807 0.000000000000
C 0.020706814644 0.080666281277 0.000000000000
O 0.237599947165 1.294682757373 0.000000000000
C 1.142719351301 -0.920517959929 0.000000000000
H 2.098417670142 -0.405516072883 0.000000000000
H 1.071934885514 -1.565049400518 0.878325435704
H 1.071934885514 -1.565049400518 -0.878325435704
H -1.431366613082 -1.534743789971 0.000000000000
H -1.919080888815 -0.070395132935 0.877121496922
H -1.919080888815 -0.070395132935 -0.877121496922
}
set scf_type df
set basis 6-311G
set reference uks
mol.update_geometry()
mol.symmetrize(1e-3)
e, wfn = optimize(‘wB97M-V’, return_wfn=True)
oeprop(wfn, “MULLIKEN_CHARGES”)
mol.print_out()
print_variables()--------------------------------------------------------------------------
Geomerty_Optimize_Neutral_wB97M_V_6_311G_Output.dat (1.2 MB)
I was upgraded and able to upload the geometry optimization file. The starting geometry was the optimized geometry from 6-311G/M06-2X. The output file shows Mulliken Charges but I use Multiwfn to generate Atomic Dipole Corrected Hershfield Charges. Those show the same trend as the output from Psi4.
Just to be absolutely clear, what symmetry are you expecting? Cs or C2v?
Let me show my ignorance and say I don’t know. Time to go back to my old textbooks for a refresher. I will also do some reading on specifying symmetry in Psi4. Thanks for providing direction.
The structure you’ve provided is staggered, which means it only has a plane of symmetry (C_s) , and thus the charges on the atoms will not be symmetric, except for the hydrogen atoms on either side of the mirror plane. If you choose an eclipsed geometry (C_2v), then you have two planes of symmetry, and atoms on either side of the plane will have identical charges.
It’s exactly as Prof. Crawford said. To phrase it differently, the connectivity being symmetric doesn’t mean the molecular geometry is symmetric, and it’s the symmetry of the molecular geometry that determines the symmetry of the charges.
Your output shows hydrogen atoms 6 and 7 have the same charge, as do 9 and 10. All hydrogen charges are similar. That’s as much as I’d expect from a staggered (Cs) structure.
If Gaussian shows more symmetry, then check which geometry Gaussian optimized to: staggered or eclipsed? If it optimized to an eclipsed geometry, then you should see high symmetry among the hydrogen charges, and I believe everything is wrapped up nicely.
Please report back with whatever you find.
Please try the following c2v geometry (b3lyp/6-311G* optimized)
C 0.000000000000 0.000000000000 0.102965814809
O 0.000000000000 0.000000000000 1.311944341861
C -1.291335016568 -0.000000000000 -0.697503321657
C 1.291335016568 0.000000000000 -0.697503321657
H -2.144911944720 -0.000000000000 -0.020976036079
H 2.144911944720 0.000000000000 -0.020976036079
H -1.343339911801 0.878498098718 -1.348861972381
H -1.343339911801 -0.878498098718 -1.348861972381
H 1.343339911801 -0.878498098718 -1.348861972381
H 1.343339911801 0.878498098718 -1.348861972381
Hi All,
Many thanks for all the instruction. I am definitely learning a lot . I tried using the molecule option “symmetry c2v” in the molecule specification and starting geometry built from WebMO but got an error stating that C2V was not part of Cs. I am still working to learn this part. If a PubChem structure is available as a starting point, is that a better option than building in WebMO? I will also try the geometry provided above.
Ultimately, my goal is to compare possible reactivities on thousands of compounds, which geometry will more accurately represent the real world behavior. Acetone is just one of a larger set of symmetrical compounds in my data set so I want to handle them properly.
The geometry supplied above does give C2v optimization and symmetrical charges. Many thanks! So it is clear that the initial geometry input is important for Psi4 to recognize the correct geometry. Can you point me to any references that I can read so I can learn how to choose/generate the most appropriate initial geometry?
E.g. CCCBDB List calculated geometries
I recommend to start with b3lyp/6-31G*, for anions I choose 6-31+G*.
For further optimizations with dft I recommend the basis sets pcseg-2 or def2-tzvpp.
Hi Dr. Pfaff, I moved to the 6-311 set because some compounds in my data set contain I. I assumed it would be best to use the same basis set for all compounds to get the best comparison. Is it OK to mix basis sets when comparing Fukui Function and Dual Descriptor data between compounds?
In general you schould use the same basis set. From the BSE (Basis Set Exchange) you can get
6-311G* or def2-tzvpp for I. I strogly recommend the use of polarization functions…
Thanks for the NIST reference! That is a huge help and also the def2-tzvpp.
Much of the literature I have read also recommends the use of polarizations functions as well. I appreciate your guidance.
Because you’re new to computational chemistry, I should clarify something:
Every geometry optimization locates a stationary point (almost always a local minimum, but in rare cases a transition state or saddle point) in coordinate space. However, there’s no guarantee that it’s the local minimum you’re interested in. In the sense of “not the minimum I want,” then yes, the geometry you found is “incorrect.” In terms of “actually a local minimum,” the geometry minimization is correct.
If you have a guess as to what the local minimum looks like, e.g., I want the C2v geometry, then indeed, choosing a geometry close to it greatly increases the chance that your optimizer will get to the minimum you want.
If your criterion for what geometry you want relies on energetics, e.g., lowest enthalpy, it’s standard to do a search for multiple conformers. I’m not familiar enough with conformer searches to recommend a way to do this with confidence. (I mostly work on method development now.)
Thank you for the clarification. I am certainly learning a lot and have made months of progress with the help provided here. I think I am getting closer to being able to make fewer bad assumptions. It is challenging when you don’t know what you don’t know but that is why we read and ask questions.