[Q-e-developers] Fwd: Makov-Payne correction
Oliviero Andreussi
oliviero.andreussi at materials.ox.ac.uk
Sun Dec 19 11:33:58 CET 2010
Dear Stefano,
I reached the exact same conclusions(the missing factor one half is due
to the neglected quadrupole term in the periodic images), and I was
indeed testing what was wrong with Figure 5. I also noticed the wrong
order of magnitude in the quadrupole correction in the figure (it is
twice as large as it should be), but when I redid the calculations with
the correct factors it also seems that the correction with the factor
one half is better. I tend to think it is just an unlucky coincidence of
error cancellation, since in other analytical systems I have tried the
original makov-payne correction performs better, as it should be. I was
waiting to have a more significant statistics, but I would also say that
it is better to go back to the original formulation.
Best,
Oliviero
On 12/19/2010 12:19 AM, Stefano de Gironcoli wrote:
> 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