[Q-e-developers] Fwd: Makov-Payne correction
Stefano de Gironcoli
degironc at sissa.it
Sat Dec 18 13:40:49 CET 2010
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