5
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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/.wfx file as input, but .pdb, .xyz, .chg, .cub and .grd files 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.

$\endgroup$
-2
$\begingroup$

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!

$\endgroup$
  • 1
    $\begingroup$ As far as I know, the standard (free) mode of AIMALL supports only a few atoms and basis functions and is in my eyes utterly useless. (It's also such a slow user interface that I would not recommend buying a professional version.) This post reads more like a commercial than an actual answer, and I disagree that QTAIM is to complex to be explained here. $\endgroup$ – Martin - マーチン Apr 5 '17 at 0:24
  • $\begingroup$ No one said it is too (not to) complex, I wrote "it would take really a lot". I would be pleased to see you explaining the atomic virial theorem and the formalism of bond/ring/cage critical points in one answer! I think I would learn a lot! Moreover I expressed my enthusiasm for AIMAll since I used it a couple of times and I think it is really a nice tool..of course if you work with many atoms it makes no sense to perform QCT. So I cannot understand what's wrong with my naswer..He ansked and I replied with a possible way to solve his problem. $\endgroup$ – Pierpy Apr 6 '17 at 6:49
  • $\begingroup$ Obviously the full derivation of QTAIM would be too (I'm very sorry for my typo) long for a post. However, the basis on what the electron density is analysed, divided and what criteria are used to obtain the partial charge are not. I can understand that you like AIMAll, however, it is proprietary software (quite expensive actually) and therefore not free. The standard edition has a limitation of 12 atoms, which is barely enough to cope with benzene. If you are well versed with QTAIM, maybe you can answer my question. $\endgroup$ – Martin - マーチン Apr 6 '17 at 7:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.