[Q-e-developers] Fwd: Makov-Payne correction
Stefano de Gironcoli
degironc at sissa.it
Sun Dec 19 01:19:48 CET 2010
Dear Nicola, dear all,
I've been looking into PRB 77 115139 (2008).
The formula we are discussing is Eq. (21), is n't it ?
I think Eq. 21 is not correct.. it is based on Eq. (11) for v0corr but
this is appropriate for a system with a single point-like counter charge
(no dipole, no quadrupole) so I think it is not consistent to then
include the quadrupole in the correction because if the quadruploe is
deemed important then the correcting potential is not Eq (11) but rather
Eq. (13) .
If you plug Eq. (13) into Eq. (20) then what you get is
Delta Ecorr = 1/2 * (alpha_0 q/L * q - 2pi q/3L^3 * Q + 4pi/L^3 p * p
- 2pi Q /3L^3 * q)
which by collecting similar terms becomes
Delta Ecorr = alpha_0 q^2/2L - 2 pi /3L^3 ( q * Q - p * p)
which is just the Makov-Payne correction as coded in PW and CP before
the last update... (there is of course an additional factor of e2 in PW).
I think we should revert to it.
Something must be wrong in fig.5 of PRB 77 115139 (2008) because if
we look at L=15 a.u. and use Q=153, q=2
then the difference between the red and the blue points should be
pi q Q/ 3 /L^3 = 3.14 * 2 * 153 /3/ 15^3 = 0.09489 Ha = 0.18978 Ry ...
which is clearly not the case.
stefano
Nicola Marzari wrote:
> ciao stefano,
>
> indeed, Oliviero is retesting in these days everything -
> when we checked it with Ismaila (fig 5 of the 2008 prb)
> it seemed correct, but we need to recheck units, and the
> role of the quadrupole.
>
> Soon we'll have an update.
>
> Nicola
>
> ----- Original message -----
>
>> dear Nicola, dear all,
>> admittedly I know nothing about your 2008 PRB with Ismaila and there
>> may be some problem of units or other glitch in th eimplementation ..
>> but on the simple tests that are contained in the cluster_example
>> directory (H20, NH4+ and atomic N) the new MP correction is clearly
>> inferior to the old one in the only case where it is not negligible, the
>> charged NH4+ system.
>>
>> stefano
>>
>> comparison follows between total energies of a few molecule in simple
>> cubic boxes with 16,20,24 a.u. sides
>> The reference results are those that can be obtained by running
>> run_example in cluster_example and correspond to Martyna Tuckeman
>> correction.
>>
>> NEW MP correction
>> results-mp_new/h2o.out-16:! Total+Makov-Payne energy =
>> -44.00269835 Ry
>> results-mp_new/h2o.out-20:! Total+Makov-Payne energy =
>> -44.00238516 Ry
>> results-mp_new/h2o.out-24:! Total+Makov-Payne energy =
>> -44.00256923 Ry
>> results-mp_new/nh4+.out-16:! Total+Makov-Payne energy =
>> -32.45635049 Ry
>> results-mp_new/nh4+.out-20:! Total+Makov-Payne energy =
>> -32.45333443 Ry
>> results-mp_new/nh4+.out-24:! Total+Makov-Payne energy =
>> -32.45292630 Ry
>> results-mp_new/n.out-16:! Total+Makov-Payne energy =
>> -27.82758316 Ry results-mp_new/n.out-20:! Total+Makov-Payne energy =
>> -27.82739183 Ry results-mp_new/n.out-24:! Total+Makov-Payne
>> energy = -27.82770125 Ry
>>
>> OLD MP correction
>> results-mp_old/h2o.out-16:! Total+Makov-Payne energy =
>> -44.00242684 Ry
>> results-mp_old/h2o.out-20:! Total+Makov-Payne energy =
>> -44.00224753 Ry
>> results-mp_old/h2o.out-24:! Total+Makov-Payne energy =
>> -44.00249007 Ry
>> results-mp_old/nh4+.out-16:! Total+Makov-Payne energy =
>> -32.45145525 Ry
>> results-mp_old/nh4+.out-20:! Total+Makov-Payne energy =
>> -32.45083466 Ry
>> results-mp_old/nh4+.out-24:! Total+Makov-Payne energy =
>> -32.45148127 Ry
>> results-mp_old/n.out-16:! Total+Makov-Payne energy =
>> -27.82758316 Ry results-mp_old/n.out-20:! Total+Makov-Payne energy =
>> -27.82739183 Ry results-mp_old/n.out-24:! Total+Makov-Payne
>> energy = -27.82770125 Ry
>>
>> Martyna Tuckerman approach
>> reference/h2o.out-16:! total energy = -44.00239103 Ry
>> reference/h2o.out-20:! total energy = -44.00224062 Ry
>> reference/h2o.out-24:! total energy = -44.00248682 Ry
>> reference/nh4+.out-16:! total energy = -32.45134332
>> Ry reference/nh4+.out-20:! total energy =
>> -32.45079847 Ry reference/nh4+.out-24:! total energy =
>> -32.45146630 Ry reference/n.out-16:! total energy =
>> -27.82757465 Ry reference/n.out-20:! total energy =
>> -27.82739133 Ry reference/n.out-24:! total energy =
>> -27.82770111 Ry
>>
>>
>>
>>
>> Nicola Marzari wrote:
>>
>>> Ciao Paolo,
>>>
>>>
>>> I spoke today with Oliviero Andreussi (that wrote the email
>>> below) and he was also going to do some tests to try out.
>>>
>>> Fig. 5 in the 2008 PRB Dabo at al was showing the discrepancy between
>>> the original M-P formula and our corrected one.
>>>
>>>
>>> nic
>>>
>>>
>>>
>>> On 12/3/10 7:04 PM, Paolo Giannozzi wrote:
>>>
>>>
>>>> It would be useful if anybody with previous Makov-Payne
>>>> experience could make some tests
>>>>
>>>> P.
>>>>
>>>>
>>>>
>>>>> Dear all,
>>>>>
>>>>> I checked the equations in the subroutine PW/makov-payne.f90 and it
>>>>> seems to me that the equation implemented is a modified version of
>>>>> the original Makov-Payne correction of 1995. The modification is
>>>>> due to the fact that in the original derivation the second
>>>>> correction was computed assuming that pointlike charge that
>>>>> compensate the molecular density was placed such as to give a zero
>>>>> net dipole moment, while the implemented version is more general
>>>>> and thus include a dipolar term. Moreover the correction has the
>>>>> correct sign (the opposite of the one in the original equation).
>>>>> This to say that the equation does not contain the factor one half
>>>>> suggested in the paper by Ismaila and Nicola.
>>>>> I would suggest the following changes to some of the lines of the
>>>>> subroutine:
>>>>>
>>>>> from
>>>>>
>>>>>
>>>>> !
>>>>> ! ... Makov-Payne correction, PRB 51, 4014 (1995)
>>>>> ! ... Note that Eq. 15 has the wrong sign for the quadrupole term
>>>>> !
>>>>> corr1 = - madelung(ibrav) / alat * qq**2
>>>>> !
>>>>> aa = quadrupole
>>>>> bb = dipole(1)**2 + dipole(2)**2 + dipole(3)**2
>>>>> !
>>>>> corr2 = ( 4.D0 / 3.D0 * pi )*( qq*aa - bb ) / alat**3
>>>>>
>>>>>
>>>>> to (I also made explicit the factor e2/2.0 to convert from Hartree
>>>>> to Rydberg, that was implicitly in the original implementation, so
>>>>> e2 should also be included in the declaration of the constants)
>>>>>
>>>>> !
>>>>> ! ... Makov-Payne correction, PRB 51, 4014 (1995)
>>>>> ! ... Note that Eq. 15 has the wrong sign for the quadrupole term
>>>>> ! ... The prefactor of the quadrupole term is half the one in Eq.
>>>>> 15 ! ... according to the results derived by Dabo et al., PRB 77,
>>>>> 115139 (2008)
>>>>> !
>>>>> corr1 = - madelung(ibrav) / alat * qq**2 / 2.0D0 * e2
>>>>> !
>>>>> aa = quadrupole
>>>>> bb = dipole(1)**2 + dipole(2)**2 + dipole(3)**2
>>>>> !
>>>>> corr2 = ( 1.D0 / 3.D0 * pi )*( qq*aa - bb ) / alat**3 * e2
>>>>>
>>>>> Similarly, for the routine CPV/makov_payne.f90
>>>>>
>>>>> !
>>>>> ! ... Makov-Payne correction, PRB 51, 43014 (1995)
>>>>> ! ... Note that Eq. 15 has the wrong sign for the quadrupole term
>>>>> !
>>>>> ! 1 / 2 Ry -> a.u.
>>>>> corr1 = - 2.8373D0 / alat * charge**2 / 2.0D0
>>>>> !
>>>>> aa = quadrupole
>>>>> bb = dipole(1)**2 + dipole(2)**2 + dipole(3)**2
>>>>> !
>>>>> ! 4 / 3 -> 2 / 3 Ry -> a.u.
>>>>> corr2 = ( 2.D0 / 3.D0 * pi )*( charge*aa - bb ) / alat**3
>>>>>
>>>>> to (in this case I assumed that e2 is defined in CP and with a
>>>>> value equal to 1.0,
>>>>> otherwise, probably is just better to forget about it)
>>>>>
>>>>> !
>>>>> ! ... Makov-Payne correction, PRB 51, 43014 (1995)
>>>>> ! ... Note that Eq. 15 has the wrong sign for the quadrupole term
>>>>> ! ... The prefactor of the quadrupole term is half the one in Eq.
>>>>> 15 ! ... according to the results derived by Dabo et al., PRB 77,
>>>>> 115139 (2008)
>>>>> !
>>>>> ! 1 / 2 Ry -> a.u.
>>>>> corr1 = - 2.8373D0 / alat * charge**2 / 2.0D0 * e2
>>>>> !
>>>>> aa = quadrupole
>>>>> bb = dipole(1)**2 + dipole(2)**2 + dipole(3)**2
>>>>> !
>>>>> ! 2/ 3 -> 1 / 3 Ry -> a.u.
>>>>> corr2 = ( 1.D0 / 3.D0 * pi )*( charge*aa - bb ) / alat**3 * e2
>>>>>
>>>>> Please, if someone else can also check these corrections before
>>>>> implementing them it would be better.
>>>>>
>>>>>
>>>> ---
>>>> Paolo Giannozzi, Democritos and University of Udine, Italy
>>>>
>>>>
>>>> _______________________________________________
>>>> Q-e-developers mailing list
>>>> Q-e-developers at qe-forge.org
>>>> http://qe-forge.org/mailman/listinfo/q-e-developers
>>>>
>>>>
>>>
>>>
>
>
>
More information about the developers
mailing list