search

JuliaMolSim / Molly.jl

Molecular simulation in Julia
Other
372 stars 51 forks source link

Fix shifted force cutoff potential calculation #53

Closed noeblassel closed 2 years ago

noeblassel commented 2 years ago

I believe there is a sign mistake in function potential_cutoff(cutoff::ShiftedForceCutoff, r2, inter, params) at cutoffs.jl at line 92, potential(inter, r2, invr2, params) - (r - rc) * fc - potential(inter, cutoff.sqdist_cutoff, cutoff.inv_sqdist_cutoff, params) should read potential(inter, r2, invr2, params) + (r - rc) * fc - potential(inter, cutoff.sqdist_cutoff, cutoff.inv_sqdist_cutoff, params). Indeed the correction should subtract the derivative of the potential, not the force. The (exagerated) effect can be seen on the following pictures (the "linear correction" curves):

Current version cutoffs_original_exagerated

Correct version cutoffs_exagerated

I also implemented a spline interpolating cutoff, I can open a PR for it if you are interested.

jgreener64 commented 2 years ago

Thanks for spotting this and describing it clearly, I pushed a fix.

Definitely interested in the spline cutoff PR if you have time.

I found your repo at https://github.com/noeblassel/MolecularDynamicsJulia and it looks interesting. Sorry if the Molly interface has been changing a lot recently, but I think the changes are better in the long term. In particular I am going to write all future simulators in a vectorised way so they work with GPU and AD, I just removed the in-place version of VelocityVerlet yesterday.

noeblassel commented 2 years ago

Thanks! I posted a pull request with the spline cutoff yesterday, didn't find the time to give a thorough explanation of the proposed changes yet.