Dear LAMMPS community,
I am new to LAMMPS and given that it has a steep learning curve, there’s still a lot of it that I don’t know how to use, so please bear with me. I’ve been reading the documentations, but still unsure how to proceed.
I want to add a force to every particle in my 3D simulation that is dependent on the particle’s X and Z position. Eventually, I want to base these forces on a complicated function that I represent in a table that I interpolate from. However, for now, I’m taking a step back and trying more simple functions.
I attempted to use the pair_style table command, though from what I can gather, it only accepts radii as input along with energy and force, while what I need to do is substitute radii for X and Z positions. I have since concluded that this is most likely not the way to go — correct?
Currently, I am using Python-style variables that feed into the fix addforce command. I was able to get a uniform (not spacially dependent) configuration going, which leads me to believe that I at least have the syntax correct; however, I am unable to modify it to get the X and Z components of the force for every particle to be dependent on their position.
At every timestep, I want to access the vector of X positions and the vector of Z positions for every particle. If I can feed these vectors to Python, I should be able to add the appropriate force to each particle. Is there a way to do this that you can point me to?
One thing that I have been doing is the code below that loops through every particle but it is not working. Any ideas you would have would be extremely helpful.
label looplabel
variable loops loop 600
variable zforce python get_zforce
python get_zforce input 2 v_loops v_fieldforce return v_zforce format iff here “”"
def get_zforce(loops, fieldforce):
This function will return fieldforce if the corresponding x position is between 1/3 and 2/3 of Lxy, otherwise it will return zero
My aim is to have zforce be equal to fieldforce if and only if the x position of the particle is between 1/3 and 2/3 of Lxy
import numpy as np
num_particles = 600
Lxy = 950e-9
iter = 0
x = np.zeros(num_particles)
with open(“incl-initial-position.txt”) as open_file:
for line in open_file:
iter += 1
spl = str.split(line, " ")
if(iter > 4 and iter < 604):
x[iter-5] = spl[2]
if(x[loops-1] > Lxy/3 and x[loops-1] < 2*Lxy/3):
return fieldforce
else:
return 0
“”"
print ${loops}
print ${zforce}
I want to accumulate all zforce values here, and then pass that onto fix addforce in place of {zforce} in the “fix efield all addforce 0.0 0.0 {zforce}” command
next loops
jump SELF looplabel
Single body potentials
fix electrode all wall/colloid zlo EDGE {wcHamaker} {LJdiameter} ${WallCut} units box
fix efield all addforce 0.0 0.0 ${zforce}
Sincerely,
Haesung