velocity setting with airebo potential

Dear all,

I am doing friction simulation setting the velocity of 0.5 (angstoms/picosecond) in x direction to simulate velocity of 50m/s.

However, the result showed atoms mixed up with increase of volume and decrease of pressure.

After I increase the velocity to 500, I could get a good result of friction simulation.

This is confusing because 500 in metal unit means 50000m/s.

I have tried…

  1. Very low timestep (0.00001, 0.000001)

  2. Nvt instead of nve and temp/rescale

  3. lower neighbor 2.0

  4. lower cutoff of airebo (2.0, 1.0)

None of them worked.

Can anybody give me a hint to drive a good result of simulation with velocity of 0.5 (50m/s) in this simulation?

Below is modified in.friction example file which is from lammps package. (red colored are critical changes)

Thank you,

Sungae Lee

2d friction simulation

dimension 2

boundary p s p

units metal

atom_style full

neighbor 10 bin

neigh_modify delay 5

create geometry

lattice hex 3

region box block 0 50 0 30 -0.25 0.25 units box

create_box 4 box

mass 1 12.01

mass 2 12.01

mass 3 12.01

mass 4 12.01

atom regions

region lo-fixed block INF INF INF 2 INF INF units box

region lo-slab block INF INF INF 10 INF INF units box

region above-lo block INF INF INF 10 INF INF side out units box

region hi-fixed block INF INF 25 INF INF INF units box

region hi-slab block INF INF 17 INF INF INF units box

region below-hi block INF INF 17 INF INF INF side out units box

region lo-asperity sphere 32 10 0 5 units box

region hi-asperity sphere 18 17 0 5 units box

region lo-half-sphere intersect 2 lo-asperity above-lo units box

region hi-half-sphere intersect 2 hi-asperity below-hi units box

region outboxes intersect 2 above-lo below-hi units box

region asperities union 2 lo-half-sphere hi-half-sphere units box

create 2 surfaces with asperities

create_atoms 1 region lo-slab

create_atoms 2 region hi-slab

create_atoms 3 region lo-half-sphere

create_atoms 4 region hi-half-sphere

LJ potentials

pair_style airebo 3.0 0 0

pair_coeff * * CH.airebo C C C C

define groups

group lo region lo-slab

group lo type 1 3

group hi region hi-slab

group hi type 2 4

group lo-fixed region lo-fixed

group hi-fixed region hi-fixed

group boundary union lo-fixed hi-fixed

group mobile subtract all boundary

group outboxes region outboxes

group asperities region asperities

group air subtract outboxes asperities

group no-air subtract all air

#set group lo-fixed type 4

#set group hi-fixed type 4

initial velocities

compute new mobile temp/partial 0 1 0

velocity mobile create 300 482748 temp new units box

velocity hi set 0.5 0.0 0.0 sum yes units box

fixes

fix 1 all nve

fix 2 boundary setforce 0.0 0.0 0.0

fix 3 mobile temp/rescale 10 300 300 10 1 units box

fix_modify 3 temp new

Comments below.

Ray

Dear all,

I am doing friction simulation setting the velocity of 0.5
(angstoms/picosecond) in x direction to simulate velocity of 50m/s.

However, the result showed atoms mixed up with increase of volume and

Do you mean "results" after only 500 steps? If yes, 500 steps is next
to nothing, and these "results" are not useful at all.

decrease of pressure.

After I increase the velocity to 500, I could get a good result of friction
simulation.

This is confusing because 500 in metal unit means 50000m/s.

Indeed. A 50 km/s velocity is not friction, but ultrasonic shock.

I have tried..

1. Very low timestep (0.00001, 0.000001)

2. Nvt instead of nve and temp/rescale

3. lower neighbor 2.0

4. lower cutoff of airebo (2.0, 1.0)

Have you verified that your structure is correct in the first place?

None of them worked.

[...]

# 2d friction simulation

dimension 2

Are you sure you want to a 2D simulation?

boundary p s p

units metal

atom_style full

You don't need full for AIREBO.

neighbor 10 bin

This 10 is definitely inappropriately long.

neigh_modify delay 5

# create geometry

lattice hex 3

Are you sure the AIREBO potential will successfully describe this 2D
hex structure?

region box block 0 50 0 30 -0.25 0.25 units box

create_box 4 box

mass 1 12.01

mass 2 12.01

mass 3 12.01

mass 4 12.01

# atom regions

region lo-fixed block INF INF INF 2 INF INF units box

region lo-slab block INF INF INF 10 INF INF units box

region above-lo block INF INF INF 10 INF INF side out units box

region hi-fixed block INF INF 25 INF INF INF units box

region hi-slab block INF INF 17 INF INF INF units box

region below-hi block INF INF 17 INF INF INF side out units box

region lo-asperity sphere 32 10 0 5 units box

region hi-asperity sphere 18 17 0 5 units box

region lo-half-sphere intersect 2 lo-asperity above-lo units box

region hi-half-sphere intersect 2 hi-asperity below-hi units box

region outboxes intersect 2 above-lo below-hi units box

region asperities union 2 lo-half-sphere hi-half-sphere units box

# create 2 surfaces with asperities

create_atoms 1 region lo-slab

create_atoms 2 region hi-slab

create_atoms 3 region lo-half-sphere

create_atoms 4 region hi-half-sphere

# LJ potentials

pair_style airebo 3.0 0 0

With this setting, you are using REBO, not AIREBO.

pair_coeff * * CH.airebo C C C C

# define groups

group lo region lo-slab

group lo type 1 3

group hi region hi-slab

group hi type 2 4

group lo-fixed region lo-fixed

group hi-fixed region hi-fixed

group boundary union lo-fixed hi-fixed

group mobile subtract all boundary

group outboxes region outboxes

group asperities region asperities

group air subtract outboxes asperities

group no-air subtract all air

#set group lo-fixed type 4

#set group hi-fixed type 4

# initial velocities

compute new mobile temp/partial 0 1 0

velocity mobile create 300 482748 temp new units box

velocity hi set 0.5 0.0 0.0 sum yes units box

# fixes

fix 1 all nve

fix 2 boundary setforce 0.0 0.0 0.0

fix 3 mobile temp/rescale 10 300 300 10 1 units box

This is a bad thermostat.

fix_modify 3 temp new

#

timestep 0.0001

thermo 10

thermo_modify lost ignore flush yes

Why ignore lost atoms, which usually indicates a problem with your simulation?

dump 1 all atom 10 dump.frictionMM

run 500

Again, 500 steps is too short to obtain any meaningful results.

Ray,

Thank you very much for your comments.

I used 5000 steps instead of 500 steps for 50m/s. The number of run stemp is set for one time contact (1 stroke) between asperity.
"lost ignore" is to avoid error. I am not sure what to do when I get an error from lost atoms. This error message doesn't give me the details.

I was doing 2D simulation to reduce the debugging time because the result failed showing similar atomic movement with my main simulation. In my main simulation, I use diamond lattice with 3.57 lattice parameter in 3D. (hex is only one I could use for 2D, diamond was not available)
Below is the input script that I am using for the main simulation. I couldn't find any error in the script, but the result shows that atoms are mixed up from first 1000 run.

If you see anything wrong or doubt, could you make some comments on my input in the same way you did earlier?
Thanks a lot.

-------input----------------
units metal
dimension 3
boundary p p s

atom_style atomic
neighbor 2.5 bin
neigh_modify delay 5

#lattice
lattice diamond 3.57

read_data Data.BigDiatipEquim

region box block 0 50.2 -0.1 50 0 40 units box

# atom regions
region hi-asperity sphere 10 25 33 10 units box
region lo-fixed block 0 50.2 -1 50 0 5 units box
region lo-slab block 0 50.2 -1 50 0 20 units box
region lo-Tstat block 0 50.2 -1 50 5 12 units box
region above-lo block 0 50.2 -1 50 0 20 side out units box
region hi-fixed block 0 50.2 -1 50 35 40 units box
region hi-slab block 0 50.2 -1 50 33 40 units box
region hi-Tstat block 0 50.2 -1 50 33 35 units box
region below-hi block 0 50.2 -1 50 22 40 side out units box
region TA union 2 hi-asperity hi-slab units box

region outbox block 0 50.2 -1 50 0 40 side out units box

# LJ potentials

pair_style airebo 3 1 0
pair_coeff * * /lustre/work/apps/lammps-6Dec12/potentials/CH.airebo C C

# define groups

group lo region lo-slab
group lo type 1
group lo-Tstat region lo-Tstat
group hi-Tstat region hi-Tstat
group hi region TA
group hi type 2
group lo-fixed region lo-fixed
group hi-fixed region hi-fixed
group Tstat union lo-Tstat hi-Tstat
group boundary union lo-fixed hi-fixed
group mobile subtract all boundary

# initial velocities

compute folds lo coord/atom 1.6
compute new Tstat temp/partial 1 0 1
velocity mobile create 300 482748 temp new
velocity hi set 0.5 0.0 -0.4 sum yes

# fixes
fix 3 all nvt temp 300 300 100
fix 2 boundary setforce 0.0 0.0 0.0

# Run

timestep 0.00005
thermo 1000
thermo_style custom step atoms temp press vol enthalpy ke pe etotal lx ly lz
thermo_modify temp new
thermo_modify lost ignore flush yes

dump 1 all atom 1000 dump.BigDiatip
dump 2 all xyz 1000 xyz.BigDiatip
write_restart restart.BigDiatip

run 400000

Sungae,

Comments below.

Ray

Ray,

Thank you very much for your comments.

I used 5000 steps instead of 500 steps for 50m/s. The number of run stemp is set for one time contact (1 stroke) between asperity.

5000 steps is still too short. Have you first minimized,
equilibrated, and/or run NVE simulations of the system, without any
velocity set?

"lost ignore" is to avoid error. I am not sure what to do when I get an error from lost atoms. This error message doesn't give me the details.

Lost atoms error generally indicates bad structure, bad dynamics,
and/or bad potentials. It is not wise to ignore the error. Does this
error also occur when you run only NVE?

I was doing 2D simulation to reduce the debugging time because the result failed showing similar atomic movement with my main simulation. In my main simulation, I use diamond lattice with 3.57 lattice parameter in 3D. (hex is only one I could use for 2D, diamond was not available)
Below is the input script that I am using for the main simulation. I couldn't find any error in the script, but the result shows that atoms are mixed up from first 1000 run.

What do you mean by "mixed up"? If chaos occurs within the first 1000
steps, then something is definitely wrong, which can be the structure,
dynamics, or potentials. Fix the lost atoms error first. Visualize
your structure, characterize it (RDF, bond length, bond angles,
energies), equilibrate it, etc.

Ray