Closed edwinb-ai closed 6 months ago
Have you checked the positions? Are you sure you don't have particles overlapping? Maybe they are overlapping at the boundaries.
To generate the initial coordinates you might be interested in this function:
https://m3g.github.io/Packmol.jl/stable/#Packmol.pack_monoatomic!
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.
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.
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:
virial
field in the EnergyandForces
type to accumulate the value.reset_output!
function.reducer
function.Thanks for the help, @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).
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 hereI 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.