Dear all I want to apply a force to a group of atoms in this way:
F_add = constant*(x-x0-vt)
where
x is the position of the particle
x0 is the position of the particle at the initial time
v is a constant velocity
t is the time from when the simulation is started
The system is periodic in x,y direction and I want to apply the force in x direction to every particles of the group. The problem is that when the particles cross the boundary then the box length are subtracted, but for my purpose this work wrong because velocity and time increasing during the simulation. How can I take into account how many time a particles cross to a boundary and then adding this quantity for computing the force?
In other word I want the rewrite the previous formula as
F_add = constant*(x+boxlength*count-x0-vt)
where count is the number of time that a particles has crossed the boundary in a given direction.
Dear all I want to apply a force to a group of atoms in this way:
F_add = constant*(x-x0-vt)
where
x is the position of the particle
x0 is the position of the particle at the initial time
v is a constant velocity
t is the time from when the simulation is started
The system is periodic in x,y direction and I want to apply the force in x
direction to every particles of the group. The problem is that when the
particles cross the boundary then the box length are subtracted, but for my
purpose this work wrong because velocity and time increasing during the
simulation. How can I take into account how many time a particles cross to
a boundary and then adding this quantity for computing the force?
just use "fix smd" with the "cvel" option in "tether" mode.
it does exactly what you want. you do give a reference point
to which the pulling force is pulling towards, but that is only
used to get the direction vector for the pulling velocity. that
should then go on through periodic boundaries without a
problem. remember to set r0 to 0.0 and forget about it.
people get confused about what it means.
do you think that apply a force to a cdm of a group (which I treat as rigid) is the same that apply a force to each atoms?
If I want that initially the spring is in equilibrium it’s needed that RX must be equal to the center of the mass?
Dear Axel, however there is a problem with the periodic system: if apply a velocity of 0.1 the system start to move from left to right and when the cdm reach the border of the box the boxlength are subtracted ok?
In order to avoid this problem and allow to smd to treating periodic system I think is it possible to do something like this in tether style
if(((xcm[0]-xcm_old[0])+boxlength/2.0) < 0) counter +=1 ;
xcm_old[0] is the position of the center of mass at the previous timestep
counter is an integer that take into account how many times the system goes through the boundary (how don’t know if LAMMPS just provide some function for this check)
Then I think to update dr with the formula
dr = r + counter*boxlength - r0 - r_old;
instead
dr = r - r0 - r_old;
in order to take into account the crossing through the boundary condition.
Another option is to use compute displace/atom which
will calculate the displacement of each atom from its
initial position including periodic box effects. Then
use this compute in an atom-style variable to calculate
the force you want to apply to each atom, i.e. the per-atom
vector produced by the compute is the x-x0 in your formula.
Then, use the atom-style variable as the fx variable (e.g.
v_myforce) in the fix addforce command.
Ok but I want apply the force to a rigid body and on the documentation is written:
IMPORTANT NOTE: If an atom is part of a rigid body (see the fix rigid command), it’s periodic image flags are altered, and the computed displacement may not reflect its true displacement. See the fix rigid command for details. Thus, to compute the displacement of rigid bodies as they cross periodic boundaries, you will need to post-process a dump file containing coordinates of the atoms in the bodies.
Of course I can’t post-process a dump file for computing the force
Dear Axel, however there is a problem with the periodic system: if apply a
velocity of 0.1 the system start to move from left to right and when the cdm
reach the border of the box the boxlength are subtracted ok?
since the fix smd module only applies a force vector, there is no
problem like you describe it. if you have proof that it doesn't work
like i claim, please provide an input that reproduces it, as that would
mean that there is a bug in the code.
In order to avoid this problem and allow to smd to treating periodic system
I think is it possible to do something like this in tether style
if(((xcm[0]-xcm_old[0])+boxlength/2.0) < 0) counter +=1 ;
xcm_old[0] is the position of the center of mass at the previous timestep
counter is an integer that take into account how many times the system goes
through the boundary (how don't know if LAMMPS just provide some function
for this check)
the implementation in fix smd is different and should not need this.
it automatically applied the minimum image convention when computing
the center of mass of the pull group with respect to the reference point,
which is displaced by the velocity vector in every step (without any regard
of the PBC).
Here it’s the example, if you go to log.lammps line 739/740 you see from the first column that cdm cross the boundary and if you keep the 5th column you can see a discontinuity in the force.
Maybe there’s something wrong in my script but I’ve checked many time.
yes, you might have problems with what I suggested
for rigid bodies. The right way to do this for rigid bodies
is probably to apply the force to the center-of-mass which
would require writing code in another fix that accessed
the COM data stored by fix rigid.
there is one important issue with your script.
you don't really pull. your spring constant is too
small while the pulling velocity is too high,
so that the atoms in fix spring actually are not
really moving. you have to compare column 5
and 6 of the fix smd output.
Here it's the example, if you go to log.lammps line 739/740 you see from the
first column that cdm cross the boundary and if you keep the 5th column you
can see a discontinuity in the force.
Maybe there's something wrong in my script but I've checked many time.
i think there are two overlapping issues. the one i mentioned
above, but if i adjust the fix smd parameters to act more reasonable,
there is indeed a bad jump when the pulled group reaches half the
box length. i'll have to dig into the code to see what is happening here.
there is one important issue with your script.
you don’t really pull. your spring constant is too
small while the pulling velocity is too high,
so that the atoms in fix spring actually are not
really moving. you have to compare column 5
and 6 of the fix smd output.
I agree with you, but currently I want to study the behavior of the top atoms as a function of the parameters.
Here it’s the example, if you go to log.lammps line 739/740 you see from the
first column that cdm cross the boundary and if you keep the 5th column you
can see a discontinuity in the force.
Maybe there’s something wrong in my script but I’ve checked many time.
i think there are two overlapping issues. the one i mentioned
above, but if i adjust the fix smd parameters to act more reasonable,
there is indeed a bad jump when the pulled group reaches half the
box length. i’ll have to dig into the code to see what is happening here.
For some parameters it’s possible to have stick-slip behavior, which parameters have you used?
there is one important issue with your script.
you don't really pull. your spring constant is too
small while the pulling velocity is too high,
so that the atoms in fix spring actually are not
really moving. you have to compare column 5
and 6 of the fix smd output.
I agree with you, but currently I want to study the behavior of the top
atoms as a function of the parameters.
but pulling is pointless, if there is no pulling happening.
> Here it's the example, if you go to log.lammps line 739/740 you see from
> the
> first column that cdm cross the boundary and if you keep the 5th column
> you
> can see a discontinuity in the force.
> Maybe there's something wrong in my script but I've checked many time.
i think there are two overlapping issues. the one i mentioned
above, but if i adjust the fix smd parameters to act more reasonable,
there is indeed a bad jump when the pulled group reaches half the
box length. i'll have to dig into the code to see what is happening here.
For some parameters it's possible to have stick-slip behavior, which
parameters have you used?
BTW: you cannot use the center of mass as tether point,
but you have to choose a point away from it.
the main problem is, that fix smd was designed for cases
where you pull a small group of atoms around to induce
structural changes. enforcing a consistent center for a periodic
layer is outside of the scope of what it was designed for.
so i have to look and find a way to allow this option without
losing any existing (normal) functionality.
there is one important issue with your script.
you don’t really pull. your spring constant is too
small while the pulling velocity is too high,
so that the atoms in fix spring actually are not
really moving. you have to compare column 5
and 6 of the fix smd output.
I agree with you, but currently I want to study the behavior of the top
atoms as a function of the parameters.
but pulling is pointless, if there is no pulling happening.
Here it’s the example, if you go to log.lammps line 739/740 you see from
the
first column that cdm cross the boundary and if you keep the 5th column
you
can see a discontinuity in the force.
Maybe there’s something wrong in my script but I’ve checked many time.
i think there are two overlapping issues. the one i mentioned
above, but if i adjust the fix smd parameters to act more reasonable,
there is indeed a bad jump when the pulled group reaches half the
box length. i’ll have to dig into the code to see what is happening here.
For some parameters it’s possible to have stick-slip behavior, which
parameters have you used?
BTW: you cannot use the center of mass as tether point,
but you have to choose a point away from it.
the main problem is, that fix smd was designed for cases
where you pull a small group of atoms around to induce
structural changes. enforcing a consistent center for a periodic
layer is outside of the scope of what it was designed for.
so i have to look and find a way to allow this option without
losing any existing (normal) functionality.
In fact I think the problem is that actually dr is computed as
dr = r - r0 - r_old;
but this does not take into account the fact that cdm can cross the boundary. So I think it’s possible to introduce a counter and modifiy dr in the following way:
dr = r + counter*boxlength - r0 - r_old;
then it’s needed to define a way in order to increase(decrease) the counter, for example like this
In fact I think the problem is that actually dr is computed as
dr = r - r0 - r_old;
but this does not take into account the fact that cdm can cross the
boundary. So I think it's possible to introduce a counter and modifiy dr in
the following way:
dr = r + counter*boxlength - r0 - r_old;
then it's needed to define a way in order to increase(decrease) the counter,
for example like this
adding boxlength/2.0 is needed because in order to avoid an accidental
increasing of the counter due to small oscillations.
you are correct. i must have mixed up the code in the current
repository with a rewrite that i started quite a while ago and never
found the time to finish. the actual solution is more complex
as it has to support non orthogonal cells.
you may be better off to write a custom fix for your purposes
based on fix spring that just moves the reference point for
fix spring with a constant velocity along a given vector and just
keeps adding to it, computes dx, dy, dz, relative to the c.o.m.
of the fix group (which can be a single atom, btw) and then
call domain->minimum_image(dx, dy, dz). that or something
similar should do the trick.
you may be better off to write a custom fix for your purposes
based on fix spring that just moves the reference point for
fix spring with a constant velocity along a given vector and just
keeps adding to it, computes dx, dy, dz, relative to the c.o.m.
of the fix group (which can be a single atom, btw) and then
call domain->minimum_image(dx, dy, dz). that or something
similar should do the trick.
i think i have figured it out. there are actually two calls
to domain->minimum_image required. however, this modification
doesn't fit well into the fix smd code and its purpose.
so i have (quickly) written a modification of fix spring
called fix sprint/pull that only works in tether mode and
does continuous pulling of a tethered object into a given
direction. this adds the pulling vector to the reference
point. than applied minimum image to i, then computed the
deltas to the pull group c.o.m. and applies minimum image
again. this way, i am able to have multiple recurrences
of the "top" group in your example and there are no problems
due to fix rigid/nvt manipulating the image flags in
unexpected ways. i am attaching the add-on sources and
a modified input with a few comments and a reference run log.
the fix has 8 data elements. the first 4 are the same as fix
spring (fx,fy,fz ftotal) and the you have the distance between
the tether point and the fix group's c.o.m. and the position of
the reference point.
let me know, if you have any questions about this.
Dear Axel, thanks for the effort; now all seems ok. I try to reproduce some results obtained with a simple MD code.
I have a few question:
-in you in-ak input file what is the mean of group pull id 220?
-in the other code the force are computed in the order: pair interaction -> external force -> langevin force. I do not know if LAMMPS use the same order but do you think that is possible to find different results if the order of the operation is not the same?
Dear Axel, thanks for the effort; now all seems ok. I try to reproduce some
results obtained with a simple MD code.
I have a few question:
-in you in-ak input file what is the mean of group pull id 220?
it is not used anywhere. the line was left over from some previous
attempts to avoid c.o.m. problems by pulling only a single atom
in the rigid layer. if you delete that line, it should make no difference.
-in the other code the force are computed in the order: pair interaction ->
external force -> langevin force. I do not know if LAMMPS use the same order
but do you think that is possible to find different results if the order of
the operation is not the same?
the external force and langevin parts are both fixes in lammps and
you can determine the order in which they are executed by the order
in which you define them. since the external force will add energy
to the system, the order should make a (small?) difference.
Dear Axel, thanks for the effort; now all seems ok. I try to reproduce some
results obtained with a simple MD code.
I have a few question:
-in you in-ak input file what is the mean of group pull id 220?
it is not used anywhere. the line was left over from some previous
attempts to avoid c.o.m. problems by pulling only a single atom
in the rigid layer. if you delete that line, it should make no difference.
-in the other code the force are computed in the order: pair interaction ->
external force -> langevin force. I do not know if LAMMPS use the same order
but do you think that is possible to find different results if the order of
the operation is not the same?
the external force and langevin parts are both fixes in lammps and
you can determine the order in which they are executed by the order
in which you define them. since the external force will add energy
to the system, the order should make a (small?) difference.
actually there is about a factor 2 in the average of the force but the langevin thermostat is implemented by using a gaussian distribution instead of a uniform as in LAMMPS.