Comments (5)
Thanks for reporting this, I'm not sure I've ever used 3 interactions on the GPU and it looks like that's the issue. One workaround might be to keep the loop but loop 2:N
where N
is the length of the tuple in the type domain. Using sum
in some arrangement would work but may have performance implications.
I'll have a more thorough look next week.
from molly.jl.
Thank you for looking into it! I would like to let you know that the reason for the reappearance of the error after modifying the force evaluation was because of the energy calculation (probably because of the same reason). On src/cuda.jl
lines 211-216, I found the same structure as in the force calculation:
pe = potential_energy(inters[1], dr, coord_i, coord_j, atoms[i], atoms[j],
boundary, special)
for inter in inters[2:end]
pe += potential_energy(inter, dr, coord_i, coord_j, atoms[i], atoms[j],
boundary, special)
end
Making the changes on both the force and energy calculations now makes the errors completely disappear for any combination of three or more pairwise potentials.
Regarding the loop over 2:N, I tried but found the error reappeared. I will also try to figure out the reason.
from molly.jl.
In particular this occurs when different interactions are used, using three of the same interactions is okay.
This works:
f = force_gpu(inters[1], dr, coord_i, coord_j, atoms[i], atoms[j], boundary, special)
f += force_gpu(inters[2], dr, coord_i, coord_j, atoms[i], atoms[j], boundary, special)
f += force_gpu(inters[3], dr, coord_i, coord_j, atoms[i], atoms[j], boundary, special)
But this errors:
f = force_gpu(inters[1], dr, coord_i, coord_j, atoms[i], atoms[j], boundary, special)
for inter_i in 2:3
f += force_gpu(inters[inter_i], dr, coord_i, coord_j, atoms[i], atoms[j], boundary, special)
end
Extracting it out to a function also errors.
A solution would be to auto-generate the code at the top. I tried metaprogramming for this but couldn't get it to work, I'm not the strongest at that though. I got Base.Cartesian.@nexprs
working with a fixed n
but couldn't make it a variable in the GPU kernel.
One workaround is to define functions like:
function addforces(inters::Tuple{<:Any, <:Any, <:Any}, dr, coord_i, coord_j, atom_i, atom_j, boundary, special)
return force_gpu(inters[1], dr, coord_i, coord_j, atom_i, atom_j, boundary, special) +
force_gpu(inters[2], dr, coord_i, coord_j, atom_i, atom_j, boundary, special) +
force_gpu(inters[3], dr, coord_i, coord_j, atom_i, atom_j, boundary, special)
end
This works but would require defining multiple functions. Maybe the function definitions could be written with metaprogramming.
I wonder if @vchuravy has any ideas as to why the first example in this comment works in a CUDA kernel and the second doesn't, or if there is an easy workaround.
MWE, for reference:
using Molly, CUDA
boundary = CubicBoundary(1.0u"nm")
coords = CuArray(place_atoms(100, boundary; min_dist=0.1u"nm"))
atoms = CuArray([Atom(σ=0.02u"nm", ϵ=0.1u"kJ * mol^-1") for _ in 1:100])
nf = DistanceNeighborFinder(eligible=CuArray(trues(100, 100)), dist_cutoff=0.2u"nm")
lj = LennardJones(use_neighbors=true)
coul = Coulomb(use_neighbors=true)
ss = SoftSphere(use_neighbors=true)
sys2 = System(coords=coords, atoms=atoms, boundary=boundary, neighbor_finder=nf, pairwise_inters=(lj, coul,))
sys3 = System(coords=coords, atoms=atoms, boundary=boundary, neighbor_finder=nf, pairwise_inters=(lj, coul, ss))
neighbors = find_neighbors(sys2)
forces(sys2, neighbors) # Works
forces(sys3, neighbors) # Errors
from molly.jl.
What is typeof(inters)
.
But the likely solution is something like:
f += sum(ntuple(Val(N)) do inter_i
@inline force_gpu(inters[inter_i], dr, coord_i, coord_j, atoms[i], atoms[j], boundary, special)
end)
Important is that N
is a compile time constant. Since inters
is a tupleN=length(inters)
may work. Essentially we are forcing the compiler to unroll this code statically.
from molly.jl.
Thanks Valentin. Indeed the force case seems to work with permutations of the above and length(inters)
seems to be available at compile time. This works without slowdown, and with Enzyme:
f_tuple = ntuple(length(inters)) do inter_type_i
force_gpu(inters[inter_type_i], dr, coord_i, coord_j, atom_i, atom_j, boundary, special)
end
f = sum(f_tuple)
typeof(inters)
is a Tuple
of 3 different concrete structs:
Tuple{
LennardJones{false, NoCutoff, Int64, Int64, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), 𝐋 𝐌 𝐍^-1 𝐓^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), 𝐋^2 𝐌 𝐍^-1 𝐓^-2, nothing}},
Coulomb{NoCutoff, Int64, Quantity{Float64, 𝐋^3 𝐌 𝐍^-1 𝐓^-2, Unitful.FreeUnits{(kJ, nm, mol^-1), 𝐋^3 𝐌 𝐍^-1 𝐓^-2, nothing}}, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), 𝐋 𝐌 𝐍^-1 𝐓^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), 𝐋^2 𝐌 𝐍^-1 𝐓^-2, nothing}},
SoftSphere{false, NoCutoff, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), 𝐋 𝐌 𝐍^-1 𝐓^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), 𝐋^2 𝐌 𝐍^-1 𝐓^-2, nothing}}
}
Strangely, though, the potential energy case
pe_tuple = ntuple(length(inters)) do inter_type_i
potential_energy(inters[inter_type_i], dr, coord_i, coord_j, atom_i, atom_j, boundary, special)
end
pe = sum(pe_tuple)
works for tuples of length 2 but fails on tuples of length 3 with a similar error:
Reason: unsupported call to an unknown function (call to ijl_get_nth_field_checked)
potential_energy
returns a "simpler" value than force_gpu
, a single value (both unitful and unitless cases fail) compared to a SVector
, so I was surprised that it didn't work but the force did. Perhaps it is because potential_energy
can be applied to more argument types, so it is struggling to infer the return type? I did try annotating the return type, adding an additional function boundary, using @inline
and using Val
but that did not seem to work.
Any idea what is going on there?
from molly.jl.
Related Issues (20)
- Links in Docs are broken HOT 2
- Boltzmann Constant in LAMMPS 'real' units does not work HOT 5
- AtomsBase not properly implemented HOT 4
- Inconsistent Units Crash Simulation HOT 5
- apply_coupling! method is not found for custom coupling function HOT 2
- Example for new MC membrane barostat HOT 12
- How to cite? HOT 6
- [feature request] velocity-dependent forces HOT 2
- Should we pass more properties (e.g. velocities, step number) to `force`/`potential_energy`? HOT 4
- Should we have functions to add and remove atoms to/from a `System`? HOT 3
- open an sdf file HOT 10
- Molly.jl: AD with Enyzme returns 0 gradient HOT 1
- Units in Nose Hoover temperature HOT 1
- GPU error with the example HOT 4
- return type for force for pairwise potential in CUDA HOT 1
- sys.coords does not work as Vector{Vector} HOT 2
- EAM implementation and dimension error HOT 5
- Boundary Conditions HOT 1
- UndefVarError: `σ6` not defined ERROR when simulating the NPT ensemble HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from molly.jl.