[Pw_forum] point charge-spin

[Giovanni La Penna]

As for point charges, here is a small recipe for isolated
molecules (see below) for the ESP/RESP procedure [A].
Reference is Bayly et al, J. Phys. Chem., 97, 10269 (1993).
Final point charges for organo-metallic molecules compare
fairly well with standard or published values.
Another possibility is to compute charges by integrating
over atomic basins built on the basis of QT-AIM (Bader) [B].

[A]
1) From PW or CP, write into the restart
directory (I guess into the charge-density.dat file)
the electrostatic potential in place of electron
density. This requires a call before writefile
in CPV/cpr.f90 (in case of CP). Paolo Giannozzi
wrote the subroutine to do this in CP. It would be
nice to have it in some contrib space.
In the input, it is recommended to use
explicit ion_radius(i), especially if you
want to compare with other programs (like CPMD,
for which the default radii are all 1.2 bohr).

2) Once the restart is generated, postprocess
it with PP as to print the electron density
into a friendly format, like a gaussian-cube file.
This file will contain the electrostatic potential
on a 3d-grid. PP can process restart produced also
by CP, provided disk_io='high'.

3) Write your own program to map the 3d-grid on
the molecular solvent-accessible surface. This will consist in reading
the 3d-grid, reading the molecular surface produced
with vdw radii and a rolling probe, and finally
interpolate the potential at the surface point
using neighbour points in the grid. For the molecular
surface, I recommend an old program called NSC (numerical surface
calculation):

Eisenhaber et al, J. Comput. Chem., 16, 273 (1995)

WARNING: this latter step for non-periodic systems!
Non-periodicity is assumed in the ESP/RESP procedure,
that consists in fitting the electrostatic surface
potential with a LIMITED sum of Coulomb interactions.

The interpolation can be very coarse, the entire method
is plenty of error sources.
The output of this procedure is a list with
x y z U
for each point of the molecular surface.
In case step 4 is performed with the aresp code
that was once part of Amber, the output for a water
molecule is:

     3 1072    0
                   -0.9982123E-32   0.0000000E+00   0.2313846E+00
                   -0.2610123E-32   0.1494187E+01  -0.9255383E+00
                   -0.1829851E-15  -0.1494187E+01  -0.9255383E+00
    -6.5264638e-01   2.0485851e+00   1.0162867e+00   1.6945994e+00
    -7.2948945e-01   2.4209246e+00   6.0282018e-01   1.7305686e+00
...

where 3 is the number of atoms, 1072 is the number of surface points,
0 is the net charge, the three following lines are the atomic coordinates
in bohr, and the following lines are x,y,z,U (bohr,kcal/mol)
for the 1072 surface points. This format was once produced
by Gaussian.

4) Finally, one of the RESP codes available can be used.
I used the original aresp code that was part of the Amber 4.1
code.
I think it is not worth to develop the RESP part in QE,
but for sure steps 1-3 could be done relatively easily.
And maybe with a minimum-image based periodicity.
I really don't know how to do it, but I am available
for testing.

[B]
This is simpler, but for reasons that will be apparent to
everybody in a few applications, it is rarely used.

1) Postprocess the restart with PP and write the electron density
on a fine 3d-grid (cube-file). The size of 10 pm associated
to a density energy-cutoff of 250 Ry is rather accurate.

2) Process this electron density cube-file with a real-space
QT-AIM code. I used:

Sanville et al., J. Comp. Chem., 28, 899 (2007).
http://theory.cm.utexas.edu/bader/

Charge are in ACF.dat.
This works also for periodic systems (maybe only orthorhombic?).
Problems occur with polar X-H bonds or in all cases
where the zero-flux of density comes too close to atoms desribed
with pseudo-potentials.

Maybe, this will help,

              Giovanni

============================================================
Giovanni La Penna - National research council (Cnr)
Institute for chemistry of organo-metallic compounds (Iccom)
via Madonna del Piano 10,
I-50019 Sesto Fiorentino, Firenze, Italy
tel.: +39 055 522-5264, fax: +39 055 522-5203
e-mail: glapenna at iccom.cnr.it - http://www.iccom.cnr.it/lapenna
skype: giovannilapenna
============================================================

On Wed, 1 Dec 2010, David Grifith wrote:

> Dear All
> 
> 
> I am trying to calculate the point charge and the point spin density on each
> atom of a unit cell, certainly "WITHOUT" considering a 3D vector which is
> ...