A kolmogorov/crespi/full LAMMPS question (from a mathematician)

Dear Colleagues,

We are developing a discrete-to-continuum procedure that will allow to model a graphene bilayer as a continuum sheet. We are able to compare continuum and discrete simulations results for fake simple systems (e.g., in-plane harmonic extensional, torsional, and dihedral springs + LJ or kolmogorov/crespi/z interactions between the layers) at zero temperature and the results look pretty good. The next step is to increase the degree of complexity; to this end we would like to investigate the kolmogorov/crespi/full potential, while for the moment continuing to use harmonic potentials to model other interactions. Even though the code below works for kc/z, the obvious modifications that I attempted to adjust these files for kc/full do not seem to work. My questions are:

(1) Is there anything in the kc/full potential that would prevent it from working with harmonic potentials in place of tersoff or airebo?

(2) If not, what modifications need to be made to the files below to make them agreeable to LAMMPS?

Note that, at this point, I do not need simulations to result in anything physically realistic hence the numbers I am using are bogus.

Many thanks in advance for your help.

Dmitry

-------------- in.graphene-----------
dimension 3
boundary s s s

atom_style full
neighbor 0.3 bin
neigh_modify delay 5

read_data generate_initial/initial_data.txt

mass * 1.0

Potentials

pair_style hybrid/overlay kolmogorov/crespi/z 12.0
pair_coeff * * none
pair_coeff 12 12 kolmogorov/crespi/z CC.KC C C

bond_style harmonic
bond_coeff * 5.E-1 1.42

angle_style harmonic
angle_coeff * 3.625E-2 120

dihedral_style harmonic
dihedral_coeff * 2.52E-2 -1 2

define groups

group lower type 2
group mobile type 1

initial velocities

compute new mobile temp
velocity mobile zero linear
fix 1 mobile nvt temp 1E-3 1E-3 1E-3
fix 2 lower setforce 0.0 0.0 0.0
fix 3 mobile viscous 0.1
#fix 3 all temp/rescale 100 0.1 0.1 0.01 1.0

timestep 0.003

thermo 1000
thermo_modify temp new

dump 1 mobile xyz 1000 graphene.xyz
dump 2 lower xyz 1000 substrate.xyz

run 1000000

print “Simulation complete”
--------------CC.KC----------

Kolmogorov-Crespi Potential

Cite as A.N. Kolmogorov & V. H. Crespi,

Registry-dependent interlayer potential for graphitic systems

Physical Review B 71, 235415 (2005)

z0 C0 C2 C4 C delta lambda A S

C C 3.34 15.071 12.029 4.0933 3.0030 0.578 3.629 10.238 1.0

No.

Impossible to say without seeing the error message you are getting.

P.S.: when adding an input with copy-n-paste, please enclose the block with triple backquotes “```” so that the forum software will not apply the regular typesetting rules which will render parts of it unreadable.

Thanks, Axel - using the files below, makes LAMMPS run without errors, but the mobile group of atoms remains completely stationary; it does deform for a kc/z simulation that uses the same parameters. Clearly, I am setting up kc/full incorrectly (e.g., neigh below is zero) but I am not sure where the mistake is.

---------------in.graphene-------------

boundary	s s s

atom_style	full
neighbor	0.3 bin
neigh_modify delay 5

read_data generate_initial/initial_data.txt

mass	* 1.0

# Potentials

pair_style hybrid/overlay kolmogorov/crespi/full 12.0 0
pair_coeff * * none
pair_coeff * * kolmogorov/crespi/full  CC.KC   C C

bond_style  	harmonic
bond_coeff  	* 5.E-1 1.42

angle_style	harmonic
angle_coeff	* 5.625E-2 120

dihedral_style	harmonic
dihedral_coeff	* 4.52E-2 -1 2

# define groups

group		lower type 2
group		mobile type 1

# initial velocities

compute	  	new mobile temp
velocity	mobile zero linear 
fix		1 mobile nvt temp 1E-3 1E-3 1E-3
fix		2 lower setforce 0.0 0.0 0.0
fix		3 mobile viscous 0.1
#fix		3 all temp/rescale 100 0.1 0.1 0.01 1.0

timestep	0.003

thermo		1000
thermo_modify	temp new

dump 		1 mobile xyz 1000 graphene.xyz
dump        	2 lower xyz 1000 substrate.xyz

run	   1000000

print  "Simulation complete"```
---------------------CC.KC-------------------
```# Kolmogorov-Crespi Potential
#
# Cite as A.N. Kolmogorov & V. H. Crespi, 
# Registry-dependent interlayer potential for graphitic systems
# Physical Review B 71, 235415 (2005)
# 
#      z0     C0     C2      C4       C    delta  lambda      A     S
C  C  3.34  15.071  12.029   4.0933   3.0030   0.578   3.629 10.238  1.0 0.0```
----------------------------------------------------
```LAMMPS (10 Feb 2021)
  using 8 OpenMP thread(s) per MPI task
# twisted bilayer graphene

dimension	3
boundary	s s s

atom_style	full
neighbor	0.3 bin
neigh_modify delay 5

read_data generate_initial/initial_data.txt
Reading data file ...
  orthogonal box = (-2000.0000 -2000.0000 -10.000000) to (2000.0000 2000.0000 10.000000)
  2 by 4 by 1 MPI processor grid
  reading atoms ...
  12864 atoms
  scanning bonds ...
  2 = max bonds/atom
  scanning angles ...
  3 = max angles/atom
  scanning dihedrals ...
  8 = max dihedrals/atom
  reading bonds ...
  5815 bonds
  reading angles ...
  11454 angles
  reading dihedrals ...
  22605 dihedrals
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     3 = max # of 1-2 neighbors
     6 = max # of 1-3 neighbors
    18 = max # of 1-4 neighbors
    18 = max # of special neighbors
  special bonds CPU = 0.002 seconds
  read_data CPU = 0.069 seconds

mass	* 1.0

# Potentials

pair_style hybrid/overlay kolmogorov/crespi/full 12.0 0
pair_coeff * * none
pair_coeff * * kolmogorov/crespi/full  CC.KC   C C

bond_style  	harmonic
bond_coeff  	* 5.E-1 1.42

angle_style	harmonic
angle_coeff	* 5.625E-2 120

dihedral_style	harmonic
dihedral_coeff	* 4.52E-2 -1 2

# define groups

group		lower type 2
8928 atoms in group lower
group		mobile type 1
3936 atoms in group mobile

# initial velocities

compute	  	new mobile temp
velocity	mobile zero linear
fix		1 mobile nvt temp 1E-3 1E-3 1E-3
fix		2 lower setforce 0.0 0.0 0.0
fix		3 mobile viscous 0.1
#fix		3 all temp/rescale 100 0.1 0.1 0.01 1.0

timestep	0.003

thermo		1000
thermo_modify	temp new
WARNING: Temperature for thermo pressure is not for group all (src/thermo.cpp:472)

dump 		1 mobile xyz 1000 graphene.xyz
dump        	2 lower xyz 1000 substrate.xyz

run	   10000

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

Your simulation uses code contributions which should be cited:

- @Article{Ouyang2018
 author = {W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod},
 title = {Nanoserpents: Graphene Nanoribbon Motion on Two-Dimensional Hexagonal Materials},
 journal = {Nano Letters},
 volume =  18,
 pages =   {6009}
 year =    2018,
}

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

Neighbor list info ...
  update every 1 steps, delay 5 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 12.3
  ghost atom cutoff = 12.3
  binsize = 6.15, bins = 26 26 1
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair kolmogorov/crespi/full, perpetual
      attributes: full, newton on, ghost
      pair build: full/bin/ghost
      stencil: full/ghost/bin/3d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 20.78 | 22.17 | 23.57 Mbytes
Step Temp E_pair E_mol TotEng Press Volume 
       0            0            0 1.2537168e-12 1.2537168e-12 -4.0493317e-10    85037.702 
    1000            0            0 3.9775833e-09 3.9775833e-09 -1.9708938e-08    85037.702 
    2000            0            0 3.8654076e-09 3.8654076e-09 -1.9571668e-08    85037.702 
    3000            0            0 3.7564848e-09 3.7564848e-09 -1.9435645e-08    85037.702 
    4000            0            0 3.6507187e-09 3.6507187e-09 -1.9300853e-08    85037.702 
    5000            0            0 3.548016e-09 3.548016e-09 -1.9167276e-08    85037.702 
    6000            0            0 3.448286e-09 3.448286e-09 -1.9034897e-08    85037.702 
    7000            0            0 3.3514411e-09 3.3514411e-09 -1.89037e-08    85037.702 
    8000            0            0 3.2573961e-09 3.2573961e-09 -1.877367e-08    85037.702 
    9000            0            0 3.1660683e-09 3.1660683e-09 -1.864479e-08    85037.702 
   10000            0            0 3.0773777e-09 3.0773777e-09 -1.8517047e-08    85037.702 
Loop time of 92.5367 on 64 procs for 10000 steps with 12864 atoms

Performance: 28010.517 tau/day, 108.065 timesteps/s
81.2% CPU use with 8 MPI tasks x 8 OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 40.457     | 63.365     | 86.846     | 282.2 | 68.48
Bond    | 0.84309    | 2.0935     | 3.353      |  86.1 |  2.26
Neigh   | 0          | 0          | 0          |   0.0 |  0.00
Comm    | 1.1131     | 25.656     | 49.445     | 464.1 | 27.73
Output  | 0.02486    | 0.024922   | 0.025005   |   0.0 |  0.03
Modify  | 0.61152    | 0.91839    | 1.4275     |  27.8 |  0.99
Other   |            | 0.4792     |            |       |  0.52

Nlocal:        1608.00 ave        1957 max        1260 min
Histogram: 4 0 0 0 0 0 0 0 0 4
Nghost:        1358.75 ave        1785 max         935 min
Histogram: 4 0 0 0 0 0 0 0 0 4
Neighs:         0.00000 ave           0 max           0 min
Histogram: 8 0 0 0 0 0 0 0 0 0
FullNghs:      423899.0 ave      584167 max      263741 min
Histogram: 4 0 0 0 0 0 0 0 0 4

Total # of neighbors = 3391190
Ave neighs/atom = 263.61863
Ave special neighs/atom = 5.3224502
Neighbor list builds = 0
Dangerous builds = 0

print  "Simulation complete"
Simulation complete
Total wall time: 0:01:32```

Please study the documentation of the pair style very carefully. My guess is that you have overlooked the note about assigning different molecule IDs to atoms in the different layers to identify which atom belongs to which layer.

Thanks again, Axel,

The last question: under documentation you mean this pair_style kolmogorov/crespi/full command — LAMMPS documentation
or there is a fuller description? If there is one, I was not able to find it.

Dmitry

This potential (ILP) is intended for interlayer interactions between two different layers of graphene. To perform a realistic simulation, this potential must be used in combination with intralayer potential, such as AIREBO or Tersoff potential. To keep the intralayer properties unaffected, the interlayer interaction within the same layers should be avoided. Hence, each atom has to have a layer identifier such that atoms residing on the same layer interact via the appropriate intralayer potential and atoms residing on different layers interact via the ILP. Here, the molecule id is chosen as the layer identifier, thus a data file with the “full” atom style is required to use this potential.

This is the first “Note” on the page and using the molecule ID is mentioned in the last sentence.

OK - thank you.

Dmitry

You were correct - assigning different molecule IDs solved the problem - thanks again for your help!

Dmitry