[Q-e-developers] Fwd: Makov-Payne correction
Nicola Marzari
nicola.marzari at materials.ox.ac.uk
Sat Dec 18 14:05:05 CET 2010
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
> > >
> >
> >
> >
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20101218/a748daa0/attachment.html>
More information about the developers
mailing list