Value of variable ID equal bound(all, zmax/min) jumps after 25 consecutive 10000 timestep runs

LAMMPS LAMMPS Mailing List Mirror
1

Dear all,

I’m trying to run mechanical annealing simulations in Lammps where I repeatedly make and break contact between two nano-sized metallic fragments. The idea is to simulate what happens during the scanning tunneling microscopy STM or mechanically controllable break junction MCBJ experiments conducted by the experimentalists in our research group. This has been successfully implemented in MODAS, a fortran based simulator. But, for instance, Lammps would give us the ability to put molecules between our metallic fragments, which is not currently possible with MODAS.

So, my simulations start with a neck-shaped metallic nano-contact that I slowly rupture into two fragments and then bring back into contact. I do this by freezing a few layers of atoms at the top and bottom of the initial neck structure, defining these as regions that extend to positive and negative infinity, respectively, which I allow to move with the move command and a “fix move variable” command that does the actual moving. The layers in between the frozen top and bottom layers are allowed to stretch and break and make contact again when moving the frozen layers back toward each other. This is equivalent to the way it is done in MODAS.

So, by allowing the structure to rupture I can guess how many timesteps I would need before reversing the motion to bring the ruptured fragments back into contact. This works fine in Lammps.

My problem arises when I want to repeat this process a few times, in order to simulate mechanical annealing. Ideally I want to make contact only to a certain depth just as is possible in STM and MCBJ, by checking say, after every 10000 timesteps how many atoms are in the narrowest cross-section of my structure when I bring two ruptured fragments back into contact. This is conveniently possible in Lammps through the fix ave/spatial command and so far has worked fine for me.

However, I realised that I had to use the “start stop” capability of the run command to ensure the structure keeps stretching or shortening after every 10000 timesteps. This works fine during the 25 cycles (of 10000 timesteps each) of stretching using a “variable a loop 25” loop, but on the second timestep (I use thermo 100) of the next loop “variable b loop 25” that brings the fragments back into contact, the inner boundaries of the top and bottom frozen layers (calculated using variable ID bound(all,zmax or zmin)) jump drastically in value (see below), as well as the potential and total energy (taken from log.lammps):

variable z1 delete # variables to define the displacement of the top (z1) and bottom (z2) frozen layers
variable z2 delete

variable z1 equal -0.0408 # I reverse these values from the initial rupture si
variable z2 equal 0.0408

undump file_xyz # I want separate rupture/formation trajectories for postprocessing
dump file_xyz all xyz 5000 pelicula2.xyz
dump_modify file_xyz format " %s %11.4E %11.4E %11.4E "

fix move_up push_up move variable NULL NULL v_dz1 NULL NULL v_z1 units box
fix move_dw push_dw move variable NULL NULL v_dz2 NULL NULL v_z2 units box
fix atoms_layer all ave/spatial 10000 1 10000 z center 2.04 ave one units box file atom.layer2

label loopb
variable b loop 25
run 10000 start 0 stop 500000
run 10000 start 0 stop 500000
Memory usage per processor = 4.23268 Mbytes
Step Temp temp_cen KinEng PotEng TotEng zlim_up_ zlim_dw_ z1 z2
160000 2.6471179 4.2 0.17929543 -1794.3314 -1794.1521 22.848 -20.808 -0.0408 0.0408
160100 2.6471179 4.2 0.17929543 -1634.4222 -1634.2429 29.38008 -27.34008 -0.0408 0.0408
160200 2.6471179 4.2 0.17929543 -1638.0523 -1637.873 29.38416 -27.34416 -0.0408 0.0408
160300 2.6471179 4.2 0.17929543 -1640.7751 -1640.5958 29.38824 -27.34824 -0.0408 0.0408
160400 2.6471179 4.2 0.17929543 -1642.903 -1642.7237 29.39232 -27.35232 -0.0408 0.0408
160500 2.6471179 4.2 0.17929543 -1644.6279 -1644.4486 29.3964 -27.3564 -0.0408 0.0408
160600 2.6471179 4.2 0.17929543 -1646.0522 -1645.8729 29.40048 -27.36048 -0.0408 0.0408
160700 2.6471179 4.2 0.17929543 -1647.2422 -1647.0629 29.40456 -27.36456 -0.0408 0.0408
160800 2.6471179 4.2 0.17929543 -1648.2616 -1648.0823 29.40864 -27.36864 -0.0408 0.0408
160900 2.6471179 4.2 0.17929543 -1649.174 -1648.9947 29.41272 -27.37272 -0.0408 0.0408
161000 2.6471179 4.2 0.17929543 -1650.0355 -1649.8562 29.4168 -27.3768 -0.0408 0.0408

I initially unfixed and refixed, but this gave worse problems, so the fixes above for contact fomation are identical to the ones higher up in my input script (follows below).

I would appreciate any help or suggestions, as I believe I’m close to getting this right and might just be missing something small.

INPUT SCRIPT:

LAMMPS - Tip on tip - lattice 001 Z axis.

Definition of units

units metal
atom_style atomic

Boundary conditions: shrink-wrapped in all three dims

boundary s s s

FCC Gold lattice oriented along the 001 crystal plane

lattice fcc 4.08 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1

Data file containing coordinates of initial structure

read_data input525.dat

Number of time steps (each step = 1fs or 0.001 ps)

variable num_steps_close equal 10000

Variables to define displacement of structure (velocity of separation and approach).

variable z1 equal 0.0408
variable z2 equal -0.0408
variable dz1 equal vdisplace(0.0,{z1}) variable dz2 equal vdisplace(0.0,{z2})

Mass of the type of atoms used

mass 1 196.96

Type of interaction between atoms (applied potential).

pair_style eam/alloy
pair_coeff * * Au_zhou_alloy.eam Au

Neighbour list (update).

neighbor 0.85 bin
neigh_modify delay 10

Variables to define the fixed layers of the structure.

variable zlim_up_push equal (bound(all,zmax)-4.08)
variable zlim_dw_push equal (bound(all,zmin)+4.08)

Groups and regions needed to define frozen and moving parts of the structure.

region upper block INF INF INF INF {zlim_up_push} INF move NULL NULL v_z1 units box region lower block INF INF INF INF INF {zlim_dw_push} move NULL NULL v_z2 units box

group push_up region upper
group push_dw region lower
group center subtract all push_dw push_up
group extremes union push_dw push_up

Establecimiento de la temperatura (en K y aleatoria). Solo al grupo central.

Los extremos permanecen inmóviles.

##velocity center create 4.20000001 15

Non-standard calculations

compute temp_center center temp

Thermodynamic output (temperature, energy, etc.).

thermo_style custom step temp c_temp_center ke pe etotal v_zlim_up_push v_zlim_dw_push v_z1 v_z2
thermo 100

Loop to eventually achieve repeated cycles of contact rupture and formation

##variable num_cycles equal 5
##label folderd
##variable d loop ${num_cycles} pad
##shell mkdir run$d
##shell cd run$d

Periodically generate images of the structure to check that simulation is proceeding correctly.

dump images all image 10000 image-*.jpg type type adiam 2 &
box no 0.0 view 90 180 zoom 2
dump_modify images pad 7 acolor 1 gold

Write trajectory to xyz file for post-processing.

dump file_xyz all xyz 5000 pelicula1.xyz
dump_modify file_xyz format " %s %11.4E %11.4E %11.4E "

Fixes to define motion of all the groups of atoms in the structure.

#fix temp_control center nvt temp 4.20000001 4.20000001 10.0
fix temp_control center temp/rescale 10 4.20000001 4.20000001 0.05 1.0
fix temp center nve
fix move_up push_up move variable NULL NULL v_dz1 NULL NULL v_z1 units box
fix move_dw push_dw move variable NULL NULL v_dz2 NULL NULL v_z2 units box
fix freeze extremes setforce 0.0 0.0 0.0
fix atoms_layer all ave/spatial 10000 1 10000 z center 2.04 ave one units box file atom.layer1

Variables needed to calculate number of atoms at narrowest cross-section of the structure

variable minimum equal min(f_atoms_layer[2])
variable pos equal length(f_atoms_layer[2])
variable penpos equal “v_pos-1”
variable firstbin equal f_atoms_layer[1][2]
variable lastbin equal f_atoms_layer[v_pos][2]
variable secondbin equal f_atoms_layer[2][2]
variable secondlast equal f_atoms_layer[v_penpos][2]

Rupture run

label loopa
variable a loop 25
run {num_steps_close} start 0 stop 500000 print "The minimum cross-secton is {minimum}"
if “({minimum} < 1) && ({firstbin} > 0) && ({secondbin} > 0) && ({secondlast} > 0) && (${lastbin} > 0)” then “jump SELF break1”
next a
jump SELF loopa
label break1

variable z1 delete
variable z2 delete

variable z1 equal -0.0408
variable z2 equal 0.0408

undump file_xyz
dump file_xyz all xyz 5000 pelicula2.xyz
dump_modify file_xyz format " %s %11.4E %11.4E %11.4E "

fix move_up push_up move variable NULL NULL v_dz1 NULL NULL v_z1 units box
fix move_dw push_dw move variable NULL NULL v_dz2 NULL NULL v_z2 units box
fix freeze extremes setforce 0.0 0.0 0.0
fix atoms_layer all ave/spatial 10000 1 10000 z center 2.04 ave one units box file atom.layer2

Formation run

label loopb
variable b loop 25
run {num_steps_close} start 0 stop 500000 print "The minimum cross-secton is {minimum}"
if “({minimum} > 5) && ({firstbin} > 10) && (${lastbin} > 10)” then “jump SELF break2”
next b
jump SELF loopb
label break2

#unfix move_up
#unfix move_dw
#unfix temp_control
#unfix temp
#unfix freeze

#undump file_xyz
#undump images

#shell cd …

#next d
#jump in.auneck folderd

2

unless you can reduces this complex input to something much simpler
and smaller, that still shows the undesired behavior, chances are very
slim that somebody will look into your problem and try to determine
whether something is incorrectly programmed in LAMMPS or between your
ears. :wink:

axel.

3

Dear Dr Kohlmeyer,

As always, thank you for your prompt reply.

I didn't mean to imply that there was something wrong with Lammps's programming, so I apologise if that's the way my question came across.

And I'm pretty convinced that something is incorrectly programmed between my ears with respect to this problem, and that what I'm trying to do is entirely possible with Lammps.

I had already spent a considerable amount of time simplifying the input script I included in my original question, which originally involved jumping back and forth between two separate scripts for contact rupture and formation, including using write_restart and read_restart commands. Obviously, the complexity of such a script made debugging a nightmare. So, I decided to start over.

Based on the script in my original question in the mailing list, the next simplifying step would be to remove the loops over shorter consecutive runs and just to run one long rupture simulation (takes around 1,600,000 time steps for the neck to break) and another long contact formation sim (takes around 2,000,000). Lammps transitions seamlessly between two such longer runs. And I don't even need to use the "start stop" option of the run command to accomplish this.

I also originally ran an energy equilibration 100,000 steps long, which converged rapidly to a stable energy. I further considered minimization between the shorter consecutive runs, and CG minimization usually needed only one step to converge. This proved to be incompatible with the Nfreq defined in fix ave/spatial, and so I had to abandon that. But because minimization happened after only one step, I didn't think it was the problem.

So, if you'll forgive the long email, I would be very grateful for any further suggestions and/or comments you might have for me.

Thank you very much in advance.

Wynand

4

Dear Dr Kohlmeyer,

As always, thank you for your prompt reply.

I didn't mean to imply that there was something wrong with Lammps's programming, so I apologise if that's the way my question came across.

neither did i have the impression that you wanted to imply this. you
misunderstood my statement. when something doesn't behave as expected,
there is *always* the option that this is a bug, especially when using
a feature that is not routinely used by a lot of people. so i was
simply making a neutral assessment. it has to be one of the two.

And I'm pretty convinced that something is incorrectly programmed between my ears with respect to this problem, and that what I'm trying to do is entirely possible with Lammps.

I had already spent a considerable amount of time simplifying the input script I included in my original question, which originally involved jumping back and forth between two separate scripts for contact rupture and formation, including using write_restart and read_restart commands. Obviously, the complexity of such a script made debugging a nightmare. So, I decided to start over.

Based on the script in my original question in the mailing list, the next simplifying step would be to remove the loops over shorter consecutive runs and just to run one long rupture simulation (takes around 1,600,000 time steps for the neck to break) and another long contact formation sim (takes around 2,000,000). Lammps transitions seamlessly between two such longer runs. And I don't even need to use the "start stop" option of the run command to accomplish this.

that is not the kind of simplification that i had in mind. you are
still trying to do a meaningful calculation. that is often not needed.
for something like you describe, the first thing to do is a very
careful examination and visualization of the trajectory. there may be
something usual happening right when the "jump" happens. a second
option would be to try reproducing the issue from a prerecorded
trajectory in a dump file using the rerun command.

I also originally ran an energy equilibration 100,000 steps long, which converged rapidly to a stable energy. I further considered minimization between the shorter consecutive runs, and CG minimization usually needed only one step to converge. This proved to be incompatible with the Nfreq defined in fix ave/spatial, and so I had to abandon that. But because minimization happened after only one step, I didn't think it was the problem.

So, if you'll forgive the long email, I would be very grateful for any further suggestions and/or comments you might have for me.

the most important thing would be to identify a reliable trigger for
the issue and make some simple schematic test, e.g. using fix move.

5

Dear Axel,

Thanks for your additional comments and suggestions. I will approach the problem the way you recommend.

By looking at the images I generate every 5000 time steps, I've already picked up that the frozen layers I define at the beginning of my input script revert back to their initial position at the moment of the jump (i.e., the top and bottom of the input structure -- see the attached images for t = 0.0 and 230,000.0 fs before the jump and t = 240,000.0 fs after it).

I will try to see if I can remedy this and let you and the mailing list know how it goes.

Wynand

image-0000000.jpg

image-0230000.jpg

image-0240000.jpg

6

Hi Axel,

I managed to solve my most immediate problem. In the end it was something small incorrectly programmed between my ears :wink:

I guess you wouldn't remember with all the questions posted on the mailing list, but I run nano-wire rupture and contact formation simulations. I've been trying to repeat consecutive cycles of rupture and contact formation to simulate mechanical annealing and ran into a snag that has taken much time and effort to isolate.

The problem was that I was using the start stop capability of the run command incorrectly. I applied the same start stop values over an entire cycle of rupture and contact formation when I should have just reset the time step to zero after rupture, and commenced contact formation from where the rupture simulation left off. I write separate trajectories for these processes anyway (for post-processing purposes) and don't need to start the formation simulation at the time step the rupture simulation has ended. So now I essentially run two different simulations from my input script, each with their own start stop values --> problem solved!

I knew it was a small problem, but my complex script made it hard to catch.

Regards
Wynand