Dear Lammps users,
I am trying to simulate polymers in spherical confinement with fixed boundary conditions and lj units. However, each time I get an error showing my bond atoms are missing. I have also fixed a wall using the fix wall/region command. Could you please suggest where I am going wrong? I am using LAMMPS 21 July 2020 version. I am attaching my script below:
units lj
boundary f f f
atom_style full
region 1r sphere 0.0 0.0 0.0 8 side in
read_data system.config
log logfile.txt
pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0 1.0 #1-all
pair_modify shift yes #truncated-shifted LJ
bond_style harmonic
bond_coeff 1 100.0 1.0
special_bonds lj 0.0 1.0 1.0
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes page 100000 one 2000
angle_style harmonic
angle_coeff * 0 180.0
velocity all create 1 79284 rot yes dist gaussian
fix 1 all wall/region 1r lj126 1.0 1.0 1.0
#fix 1 all wall/reflect xlo EDGE xhi EDGE ylo EDGE yhi EDGE zlo EDGE zhi EDGE
fix 2 all nvt temp 1.0 1.0 100
thermo_style custom step temp pe ke press spcpu
timestep 0.001
dump 1 all dcd 5000 trajafter10M.dcd
dump_modify 1 unwrap yes
minimize 1.0e-4 1.0e-6 1000000 10000
write_restart restart.min
thermo 10000
run 10000000
write_restart restartafter1M.min
print “Done”
Hi @pratyasha,
first of all it might be a good idea if possible to update your LAMMPS version since you are using a pretty old one right now. This might already remove potential issues you might end up dealing with.
You are then providing your input script but no log file and no data file which makes it impossible to guess where your simulation crashes and narrow the problem. A missing atoms in bond evaluation can come from both atoms going out of the box with fix boundary but also from atoms migrating too far from one another due to bad geometry and/or forcefield.
That being said, your pair_coeff parameters and wall/region command are a bit odd. If you want a purely repulsive interaction, you should put the cut off radius value at the energy minimum, not the \sigma distance below which the energy is positive (particles “penetrates” one another). 2.5 is generally a good cut off radius value, use 2^{\frac{1}{6}} if you want a purely repulsive model. I suggest you read articles using standard polymer model with reduced units (for example this one) to get a better understanding on how they work and what to expect. Reproducing these as a first step can give you confidence in your model and what you are doing and set you on track to making more complicated systems.
Dear Germain,
Thank You for your input. I will update the LAMMPS version. Meanwhile, I have also increased the cutoff and I am getting the same error as before. I am attaching the datafile and logfile below:
LAMMPS (21 Jul 2020)
Reading data file …
orthogonal box = (-7.1040000 -6.9505000 -6.8095000) to (6.8960000 7.0495000 7.1905000)
2 by 2 by 4 MPI processor grid
reading atoms …
250 atoms
scanning bonds …
1 = max bonds/atom
scanning angles …
1 = max angles/atom
reading bonds …
245 bonds
reading angles …
240 angles
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0.0 0.0 0.0
special bond factors coul: 0.0 0.0 0.0
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.013 seconds
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0.0 1.0 1.0
special bond factors coul: 0.0 0.0 0.0
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.000 seconds
WARNING: Using ‘neigh_modify every 1 delay 0 check yes’ setting during minimization (…/min.cpp:186)
Neighbor list info …
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 1.42246
ghost atom cutoff = 1.42246
binsize = 0.71123, bins = 20 20 20
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Setting up cg style minimization …
Unit style : lj
Current step : 0
WARNING: Communication cutoff 1.42246 is shorter than a bond length based estimate of 1.8. This may lead to errors. (…/comm.cpp:680)
Per MPI rank memory allocation (min/avg/max) = 7.843 | 7.91 | 8.111 Mbytes
Step Temp PotEng KinEng Press S/CPU
0 1 0.55103828 1.494 0.049592584 0
67 1 9.9345165e-11 1.494 0.090741565 8886.537
Loop time of 0.00754725 on 16 procs for 67 steps with 250 atoms
97.0% CPU use with 16 MPI tasks x no OpenMP threads
Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
0.551038279294794 9.93451653763446e-11 9.93451653763446e-11
Force two-norm initial, final = 351.28249926693053 0.0020474155562859864
Force max component initial, final = 68.95596715985492 0.00040999678067931225
Final line search alpha, max atom move = 0.5 0.00020499839033965612
Iterations, force evaluations = 67 194
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
Pair | 1.9707e-05 | 5.913e-05 | 0.00014541 | 0.0 | 0.78
Bond | 9.2823e-05 | 0.0005137 | 0.0011666 | 0.0 | 6.81
Neigh | 8.9834e-05 | 0.00011135 | 0.00013288 | 0.0 | 1.48
Comm | 0.0023211 | 0.0029219 | 0.0035893 | 0.7 | 38.71
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 3.7129e-05 | 0.00013454 | 0.00024476 | 0.0 | 1.78
Other | | 0.003807 | | | 50.44
Nlocal: 15.6250 ave 34 max 2 min
Histogram: 2 1 4 0 3 2 1 2 0 1
Nghost: 32.1250 ave 53 max 20 min
Histogram: 2 1 4 4 1 2 1 0 0 1
Neighs: 5.93750 ave 32 max 0 min
Histogram: 7 2 5 1 0 0 0 0 0 1
Total # of neighbors = 95
Ave neighs/atom = 0.38
Ave special neighs/atom = 5.7600000
Neighbor list builds = 3
Dangerous builds = 0
System init for write_restart …
WARNING: Communication cutoff 1.42246 is shorter than a bond length based estimate of 1.8. This may lead to errors. (…/comm.cpp:680)
Setting up Verlet run …
Unit style : lj
Current step : 67
Time step : 0.0001
WARNING: Communication cutoff 1.42246 is shorter than a bond length based estimate of 1.8. This may lead to errors. (…/comm.cpp:680)
Per MPI rank memory allocation (min/avg/max) = 6.718 | 6.785 | 6.986 Mbytes
Step Temp PotEng KinEng Press S/CPU
67 1 9.9345165e-11 1.494 0.090741565 0
10000 0.76131073 0.35662032 1.1373982 0.068260609 58936.957
20000 0.80193428 0.29602016 1.1980898 0.0022126465 65296.049
30000 0.76454199 0.35159847 1.1422257 -0.013831385 65516.375
ERROR on proc 9: Bond atoms 147 148 missing on proc 9 at step 32351 (…/ntopo_bond_all.cpp:60)
My guess is the problem lies with the pair_coeff settings. By truncating the LJ potential at 1 sigma, you are not truncating it in the minimum, but at the point where epsilon is zero. That means, you have a substantial jump in forces which can cause all kinds of problems.
You probably want something like this:
pair_style lj/cut $(2.0^(1.0/6.0))
pair_coeff * * 1.0 1.0
pair_modify shift yes
comm_modify cutoff 2.0
There should be no need to change the page and one settings for your kind of system
Dear Dr. Axel,
Thank you for your suggestion. I did implement the new cutoff and this is my current script.
scriptfile.txt (858 Bytes)
However, I am still getting the same error. I have noticed that my system does run with normal lj. However, with purely repulsive forces, I am getting errors. I am attaching the log file too. I am not able to understand why this is happening. Could you please suggest what
log (3.7 KB)
LAMMPS (21 Jul 2020)
Reading data file …
orthogonal box = (-7.1040000 -6.9505000 -6.8095000) to (6.8960000 7.0495000 7.1905000)
2 by 2 by 4 MPI processor grid
reading atoms …
250 atoms
scanning bonds …
1 = max bonds/atom
scanning angles …
1 = max angles/atom
reading bonds …
245 bonds
reading angles …
240 angles
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0.0 0.0 0.0
special bond factors coul: 0.0 0.0 0.0
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.010 seconds
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0.0 1.0 1.0
special bond factors coul: 0.0 0.0 0.0
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.000 seconds
Neighbor list info …
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 1.422462
ghost atom cutoff = 2
binsize = 0.71123102, bins = 20 20 20
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Setting up cg style minimization …
Unit style : lj
Current step : 0
Per MPI rank memory allocation (min/avg/max) = 7.848 | 7.915 | 8.116 Mbytes
Step Temp PotEng KinEng Press S/CPU
0 0.5 0.55103828 0.747 0.0042208636 0
59 0.5 8.8827844e-11 0.747 0.045371602 8827.8263
Loop time of 0.00669307 on 16 procs for 59 steps with 250 atoms
100.0% CPU use with 16 MPI tasks x no OpenMP threads
Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
0.551038279304385 8.8827844465827e-11 8.8827844465827e-11
Force two-norm initial, final = 347.69666751410654 0.0024802832289074167
Force max component initial, final = 68.95596715985492 0.0004333265086616146
Final line search alpha, max atom move = 0.5 0.0002166632543308073
Iterations, force evaluations = 59 187
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
Pair | 1.9925e-05 | 5.2014e-05 | 0.00014715 | 0.0 | 0.78
Bond | 9.6041e-05 | 0.00045895 | 0.0011975 | 0.0 | 6.86
Neigh | 9.983e-05 | 0.00011251 | 0.00012619 | 0.0 | 1.68
Comm | 0.0023069 | 0.0029668 | 0.0035301 | 0.7 | 44.33
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 4.7639e-05 | 0.00011651 | 0.00025222 | 0.0 | 1.74
Other | | 0.002986 | | | 44.62
Nlocal: 15.6250 ave 34 max 2 min
Histogram: 2 1 4 0 3 2 1 2 0 1
Nghost: 48.8125 ave 68 max 32 min
Histogram: 1 5 0 2 0 1 2 2 2 1
Neighs: 6.00000 ave 33 max 0 min
Histogram: 7 2 5 1 0 0 0 0 0 1
Total # of neighbors = 96
Ave neighs/atom = 0.384
Ave special neighs/atom = 5.7600000
Neighbor list builds = 3
Dangerous builds = 0
System init for write_restart …
Setting up Verlet run …
Unit style : lj
Current step : 59
Time step : 0.0005
Per MPI rank memory allocation (min/avg/max) = 6.723 | 6.79 | 6.991 Mbytes
Step Temp PotEng KinEng Press S/CPU
59 0.5 8.8827844e-11 0.747 0.045371602 0
ERROR on proc 1: Bond atoms 225 226 missing on proc 1 at step 7353 (…/ntopo_bond_all.cpp:60)
Here is the problem:
You have a confinement potential using a spherical region centered on 0.0 0.0 0.0 with a radius of 9 sigma, but you also have a box with fixed boundaries that is smaller than your confinement potential. So you can “loose” atoms from the boundary even though they are within the confined area.
Try increasing your box dimensions, or shrinking your confinement region or changing to shrinkwrap boundaries.
The fact that it works with a LJ potential that has the attractive branch included is not a surprise: the system will “clump together” and thus it is much less likely for atoms to get close to the boundary.
Dear Dr. Axel,
Thank You. I have tried using shrinkwrap boundary conditions. The system works well.
Because I want to see how the polymer might behave in a spherical confinement,I am interested in fixed boundary onditions. But the problem remains with fixed boundary condition. If I shrink my confinement region ,I get the error: ERROR on proc 9: Particle outside surface of region used in fix wall/region (…/fix_wall_region.cpp:300).
man, did you use a larger simulation box?
Are all atoms in the initial configuration inside the confinement region??
You (obviously) cannot shrink your confinement to be smaller than your atom positions.
Even though this is the “LAMMPS Beginners” category and thus there is some extra allowance for asking obvious questions, this is beginning to become very tedious. At some point you will have to do your own thinking and learning to resolve such (rather obvious!) problems of your own and it looks to me the point has come.
This logic makes NO sense. The boundary that you want to study is defined by your confinement and not the box. Thus there is no reason to favor fixed boundaries over shrink-wrap.
I finally understood my mistake. Thank you for pointing things out.