Dear Giovanni<br><br> Many thanks for your kind attention and very detailed answer. It has raised up our knowledge around the topic.<br><br>-- <br>Sincerely Yours<br>David G.<br>JCU<br><br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
<br>
<br>
As for point charges, here is a small recipe for isolated<br>
molecules (see below) for the ESP/RESP procedure [A].<br>
Reference is Bayly et al, J. Phys. Chem., 97, 10269 (1993).<br>
Final point charges for organo-metallic molecules compare<br>
fairly well with standard or published values.<br>
Another possibility is to compute charges by integrating<br>
over atomic basins built on the basis of QT-AIM (Bader) [B].<br>
<br>
[A]<br>
1) From PW or CP, write into the restart<br>
directory (I guess into the charge-density.dat file)<br>
the electrostatic potential in place of electron<br>
density. This requires a call before writefile<br>
in CPV/cpr.f90 (in case of CP). Paolo Giannozzi<br>
wrote the subroutine to do this in CP. It would be<br>
nice to have it in some contrib space.<br>
In the input, it is recommended to use<br>
explicit ion_radius(i), especially if you<br>
want to compare with other programs (like CPMD,<br>
for which the default radii are all 1.2 bohr).<br>
<br>
2) Once the restart is generated, postprocess<br>
it with PP as to print the electron density<br>
into a friendly format, like a gaussian-cube file.<br>
This file will contain the electrostatic potential<br>
on a 3d-grid. PP can process restart produced also<br>
by CP, provided disk_io='high'.<br>
<br>
3) Write your own program to map the 3d-grid on<br>
the molecular solvent-accessible surface. This will consist in reading<br>
the 3d-grid, reading the molecular surface produced<br>
with vdw radii and a rolling probe, and finally<br>
interpolate the potential at the surface point<br>
using neighbour points in the grid. For the molecular<br>
surface, I recommend an old program called NSC (numerical surface<br>
calculation):<br>
<br>
Eisenhaber et al, J. Comput. Chem., 16, 273 (1995)<br>
<br>
WARNING: this latter step for non-periodic systems!<br>
Non-periodicity is assumed in the ESP/RESP procedure,<br>
that consists in fitting the electrostatic surface<br>
potential with a LIMITED sum of Coulomb interactions.<br>
<br>
The interpolation can be very coarse, the entire method<br>
is plenty of error sources.<br>
The output of this procedure is a list with<br>
x y z U<br>
for each point of the molecular surface.<br>
In case step 4 is performed with the aresp code<br>
that was once part of Amber, the output for a water<br>
molecule is:<br>
<br>
3 1072 0<br>
-0.9982123E-32 0.0000000E+00 0.2313846E+00<br>
-0.2610123E-32 0.1494187E+01 -0.9255383E+00<br>
-0.1829851E-15 -0.1494187E+01 -0.9255383E+00<br>
-6.5264638e-01 2.0485851e+00 1.0162867e+00 1.6945994e+00<br>
-7.2948945e-01 2.4209246e+00 6.0282018e-01 1.7305686e+00<br>
...<br>
<br>
where 3 is the number of atoms, 1072 is the number of surface points,<br>
0 is the net charge, the three following lines are the atomic coordinates<br>
in bohr, and the following lines are x,y,z,U (bohr,kcal/mol)<br>
for the 1072 surface points. This format was once produced<br>
by Gaussian.<br>
<br>
4) Finally, one of the RESP codes available can be used.<br>
I used the original aresp code that was part of the Amber 4.1<br>
code.<br>
I think it is not worth to develop the RESP part in QE,<br>
but for sure steps 1-3 could be done relatively easily.<br>
And maybe with a minimum-image based periodicity.<br>
I really don't know how to do it, but I am available<br>
for testing.<br>
<br>
[B]<br>
This is simpler, but for reasons that will be apparent to<br>
everybody in a few applications, it is rarely used.<br>
<br>
1) Postprocess the restart with PP and write the electron density<br>
on a fine 3d-grid (cube-file). The size of 10 pm associated<br>
to a density energy-cutoff of 250 Ry is rather accurate.<br>
<br>
2) Process this electron density cube-file with a real-space<br>
QT-AIM code. I used:<br>
<br>
Sanville et al., J. Comp. Chem., 28, 899 (2007).<br>
<a href="http://theory.cm.utexas.edu/bader/" target="_blank">http://theory.cm.utexas.edu/bader/</a><br>
<br>
Charge are in ACF.dat.<br>
This works also for periodic systems (maybe only orthorhombic?).<br>
Problems occur with polar X-H bonds or in all cases<br>
where the zero-flux of density comes too close to atoms desribed<br>
with pseudo-potentials.<br>
<br>
Maybe, this will help,<br>
<br>
Giovanni<br>
<br>
============================================================<br>
Giovanni La Penna - National research council (Cnr)<br>
Institute for chemistry of organo-metallic compounds (Iccom)<br>
via Madonna del Piano 10,<br>
I-50019 Sesto Fiorentino, Firenze, Italy<br>
tel.: +39 055 522-5264, fax: +39 055 522-5203<br>
e-mail: <a href="mailto:glapenna@iccom.cnr.it">glapenna@iccom.cnr.it</a> - <a href="http://www.iccom.cnr.it/lapenna" target="_blank">http://www.iccom.cnr.it/lapenna</a><br>
skype: giovannilapenna<br>
============================================================<br>
<br>
On Wed, 1 Dec 2010, David Grifith wrote:<br>
<br>
> Dear All<br>
><br>
><br>
> I am trying to calculate the point charge and the point spin density on each<br>
> atom of a unit cell, certainly "WITHOUT" considering a 3D vector which is<br>
> ...<br>
<br>
</blockquote><br></div><br><br clear="all"><br><br><br>