[lammps-users] The number of atom types affects the simulation results

Lee_Miru · January 25, 2021, 4:32pm

Dear all

I’m designing a simulation where there are more than one types of atoms.
And I found this behavior that, given the exact same physical setup, setting different numbers of possible atom type through
“create_box” command gives me a different result.

To illustrate the problem that I’m encountering, here’s a much simpler system that still shows the problem.
I set up a simple cubic crystal system with two harmonic bonds connecting the nearest and the next nearest neighbors.
For pair_style, I’m using the following command:

pair_style hybrid lj/smooth/linear 3
pair_coeff * * none
pair_coeff 1 1 lj/smooth/linear 0 1

This command, regardless of how many atom types I have, I should see no pair interaction between atoms.

Then I call the following quantity:
compute peratom MD pe/atom
compute pe all reduce sum c_peratom

But the results of the above calculation are different when I’m setting
create_box 2 REGION bond/types 2 extra/bond/per/atom 18
create_atoms 1 box
and

create_box 1 REGION bond/types 2 extra/bond/per/atom 18

create_atoms 1 box

Here, although I’m stating that I may have two different atom types in my simulation,
I’m only using one atom type which is “1”.

The discrepancy goes away when I’m setting pair_coeff as

pair_coeff * * lj/smooth/linear 0 1

which gives me the same result as the case above where I’m stating only one atom type will be used.

Please elaborate me where the discrepancy is coming from.

Thanks in advance
Best,

akohlmey · January 25, 2021, 6:14pm

Several important pieces of information are missing.

how large is the discrepancy you are seeing?
what version of LAMMPS are you using?
proper input files to reproduce your observations rather than descriptions with individual lines of input, but apparently pieces missing.

Axel

akohlmey · January 25, 2021, 6:50pm

I see no difference when I run these with the current LAMMPS version 24 December 2020 and one CPU.

It would be much more important to see the output for the first steps and everything leading up to that. also whether you would be running in parallel etc.
to get consistent initialization of the velocities with different numbers of processors, you have to add the keyword/value “loop geom” to the “velocity” command

axel.

Lee_Miru · January 25, 2021, 6:33pm

Dear Axel

Thanks for the quick response.

My version is LAMMPS (20 Nov 2019).

For this simple setup, the two results are different by roughly 10 percent.

Regarding the system that I actually working with, the different is much larger.

Below I write down the input file and the output data.

The input file is following:

log LOG.0
shell “echo 0 $(date +%d-%b-%Y-%T) START >> JOBINFO”

units lj
dimension 3
boundary p p s
atom_style bond

CREATE INITIAL SYSTEM

lattice sc 1.0

region REGION block -0.5 7.500000000000e+00 -0.5 7.500000000000e+00 0 6
create_box 2 REGION bond/types 2 extra/bond/per/atom 18
create_atoms 1 box

variable layer atom “floor((z+0.001))”

FROZEN group: bottom layer

variable tmp atom v_layer==0
group FROZEN variable tmp
variable tmp delete

bond_style harmonic

SET LAMMPS RUNTIME LIMIT: clean exit after RUNTIME hours

timer timeout 23:00:00 every 10000

mass / md timestep

mass * 1.0
timestep 0.001

SET INTERACTIONS

pair_style hybrid lj/smooth/linear 3
pair_coeff * * none
pair_coeff * * lj/smooth/linear 0 1
comm_modify vel yes
neigh_modify every 10 delay 0 check yes

REPLACE BOND PARAMETERS WITH DESIRED ONES (!)

bond_coeff 1 10 1
bond_coeff 2 10 1.414213562373e+00

group SUB type 1 2
create_bonds many SUB SUB 1 9.900000000000e-01 1.010000000000e+00
create_bonds many SUB SUB 2 1.404213562373e+00 1.424213562373e+00

MD group: everything excluding bottom layer and beads

group MD subtract SUB FROZEN

INITIAL CONDITIONS / ASSIGN VARIABLES AND COMPUTES

assign initial velocities to substrate particles, remove any com motion / global rotation

velocity MD create 0.15 18618316 mom yes rot yes

fix MD MD nve

temperature of the MD group (standard ‘temp’ is too low, since this includes frozen particles)

compute MDTEMP MD temp
compute peratom MD pe/atom
compute pe all reduce sum c_peratom

write_data init.0

thermo 500
thermo_style custom step temp c_MDTEMP c_pe
thermo_modify norm no

run 5000

#write_data equil.0

PRODUCTION STAGE

run 100000
shell “echo 0 $(date +%d-%b-%Y-%T) FINISH >> JOBINFO”
quit

The corresponding output is (when setting the number of type = 2)

Step Temp c_MDTEMP c_pe

5000 0.067266531 0.078506892 40.152896
5500 0.067288419 0.078532437 39.913958
6000 0.066552635 0.077673702 40.45838
6500 0.069867469 0.081542451 38.228388
7000 0.070450536 0.082222949 37.884144
7500 0.068635944 0.080105136 39.053247
8000 0.067865416 0.079205851 39.781417
8500 0.07022641 0.081961371 38.002604
9000 0.068169544 0.0795608 39.397138
9500 0.068233379 0.079635301 39.403529
10000 0.070104774 0.081819409 38.194574
10500 0.0664537 0.077558234 40.321111
11000 0.069345706 0.0809335 38.531507
11500 0.067488085 0.078765467 40.137109
12000 0.070170954 0.081896649 38.187992
12500 0.069187902 0.080749326 39.085835
13000 0.069911441 0.08159377 38.623282
13500 0.071996657 0.08402743 37.099154
14000 0.073405097 0.085671223 36.223269
14500 0.071799285 0.083797077 37.312354
15000 0.076258406 0.089001325 34.49032
15500 0.071472916 0.08341617 37.257779
16000 0.072022699 0.084057823 36.857586
16500 0.071846726 0.083852446 37.003052
17000 0.067509739 0.078790739 39.740275
17500 0.071276402 0.083186819 37.311916
18000 0.071823357 0.083825171 37.313524
18500 0.073082271 0.085294452 36.207564
19000 0.073069549 0.085279605 36.153013
19500 0.072553933 0.084677828 36.606639
20000 0.071916984 0.083934444 36.781882
20500 0.072994523 0.085192041 36.31569
21000 0.072256629 0.084330843 36.784576
21500 0.074921674 0.087441222 35.13166

And when setting the number of type = 1

Step Temp c_MDTEMP c_pe

5000 0.073304528 0.085553849 35.931852
5500 0.060067827 0.070105271 44.717848
6000 0.065163119 0.076051995 41.547947
6500 0.068057504 0.079430037 39.616606
7000 0.056450572 0.065883565 46.982729
7500 0.064316954 0.075064434 41.808312
8000 0.061735318 0.072051403 43.832576
8500 0.065540897 0.076492901 41.077493
9000 0.071118584 0.083002629 37.581507
9500 0.066207165 0.077270503 40.380476
10000 0.060717219 0.070863177 44.238874
10500 0.060347656 0.07043186 44.20803
11000 0.064025711 0.074724524 42.269167
11500 0.064187221 0.074913023 41.862736
12000 0.063765906 0.074421305 42.115145
12500 0.061483035 0.071756962 43.972857
13000 0.05863742 0.06843584 45.711041
13500 0.060020829 0.070050419 44.864965
14000 0.06643381 0.077535021 40.622482
14500 0.052738521 0.061551225 49.589298
15000 0.061248188 0.071482872 43.980723
15500 0.064726856 0.075542831 41.323342
16000 0.065004547 0.075866926 41.211363
16500 0.061268131 0.071506148 43.726626
17000 0.063068894 0.073607821 42.827017
17500 0.063384651 0.073976342 42.397989
18000 0.064769322 0.075592394 41.456778
18500 0.059161239 0.06904719 45.299969
19000 0.060616103 0.070745165 44.432415
19500 0.059887791 0.06989515 44.743213
20000 0.069074152 0.080616568 38.822442
20500 0.065382365 0.076307877 41.210434
21000 0.067379641 0.078638902 39.954273
21500 0.061553908 0.071839678 43.702575

Please look at c_pe which is roughly 10% higher for the latter case.

Best,
Miru

Lee_Miru · January 26, 2021, 9:57am

Dear Axel

Thanks for the input.
Now I’m running the latest version 24 Dec 2020 and still have the same issue.

I’ve also added the “loop geom” option to the “velocity” command.
I see no improvement when using only one core too.

But I do see different initializations which might be the culprit.

when setting only one possible atom type I see non-zero potential at the 0th time step:

Step Temp c_MDTEMP c_pe
0 0.12852349 0.15 1.9652911e-22

Performance: 246110.405 tau/day, 2848.500 timesteps/s

whereas when setting two or more possible atom types I see exactly zero potential:

Step Temp c_MDTEMP c_pe
0 0.12852349 0.15 0

Performance: 87748.477 tau/day, 1015.607 timesteps/s

Plus, the simulation runs significantly slower.

Best,
Miru

akohlmey · January 26, 2021, 11:30am

I need more details to track this down. I can only resolve problems if I can reproduce them. and so far I am unable to.

what is the platform you are running on (OS, hardware)? what is the compiler and what are the compiler settings used to compile LAMMPS?
what is the output if you run LAMMPS with just the one command. “info config communication”?
what is the command line you are using to run each input.

also, please send to me personally (not to the mailing list):
exact copies of the two input files you are using
copies of the two log files that you generate with the different type settings

thanks,
axel.

akohlmey · January 26, 2021, 1:04pm

thanks for the files.

there are two main inconsistencies:

you have different random number seeds for the velocity command in the in.1 versus the in.2 input
that will lead to diverging but still valid trajectories
you are using a different force field definition. in the input you sent to the mailing list you had
pair_style hybrid lj/smooth/linear 3
pair_coeff * * none
pair_coeff * * lj/smooth/linear 0 1

while in the files you sent you have
pair_style hybrid lj/smooth/linear 3
pair_coeff * * none
pair_coeff 1 1 lj/smooth/linear 0 1

that is causing different neighbor list constructions and it seems there may be a bug in the neighborlist code because of that, which then also affects the bond creation.
the different bond creation then changes the number of actual pairs in the neighbor list, as suddenly fewer bonds are created which means that non-bonded interactions are no longer completely excluded. while the computation of those interactions will not result in forces, the loop over the pairs and computing the distances to determine whether they are within the cutoff range and then computing potential energies that are then multiplied with zero is what is causing the slower performance.

your workaround
pair_style hybrid lj/smooth/linear 3
pair_coeff * * lj/smooth/linear 0 1

is the same as the first case.

I will have to have somebody else look into whether the neighborlist situation is correct or not.
as you know, your system is highly artificial, so you have to look into what exactly is causing the issue in your real world example.
but due to the inconsistencies in what you have been sending to the mailing list and the files you sent to me, there is a good chance that you have inconsistent inputs and the behavior of LAMMPS is correct. I cannot tell for certain based on the information I have. You need to re-check everything you did with more care.

axel.

Lee_Miru · January 26, 2021, 4:57pm

Axel,

Thanks for the detailed explanation.

The first inconsistency is minor, although I better had to control it, since I have a mother code from which I derive a code with small modification.
The mother code assigns a random number to the seed whenever I create a son code.
Each run I create a new son code, and I consistently see this issue for different runs.

The second inconsistency is my mistake.
I copied the code from another code without my noticing, while the outputs that follow after are indeed from the commands with the issue:

pair_style hybrid lj/smooth/linear 3
pair_coeff * * none
pair_coeff 1 1 lj/smooth/linear 0 1

As you speculated, the issue seems to be in constructing the neighbor list.
The discrepancy in the results goes away when I build the system with the following command first

pair_style hybrid lj/smooth/linear 3
pair_coeff * * none
pair_coeff * * lj/smooth/linear 0 1

And then reconfigure the pair_coeff as I have intended after creating bonds.

This way the slowing down issue also disappear.

Anyways, this has been tremendously helpful, thank you.

Best,
Miru

akohlmey · February 11, 2021, 4:23pm

As a follow up to this discussion, there are two important comments.

We have looked at this case carefully and identified two bugs in the LAMMPS code that are now fixed in the 10 February 2021 patch release.
Your input is incorrect and was creating bonds in the second create_bonds command where it should not. This is because of one of the bugs in 1). You would need to change the special_bonds setting from the default of 0.0 0.0 0.0 to lj/coul 0.0 1.0 1.0 for the second command to work correctly. Please see the note on this subject in the updated manual for the new patch release.

Axel.

Lee_Miru · February 11, 2021, 6:10pm

Axel,

This is great!!
Thank you very much for your work and support
I also thank those who have looked into this issue and addressed it.

Best,
Miru

akohlmey · February 11, 2021, 6:49pm

The credit for figuring this out (and it must have taken considerable effort) goes to “Mr. LAMMPS” himself, Steve Plimpton.