My script is taking too long to run

Igwe · April 27, 2022, 2:50am

Hello there,
I am trying to simulate silica (SiO2) glass using the BKS potential. But my simulation is taking painfully long time to run. Sadly I can’t identify the error.
Please, how can I optimize the code in order to shorten the computation time?
I have attached my input script below.

# input file for silica glass
units		real
dimension 	3
boundary 	p p p
atom_style	charge
#create atoms
lattice custom  15.5000 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 &
                basis 0.0 0.0 0.0
region          my_box block 0 35.1 0 35.1 0 35.1            
create_box      2 my_box
create_atoms    1 random 144 12345 my_box
create_atoms    2 random 288 6789 my_box
mass 1		28.8855  #Si
mass 2		15.9994	 #O
set type	1 charge  2.4  #Si
set type	2 charge -1.2  #O	
velocity all create 100 28459 rot yes dist gaussian mom yes
#potential
include		pot.dat
#outputs
thermo		100
thermo_style	custom step temp pe press vol lx density
dump		2 all custom 1000 md.lammpstrj id type x y z
# initial minimisation
minimize	1.0e-10	1.0e-10	100000	100000
# Initial mixing
fix		1 all nvt temp 2700 2700 100
run		200000
unfix		1
fix		1 all npt temp 2700 2700 100 iso 0 0 1000
run		100000
unfix		1
# cooling at  K/ps
fix		1 all npt temp 2700 300 100 iso 0 0 1000
run		2400000
unfix		1
# final relaxation
fix		1 all npt temp 300 300 100 iso 0 0 1000
run		100000
unfix		1
write_data	S300K.dat
write_restart	s300K.rest

Thank you in advance.
Igwe.

akohlmey · April 27, 2022, 3:47am

Please also provide:

what is your LAMMPS version?
the content of the “pot.dat” file
a log of the output
the command line you are using

Igwe · April 27, 2022, 7:47am

LAMMPS version: 15 Nov. 2018

the content of the “pot.dat” file:

# BKS potential
pair_style	buck/coul/long 5.5 10.0
pair_coeff	1 1 0.0		1.0	0.0		# Si-Si
pair_coeff	1 2 414512.50	0.20520	2075.270	# Si-O
pair_coeff	2 2 31932.360	0.26231	4030.1136	# o-o

# MD parameters
kspace_style	ewald 1.0e-5
neighbor	2.0 bin
neigh_modify	every 1 check yes
timestep	1.0

a log of the output

LAMMPS (15 Nov 2018)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:87)
  using 1 OpenMP thread(s) per MPI task
# input file for silica glass

units		real
dimension 	3
boundary 	p p p
atom_style	charge

#create atoms
lattice custom  15.5000 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0                 basis 0.0 0.0 0.0
Lattice spacing in x,y,z = 15.5 15.5 15.5

region          my_box block 0 35.1 0 35.1 0 35.1
create_box      2 my_box
Created orthogonal box = (0 0 0) to (544.05 544.05 544.05)
  1 by 1 by 1 MPI processor grid

create_atoms    1 random 144 12345 my_box
Created 144 atoms
  Time spent = 0.00841308 secs
create_atoms    2 random 288 6789 my_box
Created 288 atoms
  Time spent = 0 secs

mass 1		28.8855  #Si
mass 2		15.9994	 #O

set type	1 charge  2.4  #Si
  144 settings made for charge
set type	2 charge -1.2  #O
  288 settings made for charge


velocity all create 100 28459 rot yes dist gaussian mom yes

#potential
include		pot.dat
# BKS potential
pair_style	buck/coul/long 5.5 10.0
pair_coeff	1 1 0.0		1.0	0.0		# Si-Si
pair_coeff	1 2 414512.50	0.20520	2075.270	# Si-O
pair_coeff	2 2 31932.360	0.26231	4030.1136	# o-o

# MD parameters
kspace_style	ewald 1.0e-5
neighbor	2.0 bin
neigh_modify	every 1 check yes
timestep	1.0

#outputs
thermo		100
thermo_style	custom step temp pe press vol lx density
dump		2 all custom 1000 md.lammpstrj id type x y z

# initial minimisation
minimize	1.0e-10	1.0e-10	100000	100000
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:168)
Ewald initialization ...
  using 12-bit tables for long-range coulomb (../kspace.cpp:321)
  G vector (1/distance) = 0.238709
  estimated absolute RMS force accuracy = 0.00460791
  estimated relative force accuracy = 1.38766e-05
  KSpace vectors: actual max1d max3d = 2032281 99 3940299
                  kxmax kymax kzmax  = 99 99 99
Neighbor list info ...
  update every 1 steps, delay 0 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 12
  ghost atom cutoff = 12
  binsize = 6, bins = 91 91 91
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair buck/coul/long, perpetual
      attributes: half, newton on
      pair build: half/bin/atomonly/newton
      stencil: half/bin/3d/newton
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 624.5 | 624.5 | 624.5 Mbytes
Step Temp PotEng Press Volume Lx Density 
       0          100   -1285.8304 -0.071549819 1.6103358e+08       544.05 9.040666e-05

the command line you are using:
"lmp_serial < SiO2_B.in

Thank you.

akohlmey · April 27, 2022, 8:21am

I have two more questions:

why do you create such a large box? And as a result such a low density initial system?
why do you use kspace style ewald and not pppm?

Igwe · April 27, 2022, 11:44am

I changed the simulation box to a smaller one ( 0 15.0 0 15.0 0 15.0), changed the kspace style to pppm. However, the system is running very slow and gives this error message:
ERROR on proc 0: Out of range atoms - cannot compute PPPM (…/pppm.cpp:1945).

Please what can I do to solve this problem? Below is the log of the output:

LAMMPS (15 Nov 2018)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (…/comm.cpp:87)
using 1 OpenMP thread(s) per MPI task

input file for silica glass

units real
dimension 3
boundary p p p
atom_style charge

#create atoms
lattice custom 15.5000 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0
Lattice spacing in x,y,z = 15.5 15.5 15.5

region my_box block 0 15.0 0 15.0 0 15.0
create_box 2 my_box
Created orthogonal box = (0 0 0) to (232.5 232.5 232.5)
1 by 1 by 1 MPI processor grid

create_atoms 1 random 144 12345 my_box
Created 144 atoms
Time spent = 0.00102997 secs
create_atoms 2 random 288 6789 my_box
Created 288 atoms
Time spent = 0 secs

mass 1 28.8855 #Si
mass 2 15.9994 #O

set type 1 charge 2.4 #Si
144 settings made for charge
set type 2 charge -1.2 #O
288 settings made for charge

velocity all create 100 28459 rot yes dist gaussian mom yes

#potential
include pot.dat

BKS potential

pair_style buck/coul/long 5.5 10.0
pair_coeff 1 1 0.0 1.0 0.0 # Si-Si
pair_coeff 1 2 414512.50 0.20520 2075.270 # Si-O
pair_coeff 2 2 31932.360 0.26231 4030.1136 # o-o

MD parameters

kspace_style pppm 1.0e-4
neighbor 2.0 bin
neigh_modify delay 2 every 1 check yes
timestep 1.0
#delete_atoms overlap 0.1 all all

#outputs
thermo 100
thermo_style custom step temp pe press vol lx density
dump 2 all custom 1000 md.lammpstrj id type x y z

initial minimisation

minimize 1.0e-10 1.0e-10 100000 100000
WARNING: Using ‘neigh_modify every 1 delay 0 check yes’ setting during minimization (…/min.cpp:168)
PPPM initialization …
using 12-bit tables for long-range coulomb (…/kspace.cpp:321)
G vector (1/distance) = 0.217539
grid = 72 72 72
stencil order = 5
estimated absolute RMS force accuracy = 0.0441621
estimated relative force accuracy = 0.000132993
using double precision FFTs
3d grid and FFT values/proc = 456533 373248
Neighbor list info …
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 39 39 39
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 55.19 | 55.19 | 55.19 Mbytes
Step Temp PotEng Press Volume Lx Density
0 100 -2947.8847 -4.2259964 12568078 232.5 0.0011583719
100 100 -9424.8238 -7.6864537 12568078 232.5 0.0011583719
200 100 -15072.744 -9.3869562 12568078 232.5 0.0011583719
300 100 -21663.574 -9.5092738 12568078 232.5 0.0011583719
400 100 -27441.332 -7.9772078 12568078 232.5 0.0011583719
500 100 -30780.863 -9.7996849 12568078 232.5 0.0011583719
600 100 -35035.801 -9.9700243 12568078 232.5 0.0011583719
700 100 -39165.413 -10.269136 12568078 232.5 0.0011583719
800 100 -43405.625 -10.158131 12568078 232.5 0.0011583719
900 100 -48058.149 -10.925422 12568078 232.5 0.0011583719
1000 100 -52427.961 -9.2218171 12568078 232.5 0.0011583719
1100 100 -55482.942 -8.9668087 12568078 232.5 0.0011583719
1200 100 -58534.132 -8.6989445 12568078 232.5 0.0011583719
1300 100 -62857.09 -7.3863786 12568078 232.5 0.0011583719
1400 100 -65649.199 -6.9545053 12568078 232.5 0.0011583719
1500 100 -67358.709 -7.9208115 12568078 232.5 0.0011583719
1600 100 -68935.352 -8.384078 12568078 232.5 0.0011583719
1700 100 -71855.668 -7.5914616 12568078 232.5 0.0011583719
1800 100 -73966.802 -6.2087742 12568078 232.5 0.0011583719
1900 100 -75759.098 -7.3707116 12568078 232.5 0.0011583719
2000 100 -78428.089 -5.9684084 12568078 232.5 0.0011583719
2100 100 -80146.737 -5.8971696 12568078 232.5 0.0011583719
2200 100 -81286.488 -7.8780854 12568078 232.5 0.0011583719
2300 100 -82961.738 -7.0036067 12568078 232.5 0.0011583719
2400 100 -84446.055 -5.5900274 12568078 232.5 0.0011583719
2500 100 -86407.021 -6.9769645 12568078 232.5 0.0011583719
2600 100 -88013.235 -5.118597 12568078 232.5 0.0011583719
2700 100 -89499.036 -6.8546617 12568078 232.5 0.0011583719
2800 100 -91395.314 -4.8720255 12568078 232.5 0.0011583719
2900 100 -93081.746 -4.4753371 12568078 232.5 0.0011583719
3000 100 -94620.013 -5.0686173 12568078 232.5 0.0011583719
3100 100 -95666.492 -5.7505246 12568078 232.5 0.0011583719
3200 100 -97591.92 -3.9778653 12568078 232.5 0.0011583719
3300 100 -98758.936 -3.9486505 12568078 232.5 0.0011583719
3400 100 -100146.3 -4.3258372 12568078 232.5 0.0011583719
3500 100 -101445.26 -4.7338972 12568078 232.5 0.0011583719
3600 100 -102822.93 -4.3075052 12568078 232.5 0.0011583719
3700 100 -103991.9 -3.0897534 12568078 232.5 0.0011583719
3800 100 -104889.09 -2.8390408 12568078 232.5 0.0011583719
3900 100 -105831.29 -2.5734776 12568078 232.5 0.0011583719
4000 100 -106779.8 -2.9171904 12568078 232.5 0.0011583719
4100 100 -108174.28 -2.3989923 12568078 232.5 0.0011583719
4200 100 -108837.64 -2.7518206 12568078 232.5 0.0011583719
4300 100 -109943.78 -2.6966478 12568078 232.5 0.0011583719
4400 100 -110502.98 -2.5154998 12568078 232.5 0.0011583719
4500 100 -111424.8 -2.452291 12568078 232.5 0.0011583719
4560 100 -6.3742531e+84 -6.0059497e+82 12568078 232.5 0.0011583719
Loop time of 1289.96 on 1 procs for 4560 steps with 432 atoms

99.8% CPU use with 1 MPI tasks x 1 OpenMP threads

Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-2947.88474049 -6.37425306156e+84 -6.37425306156e+84
Force two-norm initial, final = 179.588 1.84621e+99
Force max component initial, final = 31.9848 1.26649e+99
Final line search alpha, max atom move = 8.97654e-114 1.13687e-14
Iterations, force evaluations = 4560 6561

MPI task timing breakdown:

Pair | 0.29806 | 0.29806 | 0.29806 | 0.0 | 0.02
Kspace | 1288.8 | 1288.8 | 1288.8 | 0.0 | 99.91
Neigh | 0.071902 | 0.071902 | 0.071902 | 0.0 | 0.01
Comm | 0.041991 | 0.041991 | 0.041991 | 0.0 | 0.00
Output | 0.5452 | 0.5452 | 0.5452 | 0.0 | 0.04
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.183 | | | 0.01

Nlocal: 432 ave 432 max 432 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 132 ave 132 max 132 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 317 ave 317 max 317 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 317
Ave neighs/atom = 0.733796
Neighbor list builds = 350
Dangerous builds = 0

Initial mixing

fix 1 all nvt temp 2700 2700 100
run 200000
PPPM initialization …
using 12-bit tables for long-range coulomb (…/kspace.cpp:321)
G vector (1/distance) = 0.217539
grid = 72 72 72
stencil order = 5
estimated absolute RMS force accuracy = 0.0441621
estimated relative force accuracy = 0.000132993
using double precision FFTs
3d grid and FFT values/proc = 456533 373248
Per MPI rank memory allocation (min/avg/max) = 54.07 | 54.07 | 54.07 Mbytes
Step Temp PotEng Press Volume Lx Density
4560 100 -6.3742531e+84 -6.0059497e+82 12568078 232.5 0.0011583719
ERROR on proc 0: Out of range atoms - cannot compute PPPM (…/pppm.cpp:1945)
Last command: run 200000

akohlmey · April 27, 2022, 11:57am

Slow compared to what?

There is a general issue with starting from random positions and that is that you are creating a very high potential energy configuration, since you have no provisions to reduce overlaps when creating the initial geometry. In addition, you are using a pair potential, that has a known problem with close contacts. You can see this from the massive drop in potential energy during minimization. Some atoms have overcome the barrier by the Born potential and now the Coulomb part is dominant and opposite charged atoms will move on top of each other. The PPPM error is a consequence of that.

On a more general level, what kind of equilibrated system (what type of structure, what density, what thermodynamic state etc.) is it that you want to generate here? That information will be crucial to determine the best approach to design a suitable workflow.

To avoid the immediate issue, you need to use a different pair style for the initial minimization: one that doesn’t have a discontinuity a short separation (like Coulomb) and that will guarantee a suitable separation of atoms before switching to the desired force field. One option would be pair style soft.

akohlmey · April 27, 2022, 2:40pm

I have experimented with your input and it looks like your force field parameters are not correct.
It cannot run an MD for a pre-equilibrated quartz crystal.

Please note that BKS needs to be augmented with a short range repulsive potential to keep Si-Si and Si-O atoms overlapping due to the strongly attractive nature of the potential after atoms have come too close.