[Pw_forum] question on potential energy in constant pressure simulations

[Sandro Scandolo]

Dear Kostya,

an additional note on the stress tensor, perhaps trivial for you, but 
hopefully useful for other less experienced users:

the stress tensor calculated in the Q-E codes is based on the expression 
derived by Nielsen and Martin (PRL 50, 697, 1983). Nielsen and Martin's 
derivation is based on the Hellmann-Feynman theorem, and therefore 
assumes that the basis set is complete. It can be shown that the 
derivation is valid also when the basis set is not complete, as long as 
the derivative is taken in such a way that the basis set does not 
change, to first order, upon infinitesimal changes of the strain tensor.

As a consequence, you should expect the calculated stress tensor to 
coincide with the derivative of the energy with respect to the strain, 
as calculated by small finite strains, *only* if you are able to keep 
the same number of plane waves in the calculations with and without 
strain. For isotropic strains (i.e. volume changes, where V->V+dV) this 
is easily done by changing the energy cut-off by [(V+dV)/V]^2/3 . For 
anisotropic strains (i.e. for the off-diagonal elements of the stress 
tensor) it may be more tricky. This is the correct procedure for testing 
if the stress tensor is calculated correctly by the code. However, as 
you probably know, when the basis set is not complete, a better (in 
physical terms) estimate of the stress tensor is obtained by 
differentiating the total energy (or potential energy in your notation) 
calculated by keeping the wfc cut-off (ecutwfc) fixed in the finite 
difference calculation. Unfortunately this leads to discontinous 
pressure-volume  (or stress-strain) curves, as the basis set size 
changes discontinuosly upon strain changes, if the cut-off is fixed (NB: 
to be more precise, it leads to a continuous curve only if the k-point 
sampling is infinitely dense, which is never the case in practice, and 
certainly not with Gamma-point only).

The "qcutz" trick has been introduced as a partial relief of the above 
problem. By introducing a (smooth) penalty function for plane waves 
whose kinetic energy exceeds a given *fixed* cut off (ecfixed), you can 
forget about the need to keep the cut-off fixed in the construction of 
the basis set (ecutwfc), as  long as ecutwfc is above the tail of the 
smooth penalty function. All this is described in M. Bernasconi et al, 
J. Phys. Chem. Solids 56, 501 (1995). As a strightforward consequence, 
the stress tensor calculated with this trick should be much less 
sensitive to whether it is calculated at constant cut-off or at constant 
number of plane waves. Hence, Nielsen and Martin's expression, with this 
modified functional, gives you a stress tensor in much better agreement 
with the physically sounder expression obtained by keeping the cut-off 
fixed.  However, you should not expect it to coincide exactly with the 
derivative of the potential energy when the derivative is calculated by 
keeping the cut-off fixed. You should however expect it to coincide 
exactly with the derivative taken at constant number of plane waves 
(like for the stress tensor without "qcutz").

Regards,
Sandro


Konstantin Kudin wrote:

> Hi Nicola,
>
> Thanks for the info! I was under the impression that ALL bugs were
>fixed by now, but apparently - not. I do in fact use ultra-soft PSPs,
>so that would explain my results.
>
> An important consideration for me is that I did not intend just to
>test the stress, but rather to use it for some production jobs. So
>unfortunately it looks like this has to wait.
>
> Kostya
>
>
>--- Nicola Marzari <marzari at MIT.EDU> wrote:
>
>  
>
>>Dear Kostya,
>>
>>great that your are testing the constant-pressure cp.
>>
>>There is still an unresolved issue (bug...) for ultrasoft psp, in
>>the diagonal terms only - that could explain your results. On the
>>bright
>>side, a lot of other arcane issues (bugs) that had crept in the CVS
>>have been recently found  and fixed  by Carlo C. and Paolo G.
>>
>>The safest first step would be to look at the numerical derivative of
>>the ground-state energy with respect to your variable lattice
>>vector - should be equal with the calculated stress
>>(again, if you are using psp, it will not...).
>>
>>
>>
>>			nicola
>>
>>
>>Konstantin Kudin wrote:
>>    
>>
>>> Hi all,
>>>
>>> I am doing the CP dynamics with one of the lattice vectors being
>>>optimized as well.
>>>
>>> I am using a very recent CVS (yesterday's), which should have the
>>>correct stress (knocking on wood).
>>>
>>> The parameters for the modified kinetic energy functional are:
>>>ecutwfc = 30.0, ecfixed = 25.0, qcutz = 25.0, q2sigma = 3.0
>>> This is for what used to be just ecutwfc=25.0 in constant volume
>>>simulations.
>>>
>>> So what I am seeing is that while the stress becomes smaller, the
>>>potential energy of the system goes up at the same time with the
>>>decreasing stress.
>>>
>>> Should not the decrease in stress correspond to the *decreasing*
>>>potential energy? Or is there some funny business with this
>>>      
>>>
>>modified
>>    
>>
>>>kinetic energy functional such that I should disregard the
>>>      
>>>
>>potential
>>    
>>
>>>energy values and trust the stress instead ?
>>>
>>> Thanks!
>>> Kostya
>>>
>>>
>>>__________________________________________________
>>>Do You Yahoo!?
>>>Tired of spam?  Yahoo! Mail has the best spam protection around 
>>>http://mail.yahoo.com 
>>>_______________________________________________
>>>Pw_forum mailing list
>>>Pw_forum at pwscf.org
>>>http://www.democritos.it/mailman/listinfo/pw_forum
>>>      
>>>
>>-- 
>>---------------------------------------------------------------------
>>Prof Nicola Marzari   Department of Materials Science and Engineering
>>13-5066   MIT   77 Massachusetts Avenue   Cambridge MA 02139-4307 USA
>>tel 617.4522758 fax 2586534 marzari at mit.edu http://quasiamore.mit.edu
>>_______________________________________________
>>Pw_forum mailing list
>>Pw_forum at pwscf.org
>>http://www.democritos.it/mailman/listinfo/pw_forum
>>
>>    
>>
>
>
>__________________________________________________
>Do You Yahoo!?
>Tired of spam?  Yahoo! Mail has the best spam protection around 
>http://mail.yahoo.com 
>_______________________________________________
>Pw_forum mailing list
>Pw_forum at pwscf.org
>http://www.democritos.it/mailman/listinfo/pw_forum
>  
>