[Q-e-developers] recompute the total energy after scf is finished

[Paolo Giannozzi]

What is done and how and why is explained in some detail in the 
comments in electrons.f90:

        ! ... eband  = \sum_v \epsilon_v    is calculated by sum_band
        ! ... deband = - \sum_v <\psi_v | V_h + V_xc |\psi_v>
        ! ... eband + deband = \sum_v <\psi_v | T + Vion |\psi_v>
        !
        deband = delta_e()
	...
       CALL mix_rho ( rho, rhoin, mixing_beta, dr2, tr2_min, iter,&
			nmix, iunmix, conv_elec )
 	...
       IF ( .NOT. conv_elec ) THEN
	...
           CALL v_of_rho( rhoin, rho_core, rhog_core, &
                          ehart, etxc, vtxc, eth, etotefield, charge, v)
	...
           ! ... estimate correction needed to have variational energy:
           ! ... T + E_ion (eband + deband) are calculated in sum_band
           ! ... and delta_e using the output charge density rho;
           ! ... E_H (ehart) and E_xc (etxc) are calculated in v_of_rho
           ! ... above, using the mixed charge density rhoin%of_r.
           ! ... delta_escf corrects for this difference at first order
           !
           descf = delta_escf()
	...
	ELSE
	...
       ! ... note that rho is here the output, not mixed, charge density
       ! ... so correction for variational energy is no longer needed
           descf = 0._dp
	...
	ENDIF

P.

On Mon, 2013-11-04 at 12:13 +0100, Lorenzo Paulatto wrote:
> Dear all,
> I'm trying to recompute the total energy after the scf cycle is finished 
> recycling the wavefunctions that are stored to file at the end.
> 
> The worflow is like this:
> 1. call weights to recompute bands occupations, and eventually demet 
> (smearing contribution)
> 2. call sum_bands to recompute rho from the wfcs+occupations, this also 
> computes eband
> 3. call v_of_rho to recompute the potential, this also computes ehart, 
> etxc, etxcc, etc.
> 4. call deband = delta_e() = \int v_(r) rho_r(r) dr (function from 
> electrons)
> 5. recompute ewald term
> 6. etot = eband + deband + (etxc - etxcc) + ewld + ehart + demet
> 
> All the terms are identical to the SCF calculation, except for deband 
> which is slightly different.
> 
> I've noticed that the value deband computed during the scf cycle uses 
> the value of v(r) from the second-last iteration. I.e. if the 
> calculations converges in N iterations:
>   deband_N = \int v_{N-1}(r) \rho_N(r) dr
> 
> I would expect the two values to differ by something of the order
>   descf = \int (\rho_{N-1)-rho_N(r)) v_{N-1}(r) dr
> but it turns out that the difference is more of order sqrt(descf).
> 
>  From this I deduce that the form I use for the total energy may not be 
> the correct variational one. Finally, I'm confused. Anybody has the time 
> to help me find out what I may be doing wrong?
> 
> kind regards
> 
> P.S. I can provide the modified files (electrons.f90 + one new file) to 
> anybody on request, I'm not posting them here for size concerns.
> 

-- 
Paolo Giannozzi, Dept. Chemistry&Physics&Environment, 
Univ. Udine, via delle Scienze 208, 33100 Udine, Italy
Phone +39-0432-558216, fax +39-0432-558222