Hi everybody,
I’m using LAMMPS to simulate the heating effects in Zhang’s paper, pages 4-5: Phys. Rev. A 76, 012719 (2007) - Molecular-dynamics simulations of cold single-species and multispecies ion ensembles in a linear Paul trap (aps.org).
The idea is that in a system of around 100 Barium ions, at every time step, I need to give a velocity kick to a random number among these 100 ions. I want to do this by applying a force to some ions via ‘fix addforce’, but this applies the force to the entire ion ensemble. Is there a way I can pick and choose among these ions to give kicks to? Thank you so much!
How about something like
variable r equal ceil(random(0,count(all),12345))
fix r ave/time all 1 1 1 v_r # makes RNG fire once per step, not per atom
variable chosen atom id==f_r
variable forcex atom v_chosen*v_my_force_variable
fix a addforce all v_forcex 0 0
If I have my LAMMPS right, this puts your desired force on a randomly chosen atom: the atom whose ID equals r at every step will have v_chosen == 1 (and it will be 0 for all other atoms).
Thank you so much! I want to add the force in all three dimensions. How do I get the vx, vy and vz of chosen? I’m thinking about something like this:
I have wx, wy, wz which are the forces in each direction that I will multiply to each v component.
variable forcex atom v_chosen_x * wx
variable forcey atom v_chosen_y * wy
etc.
Is that the right way to go? Or is this not the right syntax to get vx and vy? Thank you!!
lines = [‘\n# Define ion-neutral heating for a species…’,
f’group {uid} type {iid}‘,
f’variable r equal ceil(random(0,count(all),12345))’,
f’fix {uid} {iid} ave/time 1 1 1 v_r’,
f’variable chosen atom id==f_r’,
f’variable forcex atom v_chosen_x*{valuewx:e}',
f’variable forcey atom v_chosen_y{valuewy:e}',
f’variable forcez atom v_chosen_z{value*wz:e}‘,
f’fix {uid} {iid} addforce all v_forcex v_forcey v_forcez\n’]
Here’s my current code (I’m using a python wrapper for LAMMPS). It’s saying that I’m “Replacing a fix, but new style != old style (src/modify.cpp:856).” How should I fix it? Thank you so much
It sounds like you have another fix that has the same name as one of the names I randomly chose. Please check the documentation for further details.
Also, please make sure you understand those few lines that I have written (including asking questions about what you don’t understand). One of the standard recommendations for posting here is that you should always start with a simple system and work your way up to something more complex. So just verify first that what I wrote even works for an x-force before trying to do y, z, etc.
Thank you! I don’t understand the second and last lines of your code, where you wrote fix r and fix a. I was looking at the LAMMPS syntax and it says that it should be “fix ID group-ID” etc, and there isn’t a place for a and r. The code didn’t work at first so I added some uid and group id to it, and the replacing a fix error pops up. Is this because there are 2 fix commands in the code?
There absolutely is. Those are the fix IDs.
OH! So the ID is actually the name we give the fix command! I finally understand it. Thank you so much! LAMMPS says that it doesn’t expect “all” so I deleted “all” from the second and last lines and I finally got it working.
Would you mind helping me set it up in 3 dimensions? Is there a way for me to get the x, y and z components of v_chosen? Thank you!
You can try something like
variable r equal ceil(random(0,count(all),12345))
fix r ave/time all 1 1 1 v_r # makes RNG fire once per step, not per atom
variable chosen atom id==f_r
variable force_x atom v_chosen*v_some_value_x
variable force_y atom v_chosen*v_some_value_y
variable force_z atom v_chosen*v_some_value_z
fix a addforce all v_force_x v_force_y v_force_z
Please try to understand what this is doing before using it in production. In particular, try adding the following lines and running for a few steps to understand what exactly the code is doing:
# diagnostics
fix printout all print 1 "$(f_r)"
dump printrandoms all custom 1 printrandoms.txt id v_chosen v_force_x v_force_y v_force_z
Thank you so much! I really appreciate you taking the time to help me!
I got it working for all 3 directions, and the force I apply is very small, around 1e-20, but for some reason the velocity for each ion in the ensemble jumps really high after only 1 timestep! I’m thinking about applying the force every 500 timesteps, so I’m rewriting the code to be
fix r ave/time all 500 1 1000 v_r
while keeping the rest the same, but now there’s an error that says Fix in variable not computed at a compatible time. The LAMMPS documentation says this means that the fix command is requesting a value on a non-allowed timestep, but I assume that the 1000th timestep should be allowed. I changed Nevery and Nfreq to other values and still got the same error. Would you mind pointing out what I might have done wrong? Thank you!
I tried to put 500 in the addforce command too
fix a addforce all v_force_x v_force_y v_force_z every 500
and the jump is still very big!
How do you know the force is very small?
Have you confirmed that the force you apply is actually as small as you think it is?
I don’t know what’s going on either, especially since the problem is very likely to come from either somewhere else in your code or something in the Python wrapper you’re using, neither of which we have access to or we are familiar with.
I think I got the code to work now! Thanks so much for providing me with the original code!