Computing virial and pressure

[edwinb-ai]

Hello, thanks for the great package!

This is my second attempt at building something with this package.

I have a very simple code for doing Brownian Dynamics using a continuous pseudo hard-sphere potential in 3D.
This is the code in a gist

However, when the energy and virial get updated in the energy_and_forces! function, which is summarized here

function energy_and_forces!(x,y,i,j,d2,cutoff,output::EnergyAndForces)
    d = sqrt(d2)
    r = y - x
    if d < cutoff
        (uij, fij) = pseudohs(d)
        sums = @. fij * r / d
        output.virial += dot(sums, r)
    end
    output.energy += uij
    output.forces[i] = @. output.forces[i] + (fij * r / d)
    output.forces[j] = @. output.forces[j] - (fij * r / d)

    return output
end

I always get very large values for energies and the virial. I have been trying to understand what the problem is here, but I don't know where to look into further. Any help is extremely appreciated.

[lmiq]

Have you checked the positions? Are you sure you don't have particles overlapping?
Maybe they are overlapping at the boundaries.

[lmiq]

To generate the initial coordinates you might be interested in this function:

https://m3g.github.io/Packmol.jl/stable/#Packmol.pack_monoatomic!

[edwinb-ai]

Hello, thanks for the quick response. I have checked my initial configuration and no overlaps exist. Furthermore, the energy seems to be accumulated throughout the integration steps and keeps increasing constantly and large values of the energy are obtained.

On the virial computation, it seems that I get the expected values inside the energy_and_forces! function but outside the function, after calling map_pairwise!, the variable that accumulates the computed values get reset to zero. I wonder if this is the effect of calling the reset_output! function after each integration step.

[edwinb-ai]

The problem of the energy calculation was solved by removing an else case I had in the definition of my interaction potential function, so it went from

function pseudohs(rij; lambda=50.0)
    b_param = lambda / (lambda - 1.0)
    a_param = lambda * b_param^(lambda - 1.0)
    ktemp = 1.0

    if rij < b_param
        uij = (a_param / ktemp) * ((1.0 / rij)^lambda - (1.0 / rij)^(lambda - 1.0))
        uij += (1.0 / ktemp)
        fij = lambda * (1.0 / rij)^(lambda + 1.0)
        fij -= (lambda - 1.0) * (1.0 / rij)^lambda
        fij *= -a_param / ktemp
    else
        uij = 0.0
        fij = 0.0
    end

    return uij, fij
end

to the following

function pseudohs(rij; lambda=50.0)
    b_param = lambda / (lambda - 1.0)
    a_param = lambda * b_param^(lambda - 1.0)
    ktemp = 1.0
    uij = 0.0
    fij = 0.0

    if rij < b_param
        uij = (a_param / ktemp) * ((1.0 / rij)^lambda - (1.0 / rij)^(lambda - 1.0))
        uij += (1.0 / ktemp)
        fij = lambda * (1.0 / rij)^(lambda + 1.0)
        fij -= (lambda - 1.0) * (1.0 / rij)^lambda
        fij *= -a_param / ktemp
    end

    return uij, fij
end

By removing the else case, the energy computed stopped increasing throughout the simulation. However, the computation of the virial has the same behavior, it keeps increasing.

[edwinb-ai]

I have fixed the issue, and updated the code in the gist in case anyone find this useful.

In the end, the fixes were the following:

• Include a new virial field in the EnergyandForces type to accumulate the value.
• Reset the value of this field in the reset_output! function.
• Include the same behavior of summing the values of the energy as for the virial in the reducer function.
• The interaction potential case in the definition, which doesn't include a case of setting the values to zero. (I don't understand completely why this would be a problem, but it worked)

Thanks for the help, @lmiq !

[lmiq]

It seems reasonable that you needed to update the energy and the viral in the reducer/reset functions as well.

I don't think the else there plays any role the solution of the problem, maybe that change got mixed with the other ones?

A small note:

Here:

    if d < cutoff
        (uij, fij) = pseudohs(d)
        sumies = @. fij * r / d

if the cutoff used is the same as the system.cutoff (which appears to be), the conditional is redundant. That function will only be called for pairs of particles within that cutoff. An additional conditional there only is useful if the cutoff for the particle interactions is smaller than the actual cutoff used to build the cell lists (which might happen, for example, if multiple time-stepping is used).