In GAMESS, a calculation of optimisation gives data about partial atomic charges (like Mulliken or Lowdin charges). Can you please let me know how to calculate the other types of partial atomic charges (like Hirshfeld, Maslen or Politzer charges)?
FYI: I have used B3LYP and MP2 for optimization.
Unfortunately GAMESS does not support any other than the named population analyses.
You will need to use a different program. For that matter I suggest MultiWFN, which comes with a variety of different population analyses and is still in active development. From the website:
Population analysis. Hirshfeld, VDD, Mulliken, Löwdin, Modified MPA (including three methods: SCPA, Stout & Politzer, Bickelhaupt), Becke, ADCH (Atomic dipole moment corrected Hirshfeld), CHELPG, Merz-Kollmann and AIM methods are supported.
The currently developed version MultiWFN 3.4 should be able to even read GAMESS output directly as stated in the update log:
GAMESS-US output file now can be used as input file (not comprehensively tested, currently only single point task at HF/DFT level is formally supported). The suffix of output file should be changed to .gms so that Multiwfn can properly recognize it
As a workaround for older versions, and a reliable way to import this data, you can export a .wfn file with GAMESS. I will give a brief summary for this, since I am still using MultiWFN 3.3.8. From the corresponding manual:
For example, the wavefunction represented by GTFs is enough for Hirshfeld population, so you can use
.fch/.molden/.31~.40/.wfn/.wfxfile as input, but.pdb,.xyz,.chg,.cuband.grdfiles do not carry any wavefunction information hence cannot be used.
As I have previously stated, .wfn files are technically outdated, especially when it comes to larger basis sets. However, for the time being, they are our only option to make it work.
I am using GAMESS version % DEC 2014 (R1) and the following input to demonstrate the procedure. I am using a Linux environment, so you might have to adjust this guide accordingly.
$CONTRL SCFTYP=RHF ! Restricted calculation.
RUNTYP=OPTIMIZE ! Geometry optimisation.
COORD=ZMT ! Z-matix specifies coordinates.
AIMPAC=.TRUE. ! Requests wfn file to be written.
$END
$BASIS GBASIS=STO NGAUSS=3 ! STO-3G minimal basis
$END
$GUESS GUESS=HUCKEL ! Start with Hückel MO.
$END
! Following group contains molecule specification.
! Blank lines are important.
$DATA
Water
Cnv 2
O
H 1 rOH
H 1 rOH 2 aHOH
rOH=1.09
aHOH=110.0
$END
Note that the indentations before the keywords are important. ! indicates a comment. Any additional output will be written to the punch file. This is dependent on how you set up GAMESS, but it normally is located in ~/scr/<rootfilename>.dat.
To extract the .wfn file you can use the following little script, or extract it manually from the .dat file.
#Script to extract wfn file from GAMESS dat (punch) file
#! /bin/bash
[ ! -z $1 ] && inputfile="$1" || exit 1
outputfile="${inputfile%.*}.wfn"
startpattern="----- TOP OF INPUT FILE FOR BADER'S AIMPAC PROGRAM -----"
endpattern="----- END OF INPUT FILE FOR BADER'S AIMPAC PROGRAM -----"
sed "/$startpattern/,/$endpattern/!d;//d" "$inputfile" > "$outputfile"
This should give you the following .wfn file for analysis:
Water GAUSSIAN 5 MOL ORBITALS 21 PRIMITIVES 3 NUCLEI O 1 (CENTRE 1) 0.00000000 0.00000000 -0.14560943 CHARGE = 8.0 H 2 (CENTRE 2) -1.43257773 0.00000000 1.05592484 CHARGE = 1.0 H 3 (CENTRE 3) 1.43257773 0.00000000 1.05592484 CHARGE = 1.0 CENTRE ASSIGNMENTS 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 3 3 CENTRE ASSIGNMENTS 3 TYPE ASSIGNMENTS 1 1 1 1 1 1 2 2 2 3 3 3 4 4 4 1 1 1 1 1 TYPE ASSIGNMENTS 1 EXPONENTS 1.3070932E+02 2.3808866E+01 6.4436083E+00 5.0331513E+00 1.1695961E+00 EXPONENTS 3.8038896E-01 5.0331513E+00 1.1695961E+00 3.8038896E-01 5.0331513E+00 EXPONENTS 1.1695961E+00 3.8038896E-01 5.0331513E+00 1.1695961E+00 3.8038896E-01 EXPONENTS 3.4252509E+00 6.2391373E-01 1.6885540E-01 3.4252509E+00 6.2391373E-01 EXPONENTS 1.6885540E-01 MO 1 OCC NO = 2.00000000 ORB. ENERGY = -20.25158070 4.22735195E+00 4.08851077E+00 1.27421022E+00 -6.18804850E-03 8.27701479E-03 6.24678653E-03 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 6.97684832E-03 4.38722945E-03 6.95010854E-04 -1.54630881E-03 -1.49552256E-03 -4.66089059E-04 -1.54630881E-03 -1.49552256E-03 -4.66089059E-04 MO 2 OCC NO = 2.00000000 ORB. ENERGY = -1.25754111 -9.93961353E-01 -9.61316148E-01 -2.99600254E-01 -2.02173727E-01 2.70423693E-01 2.04092795E-01 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 2.05800516E-01 1.29412887E-01 2.05011755E-02 4.30889918E-02 4.16737970E-02 1.29879022E-02 4.30889918E-02 4.16737970E-02 1.29879022E-02 MO 3 OCC NO = 2.00000000 ORB. ENERGY = -0.59384261 0.00000000E+00 0.00000000E+00 0.00000000E+00 -0.00000000E+00 0.00000000E+00 0.00000000E+00 1.02654522E+00 6.45519180E-01 1.02261083E-01 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 -1.24405583E-01 -1.20319664E-01 -3.74983838E-02 1.24405583E-01 1.20319664E-01 3.74983838E-02 MO 4 OCC NO = 2.00000000 ORB. ENERGY = -0.45973131 -4.42352599E-01 -4.27824176E-01 -1.33334109E-01 -1.28843153E-01 1.72338126E-01 1.30066155E-01 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 -1.26642108E+00 -7.96359559E-01 -1.26156734E-01 -8.17262213E-02 -7.90420432E-02 -2.46339524E-02 -8.17262213E-02 -7.90420432E-02 -2.46339524E-02 MO 5 OCC NO = 2.00000000 ORB. ENERGY = -0.39261686 0.00000000E+00 0.00000000E+00 0.00000000E+00 -0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 1.67545020E+00 1.05356804E+00 1.66902879E-01 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 END DATA RHF ENERGY = -74.9659012162 VIRIAL(-V/T) = 2.00600297
You can now use this file and load it with MultiWFN. After the start-up process you should find the key values printed to the screen like this (I abridged it a little bit):
Multiwfn -- A Multifunctional Wavefunction Analyzer (for Linux)
Version 3.3.8, release date: 2015-Dec-1
[...]
System energy: -74.965901216200 Hartree, Virial ratio: 2.00600297
Total/Alpha/Beta electrons: 10.0000 5.0000 5.0000
Net charge: 0.00000 Expected multiplicity: 1
The number of orbitals: 5, Atoms: 3, GTFs: 21
This is restricted close-shell single-determinant wavefunction
Title line of this file: Water
Formula: H2 O1
Molecule weight: 18.01528
Loaded water-sto3g.wfn successfully!
------------ Main function menu ------------
[...]
7 Population analysis
[...]
You may now choose 7 and enter the following menu:
============== Population analysis ==============
0 Return
1 Hirshfeld population
2 Voronoi deformation density (VDD) population
10 Becke atomic charge with atomic dipole moment correction
11 Atomic dipole corrected Hirshfeld population (ADCH)
12 CHELPG ESP fitting charge
13 Merz-Kollmann (MK) ESP fitting charge
14 AIM charge
Follow the guidance of the manual and the on-screen instructions to produce the required charges. For comparison here are the computed Hirshfeld charges, written to water-sto3g.chg with the following command sequence 7; 1; 1; y.
O 0.000000 0.000000 -0.077053 -0.296017 H -0.758088 0.000000 0.558771 0.147959 H 0.758088 0.000000 0.558771 0.147959
The element identifier is followed by the Cartesian coordinates, which are followed by the charges.
In principle all other implemented charges are equally as easy to compute. Keep in mind, that basin based charges (e.g. AIM) need more computing power and depend heavily on molecule and basis set size.
Apart from this procedure, you can obtain natural charges with the NBO6 program, which conveniently interfaces with GAMESS. Unfortunately you have to buy it.
Partial charges are tricky because they are not observables, therefore it is difficult to define them rigorously, however in same cases it can be done. I usually never use analysis of the population-based charges since they can be very arbitrary in assigning electrons to a given atom. I know with Gaussian you can obtain Hirshfeld charges (therefore fitting the ESP) but honestly in GAMESS I do not know. A more rigorous way of dividing the electron density is to use the Quantum Theory of Atoms in Molecules. I cannot go into details since it would take really a lot, however it is possible to obtain all kind of atomic properties partitioning the molecular electron density in basins of charge containing the nuclei, following rigorous topological arguments. From the atomic densities it is possible to obtain the expectation values of observables (atomic dipole/quadrupole moments, energy, magnetic properties etc...) and integrating the atomic density it is possible to know "how many" electrons are contained in the atomic basis, therefore the atomic partial charge. The free software AIMALL can definitely help you. Using GAMESS it is possible to generate wavefunction files (file.wfn) that contain the information needed to AIMALL to perform the decomposition of the molecular density. It is very easy and free, I strongly suggest you this approach.
In case you were dealing with planar molecules, there is another rigorous way to obtain information about the partial charge of atoms via the Atomic Polar Tensors, defined as $$P_{\sigma I}=\frac{\partial \mu_{\sigma}}{\partial r_{I}}$$ where sigma can be the x,y,z molecular dipole component and I is one of the $3N$ cartesian coordinates of the atoms defining the molecule. This tensor is analytically calculated by almost every Quantum Chemistry software, GAMESS included.
I know it is difficult to be absolutely exhaustive within an answer, there are many more things and details to say, however I give you some useful link where you can pick up more theory if needed
AIMALL web page: http://aim.tkgristmill.com/
Paper of Dinur of the charge in planar molecules: http://pubs.acs.org/doi/abs/10.1021/j100169a030
I really hope it can be helpful!
