Recently added pair_style harmonic/cut seems to be wrong

LAMMPS LAMMPS Development
,
1

Dear LAMMPS users,

I noticed that the pair_style harmonic/cut has been added in the latest available version of LAMMPS (lammps-17Feb2022).

It appears that in the web guide at pair_style harmonic/cut command — LAMMPS documentation, the following example:
pair_style harmonic/cut 2.5
should be replaced with:
pair_style harmonic/cut
Indeed, otherwise I get an error message presumably coming from the line 149 of the file pair_harmonic_cut.cpp:

     if (narg > 0) error->all(FLERR, "Illegal pair_style command");

The possible mistake in the source code pair_harmonic_cut.cpp I referred to in the title of this post is related to the functional form of the potential.
Lines 102 and 103 of this code read:

    const double philj = factor_lj * delta * k[itype][jtype];
    const double fpair = delta * philj / r;

while lines 299 and 300 read:

    const double philj = factor_lj * delta * k[itype][jtype];
    fforce = delta * philj / r;

If the potential form, as reported in the web page mentioned above, is E = k (r_c -r)^2 for r < r_c, and philj = U - Ecutoff and fforce = -(1/r) dU/dr, then I would expect lines 102 and 299 to be replaced by:

    const double philj = factor_lj * delta * delta * k[itype][jtype];

and lines 103 and 300 respectively by:

    const double fpair = 2.0 * philj / (r * delta);
    fforce = 2.0 * philj / (r * delta);

I simulated a bi-disperse system of particles in the bulk with this potential by using tables and the energy agrees with the one obtained from using the new pair_style harmonic/cut only if I change the mentioned lines of code in the file pair_harmonic_cut.cpp, otherwise the difference is of about two orders of magnitude.

In any case, also by changing these lines of code, the energy obtained by using pair_style harmonic/cut fluctuates much more respect to the case in which tables are used.

Am I making some mistake?
Thank you for any suggestions.

Kind regards,
Fabio

2

Thanks for your report.

The documentation had already been corrected in the current development version: Please see
https://docs.lammps.org/latest/pair_harmonic_cut.html

You are correct about energy and force. It looks like they got swapped due to misreading of the meaning of columns in the table file for pair style table.

This should be corrected in this commit Collected small changes and fixes by akohlmey · Pull Request #3177 · lammps/lammps · GitHub
which should be included in the next patch release of LAMMPS.

3

Dear Axel,

Thank you very much for your reply.

I am still having some problems with a bi-disperse system of harmonic repulsive particles.
I am just testing equilibration in the NVT and NVE ensembles in the bulk with time step dt=0.005 in lj units.

In the mono-disperse case, simulations with Tables or with the new pair_style give compatible results (temperature (T), pressure (P), energy (E) vs time).

In the case of a 50:50 bi-disperse system of harmonic repulsive particles, for NVT (T=0.0005) the new pair_style gives T~0.0007 with large fluctuations, P very low and E larger and very fluctuating respect to the case with Tables.
If dt=0.001, T~0.00053 with smaller fluctuations, but for P and E there are no changes with respect to the case with dt=0.005.

For NVE, in the case of the new pair_style T and P diverge quadratically in time in order to keep E ~ constant.

For the bi-disperse system it looks as if particle diameters are not properly taken into account.

Again, I’m not sure if I’m doing something wrong.

Thanks again for any suggestions.

Kind regards,
Fabio

4

I would have to see your input deck and your tables for any specific comments.

5

Dear Axel,

I tried to attached some file, but unfortunately being a new user I cannot do that, so I’m writing below the input LAMMPS script, log file, and initial and final part of tables.
If necessary, I’ll write the full tables (they are ~ 20 000 lines long).

Thank you for your time.

Kind regards,
Fabio

P.S.: using these tables with linear interpolation implies an error respect to the exact functional form of the potential of the same order of the functional form if using float.

Simulation NVE of bi-disperse system with pair_style

input.in:

# Binary additive repulsive harmonic in 3D

units	 	lj  
dimension	3
atom_style	atomic
boundary	p p p

read_data	config_equilibrated_T0.0005.dat

### VARIABLES ###
variable      T1	equal 0.0005	
variable      dt	equal 0.005  # LJ default = 0.005
variable      Nrun	equal 10000000
variable      Nthermo	equal ${Nrun}/1000
variable      Ndump	equal ${Nrun}/100
#################

### PARAMETERS ###
mass		* 1.0

pair_style	harmonic/cut
pair_coeff 1 1	1.0 1.0
pair_coeff 1 2	1.0 1.2
pair_coeff 2 2	1.0 1.4

##################

timestep	${dt}

fix		1 all nve

# Assigning random velocities Maxwell distributed at T=T1
velocity all create ${T1} 54321
run 0
velocity all scale ${T1}


thermo_style	custom step temp press pe
thermo		${Nthermo}

dump		1 all atom ${Ndump} pos.dat
dump_modify	1 scale no


run		${Nrun}

log.lammps:

LAMMPS (17 Feb 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
  using 1 OpenMP thread(s) per MPI task
# Binary additive repulsive harmonic in 3D

units	 	lj
dimension	3
atom_style	atomic
boundary	p p p

read_data	config_equilibrated_T0.0005.dat
Reading data file ...
  orthogonal box = (0 0 0) to (18.155025 18.155025 18.155025)
  1 by 2 by 3 MPI processor grid
  reading atoms ...
  4000 atoms
  read_data CPU = 0.006 seconds

### VARIABLES ###
variable      T1	equal 0.0005
variable      dt	equal 0.005  # LJ default = 0.005
variable      Nrun	equal 10000000
variable      Nthermo	equal ${Nrun}/1000
variable      Nthermo	equal 10000000/1000
variable      Ndump	equal ${Nrun}/100
variable      Ndump	equal 10000000/100
#################

### PARAMETERS ###
mass		* 1.0
pair_style	harmonic/cut
pair_coeff 1 1	1.0 1.0
pair_coeff 1 2	1.0 1.2
pair_coeff 2 2	1.0 1.4
##################

timestep	${dt}
timestep	0.005

fix		1 all nve

# Assigning random velocities Maxwell distributed at T=T1
velocity all create ${T1} 54321
velocity all create 0.0005 54321
run 0
  generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 1.7
  ghost atom cutoff = 1.7
  binsize = 0.85, bins = 22 22 22
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair harmonic/cut, perpetual
      attributes: half, newton on
      pair build: half/bin/atomonly/newton
      stencil: half/bin/3d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 3.117 | 3.117 | 3.119 Mbytes
Step Temp E_pair E_mol TotEng Press 
       0       0.0005 0.00044873224            0 0.0011985447  0.011035339 
Loop time of 3.81917e-06 on 6 procs for 0 steps with 4000 atoms

82.9% CPU use with 6 MPI tasks x 1 OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0          | 0          | 0          |   0.0 |  0.00
Neigh   | 0          | 0          | 0          |   0.0 |  0.00
Comm    | 0          | 0          | 0          |   0.0 |  0.00
Output  | 0          | 0          | 0          |   0.0 |  0.00
Modify  | 0          | 0          | 0          |   0.0 |  0.00
Other   |            | 3.819e-06  |            |       |100.00

Nlocal:        666.667 ave         687 max         652 min
Histogram: 1 0 2 0 2 0 0 0 0 1
Nghost:        1032.33 ave        1060 max        1001 min
Histogram: 1 1 0 0 0 2 0 0 0 2
Neighs:        3793.33 ave        3925 max        3712 min
Histogram: 1 1 1 1 0 1 0 0 0 1

Total # of neighbors = 22760
Ave neighs/atom = 5.69
Neighbor list builds = 0
Dangerous builds = 0
velocity all scale ${T1}
velocity all scale 0.0005


thermo_style	custom step temp press pe
thermo		${Nthermo}
thermo		10000

dump		1 all atom ${Ndump} pos.dat
dump		1 all atom 100000 pos.dat
dump_modify	1 scale no


run		${Nrun}
run		10000000
  generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 4.143 | 4.144 | 4.145 Mbytes
Step Temp Press PotEng 
       0       0.0005  0.011035339 0.00044873224 
   10000    19.658166    13.302539   0.20109979 
   20000    75.966326    50.921819   0.18321388 
   30000    137.55744    92.078141   0.17662565 
   40000    201.19151    134.59502   0.16336673 
   50000    266.65984    178.34571   0.16396462 
   60000    333.79777    223.21087   0.15879781 
   70000    402.61639     269.1974   0.15421099 
   80000    472.86276    316.14078   0.15301772 
   90000    544.57066     364.0559   0.15114932 
  100000    617.55304    412.83028   0.14550086 
  110000    691.64192    462.33925   0.14361545 
................
 3100000    66502.333    44442.559   0.18280083 
 3110000     66879.59    44694.681   0.19483811 
 3120000    67258.712    44948.033   0.17833582 
 3130000    67637.775    45201.358   0.18770853 
 3140000    68018.381     45455.71   0.18850005 
 3150000    68399.488    45710.394   0.18461318 
ERROR: Lost atoms: original 4000 current 3997 (src/thermo.cpp:435)
Last command: run		${Nrun}

Simulation NVE of bi-disperse system with tables

input.in:

# Binary additive repulsive harmonic in 3D

units	 	lj  
dimension	3
atom_style	atomic
boundary	p p p

read_data	config_equilibrated_T0.0005.dat

### VARIABLES ###
variable      T1	equal 0.0005	
variable      dt	equal 0.005  # LJ default = 0.005
variable      Nrun	equal 10000000
variable      Nthermo	equal ${Nrun}/1000
variable      Ndump	equal ${Nrun}/100
#################

### PARAMETERS ###
mass		* 1.0
pair_style	table linear 20000
pair_coeff 1 1	Table_lammps-potential_epsilon1_N20k_11.dat AA 1.0
pair_coeff 1 2	Table_lammps-potential_epsilon1_N20k_12.dat AB 1.2
pair_coeff 2 2	Table_lammps-potential_epsilon1_N20k_22.dat BB 1.4
##################

timestep	${dt}

fix		1 all nve


# Assigning random velocities Maxwell distributed at T=T1
velocity all create ${T1} 54321
run 0
velocity all scale ${T1}


thermo_style	custom step temp press pe
thermo		${Nthermo}

dump		1 all atom ${Ndump} pos.dat
dump_modify	1 scale no


run		${Nrun}

log.lammps:

LAMMPS (17 Feb 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
  using 1 OpenMP thread(s) per MPI task
# Binary additive repulsive harmonic in 3D

units	 	lj
dimension	3
atom_style	atomic
boundary	p p p

read_data	config_equilibrated_T0.0005.dat
Reading data file ...
  orthogonal box = (0 0 0) to (18.155025 18.155025 18.155025)
  1 by 2 by 3 MPI processor grid
  reading atoms ...
  4000 atoms
  read_data CPU = 0.005 seconds

### VARIABLES ###
variable      T1	equal 0.0005
variable      dt	equal 0.005  # LJ default = 0.005
variable      Nrun	equal 10000000
variable      Nthermo	equal ${Nrun}/1000
variable      Nthermo	equal 10000000/1000
variable      Ndump	equal ${Nrun}/100
variable      Ndump	equal 10000000/100
#################

### PARAMETERS ###
mass		* 1.0

pair_style	table linear 20000
pair_coeff 1 1	Table_lammps-potential_epsilon1_N20k_11.dat AA 1.0
WARNING: 9982 of 20000 force values in table AA are inconsistent with -dE/dr.
WARNING:  Should only be flagged at inflection points (src/pair_table.cpp:465)
pair_coeff 1 2	Table_lammps-potential_epsilon1_N20k_12.dat AB 1.2
WARNING: 9962 of 20000 force values in table AB are inconsistent with -dE/dr.
WARNING:  Should only be flagged at inflection points (src/pair_table.cpp:465)
pair_coeff 2 2	Table_lammps-potential_epsilon1_N20k_22.dat BB 1.4
WARNING: 9974 of 20000 force values in table BB are inconsistent with -dE/dr.
WARNING:  Should only be flagged at inflection points (src/pair_table.cpp:465)

##################

timestep	${dt}
timestep	0.005

fix		1 all nve


# Assigning random velocities Maxwell distributed at T=T1
velocity all create ${T1} 54321
velocity all create 0.0005 54321
run 0
  generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 1.7
  ghost atom cutoff = 1.7
  binsize = 0.85, bins = 22 22 22
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair table, perpetual
      attributes: half, newton on
      pair build: half/bin/atomonly/newton
      stencil: half/bin/3d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 3.117 | 3.117 | 3.119 Mbytes
Step Temp E_pair E_mol TotEng Press 
       0       0.0005 0.0005783626            0 0.0013281751   0.01407435 
Loop time of 3.234e-06 on 6 procs for 0 steps with 4000 atoms

103.1% CPU use with 6 MPI tasks x 1 OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0          | 0          | 0          |   0.0 |  0.00
Neigh   | 0          | 0          | 0          |   0.0 |  0.00
Comm    | 0          | 0          | 0          |   0.0 |  0.00
Output  | 0          | 0          | 0          |   0.0 |  0.00
Modify  | 0          | 0          | 0          |   0.0 |  0.00
Other   |            | 3.234e-06  |            |       |100.00

Nlocal:        666.667 ave         687 max         652 min
Histogram: 1 0 2 0 2 0 0 0 0 1
Nghost:        1032.33 ave        1060 max        1001 min
Histogram: 1 1 0 0 0 2 0 0 0 2
Neighs:        3793.33 ave        3925 max        3712 min
Histogram: 1 1 1 1 0 1 0 0 0 1

Total # of neighbors = 22760
Ave neighs/atom = 5.69
Neighbor list builds = 0
Dangerous builds = 0
velocity all scale ${T1}
velocity all scale 0.0005


thermo_style	custom step temp press pe
thermo		${Nthermo}
thermo		10000

dump		1 all atom ${Ndump} pos.dat
dump		1 all atom 100000 pos.dat
dump_modify	1 scale no


run		${Nrun}
run		10000000
  generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 4.143 | 4.144 | 4.145 Mbytes
Step Temp Press PotEng 
       0       0.0005   0.01407435 0.0005783626 
   10000 0.00050282195  0.014172699 0.00057412447 
   20000  0.000499744  0.014163817 0.00057873735 
   30000 0.00049285284  0.014474179 0.00058907275 
   40000 0.00049964862  0.014226466 0.00057888114 
   50000 0.00049075864  0.014294125 0.00059221687 
   60000 0.00049603576    0.0141714 0.00058429815 
   70000 0.00050257565  0.014151142 0.00057449748 
   80000  0.000495562  0.014338538 0.00058500129 
   90000 0.00049429613    0.0142483 0.00058691007 
  100000 0.00050150504  0.014219893 0.00057610701 
  110000 0.0004958415  0.014308043 0.00058458982 
.............................
 9950000 0.00050087683  0.014307953 0.0005770781 
 9960000 0.00049839937  0.014213436 0.00058079255 
 9970000 0.0004997317  0.014119572 0.00057879861 
 9980000 0.00049935678   0.01436484  0.000579349 
 9990000 0.00049854215  0.014239987 0.00058056629 
10000000 0.00050358524  0.014076311 0.00057300705 
Loop time of 471.286 on 6 procs for 10000000 steps with 4000 atoms

Performance: 9166403.191 tau/day, 21218.526 timesteps/s
89.6% CPU use with 6 MPI tasks x 1 OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 166.38     | 178.43     | 193.44     |  73.8 | 37.86
Neigh   | 7.8755     | 7.9921     | 8.2227     |   3.9 |  1.70
Comm    | 205.96     | 218.03     | 228.16     |  60.3 | 46.26
Output  | 0.051892   | 0.054067   | 0.063472   |   1.8 |  0.01
Modify  | 34.573     | 35.134     | 35.587     |   5.9 |  7.45
Other   |            | 31.65      |            |       |  6.72

Nlocal:        666.667 ave         685 max         659 min
Histogram: 1 2 2 0 0 0 0 0 0 1
Nghost:        1035.83 ave        1060 max        1022 min
Histogram: 2 1 0 0 1 1 0 0 0 1
Neighs:        3799.83 ave        3895 max        3758 min
Histogram: 3 1 0 0 0 0 1 0 0 1

Total # of neighbors = 22799
Ave neighs/atom = 5.69975
Neighbor list builds = 28749
Dangerous builds = 0
Total wall time: 0:07:51

TABLES:

file: Table_lammps-potential_epsilon1_N20k_11.dat

# Table for repulsive harmonic potential. sigma_AA=1.0; sigma_AB=1.2; sigma_BB=1.4

AA

N 20000

1 0.100045 0.809919 1.799910
2 0.100090 0.809838 1.799820
3 0.100135 0.809757 1.799730
4 0.100180 0.809676 1.799640
5 0.100225 0.809595 1.799550
6 0.100270 0.809514 1.799460
7 0.100315 0.809433 1.799370
8 0.100360 0.809352 1.799280
9 0.100405 0.809271 1.799190
10 0.100450 0.809190 1.799100
......................
19991 0.999595 0.000000 0.000810
19992 0.999640 0.000000 0.000720
19993 0.999685 0.000000 0.000630
19994 0.999730 0.000000 0.000540
19995 0.999775 0.000000 0.000450
19996 0.999820 0.000000 0.000360
19997 0.999865 0.000000 0.000270
19998 0.999910 0.000000 0.000180
19999 0.999955 0.000000 0.000090
20000 1.000000 0.000000 0.000000

file: Table_lammps-potential_epsilon1_N20k_12.dat

# Table for repulsive harmonic potential. sigma_AA=1.0; sigma_AB=1.2; sigma_BB=1.4

AB

N 20000

1 0.100055 1.209879 2.199890
2 0.100110 1.209758 2.199780
3 0.100165 1.209637 2.199670
4 0.100220 1.209516 2.199560
5 0.100275 1.209395 2.199450
6 0.100330 1.209274 2.199340
7 0.100385 1.209153 2.199230
8 0.100440 1.209032 2.199120
9 0.100495 1.208911 2.199010
10 0.100550 1.208790 2.198900
.......................
19991 1.199505 0.000000 0.000990
19992 1.199560 0.000000 0.000880
19993 1.199615 0.000000 0.000770
19994 1.199670 0.000000 0.000660
19995 1.199725 0.000000 0.000550
19996 1.199780 0.000000 0.000440
19997 1.199835 0.000000 0.000330
19998 1.199890 0.000000 0.000220
19999 1.199945 0.000000 0.000110
20000 1.200000 0.000000 0.000000

file: Table_lammps-potential_epsilon1_N20k_22.dat

# Table for repulsive harmonic potential. sigma_AA=1.0; sigma_AB=1.2; sigma_BB=1.4

BB

N 20000

1 0.100065 1.689831 2.599870
2 0.100130 1.689662 2.599740
3 0.100195 1.689493 2.599610
4 0.100260 1.689324 2.599480
5 0.100325 1.689155 2.599350
6 0.100390 1.688986 2.599220
7 0.100455 1.688817 2.599090
8 0.100520 1.688648 2.598960
9 0.100585 1.688479 2.598830
10 0.100650 1.688311 2.598700
.......................
19991 1.399415 0.000000 0.001170
19992 1.399480 0.000000 0.001040
19993 1.399545 0.000000 0.000910
19994 1.399610 0.000000 0.000780
19995 1.399675 0.000000 0.000650
19996 1.399740 0.000000 0.000520
19997 1.399805 0.000000 0.000390
19998 1.399870 0.000000 0.000260
19999 1.399935 0.000000 0.000130
20000 1.400000 0.000000 0.000000
6

Thanks, please try to apply the following change:

  diff --git a/src/EXTRA-PAIR/pair_harmonic_cut.cpp b/src/EXTRA-PAIR/pair_harmonic_cut.cpp
  index f440909964..0cd9261aee 100644
  --- a/src/EXTRA-PAIR/pair_harmonic_cut.cpp
  +++ b/src/EXTRA-PAIR/pair_harmonic_cut.cpp
  @@ -189,10 +189,10 @@ double PairHarmonicCut::init_one(int i, int j)
   {
     if (setflag[i][j] == 0) {
       cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
  -    cut[j][i] = cut[i][j];
       k[i][j] = mix_energy(k[i][i], k[j][j], cut[i][i], cut[j][j]);
  -    k[j][i] = k[i][j];
     }
  +  k[j][i] = k[i][j];
  +  cut[j][i] = cut[i][j];
     return cut[i][j];
   }
7

Dear Axel,

That’s great, now it works perfectly fine also for the bi-disperse case.
Thank you so much.

Kind regards,
Fabio

1 Like
8

Excellent! Thanks for checking and reporting. Users like you, that are careful and don’t mind reporting inconsistencies, are crucial for the success and reliability of a large and collaborative software package like LAMMPS. I wish there were more.

The corresponding changes are already merged into the development branch and will soon be released in the next patch release of LAMMPS (probably next week).

9

The corresponding changes are already merged into the development branch and will soon be released in the next patch release of LAMMPS (probably next week).

Cool.

Thank you very much to be very responsive and very efficient in solving these problems.