Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Map forces onto beads #80

Open
wants to merge 13 commits into
base: main
Choose a base branch
from
Open

Conversation

StephMcCallum
Copy link

Adding comments to start project

Copy link

@erjank erjank left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Congrats on your first pull request! Here's to many more.

@chrisjonesBSU
Copy link
Member

@StephMcCallum here are some hints to for when you're working on this at the workshop.

  1. GRiTS works with gsd files (trajectories of simulations) and the python package is here

  2. You're mapping (i.e. averaging) properties of a group of atoms onto a single particle, you can see what is already done for mass and position around the lines where we left comments in this PR.

  3. You can think of a GSD file as a sequence of frames, where each frame is a snapshot of a simulation that has all the particle data. You can access frames like:

import gsd.hoomd

trajectory = gsd.hoomd.open("put-the-path-to-a-gsd-file-here.gsd")
first_frame = trajectory[0] # Using indexing to access a frame
last_frame = trajectory[-1]

for frame in trajectory: # Use a for loop to iterate through all of the frames
    particle_velocities = frame.particles.velocity
  1. Accessing particle velocities is shown in the example above. I think a good place to start is mapping the velocities to the coarse-grain bead. I know the ultimate goal is to map forces, but I think mapping velocities can be helpful too

  2. Accessing particle forces from a gsd file is a little trickier than the example shown above. So, start with velocities, and I'll leave another comment about how to do that. I'll also send you a gsd file that has the particle forces logged.

  3. Remember that PRs are a place to talk about the code you're trying to make, so if you have questions, this is a great place to leave them, rather than emailing back-and-forth, for example.

@marjanalbooyeh
Copy link
Contributor

marjanalbooyeh commented Jun 24, 2024

Hi @StephMcCallum,
Here's a little code snippet on how to access particle forces from a gsd file (number 5 in @chrisjonesBSU's comment above). I will share some gsd files that have the forces logged in them with you on Slack

Once you load a trajectory and iterate through the frames, you can check what kind of stuff is logged for each frame with:

print(frame.log.keys())

For example, for a gsd file that I'm currently looking at right now, the result looks like this:

dict_keys(['flowermd/base/simulation/Simulation/timestep',
'flowermd/base/simulation/Simulation/tps',  
'md/compute/ThermodynamicQuantities/kinetic_temperature', 
'md/compute/ThermodynamicQuantities/potential_energy',
'md/compute/ThermodynamicQuantities/kinetic_energy', 
'md/compute/ThermodynamicQuantities/volume', 
'md/compute/ThermodynamicQuantities/pressure', 
'md/compute/ThermodynamicQuantities/pressure_tensor', 
'md/pair/LJ/energy', 'particles/md/pair/LJ/forces', 
'particles/md/pair/LJ/energies', 'md/special_pair/LJ/energy', 'particles/md/special_pair/LJ/forces', 
'particles/md/special_pair/LJ/energies', 'md/bond/Harmonic/energy', 'particles/md/bond/Harmonic/forces', 
'particles/md/bond/Harmonic/energies', 'md/angle/Harmonic/energy', 'particles/md/angle/Harmonic/forces', 
'particles/md/angle/Harmonic/energies', 'md/dihedral/OPLS/energy', 'particles/md/dihedral/OPLS/forces', 
'particles/md/dihedral/OPLS/energies'])

The key names are a bit weird but they have all the info that you need. As you can see the per particle forces are stored based on the force type i.e. LJ pair forces, Harmonic bond forces, Harmonic angle forces , etc and their key name starts with particles/.

For example, I have a gsd file that I will share with you that contains 300 PPS monomers (each monomer has 7 atoms, therefore 2100 total particles) and the number of frames in this gsd file is 101. If you want to access the forces generated from the Lennard-Jones interactions for each particle in each frame of this trajectory, you can do something like this:

import numpy as np
traj_lj_forces = [] #initiate an empty list to store all frame forces
for frame in traj:
    lj_forces = frame.log['particles/md/pair/LJ/forces'] 
    traj_lj_forces.append(lj_forces)  #for each frame, append particle forces to the list
traj_lj_forces = np.asarray(traj_lj_forces) #convert the list to a numpy array, so that array operations could be done easily
print(traj_lj_forces.shape)

The traj_lj_forces numpy array now contains all LJ force vectors for this gsd file, which is of size (101, 2100, 3)
You can basically do the same for other force types.

bond forces:'particles/md/bond/Harmonic/forces'
angle forces: 'particles/md/angle/Harmonic/forces'
dihedral forces: 'particles/md/dihedral/OPLS/forces'

I hope this makes it more clear. please let me know if you have any questions :)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the additons so far look good. Now, you'll have to add the coarse-grained velocities to the new gsd frame that GRiTS is making.

You can see what is done for other properties around lines 780:

                new_snap.particles.N = len(typeid)
                new_snap.particles.position = position
                new_snap.particles.image = images
                new_snap.particles.typeid = typeid.astype(int)
                new_snap.particles.types = types
                new_snap.particles.mass = mass

So, you can define the new_snap.particles.velocity property with the velocity list being populated.

Copy link

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants