Has anyone used LAMMPS to do this?

LAMMPS LAMMPS Beginners
,
1

Hi,
I am trying to modify the lammps demo example of 2D ellipsoid to create a 3D system of ellipsoidal particles. The code stops saying “Lost atoms: original 1000 current 7”. Can anyone help me figure out this problem? My input script is as follows:

---------------------------------------------------------------

SRD diffusion demo - ellipsoids

units lj
atom_style ellipsoid
atom_modify first big
dimension 3
boundary p p p

create big ellipsoidal particles

lattice sc 0.4
region box block 0 10 0 10 0 10
create_box 1 box
create_atoms 1 region box

set type 1 mass 1.0
set type 1 shape 3.0 1.0 1.0
group big type 1
set group big quat/random 29898

velocity big create 1.44 87287 loop geom

equilibrate big particles

pair_style gayberne 1.0 3.0 1.0 4.0
pair_coeff 1 1 1.0 1.0 1 1 1 1 1 1

minimize 0.0 1.0e-8 1000 100000

neighbor 0.3 bin
neigh_modify delay 0 every 1 check yes

fix 1 big nve/asphere

compute rot big temp/asphere
compute diameter all property/atom shapex shapey shapez
compute q all property/atom quatw quati quatj quatk

dump 1 all custom 10 dump.ellipsoid.equil id type x y z c_diameter[1] c_diameter[2] c_diameter[3] c_q[1] c_q[2] c_q[3] c_q[4]
dump_modify 1 colname c_q[1] quatw colname c_q[2] quati colname c_q[3] quatj colname c_q[4] quatk
dump_modify 1 colname c_diameter[1] shapex colname c_diameter[2] shapey colname c_diameter[3] shapez

thermo_style custom step temp c_rot epair etotal press
thermo 100

run 1000

undump 1
unfix 1

timestep 0.0005

fix 1 big nve/asphere

#thermo_modify temp tbig
thermo 1000

dump 1 all custom 1000 dump.ellipsoid id type x y z c_diameter[1] c_diameter[2] c_diameter[3] c_q[1] c_q[2] c_q[3] c_q[4]
dump_modify 1 colname c_q[1] quatw colname c_q[2] quati colname c_q[3] quatj colname c_q[4] quatk
dump_modify 1 colname c_diameter[1] shapex colname c_diameter[2] shapey colname c_diameter[3] shapez

run 100000

--------------------------------------------------------------------------------

My output error is:

-------------------------------------------------------------------------------

LAMMPS (23 Jun 2022 - Update 2)
Lattice spacing in x,y,z = 1.3572088 1.3572088 1.3572088
Created orthogonal box = (0 0 0) to (13.572088 13.572088 13.572088)
1 by 2 by 2 MPI processor grid
Created 1000 atoms
using lattice units in orthogonal box = (0 0 0) to (13.572088 13.572088 13.572088)
create_atoms CPU = 0.004 seconds
Setting atom values …
1000 settings made for mass
Setting atom values …
1000 settings made for shape
1000 atoms in group big
Setting atom values …
1000 settings made for quat/random

CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE

Your simulation uses code contributions which should be cited:

  • pair gayberne command:
    The log file lists these citations in BibTeX format.

CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE

WARNING: Using ‘neigh_modify every 1 delay 0 check yes’ setting during minimization (…/min.cpp:187)
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info …
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.3
ghost atom cutoff = 4.3
binsize = 2.15, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair gayberne, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Setting up cg style minimization …
Unit style : lj
Current step : 0
Per MPI rank memory allocation (min/avg/max) = 6.031 | 6.031 | 6.031 Mbytes
Step Temp E_pair E_mol TotEng Press
0 1.44 1.9779748e+43 0 1.9779748e+43 -3.4881945e+47
1000 1.44 57.545011 0 59.702851 155.89971
Loop time of 75.7086 on 4 procs for 1000 steps with 1000 atoms

94.3% CPU use with 4 MPI tasks x no OpenMP threads

Minimization stats:
Stopping criterion = max iterations
Energy initial, next-to-last, final =
1.97797480969295e+43 57.5450105571295 57.5450105571295
Force two-norm initial, final = 2.4602815e+51 13261.809
Force max component initial, final = 1.5275392e+51 8768.2891
Final line search alpha, max atom move = 3.5639792e-07 0.003125
Iterations, force evaluations = 1000 7506

MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total

Pair | 62.329 | 64.663 | 66.906 | 27.7 | 85.41
Neigh | 0.11416 | 0.12212 | 0.13509 | 2.4 | 0.16
Comm | 8.2391 | 10.526 | 12.836 | 68.8 | 13.90
Output | 2.1197e-05 | 2.4109e-05 | 3.2686e-05 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.3976 | | | 0.53

Nlocal: 250 ave 261 max 244 min
Histogram: 2 0 0 1 0 0 0 0 0 1
Nghost: 2008.75 ave 2029 max 1985 min
Histogram: 1 0 0 1 0 0 0 1 0 1
Neighs: 16485 ave 17202 max 16026 min
Histogram: 2 0 0 0 1 0 0 0 0 1

Total # of neighbors = 65940
Ave neighs/atom = 65.94
Neighbor list builds = 201
Dangerous builds = 0
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Setting up Verlet run …
Unit style : lj
Current step : 1000
Time step : 0.005
Per MPI rank memory allocation (min/avg/max) = 7.055 | 7.055 | 7.055 Mbytes
Step Temp c_rot E_pair TotEng Press
1000 1.44 0.71963982 57.545011 59.702851 155.89971
ERROR: Lost atoms: original 1000 current 7 (…/thermo.cpp:481)
Last command: run 1000

----------------------------------------------------------------------

Regards,
Jenis

2

Typically lost atoms is due to overlapping atoms or bad parameters in the force field. There are various ways to push apart overlapping atoms frequently discussed here on the forum.

3

Your input uses

but your output shows

5

In the beginning, the system seems to have very large E_pair, TotEng and Press. I have manage to bring it down using minimize command. But it still seems insufficient. Is there anything that can be done to further bring it down?
Also, for such a system, should E_mol be always zero?

6

yes that too! It is confusing as to where the timestep of 0.005 comes from. Would you be able to correct me in this?

7

My bad, I missed the run 1000 command that happens before the timestep 0.0005 command. Since there is no timestep command before your first run, it uses the default of 0.005. You should move the timestep command earlier in the script.

You can minimize for longer. Note that your minimization stops because it reaches the number of steps you specified (1000), not because it reaches a force or energy minimum.

Your system has no bonds, angles, dihedrals, or impropers, so there is no energy contribution from them. You can read up on the meaning of thermo keywords in the thermo_style command.