# Average squared energy & co

Dear all,

I need to compute the average of the SQUARED energy
(Kinetic+Potential) over a certain interval of time in a dynamic
simulations:
<E_total^2>
I have seen how to compute the average of the kinetic and potential
energy independently.
But I couldn't find on the manual how to define variables on their sum
and squares.
Is it possible to define a quantity
E_tot = kinetic+potential
and then
E_square = E_tot*E_tot
and then use ave and reduce
?
What precise voice of the manual should I look at ?
How is the syntax ?
Any suggestions ?

In addition, how can I print the average of the distance between all
couples of atoms <r_ij> ?
And its square <r_ij^2> ?
Is it possible with ave and reduce ?
I have found just the command for the displacement from its original coordinates
compute displace/atom, but that is not what I need...

Bests,
Daniele Scopece

The first part is easy: You can define variables that do this, e.g.:

variable U equal pe

variable K equal ke

variable E equal “v_U + v_K”

variable E2 equal “v_E*v_E”

etc…

These can then be averaged with fix ave/time.

The second part is more difficult, and I am not sure if it is possible without modifications. There is a compute pair/local which gives distances between atoms involved in pairwise interactions, which you can dump or compute with dump/local and compute/local. Constructing the full NxN matrix of distances is more difficult though, and to me it seems that this is way easier in post-processing.

Compute reduce can take the output of compute pair/local (for distance)

and average it to a single value. You could also use fix ave/histo with

output from pair/local to calc a histogram of distances.

Neither can do those operations on the square of distances, and

there is no variable option that works for “local” values. I guess

either compute pair/local or compute reduce would need a new option

to create a squared quantity.

Steve

Keep in mind your r and r^2 descriptors for point 2 may turn out to scale with system size. At least that is the case of a homogeneous system (constant number density).

= \int r * P® dr = \int r * 4\pi*r^2/V_tot dr = 4\pi/V_tot * \int r^3 dr
using V_tot ~ r_max^3

one gets for ~ r_max a linear dependency with r_max being the maximum integration radius.

Carlos

PS: you may want to check my “fast food” approach to the calculation

Added a sumsq and avesq option to compute reduce.

Will be in the next patch.

Steve