[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