I am trying to stretch a cylindrical nanopillar of a metal along the cylinder axis. I manged somehow to stretch the pillar using the fixes below and balance the forces in two slabs in the top and bottom of the pillar along the cylinder axis (or y-axis). But I do not know how to prevent the rotation of the pillar due to subsequent activated slip systems in the pillar. Any suggestion is greatly appreciated, since I need to prevent the rotation of the pillar as in the experiment.
Thanks in advance,
boundary s s s
fix 1 UpperFixed store/force
fix 2 LowerFixed store/force
compute ufyc UpperFixed reduce ave f_1
compute lfyc LowerFixed reduce ave f_2
variable avg_velocityY_up equal vcm(UpperFixed,y)
variable avg_velocityY_low equal vcm(LowerFixed,y)
variable vel equal $((10e8*(ly-100)/1e12/2))
variable vy_init_upper equal v_vel-v_avg_velocityY_up
variable vy_init_lower equal -v_vel-v_avg_velocityY_low
Define Velocity on Boundary
velocity UpperFixed set 0 v_vy_init_upper 0 sum yes units box
velocity LowerFixed set 0 v_vy_init_lower 0 sum yes units box
variable appforceY1 equal -c_ufyc
variable appforceY2 equal -c_lfyc
Define Boundary Conditions
fix 3 UpperFixed addforce 0 v_appforceY1 0
fix 4 LowerFixed addforce 0 v_appforceY2 0
fix 5 NVT nvt temp 300 300 $(100.0 * dt) tchain 5 tloop 100 nreset 1000
fix_modify 5 temp TempNVT
fix 6 NVE nve
fix 7 all dt/reset 1 NULL NULL 0.02 units box
fix 8 all balance 2500 1.0 shift xyz 10 1.0 weight time 0.8
fix 9 all momentum 1000 linear 1 1 1
run 100000 post no
For what it’s worth, when I ran a similar simulation, I used “fix addforce” to apply equal and opposite forces to two different atoms on opposite sides of the object (in the +X and -X directions, for example). This would prevent the object from spinning around the Y axis (assuming the applied force was strong enough), while allowing it to move along the Y axis. This works very well. However the object’s orientation will fluctuate a little bit around the Y axis due to thermal motion. If you cannot tolerate any rotational motion at all, I suppose you could try using “fix momentum” with the “angular” keyword (although this will prevent all rotation, not only rotation around the Y axis). I have not tried that. (If somebody has an alternative suggestion, perhaps they will chime in.)
However, I advocate using the simplest approach whenever possible. I like using “fix addforce” because it is numerically stable (for reasonable force magnitudes) and the results are easy to interpret.
I should probably terminate my message here, because I am unsure about what I am about to say in the remaining portion of this message. Here are some other comments about your input file:
I vaguely recall that the use of “fix recenter” can cause numerical instability if the atoms you are applying the fix to are touching other atoms which are not part of that group. Except in very special circumstances, I’ve had trouble using “fix recenter”. Fortunately, you are using “fix momentum”, not “fix recenter”, so I don’t think this will cause numerical stability issues. (I am just always nervous about using fixes that change the way the equations of motion are integrated.)
I don’t know if it’s safe to thermostat some of the particles using fix nvt (Nose-Hoover thermostat), and the remaining particles with fix nve. It’s probably okay, but you might want to check this before proceeding. I just don’t know enough about how the Nose-Hoover thermostat works. But I do know that it’s definitely safe to do this with fix langevin. (IE. You could apply fix langevin to some of the particles, and also apply fix nve to all of them, since fix langevin works together with fix nve.)
Again, I don’t know if there is a problem with either 1) or 2) (ie fix momentum, fix nvt + fix nve), but if you are having numerical instability, these issues may be worth looking at. (Again, in my experience, the use of fix addforce is comparatively safe.)
Thanks a lot for the feedback!
Since I am manipulating the forces in two slabs only at the top and bottom of the pillar, I am not sure balancing the forces along other directions, i.e., x and z, can generate moments. So, if moments Mx and Mz are going to be generated, then rotations should be there as well. Am I wrong, you think?
Using fix momentum with angular keyword may work, but I do not want to mess up with the atoms’ trajectories further, especially the ones not located in the boundary layers. Using a periodic boundary condition along the cylinder can resolve this issue of rotation, but I do not know how exactly it affects the atoms’ trajectories along the periodic direction. Also, having periodic boundary condition can provide us with the luxury of using the “fix deform” in that direction, but again, I am not sure interfering with the atoms’ trajectories by using “remap x” keyword how will affect the final results. LAMMPS developers certainly know these stuff much better, if this seems to be nonsense or has negligible effects.
Per Axel on the LAMMPS forum, the system may suffer from flying ice cube syndrome because of using fix recenter. He may have other comments, if I am wrong.
Per Axel again, since we have free surfaces everywhere, and no periodic boundary condition exists to trap the heat in any direction, thermostating the whole system may not make so much sense here (he may have different opinion if I am quoting him wrong here). So, thermostating only some parts of the system, for example in the boundary layers or slabs, can help thermostating the whole pillar, of course using NVE integration in non-boundary atoms. As you said, any kind of thermostating should work fine for this purpose.