[Q-e-developers] Fwd: Makov-Payne correction
Ismaïla Dabo
daboi at cermics.enpc.fr
Mon Dec 20 04:37:47 CET 2010
Dear all,
Thank you very much for the feedback.
I attach the code I used to calculate the results presented in Fig. 5,
hoping it is indeed the final version. Unfortunately, I don't have
numerical python libraries installed on my laptop. Consequently, I
cannot run computational verifications. Nevetheless, it's probable
that the results are actually in Ha and not in Ry, which should not
alter the comparison in this specific case.
As for the Eq. (21), I still believe the expression is correct "in the
case of a single point countercharge".
For corrections involving more countercharges, I agree with the result
nicely derived by Stefano de Gironcoli and I admit I didn't see the
connection with Eq. (13). Therefore, we should revert to the original
equation.
Now, I have to read the Makov-Payne paper to understand to which
situation Eq. (15) of that paper refers. If Makov and Payne indeed
refer to the multiple-countercharge case when writing their equation,
my statement at the top of the 2nd column of page 115139-5 is wrong.
In which case, I should write an erratum.
Could I ask that someone circulates the Makov-Payne paper
(unfortunately, I do not have access to PROLA in my current location)?
Many thanks,
Ismaila
On Sun, Dec 19, 2010 at 12:33 AM, Oliviero Andreussi
<oliviero.andreussi at materials.ox.ac.uk> wrote:
> 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
>>>>>>
>>>>>>
>>>>>
>>>
>>>
>
>
--
Ismaila DABO
http://cermics.enpc.fr/~daboi
-------------- next part --------------
A non-text attachment was scrubbed...
Name: graph.py
Type: application/octet-stream
Size: 4363 bytes
Desc: not available
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20101219/818c49f2/attachment.obj>
More information about the developers
mailing list