tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3125 from oywg11/ilp-for-tmds

Browse Source 
Anisotropic Interlayer potential for TMDs and Metal Surfaces
This commit is contained in:
Axel Kohlmeyer
2022-02-17 11:42:15 -05:00
committed by GitHub
parent 6d2f4343d9 8b627f92f0
commit 45b92e519d
31 changed files with 4298 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_pair.rst
View File

@ -124,6 +124,7 @@ OPT.
* :doc:`hbond/dreiding/morse (o) <pair_hbond_dreiding>`
* :doc:`hdnnp <pair_hdnnp>`
* :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>`
* :doc:`ilp/tmd <pair_ilp_tmd>`
* :doc:`kolmogorov/crespi/full <pair_kolmogorov_crespi_full>`
* :doc:`kolmogorov/crespi/z <pair_kolmogorov_crespi_z>`
* :doc:`lcbop <pair_lcbop>`
@ -241,6 +242,7 @@ OPT.
* :doc:`reaxff (ko) <pair_reaxff>`
* :doc:`rebo (io) <pair_airebo>`
* :doc:`resquared (go) <pair_resquared>`
* :doc:`saip/metal <pair_saip_metal>`
* :doc:`sdpd/taitwater/isothermal <pair_sdpd_taitwater_isothermal>`
* :doc:`smd/hertz <pair_smd_hertz>`
* :doc:`smd/tlsph <pair_smd_tlsph>`

2
doc/src/pair_ilp_graphene_hbn.rst
View File

@ -159,6 +159,8 @@ Related commands
:doc:`pair_none <pair_none>`,
:doc:`pair_style hybrid/overlay <pair_hybrid>`,
:doc:`pair_style drip <pair_drip>`,
:doc:`pair_style ilp_tmd <pair_ilp_tmd>`,
:doc:`pair_style saip_metal <pair_saip_metal>`,
:doc:`pair_style pair_kolmogorov_crespi_z <pair_kolmogorov_crespi_z>`,
:doc:`pair_style pair_kolmogorov_crespi_full <pair_kolmogorov_crespi_full>`,
:doc:`pair_style pair_lebedeva_z <pair_lebedeva_z>`,

157
doc/src/pair_ilp_tmd.rst Normal file
View File

@ -0,0 +1,157 @@
.. index:: pair_style ilp/tmd
pair_style ilp/tmd command
===================================
Syntax
""""""
.. code-block:: LAMMPS
pair_style [hybrid/overlay ...] ilp/tmd cutoff tap_flag
* cutoff = global cutoff (distance units)
* tap_flag = 0/1 to turn off/on the taper function
Examples
""""""""
.. code-block:: LAMMPS
pair_style hybrid/overlay ilp/tmd 16.0 1
pair_coeff * * ilp/tmd TMD.ILP Mo S S
pair_style hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
pair_coeff * * sw/mod 1 tmd.sw.mod Mo S S NULL NULL NULL
pair_coeff * * sw/mod 2 tmd.sw.mod NULL NULL NULL Mo S S
pair_coeff * * ilp/tmd TMD.ILP Mo S S Mo S S
Description
"""""""""""
The *ilp/tmd* style computes the registry-dependent interlayer
potential (ILP) potential for transition metal dichalcogenides (TMD)
as described in :ref:`(Ouyang7) <Ouyang7>`.
.. math::
E = & \frac{1}{2} \sum_i \sum_{j \neq i} V_{ij} \\
V_{ij} = & {\rm Tap}(r_{ij})\left \{ e^{-\alpha (r_{ij}/\beta -1)}
\left [ \epsilon + f(\rho_{ij}) + f(\rho_{ji})\right ] -
\frac{1}{1+e^{-d\left [ \left ( r_{ij}/\left (s_R \cdot r^{eff} \right ) \right )-1 \right ]}}
\cdot \frac{C_6}{r^6_{ij}} \right \}\\
\rho_{ij}^2 = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_i)^2 \\
\rho_{ji}^2 = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_j)^2 \\
f(\rho) = & C e^{ -( \rho / \delta )^2 } \\
{\rm Tap}(r_{ij}) = & 20\left ( \frac{r_{ij}}{R_{cut}} \right )^7 -
70\left ( \frac{r_{ij}}{R_{cut}} \right )^6 +
84\left ( \frac{r_{ij}}{R_{cut}} \right )^5 -
35\left ( \frac{r_{ij}}{R_{cut}} \right )^4 + 1
Where :math:`\mathrm{Tap}(r_{ij})` is the taper function which provides
a continuous cutoff (up to third derivative) for interatomic separations
larger than :math:`r_c` :doc:`pair_style ilp_graphene_hbn <pair_ilp_graphene_hbn>`.
It is important to include all the pairs to build the neighbor list for
calculating the normals.
.. note::
Since each MX2 (M = Mo, W and X = S, Se Te) layer contains two
sub-layers of X atoms and one sub-layer of M atoms, the definition of the
normal vectors used for graphene and h-BN is no longer valid for TMDs.
In :ref:`(Ouyang7) <Ouyang7>`, a new definition is proposed, where for
each atom `i`, its six nearest neighboring atoms belonging to the same
sub-layer are chosen to define the normal vector `{\bf n}_i`.
The parameter file (e.g. TMD.ILP), is intended for use with *metal*
:doc:`units <units>`, with energies in meV. Two additional parameters,
*S*, and *rcut* are included in the parameter file. *S* is designed to
facilitate scaling of energies. *rcut* is designed to build the neighbor
list for calculating the normals for each atom pair.
.. note::
The parameters presented in the parameter file (e.g. TMD.ILP),
are fitted with taper function by setting the cutoff equal to 16.0
Angstrom. Using different cutoff or taper function should be careful.
These parameters provide a good description in both short- and long-range
interaction regimes. This feature is essential for simulations in high pressure
regime (i.e., the interlayer distance is smaller than the equilibrium
distance). The benchmark tests and comparison of these parameters can
be found in :ref:`(Ouyang7) <Ouyang7>`.
This potential must be used in combination with hybrid/overlay.
Other interactions can be set to zero using pair_style *none*\ .
This pair style tallies a breakdown of the total interlayer potential
energy into sub-categories, which can be accessed via the :doc:`compute pair <compute_pair>` command as a vector of values of length 2.
The 2 values correspond to the following sub-categories:
1. *E_vdW* = vdW (attractive) energy
2. *E_Rep* = Repulsive energy
To print these quantities to the log file (with descriptive column
headings) the following commands could be included in an input script:
.. code-block:: LAMMPS
compute 0 all pair ilp/tmd
variable Evdw equal c_0[1]
variable Erep equal c_0[2]
thermo_style custom step temp epair v_Erep v_Evdw
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This pair style does not support the pair_modify mix, shift, table, and
tail options.
This pair style does not write their information to binary restart
files, since it is stored in potential files. Thus, you need to
re-specify the pair_style and pair_coeff commands in an input script
that reads a restart file.
Restrictions
""""""""""""
This pair style is part of the INTERLAYER package. It is only enabled
if LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
This pair style requires the newton setting to be *on* for pair
interactions.
The TMD.ILP potential file provided with LAMMPS (see the potentials
directory) are parameterized for *metal* units. You can use this
potential with any LAMMPS units, but you would need to create your
BNCH.ILP potential file with coefficients listed in the appropriate
units, if your simulation does not use *metal* units.
Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`,
:doc:`pair_none <pair_none>`,
:doc:`pair_style hybrid/overlay <pair_hybrid>`,
:doc:`pair_style drip <pair_drip>`,
:doc:`pair_style saip_metal <pair_saip_metal>`,
:doc:`pair_style ilp_graphene_hbn <pair_ilp_graphene_hbn>`,
:doc:`pair_style pair_kolmogorov_crespi_z <pair_kolmogorov_crespi_z>`,
:doc:`pair_style pair_kolmogorov_crespi_full <pair_kolmogorov_crespi_full>`,
:doc:`pair_style pair_lebedeva_z <pair_lebedeva_z>`,
:doc:`pair_style pair_coul_shield <pair_coul_shield>`.
Default
"""""""
tap_flag = 1
----------
.. _Ouyang7:
**(Ouyang7)** W. Ouyang, et al., J. Chem. Theory Comput. 17, 7237 (2021).

156
doc/src/pair_saip_metal.rst Normal file
View File

@ -0,0 +1,156 @@
.. index:: pair_style saip/metal
pair_style saip/metal command
===================================
Syntax
""""""
.. code-block:: LAMMPS
pair_style [hybrid/overlay ...] saip/metal cutoff tap_flag
* cutoff = global cutoff (distance units)
* tap_flag = 0/1 to turn off/on the taper function
Examples
""""""""
.. code-block:: LAMMPS
pair_style hybrid/overlay saip/metal 16.0 1
pair_coeff * * saip/metal CHAu.ILP Au C H
pair_style hybrid/overlay eam rebo saip/metal 16.0
pair_coeff 1 1 eam Au_u3.eam Au NULL NULL
pair_coeff * * rebo CH.rebo NULL C H
pair_coeff * * saip/metal CHAu.ILP Au C H
Description
"""""""""""
The *saip/metal* style computes the registry-dependent interlayer
potential (ILP) potential for hetero-junctions formed with hexagonal
2D materials and metal surfaces, as described in :ref:`(Ouyang6) <Ouyang6>`.
.. math::
E = & \frac{1}{2} \sum_i \sum_{j \neq i} V_{ij} \\
V_{ij} = & {\rm Tap}(r_{ij})\left \{ e^{-\alpha (r_{ij}/\beta -1)}
\left [ \epsilon + f(\rho_{ij}) + f(\rho_{ji})\right ] -
\frac{1}{1+e^{-d\left [ \left ( r_{ij}/\left (s_R \cdot r^{eff} \right ) \right )-1 \right ]}}
\cdot \frac{C_6}{r^6_{ij}} \right \}\\
\rho_{ij}^2 = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_i)^2 \\
\rho_{ji}^2 = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_j)^2 \\
f(\rho) = & C e^{ -( \rho / \delta )^2 } \\
{\rm Tap}(r_{ij}) = & 20\left ( \frac{r_{ij}}{R_{cut}} \right )^7 -
70\left ( \frac{r_{ij}}{R_{cut}} \right )^6 +
84\left ( \frac{r_{ij}}{R_{cut}} \right )^5 -
35\left ( \frac{r_{ij}}{R_{cut}} \right )^4 + 1
Where :math:`\mathrm{Tap}(r_{ij})` is the taper function which provides
a continuous cutoff (up to third derivative) for interatomic separations
larger than :math:`r_c` :doc:`pair_style ilp_graphene_hbn <pair_ilp_graphene_hbn>`.
It is important to include all the pairs to build the neighbor list for
calculating the normals.
.. note::
To account for the isotropic nature of the isolated gold atom
electron cloud, their corresponding normal vectors (`{\bf n}_i`) are
assumed to lie along the interatomic vector `{\bf r}_ij`. Notably, this
assumption is suitable for many bulk material surfaces, for
example, for systems possessing s-type valence orbitals or
metallic surfaces, whose valence electrons are mostly
delocalized, such that their Pauli repulsion with the electrons
of adjacent surfaces are isotropic. Caution should be used in
the case of very small gold contacts, for example, nano-clusters,
where edge effects may become relevant.
The parameter file (e.g. CHAu.ILP), is intended for use with *metal*
:doc:`units <units>`, with energies in meV. Two additional parameters,
*S*, and *rcut* are included in the parameter file. *S* is designed to
facilitate scaling of energies. *rcut* is designed to build the neighbor
list for calculating the normals for each atom pair.
.. note::
The parameters presented in the parameter file (e.g. BNCH.ILP),
are fitted with taper function by setting the cutoff equal to 16.0
Angstrom. Using different cutoff or taper function should be careful.
This potential must be used in combination with hybrid/overlay.
Other interactions can be set to zero using pair_style *none*\ .
This pair style tallies a breakdown of the total interlayer potential
energy into sub-categories, which can be accessed via the :doc:`compute pair <compute_pair>` command as a vector of values of length 2.
The 2 values correspond to the following sub-categories:
1. *E_vdW* = vdW (attractive) energy
2. *E_Rep* = Repulsive energy
To print these quantities to the log file (with descriptive column
headings) the following commands could be included in an input script:
.. code-block:: LAMMPS
compute 0 all pair saip/metal
variable Evdw equal c_0[1]
variable Erep equal c_0[2]
thermo_style custom step temp epair v_Erep v_Evdw
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This pair style does not support the pair_modify mix, shift, table, and
tail options.
This pair style does not write their information to binary restart
files, since it is stored in potential files. Thus, you need to
re-specify the pair_style and pair_coeff commands in an input script
that reads a restart file.
Restrictions
""""""""""""
This pair style is part of the INTERLAYER package. It is only enabled
if LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
This pair style requires the newton setting to be *on* for pair
interactions.
The CHAu.ILP potential file provided with LAMMPS (see the potentials
directory) are parameterized for *metal* units. You can use this
potential with any LAMMPS units, but you would need to create your
BNCH.ILP potential file with coefficients listed in the appropriate
units, if your simulation does not use *metal* units.
Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`,
:doc:`pair_none <pair_none>`,
:doc:`pair_style hybrid/overlay <pair_hybrid>`,
:doc:`pair_style drip <pair_drip>`,
:doc:`pair_style ilp_tmd <pair_ilp_tmd>`,
:doc:`pair_style ilp_graphene_hbn <pair_ilp_graphene_hbn>`,
:doc:`pair_style pair_kolmogorov_crespi_z <pair_kolmogorov_crespi_z>`,
:doc:`pair_style pair_kolmogorov_crespi_full <pair_kolmogorov_crespi_full>`,
:doc:`pair_style pair_lebedeva_z <pair_lebedeva_z>`,
:doc:`pair_style pair_coul_shield <pair_coul_shield>`.
Default
"""""""
tap_flag = 1
----------
.. _Ouyang6:
**(Ouyang6)** W. Ouyang, O. Hod, and R. Guerra, J. Chem. Theory Comput. 17, 7215 (2021).

2
doc/src/pair_style.rst
View File

@ -188,6 +188,7 @@ accelerated styles exist.
* :doc:`hbond/dreiding/morse <pair_hbond_dreiding>` - DREIDING hydrogen bonding Morse potential
* :doc:`hdnnp <pair_hdnnp>` - High-dimensional neural network potential
* :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>` - registry-dependent interlayer potential (ILP)
* :doc:`ilp/tmd <pair_ilp_tmd>` - interlayer potential (ILP) potential for transition metal dichalcogenides (TMD)
* :doc:`kim <pair_kim>` - interface to potentials provided by KIM project
* :doc:`kolmogorov/crespi/full <pair_kolmogorov_crespi_full>` - Kolmogorov-Crespi (KC) potential with no simplifications
* :doc:`kolmogorov/crespi/z <pair_kolmogorov_crespi_z>` - Kolmogorov-Crespi (KC) potential with normals along z-axis
@ -305,6 +306,7 @@ accelerated styles exist.
* :doc:`reaxff <pair_reaxff>` - ReaxFF potential
* :doc:`rebo <pair_airebo>` - second generation REBO potential of Brenner
* :doc:`resquared <pair_resquared>` - Everaers RE-Squared ellipsoidal potential
* :doc:`saip/metal <pair_saip_metal>` - interlayer potential for hetero-junctions formed with hexagonal 2D materials and metal surfaces
* :doc:`sdpd/taitwater/isothermal <pair_sdpd_taitwater_isothermal>` - smoothed dissipative particle dynamics for water at isothermal conditions
* :doc:`smd/hertz <pair_smd_hertz>` -
* :doc:`smd/tlsph <pair_smd_tlsph>` -

8
doc/src/pair_sw.rst
View File

@ -29,8 +29,8 @@ Syntax
.. parsed-literal::
*maxdelcs* value = delta1 delta2 (optional)
delta1 = The minimum thershold for cosine of three-body angle
delta2 = The maximum threshold for cosine of three-body angle
delta1 = The minimum thershold for the variation of cosine of three-body angle
delta2 = The maximum threshold for the variation of cosine of three-body angle
Examples
""""""""
@ -69,7 +69,7 @@ and K of atom I within a cutoff distance :math:`a `\sigma`.
The *sw/mod* style is designed for simulations of materials when
distinguishing three-body angles are necessary, such as borophene
and transition metal dichalcogenide, which cannot be described
and transition metal dichalcogenides, which cannot be described
by the original code for the Stillinger-Weber potential.
For instance, there are several types of angles around each Mo atom in `MoS_2`,
and some unnecessary angle types should be excluded in the three-body interaction.
@ -99,7 +99,7 @@ This smoothly turns off the energy and force contributions for :math:`\left| \de
It is suggested that :math:`\delta 1` and :math:`\delta_2` to be the value around
:math:`0.5 \left| \cos \theta_1 - \cos \theta_2 \right|`, with
:math:`\theta_1` and :math:`\theta_2` as the different types of angles around an atom.
For borophene and transition metal dichalcogenide, :math:`\delta_1 = 0.25` and :math:`\delta_2 = 0.35`.
For borophene and transition metal dichalcogenides, :math:`\delta_1 = 0.25` and :math:`\delta_2 = 0.35`.
This value enables the cut-off function to exclude unnecessary angles in the three-body SW terms.
.. note::

3
doc/utils/sphinx-config/false_positives.txt
View File

@ -692,7 +692,7 @@ diagonalizers
diagonalizing
Diallo
diblock
dichalcogenide
dichalcogenides
Dickel
diel
Dietz
@ -2954,6 +2954,7 @@ safezone
Safran
Sagui
Saidi
saip
Salanne
Salles
sandia

1
examples/PACKAGES/interlayer/ilp_graphene_hbn/BNCH.ILP Symbolic link
View File

@ -0,0 +1 @@
../../../../potentials/BNCH.ILP

1
examples/PACKAGES/interlayer/ilp_graphene_hbn/CH.rebo Symbolic link
View File

@ -0,0 +1 @@
../../../../potentials/CH.rebo

1
examples/PACKAGES/interlayer/ilp_tmds/MoS2.ILP Symbolic link
View File

@ -0,0 +1 @@
../../../../potentials/MoS2.ILP

1549
examples/PACKAGES/interlayer/ilp_tmds/bilayer_MoS2_AAprime_stacking.data Normal file
View File

File diff suppressed because it is too large Load Diff

36
examples/PACKAGES/interlayer/ilp_tmds/in.mos2 Normal file
View File

@ -0,0 +1,36 @@
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./bilayer_MoS2_AAprime_stacking.data
mass * 32.065 # mass of sulphur atom , uint: a.u.=1.66X10^(-27)kg
mass 1 95.94 # mass of molebdenum atom , uint: a.u.=1.66X10^(-27)kg
mass 4 95.94
# Define potentials
pair_style hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
pair_coeff * * sw/mod 1 tmd.sw.mod Mo S S NULL NULL NULL
pair_coeff * * sw/mod 2 tmd.sw.mod NULL NULL NULL Mo S S
pair_coeff * * ilp/tmd MoS2.ILP Mo S S Mo S S
# Calculate the pair potential
compute 0 all pair ilp/tmd
compute 1 all pair sw/mod 1
compute 2 all pair sw/mod 2
variable SW1 equal c_1
variable SW2 equal c_2
variable ILP equal c_0
variable Eatt equal c_0[1]
variable Erep equal c_0[2]
thermo_style custom step etotal pe ke v_SW1 v_SW2 v_ILP v_Erep v_Eatt temp
thermo 100
thermo_modify lost error
timestep 1e-3
velocity all create 300.0 12345
fix intall all nve
run 1000

154
examples/PACKAGES/interlayer/ilp_tmds/log.7Jan22.mos2.g++.1 Normal file
View File

@ -0,0 +1,154 @@
LAMMPS (7 Jan 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./bilayer_MoS2_AAprime_stacking.data
Reading data file ...
triclinic box = (0 0 -100) to (51.15232 44.299209 100) with tilt (25.57616 0 0)
1 by 1 by 1 MPI processor grid
reading atoms ...
1536 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.007 seconds
mass * 32.065 # mass of sulphur atom , uint: a.u.=1.66X10^(-27)kg
mass 1 95.94 # mass of molebdenum atom , uint: a.u.=1.66X10^(-27)kg
mass 4 95.94
# Define potentials
pair_style hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
pair_coeff * * sw/mod 1 tmd.sw.mod Mo S S NULL NULL NULL
Reading sw potential file tmd.sw.mod with DATE: 2018-03-26
pair_coeff * * sw/mod 2 tmd.sw.mod NULL NULL NULL Mo S S
Reading sw potential file tmd.sw.mod with DATE: 2018-03-26
pair_coeff * * ilp/tmd MoS2.ILP Mo S S Mo S S
Reading ilp/tmd potential file MoS2.ILP with DATE: 2021-12-02
# Calculate the pair potential
compute 0 all pair ilp/tmd
compute 1 all pair sw/mod 1
compute 2 all pair sw/mod 2
variable SW1 equal c_1
variable SW2 equal c_2
variable ILP equal c_0
variable Eatt equal c_0[1]
variable Erep equal c_0[2]
thermo_style custom step etotal pe ke v_SW1 v_SW2 v_ILP v_Erep v_Eatt temp
thermo 100
thermo_modify lost error
timestep 1e-3
velocity all create 300.0 12345
fix intall all nve
run 1000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- ilp/graphene/hbn potential doi:10.1021/acs.nanolett.8b02848
@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,
}
- ilp/tmd potential doi/10.1021/acs.jctc.1c00782
@Article{Ouyang2021
author = {W. Ouyang, R. Sofer, X. Gao, J. Hermann, A. Tkatchenko, L. Kronik, M. Urbakh, and O. Hod},
title = {Anisotropic Interlayer Force Field for Transition Metal Dichalcogenides: The Case of Molybdenum Disulfide},
journal = {J. Chem. Theory Comput.},
volume = 17,
pages = {7237–7245}
year = 2021,
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 18
ghost atom cutoff = 18
binsize = 9, bins = 9 5 23
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair sw/mod, perpetual, skip from (4)
attributes: full, newton on
pair build: skip
stencil: none
bin: none
(2) pair sw/mod, perpetual, skip from (4)
attributes: full, newton on
pair build: skip
stencil: none
bin: none
(3) pair ilp/tmd, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(4) neighbor class addition, perpetual, copy from (3)
attributes: full, newton on
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 30.29 | 30.29 | 30.29 Mbytes
Step TotEng PotEng KinEng v_SW1 v_SW2 v_ILP v_Erep v_Eatt Temp
0 -1834.4469 -1893.9712 59.524297 -929.02881 -929.02881 -35.913549 63.00343 -98.916979 300
100 -1834.4497 -1883.3105 48.860775 -924.84264 -925.08505 -33.382796 56.58255 -89.965346 246.25629
200 -1834.4381 -1873.7072 39.269085 -922.29961 -922.52535 -28.882252 50.08277 -78.965022 197.91457
300 -1834.4483 -1881.1263 46.678028 -923.39264 -923.65627 -34.077402 51.011967 -85.089369 235.25534
400 -1834.431 -1868.0728 33.64182 -916.85743 -916.27044 -34.944916 50.414038 -85.358954 169.55338
500 -1834.4499 -1881.9059 47.456 -925.22919 -924.29582 -32.380877 44.755168 -77.136045 239.17628
600 -1834.4343 -1869.8642 35.429976 -920.97805 -919.60784 -29.278358 43.270241 -72.548599 178.56562
700 -1834.443 -1878.144 43.700934 -921.8051 -921.55671 -34.782141 49.959943 -84.742084 220.2509
800 -1834.4327 -1869.824 35.391298 -917.19193 -917.91383 -34.718247 55.349728 -90.067976 178.37068
900 -1834.4465 -1877.6741 43.227638 -923.33877 -922.50874 -31.826599 53.965592 -85.792191 217.86551
1000 -1834.4412 -1876.1808 41.739587 -923.17282 -923.49367 -29.514347 55.454643 -84.96899 210.3658
Loop time of 72.8218 on 1 procs for 1000 steps with 1536 atoms
Performance: 1.186 ns/day, 20.228 hours/ns, 13.732 timesteps/s
99.6% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 72.781 | 72.781 | 72.781 | 0.0 | 99.94
Bond | 0.00014503 | 0.00014503 | 0.00014503 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.016691 | 0.016691 | 0.016691 | 0.0 | 0.02
Output | 0.00057989 | 0.00057989 | 0.00057989 | 0.0 | 0.00
Modify | 0.013405 | 0.013405 | 0.013405 | 0.0 | 0.02
Other | | 0.01044 | | | 0.01
Nlocal: 1536 ave 1536 max 1536 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 3510 ave 3510 max 3510 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 992256 ave 992256 max 992256 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 992256
Ave neighs/atom = 646
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:01:12

154
examples/PACKAGES/interlayer/ilp_tmds/log.7Jan22.mos2.g++.4 Normal file
View File

@ -0,0 +1,154 @@
LAMMPS (7 Jan 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./bilayer_MoS2_AAprime_stacking.data
Reading data file ...
triclinic box = (0 0 -100) to (51.15232 44.299209 100) with tilt (25.57616 0 0)
2 by 2 by 1 MPI processor grid
reading atoms ...
1536 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.007 seconds
mass * 32.065 # mass of sulphur atom , uint: a.u.=1.66X10^(-27)kg
mass 1 95.94 # mass of molebdenum atom , uint: a.u.=1.66X10^(-27)kg
mass 4 95.94
# Define potentials
pair_style hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
pair_coeff * * sw/mod 1 tmd.sw.mod Mo S S NULL NULL NULL
Reading sw potential file tmd.sw.mod with DATE: 2018-03-26
pair_coeff * * sw/mod 2 tmd.sw.mod NULL NULL NULL Mo S S
Reading sw potential file tmd.sw.mod with DATE: 2018-03-26
pair_coeff * * ilp/tmd MoS2.ILP Mo S S Mo S S
Reading ilp/tmd potential file MoS2.ILP with DATE: 2021-12-02
# Calculate the pair potential
compute 0 all pair ilp/tmd
compute 1 all pair sw/mod 1
compute 2 all pair sw/mod 2
variable SW1 equal c_1
variable SW2 equal c_2
variable ILP equal c_0
variable Eatt equal c_0[1]
variable Erep equal c_0[2]
thermo_style custom step etotal pe ke v_SW1 v_SW2 v_ILP v_Erep v_Eatt temp
thermo 100
thermo_modify lost error
timestep 1e-3
velocity all create 300.0 12345
fix intall all nve
run 1000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- ilp/graphene/hbn potential doi:10.1021/acs.nanolett.8b02848
@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,
}
- ilp/tmd potential doi/10.1021/acs.jctc.1c00782
@Article{Ouyang2021
author = {W. Ouyang, R. Sofer, X. Gao, J. Hermann, A. Tkatchenko, L. Kronik, M. Urbakh, and O. Hod},
title = {Anisotropic Interlayer Force Field for Transition Metal Dichalcogenides: The Case of Molybdenum Disulfide},
journal = {J. Chem. Theory Comput.},
volume = 17,
pages = {7237–7245}
year = 2021,
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 18
ghost atom cutoff = 18
binsize = 9, bins = 9 5 23
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair sw/mod, perpetual, skip from (4)
attributes: full, newton on
pair build: skip
stencil: none
bin: none
(2) pair sw/mod, perpetual, skip from (4)
attributes: full, newton on
pair build: skip
stencil: none
bin: none
(3) pair ilp/tmd, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(4) neighbor class addition, perpetual, copy from (3)
attributes: full, newton on
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 17.98 | 17.98 | 17.98 Mbytes
Step TotEng PotEng KinEng v_SW1 v_SW2 v_ILP v_Erep v_Eatt Temp
0 -1834.4469 -1893.9712 59.524297 -929.02881 -929.02881 -35.913549 63.00343 -98.916979 300
100 -1834.4497 -1883.3105 48.860775 -924.84264 -925.08505 -33.382796 56.58255 -89.965346 246.25629
200 -1834.4381 -1873.7072 39.269085 -922.29961 -922.52535 -28.882252 50.08277 -78.965022 197.91457
300 -1834.4483 -1881.1263 46.678028 -923.39264 -923.65627 -34.077402 51.011967 -85.089369 235.25534
400 -1834.431 -1868.0728 33.64182 -916.85743 -916.27044 -34.944916 50.414038 -85.358954 169.55338
500 -1834.4499 -1881.9059 47.456 -925.22919 -924.29582 -32.380877 44.755168 -77.136045 239.17628
600 -1834.4343 -1869.8642 35.429976 -920.97805 -919.60784 -29.278358 43.270241 -72.548599 178.56562
700 -1834.443 -1878.144 43.700934 -921.8051 -921.55671 -34.782141 49.959943 -84.742084 220.2509
800 -1834.4327 -1869.824 35.391298 -917.19193 -917.91383 -34.718247 55.349728 -90.067976 178.37068
900 -1834.4465 -1877.6741 43.227638 -923.33877 -922.50874 -31.826599 53.965592 -85.792191 217.86551
1000 -1834.4412 -1876.1808 41.739587 -923.17282 -923.49367 -29.514347 55.454643 -84.96899 210.3658
Loop time of 24.6309 on 4 procs for 1000 steps with 1536 atoms
Performance: 3.508 ns/day, 6.842 hours/ns, 40.599 timesteps/s
99.7% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 22.906 | 23.627 | 24.36 | 14.0 | 95.92
Bond | 0.00010889 | 0.00011572 | 0.00012719 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.25886 | 0.992 | 1.7126 | 68.3 | 4.03
Output | 0.00050901 | 0.00054202 | 0.00062029 | 0.0 | 0.00
Modify | 0.0030735 | 0.0033236 | 0.0035581 | 0.3 | 0.01
Other | | 0.008194 | | | 0.03
Nlocal: 384 ave 387 max 381 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Nghost: 2262 ave 2265 max 2259 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 248064 ave 250002 max 246126 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Total # of neighbors = 992256
Ave neighs/atom = 646
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:24

143
examples/PACKAGES/interlayer/ilp_tmds/tmd.sw.mod Normal file
View File

@ -0,0 +1,143 @@
# DATE: 2018-03-26 UNITS: metal CONTRIBUTOR: Jin-Wu Jiang jwjiang5918@hotmail.com
# CITATION: J.-W. Jiang, Acta Mech. Solida. Sin 32, 17 (2019).
# The Stillinger-Weber parameters, for transition-metal dichalcogenide (TMD) lateral heterostructures.
# M = Mo, W; X = S, Se, Te
# these entries are in LAMMPS "metal" units:
# epsilon = eV; sigma = Angstroms
# other quantities are unitless
# format of a single entry (one or more lines):
# element 1, element 2, element 3,
# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol
# M-X-X terms
Mo S S 1.000 1.252 2.523 67.883 1.000 0.143 6.918 7.223 4 0 0.0
Mo Se Se 1.000 0.913 3.672 32.526 1.000 0.143 5.737 27.084 4 0 0.0
Mo Te Te 1.000 0.880 4.097 23.705 1.000 0.143 5.086 40.810 4 0 0.0
W S S 1.000 0.889 3.558 37.687 1.000 0.143 5.664 24.525 4 0 0.0
W Se Se 1.000 0.706 4.689 25.607 1.000 0.143 5.476 65.662 4 0 0.0
W Te Te 1.000 0.778 4.632 21.313 1.000 0.143 4.326 62.148 4 0 0.0
# X-M-M terms
S Mo Mo 1.000 1.252 2.523 62.449 1.000 0.143 6.918 7.223 4 0 0.0
S W W 1.000 0.889 3.558 33.553 1.000 0.143 5.664 24.525 4 0 0.0
Se Mo Mo 1.000 0.913 3.672 27.079 1.000 0.143 5.737 27.084 4 0 0.0
Se W W 1.000 0.706 4.689 23.218 1.000 0.143 5.476 65.662 4 0 0.0
Te Mo Mo 1.000 0.880 4.097 20.029 1.000 0.143 5.086 40.810 4 0 0.0
Te W W 1.000 0.778 4.632 17.370 1.000 0.143 4.326 62.148 4 0 0.0
# M-X1-X2 terms
Mo S Se 1.000 0.000 0.000 46.989 1.000 0.143 0.000 0.000 4 0 0.0
Mo S Te 1.000 0.000 0.000 40.114 1.000 0.143 0.000 0.000 4 0 0.0
Mo Se S 1.000 0.000 0.000 46.989 1.000 0.143 0.000 0.000 4 0 0.0
Mo Se Te 1.000 0.000 0.000 27.767 1.000 0.143 0.000 0.000 4 0 0.0
Mo Te S 1.000 0.000 0.000 40.114 1.000 0.143 0.000 0.000 4 0 0.0
Mo Te Se 1.000 0.000 0.000 27.767 1.000 0.143 0.000 0.000 4 0 0.0
W S Se 1.000 0.000 0.000 31.065 1.000 0.143 0.000 0.000 4 0 0.0
W S Te 1.000 0.000 0.000 28.341 1.000 0.143 0.000 0.000 4 0 0.0
W Se S 1.000 0.000 0.000 31.065 1.000 0.143 0.000 0.000 4 0 0.0
W Se Te 1.000 0.000 0.000 23.362 1.000 0.143 0.000 0.000 4 0 0.0
W Te S 1.000 0.000 0.000 28.341 1.000 0.143 0.000 0.000 4 0 0.0
W Te Se 1.000 0.000 0.000 23.362 1.000 0.143 0.000 0.000 4 0 0.0
# X-M1-M2 terms
S Mo W 1.000 0.000 0.000 45.775 1.000 0.143 0.000 0.000 4 0 0.0
S W Mo 1.000 0.000 0.000 45.775 1.000 0.143 0.000 0.000 4 0 0.0
Se Mo W 1.000 0.000 0.000 25.074 1.000 0.143 0.000 0.000 4 0 0.0
Se W Mo 1.000 0.000 0.000 25.074 1.000 0.143 0.000 0.000 4 0 0.0
Te Mo W 1.000 0.000 0.000 18.652 1.000 0.143 0.000 0.000 4 0 0.0
Te W Mo 1.000 0.000 0.000 18.652 1.000 0.143 0.000 0.000 4 0 0.0
# zero terms
Mo Mo Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Mo W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Mo S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Mo Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Mo Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo W Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo W W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo W S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo W Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo W Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo S Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo S W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Se Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Se W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Te Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Te W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Mo Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Mo W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Mo S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Mo Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Mo Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W W Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W W W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W W S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W W Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W W Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W S Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W S W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Se Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Se W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Te Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Te W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Mo S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Mo Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Mo Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S W S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S W Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S W Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S S Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S S W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S S S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S S Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S S Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Se Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Se W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Se S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Se Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Se Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Te Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Te W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Te S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Te Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Te Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Mo S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Mo Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Mo Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se W S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se W Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se W Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se S Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se S W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se S S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se S Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se S Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Se Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Se W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Se S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Se Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Se Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Te Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Te W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Te S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Te Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Te Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Mo S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Mo Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Mo Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te W S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te W Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te W Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te S Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te S W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te S S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te S Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te S Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Se Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Se W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Se S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Se Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Se Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Te Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Te W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Te S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Te Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Te Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0

218
examples/PACKAGES/interlayer/saip_metal/3Lgold_1Lgr_atop_sliding.data Normal file
View File

@ -0,0 +1,218 @@
Gold/graphene heterjunction atop2
206 atoms
2 atom types
0 17.216640 xlo xhi
0 14.910048 ylo yhi
-30 30.000000 zlo zhi
8.60832 0 0 xy xz yz
Atoms
1 1 1 0 -8.60832 -4.97002 3.5
2 1 1 0 -7.1736 -4.14166 5.85461
3 1 1 0 -5.73888 -3.31334 8.20922
4 1 1 0 -7.1736 -2.48501 3.5
5 1 1 0 -5.73888 -1.65666 5.85461
6 1 1 0 -4.30416 -0.828336 8.20922
7 1 1 0 -5.73888 0 3.5
8 1 1 0 -4.30416 0.828352 5.85461
9 1 1 0 -2.86944 1.65667 8.20922
10 1 1 0 -4.30416 2.48501 3.5
11 1 1 0 -2.86944 3.31336 5.85461
12 1 1 0 -1.43472 4.14168 8.20922
13 1 1 0 -2.86944 4.97002 3.5
14 1 1 0 -1.43472 5.79837 5.85461
15 1 1 0 0 6.62669 8.20922
16 1 1 0 -1.43472 7.45502 3.5
17 1 1 0 0 8.28337 5.85461
18 1 1 0 1.43472 9.11167 8.20922
19 1 1 0 -5.73888 -4.97002 3.5
20 1 1 0 -4.30416 -4.14166 5.85461
21 1 1 0 -2.86944 -3.31334 8.20922
22 1 1 0 -4.30416 -2.48501 3.5
23 1 1 0 -2.86944 -1.65666 5.85461
24 1 1 0 -1.43472 -0.828336 8.20922
25 1 1 0 -2.86944 0 3.5
26 1 1 0 -1.43472 0.828352 5.85461
27 1 1 0 0 1.65667 8.20922
28 1 1 0 -1.43472 2.48501 3.5
29 1 1 0 0 3.31336 5.85461
30 1 1 0 1.43472 4.14168 8.20922
31 1 1 0 0 4.97002 3.5
32 1 1 0 1.43472 5.79837 5.85461
33 1 1 0 2.86944 6.62669 8.20922
34 1 1 0 1.43472 7.45502 3.5
35 1 1 0 2.86944 8.28337 5.85461
36 1 1 0 4.30416 9.11167 8.20922
37 1 1 0 -2.86944 -4.97002 3.5
38 1 1 0 -1.43472 -4.14166 5.85461
39 1 1 0 0 -3.31334 8.20922
40 1 1 0 -1.43472 -2.48501 3.5
41 1 1 0 0 -1.65666 5.85461
42 1 1 0 1.43472 -0.828336 8.20922
43 1 1 0 0 0 3.5
44 1 1 0 1.43472 0.828352 5.85461
45 1 1 0 2.86944 1.65667 8.20922
46 1 1 0 1.43472 2.48501 3.5
47 1 1 0 2.86944 3.31336 5.85461
48 1 1 0 4.30416 4.14168 8.20922
49 1 1 0 2.86944 4.97002 3.5
50 1 1 0 4.30416 5.79837 5.85461
51 1 1 0 5.73888 6.62669 8.20922
52 1 1 0 4.30416 7.45502 3.5
53 1 1 0 5.73888 8.28337 5.85461
54 1 1 0 7.1736 9.11167 8.20922
55 1 1 0 0 -4.97002 3.5
56 1 1 0 1.43472 -4.14166 5.85461
57 1 1 0 2.86944 -3.31334 8.20922
58 1 1 0 1.43472 -2.48501 3.5
59 1 1 0 2.86944 -1.65666 5.85461
60 1 1 0 4.30416 -0.828336 8.20922
61 1 1 0 2.86944 0 3.5
62 1 1 0 4.30416 0.828352 5.85461
63 1 1 0 5.73888 1.65667 8.20922
64 1 1 0 4.30416 2.48501 3.5
65 1 1 0 5.73888 3.31336 5.85461
66 1 1 0 7.1736 4.14168 8.20922
67 1 1 0 5.73888 4.97002 3.5
68 1 1 0 7.1736 5.79837 5.85461
69 1 1 0 8.60832 6.62669 8.20922
70 1 1 0 7.1736 7.45502 3.5
71 1 1 0 8.60832 8.28337 5.85461
72 1 1 0 10.0431 9.11167 8.20922
73 1 1 0 2.86944 -4.97002 3.5
74 1 1 0 4.30416 -4.14166 5.85461
75 1 1 0 5.73888 -3.31334 8.20922
76 1 1 0 4.30416 -2.48501 3.5
77 1 1 0 5.73888 -1.65666 5.85461
78 1 1 0 7.1736 -0.828336 8.20922
79 1 1 0 5.73888 0 3.5
80 1 1 0 7.1736 0.828352 5.85461
81 1 1 0 8.60832 1.65667 8.20922
82 1 1 0 7.1736 2.48501 3.5
83 1 1 0 8.60832 3.31336 5.85461
84 1 1 0 10.0431 4.14168 8.20922
85 1 1 0 8.60832 4.97002 3.5
86 1 1 0 10.0431 5.79837 5.85461
87 1 1 0 11.4777 6.62669 8.20922
88 1 1 0 10.0431 7.45502 3.5
89 1 1 0 11.4777 8.28337 5.85461
90 1 1 0 12.9124 9.11167 8.20922
91 1 1 0 5.73888 -4.97002 3.5
92 1 1 0 7.1736 -4.14166 5.85461
93 1 1 0 8.60832 -3.31334 8.20922
94 1 1 0 7.1736 -2.48501 3.5
95 1 1 0 8.60832 -1.65666 5.85461
96 1 1 0 10.0431 -0.828336 8.20922
97 1 1 0 8.60832 0 3.5
98 1 1 0 10.0431 0.828352 5.85461
99 1 1 0 11.4777 1.65667 8.20922
100 1 1 0 10.0431 2.48501 3.5
101 1 1 0 11.4777 3.31336 5.85461
102 1 1 0 12.9124 4.14168 8.20922
103 1 1 0 11.4777 4.97002 3.5
104 1 1 0 12.9124 5.79837 5.85461
105 1 1 0 14.3472 6.62669 8.20922
106 1 1 0 12.9124 7.45502 3.5
107 1 1 0 14.3472 8.28337 5.85461
108 1 1 0 15.7819 9.11167 8.20922
109 2 2 0 -11.0678 -6.39002 0
110 2 2 0 -9.83808 -5.68001 0
111 2 2 0 -9.83808 -4.26001 0
112 2 2 0 -8.60832 -3.55001 0
113 2 2 0 -8.60832 -2.13001 0
114 2 2 0 -7.37856 -1.42 0
115 2 2 0 -7.37856 0 0
116 2 2 0 -6.1488 0.710007 0
117 2 2 0 -6.1488 2.13001 0
118 2 2 0 -4.91904 2.84001 0
119 2 2 0 -4.91904 4.26001 0
120 2 2 0 -3.68928 4.97002 0
121 2 2 0 -3.68928 6.39002 0
122 2 2 0 -2.45952 7.10003 0
123 2 2 0 -8.60832 -6.39002 0
124 2 2 0 -7.37856 -5.68001 0
125 2 2 0 -7.37856 -4.26001 0
126 2 2 0 -6.1488 -3.55001 0
127 2 2 0 -6.1488 -2.13001 0
128 2 2 0 -4.91904 -1.42 0
129 2 2 0 -4.91904 0 0
130 2 2 0 -3.68928 0.710007 0
131 2 2 0 -3.68928 2.13001 0
132 2 2 0 -2.45952 2.84001 0
133 2 2 0 -2.45952 4.26001 0
134 2 2 0 -1.22976 4.97002 0
135 2 2 0 -1.22976 6.39002 0
136 2 2 0 0 7.10003 0
137 2 2 0 -6.1488 -6.39002 0
138 2 2 0 -4.91904 -5.68001 0
139 2 2 0 -4.91904 -4.26001 0
140 2 2 0 -3.68928 -3.55001 0
141 2 2 0 -3.68928 -2.13001 0
142 2 2 0 -2.45952 -1.42 0
143 2 2 0 -2.45952 0 0
144 2 2 0 -1.22976 0.710007 0
145 2 2 0 -1.22976 2.13001 0
146 2 2 0 0 2.84001 0
147 2 2 0 0 4.26001 0
148 2 2 0 1.22976 4.97002 0
149 2 2 0 1.22976 6.39002 0
150 2 2 0 2.45952 7.10003 0
151 2 2 0 -3.68928 -6.39002 0
152 2 2 0 -2.45952 -5.68001 0
153 2 2 0 -2.45952 -4.26001 0
154 2 2 0 -1.22976 -3.55001 0
155 2 2 0 -1.22976 -2.13001 0
156 2 2 0 0 -1.42 0
157 2 2 0 0 0 0
158 2 2 0 1.22976 0.710007 0
159 2 2 0 1.22976 2.13001 0
160 2 2 0 2.45952 2.84001 0
161 2 2 0 2.45952 4.26001 0
162 2 2 0 3.68928 4.97002 0
163 2 2 0 3.68928 6.39002 0
164 2 2 0 4.91904 7.10003 0
165 2 2 0 -1.22976 -6.39002 0
166 2 2 0 0 -5.68001 0
167 2 2 0 0 -4.26001 0
168 2 2 0 1.22976 -3.55001 0
169 2 2 0 1.22976 -2.13001 0
170 2 2 0 2.45952 -1.42 0
171 2 2 0 2.45952 0 0
172 2 2 0 3.68928 0.710007 0
173 2 2 0 3.68928 2.13001 0
174 2 2 0 4.91904 2.84001 0
175 2 2 0 4.91904 4.26001 0
176 2 2 0 6.1488 4.97002 0
177 2 2 0 6.1488 6.39002 0
178 2 2 0 7.37856 7.10003 0
179 2 2 0 1.22976 -6.39002 0
180 2 2 0 2.45952 -5.68001 0
181 2 2 0 2.45952 -4.26001 0
182 2 2 0 3.68928 -3.55001 0
183 2 2 0 3.68928 -2.13001 0
184 2 2 0 4.91904 -1.42 0
185 2 2 0 4.91904 0 0
186 2 2 0 6.1488 0.710007 0
187 2 2 0 6.1488 2.13001 0
188 2 2 0 7.37856 2.84001 0
189 2 2 0 7.37856 4.26001 0
190 2 2 0 8.60832 4.97002 0
191 2 2 0 8.60832 6.39002 0
192 2 2 0 9.83808 7.10003 0
193 2 2 0 3.68928 -6.39002 0
194 2 2 0 4.91904 -5.68001 0
195 2 2 0 4.91904 -4.26001 0
196 2 2 0 6.1488 -3.55001 0
197 2 2 0 6.1488 -2.13001 0
198 2 2 0 7.37856 -1.42 0
199 2 2 0 7.37856 0 0
200 2 2 0 8.60832 0.710007 0
201 2 2 0 8.60832 2.13001 0
202 2 2 0 9.83808 2.84001 0
203 2 2 0 9.83808 4.26001 0
204 2 2 0 11.0678 4.97002 0
205 2 2 0 11.0678 6.39002 0
206 2 2 0 12.2976 7.10003 0

1
examples/PACKAGES/interlayer/saip_metal/Au_u3.eam Symbolic link
View File

@ -0,0 +1 @@
../../../../potentials/Au_u3.eam

1
examples/PACKAGES/interlayer/saip_metal/CH.rebo Symbolic link
View File

@ -0,0 +1 @@
../../../../potentials/CH.rebo

1
examples/PACKAGES/interlayer/saip_metal/CHAu.ILP Symbolic link
View File

@ -0,0 +1 @@
../../../../potentials/CHAu.ILP

44
examples/PACKAGES/interlayer/saip_metal/in.gold_gr Normal file
View File

@ -0,0 +1,44 @@
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./3Lgold_1Lgr_atop_sliding.data
# global group definition
group gold type 1
group gra type 2
# Define mass
mass * 12.0107 # mass of carbon atom , uint: a.u.=1.66X10^(-27)kg
mass 1 196.96657 # mass of gold atom , uint: a.u.=1.66X10^(-27)kg
# Define potentials
pair_style hybrid/overlay eam rebo saip/metal 16.0
pair_coeff 1 1 eam ./Au_u3.eam
pair_coeff * * rebo ./CH.rebo NULL C
pair_coeff * * saip/metal ./CHAu.ILP Au C
# compute energy
compute 0 all pair rebo
compute 1 all pair eam
compute 2 all pair saip/metal
variable REBO equal c_0
variable EAM equal c_1
variable ILP equal c_2
thermo_style custom step etotal pe ke v_REBO v_ILP temp
thermo 100
thermo_modify lost error
# Creat initial velocity
velocity all set 0.0 0.0 0.0
velocity gra create 300.0 4928459 mom yes rot yes dist gaussian
velocity gold create 300.0 4928459 mom yes rot yes dist gaussian
# Integration
fix intsub gold nve
fix intrib gra nve
timestep 1e-3
run 1000

164
examples/PACKAGES/interlayer/saip_metal/log.7Jan22.gold_gr.g++.1 Normal file
View File

@ -0,0 +1,164 @@
LAMMPS (7 Jan 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./3Lgold_1Lgr_atop_sliding.data
Reading data file ...
triclinic box = (0 0 -30) to (17.21664 14.910048 30) with tilt (8.60832 0 0)
1 by 1 by 1 MPI processor grid
reading atoms ...
206 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.003 seconds
# global group definition
group gold type 1
108 atoms in group gold
group gra type 2
98 atoms in group gra
# Define mass
mass * 12.0107 # mass of carbon atom , uint: a.u.=1.66X10^(-27)kg
mass 1 196.96657 # mass of gold atom , uint: a.u.=1.66X10^(-27)kg
# Define potentials
pair_style hybrid/overlay eam rebo saip/metal 16.0
pair_coeff 1 1 eam ./Au_u3.eam
Reading eam potential file ./Au_u3.eam with DATE: 2007-06-11
pair_coeff * * rebo ./CH.rebo NULL C
Reading rebo potential file ./CH.rebo with DATE: 2018-7-3
pair_coeff * * saip/metal ./CHAu.ILP Au C
Reading saip/metal potential file ./CHAu.ILP with DATE: 2021-12-02
# compute energy
compute 0 all pair rebo
compute 1 all pair eam
compute 2 all pair saip/metal
variable REBO equal c_0
variable EAM equal c_1
variable ILP equal c_2
thermo_style custom step etotal pe ke v_REBO v_ILP temp
thermo 100
thermo_modify lost error
# Creat initial velocity
velocity all set 0.0 0.0 0.0
velocity gra create 300.0 4928459 mom yes rot yes dist gaussian
velocity gold create 300.0 4928459 mom yes rot yes dist gaussian
# Integration
fix intsub gold nve
fix intrib gra nve
timestep 1e-3
run 1000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- ilp/graphene/hbn potential doi:10.1021/acs.nanolett.8b02848
@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,
}
- saip/metal potential doi.org/10.1021/acs.jctc.1c00622
@Article{Ouyang2021
author = {W. Ouyang, O. Hod, and R. Guerra},
title = {Registry-Dependent Potential for Interfaces of Gold with Graphitic Systems},
journal = {J. Chem. Theory Comput.},
volume = 17,
pages = {7215-7223}
year = 2021,
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 18
ghost atom cutoff = 18
binsize = 9, bins = 3 2 7
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair eam, perpetual, skip from (4)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
(2) pair rebo, perpetual, skip from (3)
attributes: full, newton on, ghost
pair build: skip/ghost
stencil: none
bin: none
(3) pair saip/metal, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(4) neighbor class addition, perpetual
attributes: half, newton on
pair build: half/bin/newton/tri
stencil: half/bin/3d/tri
bin: standard
Per MPI rank memory allocation (min/avg/max) = 9.747 | 9.747 | 9.747 Mbytes
Step TotEng PotEng KinEng v_REBO v_ILP Temp
0 -1121.4621 -1129.3728 7.9107209 -724.70925 -6.0302289 298.53659
100 -1121.4483 -1124.9731 3.5248501 -723.03272 -5.9765533 133.02159
200 -1121.4403 -1125.2912 3.8509646 -722.66784 -6.0468507 145.32858
300 -1121.4424 -1126.4665 5.0240531 -722.72787 -6.0350568 189.59886
400 -1121.4419 -1125.1443 3.7023978 -722.59976 -5.8294548 139.72193
500 -1121.4413 -1125.2711 3.8297939 -722.5342 -6.0189944 144.52963
600 -1121.4449 -1125.8808 4.4359049 -722.77707 -5.8685221 167.40319
700 -1121.4489 -1126.1966 4.747709 -723.13681 -5.8666379 179.17012
800 -1121.4443 -1125.8469 4.402607 -722.94487 -6.0094567 166.14658
900 -1121.4424 -1125.8437 4.4013317 -722.94918 -5.8699702 166.09846
1000 -1121.4467 -1125.7453 4.2986881 -722.66682 -6.0651168 162.22487
Loop time of 6.43246 on 1 procs for 1000 steps with 206 atoms
Performance: 13.432 ns/day, 1.787 hours/ns, 155.462 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 6.4201 | 6.4201 | 6.4201 | 0.0 | 99.81
Bond | 8.9059e-05 | 8.9059e-05 | 8.9059e-05 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0071852 | 0.0071852 | 0.0071852 | 0.0 | 0.11
Output | 0.00026031 | 0.00026031 | 0.00026031 | 0.0 | 0.00
Modify | 0.0019433 | 0.0019433 | 0.0019433 | 0.0 | 0.03
Other | | 0.002875 | | | 0.04
Nlocal: 206 ave 206 max 206 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 2187 ave 2187 max 2187 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 158548 ave 158548 max 158548 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 158548
Ave neighs/atom = 769.65049
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:06

164
examples/PACKAGES/interlayer/saip_metal/log.7Jan22.gold_gr.g++.4 Normal file
View File

@ -0,0 +1,164 @@
LAMMPS (7 Jan 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./3Lgold_1Lgr_atop_sliding.data
Reading data file ...
triclinic box = (0 0 -30) to (17.21664 14.910048 30) with tilt (8.60832 0 0)
2 by 2 by 1 MPI processor grid
reading atoms ...
206 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.003 seconds
# global group definition
group gold type 1
108 atoms in group gold
group gra type 2
98 atoms in group gra
# Define mass
mass * 12.0107 # mass of carbon atom , uint: a.u.=1.66X10^(-27)kg
mass 1 196.96657 # mass of gold atom , uint: a.u.=1.66X10^(-27)kg
# Define potentials
pair_style hybrid/overlay eam rebo saip/metal 16.0
pair_coeff 1 1 eam ./Au_u3.eam
Reading eam potential file ./Au_u3.eam with DATE: 2007-06-11
pair_coeff * * rebo ./CH.rebo NULL C
Reading rebo potential file ./CH.rebo with DATE: 2018-7-3
pair_coeff * * saip/metal ./CHAu.ILP Au C
Reading saip/metal potential file ./CHAu.ILP with DATE: 2021-12-02
# compute energy
compute 0 all pair rebo
compute 1 all pair eam
compute 2 all pair saip/metal
variable REBO equal c_0
variable EAM equal c_1
variable ILP equal c_2
thermo_style custom step etotal pe ke v_REBO v_ILP temp
thermo 100
thermo_modify lost error
# Creat initial velocity
velocity all set 0.0 0.0 0.0
velocity gra create 300.0 4928459 mom yes rot yes dist gaussian
velocity gold create 300.0 4928459 mom yes rot yes dist gaussian
# Integration
fix intsub gold nve
fix intrib gra nve
timestep 1e-3
run 1000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- ilp/graphene/hbn potential doi:10.1021/acs.nanolett.8b02848
@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,
}
- saip/metal potential doi.org/10.1021/acs.jctc.1c00622
@Article{Ouyang2021
author = {W. Ouyang, O. Hod, and R. Guerra},
title = {Registry-Dependent Potential for Interfaces of Gold with Graphitic Systems},
journal = {J. Chem. Theory Comput.},
volume = 17,
pages = {7215-7223}
year = 2021,
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 18
ghost atom cutoff = 18
binsize = 9, bins = 3 2 7
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair eam, perpetual, skip from (4)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
(2) pair rebo, perpetual, skip from (3)
attributes: full, newton on, ghost
pair build: skip/ghost
stencil: none
bin: none
(3) pair saip/metal, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(4) neighbor class addition, perpetual
attributes: half, newton on
pair build: half/bin/newton/tri
stencil: half/bin/3d/tri
bin: standard
Per MPI rank memory allocation (min/avg/max) = 9.299 | 9.299 | 9.3 Mbytes
Step TotEng PotEng KinEng v_REBO v_ILP Temp
0 -1121.4621 -1129.3728 7.9107209 -724.70925 -6.0302289 298.53659
100 -1121.4483 -1124.9731 3.5248501 -723.03272 -5.9765533 133.02159
200 -1121.4403 -1125.2912 3.8509646 -722.66784 -6.0468507 145.32858
300 -1121.4424 -1126.4665 5.0240531 -722.72787 -6.0350568 189.59886
400 -1121.4419 -1125.1443 3.7023978 -722.59976 -5.8294548 139.72193
500 -1121.4413 -1125.2711 3.8297939 -722.5342 -6.0189944 144.52963
600 -1121.4449 -1125.8808 4.4359049 -722.77707 -5.8685221 167.40319
700 -1121.4489 -1126.1966 4.747709 -723.13681 -5.8666379 179.17012
800 -1121.4443 -1125.8469 4.402607 -722.94487 -6.0094567 166.14658
900 -1121.4424 -1125.8437 4.4013317 -722.94918 -5.8699702 166.09846
1000 -1121.4467 -1125.7453 4.2986881 -722.66682 -6.0651168 162.22487
Loop time of 2.44776 on 4 procs for 1000 steps with 206 atoms
Performance: 35.298 ns/day, 0.680 hours/ns, 408.537 timesteps/s
99.7% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.6906 | 2.0172 | 2.3939 | 19.2 | 82.41
Bond | 5.4278e-05 | 6.3818e-05 | 7.3153e-05 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.049363 | 0.42514 | 0.75161 | 41.7 | 17.37
Output | 0.00018596 | 0.00020656 | 0.00024754 | 0.0 | 0.01
Modify | 0.00063656 | 0.00074701 | 0.00089514 | 0.0 | 0.03
Other | | 0.004362 | | | 0.18
Nlocal: 51.5 ave 61 max 44 min
Histogram: 1 1 0 0 0 0 1 0 0 1
Nghost: 1670.75 ave 1699 max 1641 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 39637 ave 47080 max 33737 min
Histogram: 1 1 0 0 0 0 1 0 0 1
Total # of neighbors = 158548
Ave neighs/atom = 769.65049
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:02

18
potentials/CHAu.ILP Normal file
View File

@ -0,0 +1,18 @@
# DATE: 2021-12-02 UNITS: metal CONTRIBUTOR: Wengen Ouyang w.g.ouyang@gmail.com
# CITATION: W. Ouyang, O. Hod, and R. Guerra, J. Chem. Theory Comput. 17, 7215 (2021).
# Semi-anisotropic interlayer Potential (SAIP) for graphene/gold, benzene/gold heterojunctions
# The parameters below are fitted against the PBE + D3 DFT reference data.
# The parameters for C-C, C-H and H-H are taken from Nano Letters 18, 6009-6016 (2018).
#
# --------------- Repulsion Potential ----------------++++++++++++++ Vdw Potential ++++++++++++++++*****
# beta(A) alpha delta(A) epsilon(meV) C(meV) d sR reff(A) C6(meV*A^6) S rcut
C C 3.205843 7.511126 1.235334 1.528338E-5 37.530428 15.499947 0.7954443 3.681440 25.714535E3 1.0 2.0
H H 3.974540 6.53799 1.080633 0.6700556 0.8333833 15.022371 0.7490632 2.767223 1.6159581E3 1.0 1.2
C H 2.642950 12.91410 1.020257 0.9750012 25.340996 15.222927 0.8115998 3.887324 5.6874617E3 1.0 1.5
Au C 3.6913278482 13.5655648421 1.0175514400 0.0070964784 -0.0010368264 11.0586486772 1.0635582839 3.7552608806 81.5847131142 1000.0 1.0
Au H 3.7899616131 4.3065009639 10.6811825156 0.2250887843 -0.1116891520 18.6149213872 0.9833194192 3.3507558896 70.6865381780 1000.0 1.0
Au Au 3.6671967387 12.8109735143 1.0353581041 0.0000000000 0.0000000000 10.1628585345 1.0642897301 3.7372959779 0.0000000000 1000.0 1.0
# Symmetric part
C Au 3.6913278482 13.5655648421 1.0175514400 0.0070964784 -0.0010368264 11.0586486772 1.0635582839 3.7552608806 81.5847131142 1000.0 2.0
H Au 3.7899616131 4.3065009639 10.6811825156 0.2250887843 -0.1116891520 18.6149213872 0.9833194192 3.3507558896 70.6865381780 1000.0 1.5
H C 2.642950 12.91410 1.020257 0.9750012 25.340996 15.222927 0.8115998 3.887324 5.6874617E3 1.0 1.5

12
potentials/MoS2.ILP Normal file
View File

@ -0,0 +1,12 @@
# DATE: 2021-12-02 UNITS: metal CONTRIBUTOR: Wengen Ouyang w.g.ouyang@gmail.com
# CITATION: W. Ouyang, et al., J. Chem. Theory Comput. 17, 7237 (2021).
# Interlayer Potential (ILP) for bilayer and bulk Molybdenum Disulfide.
# The parameters below are fitted against the HSE + MBD-NL DFT reference data.
#
# --------------- Repulsion Potential ----------------++++++++++++++ Vdw Potential ++++++++++++++++*****
# beta(A) alpha delta(A) epsilon(meV) C(meV) d sR reff(A) C6(meV*A^6) S rcut
Mo Mo 5.5794503728 9.3776624643 2.0272224617 144.1517754265 97.9785701736 89.4375967588 2.0590305447 5.1220549022 491850.3161946840 1.0000 4.0
S S 3.1614021873 8.0932627454 1.9531402037 4.5867641760 118.0654664611 58.8094158385 0.2153665787 4.2996001250 148811.2434089213 1.0000 4.0
Mo S 3.6271521213 19.9713750996 7.5850312329 76.1019310047 3.3174955929 45.7203284638 0.9474703731 4.4104250318 150597.8577162951 1.0000 4.0
# symmetric part
S Mo 3.6271521213 19.9713750996 7.5850312329 76.1019310047 3.3174955929 45.7203284638 0.9474703731 4.4104250318 150597.8577162951 1.0000 4.0

38
src/INTERLAYER/pair_ilp_graphene_hbn.cpp
View File

@ -37,6 +37,7 @@
#include <cmath>
#include <cstring>
#include <map>
using namespace LAMMPS_NS;
using namespace InterLayer;
@ -56,9 +57,15 @@ static const char cite_ilp[] =
" year = 2018,\n"
"}\n\n";
// to indicate which potential style was used in outputs
static std::map<int, std::string> variant_map = {
{PairILPGrapheneHBN::ILP_GrhBN, "ilp/graphene/hbn"},
{PairILPGrapheneHBN::ILP_TMD, "ilp/tmd"},
{PairILPGrapheneHBN::SAIP_METAL, "saip/metal"}};
/* ---------------------------------------------------------------------- */
PairILPGrapheneHBN::PairILPGrapheneHBN(LAMMPS *lmp) : Pair(lmp)
PairILPGrapheneHBN::PairILPGrapheneHBN(LAMMPS *lmp) : Pair(lmp), variant(ILP_GrhBN)
{
restartinfo = 0;
one_coeff = 1;
@ -86,6 +93,14 @@ PairILPGrapheneHBN::PairILPGrapheneHBN(LAMMPS *lmp) : Pair(lmp)
dnormal = nullptr;
dnormdri = nullptr;
// for ilp/tmd
dnn = nullptr;
vect = nullptr;
pvet = nullptr;
dpvet1 = nullptr;
dpvet2 = nullptr;
dNave = nullptr;
// always compute energy offset
offset_flag = 1;
@ -104,6 +119,13 @@ PairILPGrapheneHBN::~PairILPGrapheneHBN()
memory->destroy(normal);
memory->destroy(dnormal);
memory->destroy(dnormdri);
// adds for ilp/tmd
memory->destroy(dnn);
memory->destroy(vect);
memory->destroy(pvet);
memory->destroy(dpvet1);
memory->destroy(dpvet2);
memory->destroy(dNave);
if (allocated) {
memory->destroy(setflag);
@ -194,7 +216,7 @@ void PairILPGrapheneHBN::read_file(char *filename)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, filename, "ilp/graphene/hbn", unit_convert_flag);
PotentialFileReader reader(lmp, filename, variant_map[variant], unit_convert_flag);
char *line;
// transparently convert units for supported conversions
@ -292,11 +314,15 @@ void PairILPGrapheneHBN::read_file(char *filename)
int n = -1;
for (int m = 0; m < nparams; m++) {
if (i == params[m].ielement && j == params[m].jelement) {
if (n >= 0) error->all(FLERR, "ILP potential file has duplicate entry");
if (n >= 0)
error->all(FLERR, "{} potential file {} has a duplicate entry", variant_map[variant],
filename);
n = m;
}
}
if (n < 0) error->all(FLERR, "Potential file is missing an entry");
if (n < 0)
error->all(FLERR, "{} potential file {} is missing an entry", variant_map[variant],
filename);
elem2param[i][j] = n;
cutILPsq[i][j] = params[n].rcut * params[n].rcut;
}
@ -310,9 +336,9 @@ void PairILPGrapheneHBN::read_file(char *filename)
void PairILPGrapheneHBN::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR, "Pair style ilp/graphene/hbn requires newton pair on");
error->all(FLERR, "Pair style {} requires newton pair on", variant_map[variant]);
if (!atom->molecule_flag)
error->all(FLERR, "Pair style ilp/graphene/hbn requires atom attribute molecule");
error->all(FLERR, "Pair style {} requires atom attribute molecule", variant_map[variant]);
// need a full neighbor list, including neighbors of ghosts

18
src/INTERLAYER/pair_ilp_graphene_hbn.h
View File

@ -34,16 +34,16 @@ class PairILPGrapheneHBN : public Pair {
void coeff(int, char **) override;
double init_one(int, int) override;
void init_style() override;
void ILP_neigh();
void calc_normal();
void calc_FRep(int, int);
void calc_FvdW(int, int);
double single(int, int, int, int, double, double, double, double &) override;
static constexpr int NPARAMS_PER_LINE = 13;
enum { ILP_GrhBN, ILP_TMD, SAIP_METAL }; // for telling class variants apart in shared code
protected:
int me;
int variant;
int maxlocal; // size of numneigh, firstneigh arrays
int pgsize; // size of neighbor page
int oneatom; // max # of neighbors for one atom
@ -68,6 +68,18 @@ class PairILPGrapheneHBN : public Pair {
double ***dnormdri;
double ****dnormal;
// adds for ilp/tmd
int Nnei; // max # of nearest neighbors for one atom
double **dnn;
double **vect;
double **pvet;
double ***dpvet1;
double ***dpvet2;
double ***dNave;
virtual void ILP_neigh();
virtual void calc_normal();
virtual void calc_FRep(int, int);
void read_file(char *);
void allocate();
};

885
src/INTERLAYER/pair_ilp_tmd.cpp Normal file
View File

@ -0,0 +1,885 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Wengen Ouyang (Wuhan University)
e-mail: w.g.ouyang at gmail dot com
This is a full version of the potential described in
[Ouyang et al, J. Chem. Theory Comput. 17, 7237 (2021).]
------------------------------------------------------------------------- */
#include "pair_ilp_tmd.h"
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "interlayer_taper.h"
#include "memory.h"
#include "my_page.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace InterLayer;
#define MAXLINE 1024
#define DELTA 4
#define PGDELTA 1
static const char cite_ilp_tmd[] = "ilp/tmd potential doi/10.1021/acs.jctc.1c00782\n"
"@Article{Ouyang2021\n"
" author = {W. Ouyang, R. Sofer, X. Gao, J. Hermann, A. "
"Tkatchenko, L. Kronik, M. Urbakh, and O. Hod},\n"
" title = {Anisotropic Interlayer Force Field for Transition "
"Metal Dichalcogenides: The Case of Molybdenum Disulfide},\n"
" journal = {J. Chem. Theory Comput.},\n"
" volume = 17,\n"
" pages = {7237–7245}\n"
" year = 2021,\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
PairILPTMD::PairILPTMD(LAMMPS *lmp) : PairILPGrapheneHBN(lmp)
{
variant = ILP_TMD;
single_enable = 0;
// for TMD, each atom have six neighbors
Nnei = 6;
if (lmp->citeme) lmp->citeme->add(cite_ilp_tmd);
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairILPTMD::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command");
if (!utils::strmatch(force->pair_style, "^hybrid/overlay"))
error->all(FLERR, "Pair style ilp/tmd must be used as sub-style with hybrid/overlay");
cut_global = utils::numeric(FLERR, arg[0], false, lmp);
if (narg == 2) tap_flag = utils::numeric(FLERR, arg[1], false, lmp);
}
/* ----------------------------------------------------------------------
Repulsive forces and energy
------------------------------------------------------------------------- */
void PairILPTMD::calc_FRep(int eflag, int /* vflag */)
{
int i, j, ii, jj, inum, jnum, itype, jtype, k, kk;
double prodnorm1, fkcx, fkcy, fkcz;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, fpair1;
double rsq, r, Rcut, rhosq1, exp0, exp1, r2inv, r6inv, r8inv, Tap, dTap, Vilp;
double frho1, TSvdw, TSvdw2inv, Erep, fsum, rdsq1;
int *ilist, *jlist, *numneigh, **firstneigh;
int *ILP_neighs_i;
evdwl = 0.0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double dprodnorm1[3] = {0.0, 0.0, 0.0};
double fp1[3] = {0.0, 0.0, 0.0};
double fprod1[3] = {0.0, 0.0, 0.0};
double delki[3] = {0.0, 0.0, 0.0};
double fk[3] = {0.0, 0.0, 0.0};
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
//calculate exp(-lambda*(r-z0))*[epsilon/2 + f(rho_ij)]
// loop over neighbors of owned atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
// only include the interation between different layers
if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) {
int iparam_ij = elem2param[map[itype]][map[jtype]];
Param &p = params[iparam_ij];
r = sqrt(rsq);
// turn on/off taper function
if (tap_flag) {
Rcut = sqrt(cutsq[itype][jtype]);
Tap = calc_Tap(r, Rcut);
dTap = calc_dTap(r, Rcut);
} else {
Tap = 1.0;
dTap = 0.0;
}
// Calculate the transverse distance
// note that rho_ij does not equal to rho_ji except when normals are all along z
prodnorm1 = normal[i][0] * delx + normal[i][1] * dely + normal[i][2] * delz;
rhosq1 = rsq - prodnorm1 * prodnorm1; // rho_ij
rdsq1 = rhosq1 * p.delta2inv; // (rho_ij/delta)^2
// store exponents
exp0 = exp(-p.lambda * (r - p.z0));
exp1 = exp(-rdsq1);
frho1 = exp1 * p.C;
Erep = 0.5 * p.epsilon + frho1;
Vilp = exp0 * Erep;
// derivatives
fpair = p.lambda * exp0 / r * Erep;
fpair1 = 2.0 * exp0 * frho1 * p.delta2inv;
fsum = fpair + fpair1;
// derivatives of the product of rij and ni, the result is a vector
dprodnorm1[0] =
dnormdri[i][0][0] * delx + dnormdri[i][1][0] * dely + dnormdri[i][2][0] * delz;
dprodnorm1[1] =
dnormdri[i][0][1] * delx + dnormdri[i][1][1] * dely + dnormdri[i][2][1] * delz;
dprodnorm1[2] =
dnormdri[i][0][2] * delx + dnormdri[i][1][2] * dely + dnormdri[i][2][2] * delz;
fp1[0] = prodnorm1 * normal[i][0] * fpair1;
fp1[1] = prodnorm1 * normal[i][1] * fpair1;
fp1[2] = prodnorm1 * normal[i][2] * fpair1;
fprod1[0] = prodnorm1 * dprodnorm1[0] * fpair1;
fprod1[1] = prodnorm1 * dprodnorm1[1] * fpair1;
fprod1[2] = prodnorm1 * dprodnorm1[2] * fpair1;
fkcx = (delx * fsum - fp1[0]) * Tap - Vilp * dTap * delx / r;
fkcy = (dely * fsum - fp1[1]) * Tap - Vilp * dTap * dely / r;
fkcz = (delz * fsum - fp1[2]) * Tap - Vilp * dTap * delz / r;
f[i][0] += fkcx - fprod1[0] * Tap;
f[i][1] += fkcy - fprod1[1] * Tap;
f[i][2] += fkcz - fprod1[2] * Tap;
f[j][0] -= fkcx;
f[j][1] -= fkcy;
f[j][2] -= fkcz;
// calculate the forces acted on the neighbors of atom i from atom j
ILP_neighs_i = ILP_firstneigh[i];
for (kk = 0; kk < ILP_numneigh[i]; kk++) {
k = ILP_neighs_i[kk];
if (k == i) continue;
// derivatives of the product of rij and ni respect to rk, k=0,1,2, where atom k is the neighbors of atom i
dprodnorm1[0] = dnormal[i][0][kk][0] * delx + dnormal[i][1][kk][0] * dely +
dnormal[i][2][kk][0] * delz;
dprodnorm1[1] = dnormal[i][0][kk][1] * delx + dnormal[i][1][kk][1] * dely +
dnormal[i][2][kk][1] * delz;
dprodnorm1[2] = dnormal[i][0][kk][2] * delx + dnormal[i][1][kk][2] * dely +
dnormal[i][2][kk][2] * delz;
fk[0] = (-prodnorm1 * dprodnorm1[0] * fpair1) * Tap;
fk[1] = (-prodnorm1 * dprodnorm1[1] * fpair1) * Tap;
fk[2] = (-prodnorm1 * dprodnorm1[2] * fpair1) * Tap;
f[k][0] += fk[0];
f[k][1] += fk[1];
f[k][2] += fk[2];
delki[0] = x[k][0] - x[i][0];
delki[1] = x[k][1] - x[i][1];
delki[2] = x[k][2] - x[i][2];
if (evflag)
ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0],
delki[1], delki[2]);
}
if (eflag) pvector[1] += evdwl = Tap * Vilp;
if (evflag)
ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fkcx, fkcy, fkcz, delx, dely, delz);
}
}
}
}
/* ----------------------------------------------------------------------
create ILP neighbor list from main neighbor list to calculate normals
------------------------------------------------------------------------- */
void PairILPTMD::ILP_neigh()
{
int i, j, l, ii, jj, ll, n, allnum, jnum, itype, jtype, ltype, imol, jmol, count;
double xtmp, ytmp, ztmp, delx, dely, delz, deljx, deljy, deljz, rsq, rsqlj;
int *ilist, *jlist, *numneigh, **firstneigh;
int *neighsort;
int neighptr[10], check[10];
double **x = atom->x;
int *type = atom->type;
if (atom->nmax > maxlocal) {
maxlocal = atom->nmax;
memory->destroy(ILP_numneigh);
memory->sfree(ILP_firstneigh);
memory->create(ILP_numneigh, maxlocal, "ILPTMD:numneigh");
ILP_firstneigh = (int **) memory->smalloc(maxlocal * sizeof(int *), "ILPTMD:firstneigh");
}
allnum = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// store all ILP neighs of owned and ghost atoms
// scan full neighbor list of I
ipage->reset();
for (ii = 0; ii < allnum; ii++) {
i = ilist[ii];
//initialize varibles
n = 0;
neighsort = ipage->vget();
for (ll = 0; ll < 10; ll++) {
neighptr[ll] = -1;
check[ll] = -1;
}
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = map[type[i]];
imol = atom->molecule[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = map[type[j]];
jmol = atom->molecule[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
// check if the atom i is TMD, i.e., Mo/S/W/Se
if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 ||
strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0 ||
strcmp(elements[itype], "Te") == 0) {
if (rsq != 0 && rsq < cutILPsq[itype][jtype] && imol == jmol && type[i] == type[j]) {
neighptr[n++] = j;
}
} else { // atom i is C, B, N or H.
if (rsq != 0 && rsq < cutILPsq[itype][jtype] && imol == jmol) { neighptr[n++] = j; }
}
} // loop over jj
// if atom i is Mo/W/S/Se/Te, then sorting the orders of neighbors
if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 ||
strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0 ||
strcmp(elements[itype], "Te") == 0) {
// initialize neighsort
for (ll = 0; ll < n; ll++) {
neighsort[ll] = neighptr[ll];
check[ll] = neighptr[ll];
}
// select the first neighbor of atomi
if (n == Nnei) {
neighsort[0] = neighptr[0];
check[0] = -1;
} else if (n < Nnei && n > 0) {
for (jj = 0; jj < n; jj++) { //identify the first neighbor
j = neighptr[jj];
jtype = map[type[j]];
count = 0;
for (ll = 0; ll < n; ll++) {
l = neighptr[ll];
ltype = map[type[l]];
if (l == j) continue;
deljx = x[l][0] - x[j][0];
deljy = x[l][1] - x[j][1];
deljz = x[l][2] - x[j][2];
rsqlj = deljx * deljx + deljy * deljy + deljz * deljz;
if (rsqlj != 0 && rsqlj < cutILPsq[ltype][jtype]) { count++; }
}
if (count == 1) {
neighsort[0] = neighptr[jj];
check[jj] = -1;
break;
}
} // end of idenfying the first neighbor
} else if (n > Nnei) {
fprintf(screen, "Molecule ID = %d\n", imol);
fprintf(screen, "Atom Type = %d\n", type[i]);
fprintf(screen, "Neinum = %d\n", n);
error->one(FLERR,
"There are too many neighbors for TMD atoms, please check your configuration");
}
// sort the order of neighbors of atomi
for (jj = 0; jj < n; jj++) {
j = neighsort[jj];
jtype = map[type[j]];
ll = 0;
while (ll < n) {
l = neighptr[ll];
if (check[ll] == -1) {
ll++;
continue;
}
ltype = map[type[l]];
deljx = x[l][0] - x[j][0];
deljy = x[l][1] - x[j][1];
deljz = x[l][2] - x[j][2];
rsqlj = deljx * deljx + deljy * deljy + deljz * deljz;
if (rsqlj != 0 && rsqlj < cutILPsq[ltype][jtype]) {
neighsort[jj + 1] = l;
check[ll] = -1;
break;
}
ll++;
}
} // end of sorting the order of neighbors
} else { // for B/N/C/H atoms
if (n > 3)
error->one(
FLERR,
"There are too many neighbors for B/N/C/H atoms, please check your configuration");
for (ll = 0; ll < n; ll++) { neighsort[ll] = neighptr[ll]; }
}
ILP_firstneigh[i] = neighsort;
ILP_numneigh[i] = n;
ipage->vgot(n);
if (ipage->status()) error->one(FLERR, "Neighbor list overflow, boost neigh_modify one");
}
}
/* ----------------------------------------------------------------------
Calculate the normals for each atom
------------------------------------------------------------------------- */
void PairILPTMD::calc_normal()
{
int i, j, ii, jj, inum, jnum;
int cont, id, ip, m, k, itype;
double nn, xtp, ytp, ztp, delx, dely, delz, nn2;
int *ilist, *jlist;
double Nave[3], dni[3], dpvdri[3][3];
double **x = atom->x;
int *type = atom->type;
memory->destroy(dnn);
memory->destroy(vect);
memory->destroy(pvet);
memory->destroy(dpvet1);
memory->destroy(dpvet2);
memory->destroy(dNave);
memory->create(dnn, Nnei, 3, "ILPTMD:dnn");
memory->create(vect, Nnei, 3, "ILPTMD:vect");
memory->create(pvet, Nnei, 3, "ILPTMD:pvet");
memory->create(dpvet1, Nnei, 3, 3, "ILPTMD:dpvet1");
memory->create(dpvet2, Nnei, 3, 3, "ILPTMD:dpvet2");
memory->create(dNave, 3, Nnei, 3, "ILPTMD:dNave");
// grow normal array if necessary
if (atom->nmax > nmax) {
memory->destroy(normal);
memory->destroy(dnormal);
memory->destroy(dnormdri);
nmax = atom->nmax;
memory->create(normal, nmax, 3, "ILPTMD:normal");
memory->create(dnormdri, nmax, 3, 3, "ILPTMD:dnormdri");
memory->create(dnormal, nmax, 3, Nnei, 3, "ILPTMD:dnormal");
}
inum = list->inum;
ilist = list->ilist;
//Calculate normals
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = map[type[i]];
xtp = x[i][0];
ytp = x[i][1];
ztp = x[i][2];
// Initialize the arrays
for (id = 0; id < 3; id++) {
Nave[id] = 0.0;
dni[id] = 0.0;
normal[i][id] = 0.0;
for (ip = 0; ip < 3; ip++) {
dpvdri[ip][id] = 0.0;
dnormdri[i][ip][id] = 0.0;
for (m = 0; m < Nnei; m++) {
dnn[m][id] = 0.0;
vect[m][id] = 0.0;
pvet[m][id] = 0.0;
dpvet1[m][ip][id] = 0.0;
dpvet2[m][ip][id] = 0.0;
dNave[id][m][ip] = 0.0;
dnormal[i][id][m][ip] = 0.0;
}
}
}
cont = 0;
jlist = ILP_firstneigh[i];
jnum = ILP_numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = x[j][0] - xtp;
dely = x[j][1] - ytp;
delz = x[j][2] - ztp;
vect[cont][0] = delx;
vect[cont][1] = dely;
vect[cont][2] = delz;
cont++;
}
//############################ For the dangling atoms ############################
if (cont <= 1) {
normal[i][0] = 0.0;
normal[i][1] = 0.0;
normal[i][2] = 1.0;
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormdri[i][id][ip] = 0.0;
for (m = 0; m < Nnei; m++) { dnormal[i][id][m][ip] = 0.0; }
}
}
}
//############################ For the edge atoms of TMD ################################
else if (cont > 1 && cont < Nnei) {
if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 ||
strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0) {
// derivatives of Ni[l] respect to the cont neighbors
for (k = 0; k < cont - 1; k++) {
for (ip = 0; ip < 3; ip++) {
pvet[k][ip] = vect[k][modulo(ip + 1, 3)] * vect[k + 1][modulo(ip + 2, 3)] -
vect[k][modulo(ip + 2, 3)] * vect[k + 1][modulo(ip + 1, 3)];
}
// dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l
// derivatives respect to atom l
// dNik,x/drl
dpvet1[k][0][0] = 0.0;
dpvet1[k][0][1] = vect[modulo(k + 1, Nnei)][2];
dpvet1[k][0][2] = -vect[modulo(k + 1, Nnei)][1];
// dNik,y/drl
dpvet1[k][1][0] = -vect[modulo(k + 1, Nnei)][2];
dpvet1[k][1][1] = 0.0;
dpvet1[k][1][2] = vect[modulo(k + 1, Nnei)][0];
// dNik,z/drl
dpvet1[k][2][0] = vect[modulo(k + 1, Nnei)][1];
dpvet1[k][2][1] = -vect[modulo(k + 1, Nnei)][0];
dpvet1[k][2][2] = 0.0;
// dpvet2[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l+1
// derivatives respect to atom l+1
// dNik,x/drl+1
dpvet2[k][0][0] = 0.0;
dpvet2[k][0][1] = -vect[modulo(k, Nnei)][2];
dpvet2[k][0][2] = vect[modulo(k, Nnei)][1];
// dNik,y/drl+1
dpvet2[k][1][0] = vect[modulo(k, Nnei)][2];
dpvet2[k][1][1] = 0.0;
dpvet2[k][1][2] = -vect[modulo(k, Nnei)][0];
// dNik,z/drl+1
dpvet2[k][2][0] = -vect[modulo(k, Nnei)][1];
dpvet2[k][2][1] = vect[modulo(k, Nnei)][0];
dpvet2[k][2][2] = 0.0;
}
// average the normal vectors by using the Nnei neighboring planes
for (ip = 0; ip < 3; ip++) {
Nave[ip] = 0.0;
for (k = 0; k < cont - 1; k++) { Nave[ip] += pvet[k][ip]; }
Nave[ip] /= (cont - 1);
}
// the magnitude of the normal vector
nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2];
nn = sqrt(nn2);
if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
// the unit normal vector
normal[i][0] = Nave[0] / nn;
normal[i][1] = Nave[1] / nn;
normal[i][2] = Nave[2] / nn;
// derivatives of non-normalized normal vector, dNave:3xcontx3 array
// dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
for (m = 0; m < cont; m++) {
if (m == 0) {
dNave[id][m][ip] = dpvet1[m][id][ip] / (cont - 1);
} else if (m == cont - 1) {
dNave[id][m][ip] = dpvet2[m - 1][id][ip] / (cont - 1);
} else { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m
dNave[id][m][ip] = (dpvet1[m][id][ip] + dpvet2[m - 1][id][ip]) / (cont - 1);
}
}
}
}
// derivatives of nn, dnn:contx3 vector
// dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...Nnei-1; id=0,1,2
// r[m][id]: the id's component of atom m
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] +
Nave[2] * dNave[2][m][id]) /
nn;
}
}
// dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2.
// for atom m, which is a neighbor atom of atom i, m = 0,...,Nnei-1
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2;
}
}
}
// Calculte dNave/dri, defined as dpvdri
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dpvdri[id][ip] = 0.0;
for (k = 0; k < cont; k++) { dpvdri[id][ip] -= dNave[id][k][ip]; }
}
}
// derivatives of nn, dnn:3x1 vector
dni[0] = (Nave[0] * dpvdri[0][0] + Nave[1] * dpvdri[1][0] + Nave[2] * dpvdri[2][0]) / nn;
dni[1] = (Nave[0] * dpvdri[0][1] + Nave[1] * dpvdri[1][1] + Nave[2] * dpvdri[2][1]) / nn;
dni[2] = (Nave[0] * dpvdri[0][2] + Nave[1] * dpvdri[1][2] + Nave[2] * dpvdri[2][2]) / nn;
// derivatives of unit vector ni respect to ri, the result is 3x3 matrix
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormdri[i][id][ip] = dpvdri[id][ip] / nn - Nave[id] * dni[ip] / nn2;
}
}
} // for TMD
//############################ For the edge & bulk atoms of GrhBN ################################
else {
if (cont == 2) {
for (ip = 0; ip < 3; ip++) {
pvet[0][ip] = vect[0][modulo(ip + 1, 3)] * vect[1][modulo(ip + 2, 3)] -
vect[0][modulo(ip + 2, 3)] * vect[1][modulo(ip + 1, 3)];
}
// dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l
// derivatives respect to atom l
// dNik,x/drl
dpvet1[0][0][0] = 0.0;
dpvet1[0][0][1] = vect[1][2];
dpvet1[0][0][2] = -vect[1][1];
// dNi0,y/drl
dpvet1[0][1][0] = -vect[1][2];
dpvet1[0][1][1] = 0.0;
dpvet1[0][1][2] = vect[1][0];
// dNi0,z/drl
dpvet1[0][2][0] = vect[1][1];
dpvet1[0][2][1] = -vect[1][0];
dpvet1[0][2][2] = 0.0;
// dpvet2[0][l][ip]: the derivatve of the 0 (=0,...cont-1)th Ni0 respect to the ip component of atom l+1
// derivatives respect to atom l+1
// dNi0,x/drl+1
dpvet2[0][0][0] = 0.0;
dpvet2[0][0][1] = -vect[0][2];
dpvet2[0][0][2] = vect[0][1];
// dNi0,y/drl+1
dpvet2[0][1][0] = vect[0][2];
dpvet2[0][1][1] = 0.0;
dpvet2[0][1][2] = -vect[0][0];
// dNi0,z/drl+1
dpvet2[0][2][0] = -vect[0][1];
dpvet2[0][2][1] = vect[0][0];
dpvet2[0][2][2] = 0.0;
for (ip = 0; ip < 3; ip++) { Nave[ip] += pvet[0][ip]; }
// the magnitude of the normal vector
nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2];
nn = sqrt(nn2);
if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
// the unit normal vector
normal[i][0] = Nave[0] / nn;
normal[i][1] = Nave[1] / nn;
normal[i][2] = Nave[2] / nn;
// derivatives of non-normalized normal vector, dNave:3xcontx3 array
// dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
for (m = 0; m < cont; m++) {
if (m == 0) {
dNave[id][m][ip] = dpvet1[m][id][ip] / (cont - 1);
} else if (m == cont - 1) {
dNave[id][m][ip] = dpvet2[m - 1][id][ip] / (cont - 1);
} else { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m
dNave[id][m][ip] = (dpvet1[m][id][ip] + dpvet2[m - 1][id][ip]) / (cont - 1);
}
}
}
}
// derivatives of nn, dnn:contx3 vector
// dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...Nnei-1; id=0,1,2
// r[m][id]: the id's component of atom m
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] +
Nave[2] * dNave[2][m][id]) /
nn;
}
}
// dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2.
// for atom m, which is a neighbor atom of atom i, m = 0,...,Nnei-1
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2;
}
}
}
// Calculte dNave/dri, defined as dpvdri
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dpvdri[id][ip] = 0.0;
for (k = 0; k < cont; k++) { dpvdri[id][ip] -= dNave[id][k][ip]; }
}
}
// derivatives of nn, dnn:3x1 vector
dni[0] = (Nave[0] * dpvdri[0][0] + Nave[1] * dpvdri[1][0] + Nave[2] * dpvdri[2][0]) / nn;
dni[1] = (Nave[0] * dpvdri[0][1] + Nave[1] * dpvdri[1][1] + Nave[2] * dpvdri[2][1]) / nn;
dni[2] = (Nave[0] * dpvdri[0][2] + Nave[1] * dpvdri[1][2] + Nave[2] * dpvdri[2][2]) / nn;
// derivatives of unit vector ni respect to ri, the result is 3x3 matrix
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormdri[i][id][ip] = dpvdri[id][ip] / nn - Nave[id] * dni[ip] / nn2;
}
}
} // end of cont == 2
else if (cont == 3) {
// derivatives of Ni[l] respect to the 3 neighbors
for (k = 0; k < 3; k++) {
for (ip = 0; ip < 3; ip++) {
pvet[k][ip] = vect[modulo(k, 3)][modulo(ip + 1, 3)] *
vect[modulo(k + 1, 3)][modulo(ip + 2, 3)] -
vect[modulo(k, 3)][modulo(ip + 2, 3)] * vect[modulo(k + 1, 3)][modulo(ip + 1, 3)];
}
// dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l
// derivatives respect to atom l
// dNik,x/drl
dpvet1[k][0][0] = 0.0;
dpvet1[k][0][1] = vect[modulo(k + 1, 3)][2];
dpvet1[k][0][2] = -vect[modulo(k + 1, 3)][1];
// dNik,y/drl
dpvet1[k][1][0] = -vect[modulo(k + 1, 3)][2];
dpvet1[k][1][1] = 0.0;
dpvet1[k][1][2] = vect[modulo(k + 1, 3)][0];
// dNik,z/drl
dpvet1[k][2][0] = vect[modulo(k + 1, 3)][1];
dpvet1[k][2][1] = -vect[modulo(k + 1, 3)][0];
dpvet1[k][2][2] = 0.0;
// dpvet2[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l+1
// derivatives respect to atom l+1
// dNik,x/drl+1
dpvet2[k][0][0] = 0.0;
dpvet2[k][0][1] = -vect[modulo(k, 3)][2];
dpvet2[k][0][2] = vect[modulo(k, 3)][1];
// dNik,y/drl+1
dpvet2[k][1][0] = vect[modulo(k, 3)][2];
dpvet2[k][1][1] = 0.0;
dpvet2[k][1][2] = -vect[modulo(k, 3)][0];
// dNik,z/drl+1
dpvet2[k][2][0] = -vect[modulo(k, 3)][1];
dpvet2[k][2][1] = vect[modulo(k, 3)][0];
dpvet2[k][2][2] = 0.0;
}
// average the normal vectors by using the 3 neighboring planes
for (ip = 0; ip < 3; ip++) {
Nave[ip] = 0.0;
for (k = 0; k < 3; k++) { Nave[ip] += pvet[k][ip]; }
Nave[ip] /= 3;
}
// the magnitude of the normal vector
nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2];
nn = sqrt(nn2);
if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
// the unit normal vector
normal[i][0] = Nave[0] / nn;
normal[i][1] = Nave[1] / nn;
normal[i][2] = Nave[2] / nn;
// for the central atoms, dnormdri is always zero
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) { dnormdri[i][id][ip] = 0.0; }
}
// derivatives of non-normalized normal vector, dNave:3x3x3 array
// dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
for (
m = 0; m < 3;
m++) { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m
dNave[id][m][ip] =
(dpvet1[modulo(m, 3)][id][ip] + dpvet2[modulo(m - 1, 3)][id][ip]) / 3;
}
}
}
// derivatives of nn, dnn:3x3 vector
// dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...3-1; id=0,1,2
// r[m][id]: the id's component of atom m
for (m = 0; m < 3; m++) {
for (id = 0; id < 3; id++) {
dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] +
Nave[2] * dNave[2][m][id]) /
nn;
}
}
// dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2.
// for atom m, which is a neighbor atom of atom i, m = 0,...,3-1
for (m = 0; m < 3; m++) {
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2;
}
}
}
} // end of cont == 3
else
error->one(FLERR,
"There are too many neighbors for calculating normals of B/N/C/H atoms");
} // for B/N/C/H
} // end of if(cont<Nnei)
//############################ For the bulk atoms ################################
else if (cont == Nnei) {
// derivatives of Ni[l] respect to the Nnei neighbors
for (k = 0; k < Nnei; k++) {
for (ip = 0; ip < 3; ip++) {
pvet[k][ip] = vect[modulo(k, Nnei)][modulo(ip + 1, 3)] *
vect[modulo(k + 1, Nnei)][modulo(ip + 2, 3)] -
vect[modulo(k, Nnei)][modulo(ip + 2, 3)] *
vect[modulo(k + 1, Nnei)][modulo(ip + 1, 3)];
}
// dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l
// derivatives respect to atom l
// dNik,x/drl
dpvet1[k][0][0] = 0.0;
dpvet1[k][0][1] = vect[modulo(k + 1, Nnei)][2];
dpvet1[k][0][2] = -vect[modulo(k + 1, Nnei)][1];
// dNik,y/drl
dpvet1[k][1][0] = -vect[modulo(k + 1, Nnei)][2];
dpvet1[k][1][1] = 0.0;
dpvet1[k][1][2] = vect[modulo(k + 1, Nnei)][0];
// dNik,z/drl
dpvet1[k][2][0] = vect[modulo(k + 1, Nnei)][1];
dpvet1[k][2][1] = -vect[modulo(k + 1, Nnei)][0];
dpvet1[k][2][2] = 0.0;
// dpvet2[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l+1
// derivatives respect to atom l+1
// dNik,x/drl+1
dpvet2[k][0][0] = 0.0;
dpvet2[k][0][1] = -vect[modulo(k, Nnei)][2];
dpvet2[k][0][2] = vect[modulo(k, Nnei)][1];
// dNik,y/drl+1
dpvet2[k][1][0] = vect[modulo(k, Nnei)][2];
dpvet2[k][1][1] = 0.0;
dpvet2[k][1][2] = -vect[modulo(k, Nnei)][0];
// dNik,z/drl+1
dpvet2[k][2][0] = -vect[modulo(k, Nnei)][1];
dpvet2[k][2][1] = vect[modulo(k, Nnei)][0];
dpvet2[k][2][2] = 0.0;
}
// average the normal vectors by using the Nnei neighboring planes
for (ip = 0; ip < 3; ip++) {
Nave[ip] = 0.0;
for (k = 0; k < Nnei; k++) { Nave[ip] += pvet[k][ip]; }
Nave[ip] /= Nnei;
}
// the magnitude of the normal vector
nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2];
nn = sqrt(nn2);
if (nn == 0.0) error->one(FLERR, "The magnitude of the normal vector is zero");
// the unit normal vector
normal[i][0] = Nave[0] / nn;
normal[i][1] = Nave[1] / nn;
normal[i][2] = Nave[2] / nn;
// for the central atoms, dnormdri is always zero
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) { dnormdri[i][id][ip] = 0.0; }
}
// derivatives of non-normalized normal vector, dNave:3xNneix3 array
// dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
for (
m = 0; m < Nnei;
m++) { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m
dNave[id][m][ip] =
(dpvet1[modulo(m, Nnei)][id][ip] + dpvet2[modulo(m - 1, Nnei)][id][ip]) / Nnei;
}
}
}
// derivatives of nn, dnn:Nneix3 vector
// dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...Nnei-1; id=0,1,2
// r[m][id]: the id's component of atom m
for (m = 0; m < Nnei; m++) {
for (id = 0; id < 3; id++) {
dnn[m][id] =
(Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] + Nave[2] * dNave[2][m][id]) /
nn;
}
}
// dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2.
// for atom m, which is a neighbor atom of atom i, m = 0,...,Nnei-1
for (m = 0; m < Nnei; m++) {
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2;
}
}
}
} else {
error->one(FLERR, "There are too many neighbors for calculating normals of TMD atoms");
} // end of four cases of cont
} // end of i loop
}

70
src/INTERLAYER/pair_ilp_tmd.h Normal file
View File

@ -0,0 +1,70 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(ilp/tmd,PairILPTMD);
// clang-format on
#else
#ifndef LMP_PAIR_ILP_TMD_H
#define LMP_PAIR_ILP_TMD_H
#include "pair_ilp_graphene_hbn.h"
namespace LAMMPS_NS {
class PairILPTMD : public PairILPGrapheneHBN {
public:
PairILPTMD(class LAMMPS *);
protected:
void settings(int, char **) override;
void ILP_neigh() override;
void calc_normal() override;
void calc_FRep(int, int) override;
/**************************************************************/
/* modulo operation with cycling around range */
inline int modulo(int k, int range)
{
if (k < 0) k += range;
return k % range;
}
/**************************************************************/
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: All pair coeffs are not set
All pair coefficients must be set in the data file or by the
pair_coeff command before running a simulation.
*/

12
src/INTERLAYER/pair_kolmogorov_crespi_z.cpp
View File

@ -28,8 +28,8 @@
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "potential_file_reader.h"
#include <cmath>
@ -160,15 +160,15 @@ void PairKolmogorovCrespiZ::compute(int eflag, int vflag)
if (evflag) {
ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
if (vflag_either) {
double fi[3],fj[3];
double fi[3], fj[3];
fi[0] = delx * fpair1;
fi[1] = dely * fpair1;
fi[2] = 0;
fj[0] = -delx * fpair1;
fj[1] = -dely * fpair1;
fj[2] = 0;
v_tally2_newton(i,fi,x[i]);
v_tally2_newton(j,fj,x[j]);
v_tally2_newton(i, fi, x[i]);
v_tally2_newton(j, fj, x[j]);
}
}
}
@ -246,9 +246,9 @@ void PairKolmogorovCrespiZ::coeff(int narg, char **arg)
void PairKolmogorovCrespiZ::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style kolmogorov/crespi/z requires newton pair on");
error->all(FLERR, "Pair style kolmogorov/crespi/z requires newton pair on");
neighbor->request(this,instance_me);
neighbor->request(this, instance_me);
}
/* ----------------------------------------------------------------------

245
src/INTERLAYER/pair_saip_metal.cpp Normal file
View File

@ -0,0 +1,245 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Wengen Ouyang (Wuhan University)
e-mail: w.g.ouyang at gmail dot com
This is a full version of the potential described in
[Ouyang et al, J. Chem. Theory Comput. 17, 7215-7223 (2021)]
------------------------------------------------------------------------- */
#include "pair_saip_metal.h"
#include "atom.h"
#include "citeme.h"
#include "error.h"
#include "force.h"
#include "interlayer_taper.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace InterLayer;
#define MAXLINE 1024
#define DELTA 4
#define PGDELTA 1
static const char cite_saip[] =
"saip/metal potential doi.org/10.1021/acs.jctc.1c00622\n"
"@Article{Ouyang2021\n"
" author = {W. Ouyang, O. Hod, and R. Guerra},\n"
" title = {Registry-Dependent Potential for Interfaces of Gold with Graphitic Systems},\n"
" journal = {J. Chem. Theory Comput.},\n"
" volume = 17,\n"
" pages = {7215-7223}\n"
" year = 2021,\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
PairSAIPMETAL::PairSAIPMETAL(LAMMPS *lmp) : PairILPGrapheneHBN(lmp)
{
variant = SAIP_METAL;
single_enable = 0;
if (lmp->citeme) lmp->citeme->add(cite_saip);
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairSAIPMETAL::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command");
if (!utils::strmatch(force->pair_style, "^hybrid/overlay"))
error->all(FLERR, "Pair style saip/metal must be used as sub-style with hybrid/overlay");
cut_global = utils::numeric(FLERR, arg[0], false, lmp);
if (narg == 2) tap_flag = utils::numeric(FLERR, arg[1], false, lmp);
}
/* ----------------------------------------------------------------------
Repulsive forces and energy
------------------------------------------------------------------------- */
void PairSAIPMETAL::calc_FRep(int eflag, int /* vflag */)
{
int i, j, ii, jj, inum, jnum, itype, jtype, k, kk;
double prodnorm1, fkcx, fkcy, fkcz, filp;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, fpair1;
double rsq, r, Rcut, rhosq1, exp0, exp1, Tap, dTap, Vilp;
double frho1, Erep, fsum, rdsq1;
int *ilist, *jlist, *numneigh, **firstneigh;
int *ILP_neighs_i;
evdwl = 0.0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double dprodnorm1[3] = {0.0, 0.0, 0.0};
double fp1[3] = {0.0, 0.0, 0.0};
double fprod1[3] = {0.0, 0.0, 0.0};
double delki[3] = {0.0, 0.0, 0.0};
double fk[3] = {0.0, 0.0, 0.0};
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
//calculate exp(-lambda*(r-z0))*[epsilon/2 + f(rho_ij)]
// loop over neighbors of owned atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
// only include the interaction between different layers
if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) {
int iparam_ij = elem2param[map[itype]][map[jtype]];
Param &p = params[iparam_ij];
r = sqrt(rsq);
// turn on/off taper function
if (tap_flag) {
Rcut = sqrt(cutsq[itype][jtype]);
Tap = calc_Tap(r, Rcut);
dTap = calc_dTap(r, Rcut);
} else {
Tap = 1.0;
dTap = 0.0;
}
// for atoms in bulk materials
if (strcmp(elements[map[itype]], "C") != 0 && strcmp(elements[map[itype]], "H") != 0 &&
strcmp(elements[map[itype]], "B") != 0 && strcmp(elements[map[itype]], "N") != 0) {
// Set ni along rij
exp0 = exp(-p.lambda * (r - p.z0));
frho1 = 1.0 * p.C;
Erep = 0.5 * p.epsilon + frho1;
Vilp = exp0 * Erep;
// derivatives
fpair = p.lambda * exp0 / r * Erep;
filp = fpair * Tap - Vilp * dTap / r;
fkcx = delx * filp;
fkcy = dely * filp;
fkcz = delz * filp;
f[i][0] += fkcx;
f[i][1] += fkcy;
f[i][2] += fkcz;
f[j][0] -= fkcx;
f[j][1] -= fkcy;
f[j][2] -= fkcz;
if (eflag) pvector[1] += evdwl = Tap * Vilp;
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, filp, delx, dely, delz);
} else { // for atoms in 2D materials
// Calculate the transverse distance
prodnorm1 = normal[i][0] * delx + normal[i][1] * dely + normal[i][2] * delz;
rhosq1 = rsq - prodnorm1 * prodnorm1; // rho_ij
rdsq1 = rhosq1 * p.delta2inv; // (rho_ij/delta)^2
// store exponents
exp0 = exp(-p.lambda * (r - p.z0));
exp1 = exp(-rdsq1);
frho1 = exp1 * p.C;
Erep = 0.5 * p.epsilon + frho1;
Vilp = exp0 * Erep;
// derivatives
fpair = p.lambda * exp0 / r * Erep;
fpair1 = 2.0 * exp0 * frho1 * p.delta2inv;
fsum = fpair + fpair1;
// derivatives of the product of rij and ni, the result is a vector
dprodnorm1[0] =
dnormdri[0][0][i] * delx + dnormdri[1][0][i] * dely + dnormdri[2][0][i] * delz;
dprodnorm1[1] =
dnormdri[0][1][i] * delx + dnormdri[1][1][i] * dely + dnormdri[2][1][i] * delz;
dprodnorm1[2] =
dnormdri[0][2][i] * delx + dnormdri[1][2][i] * dely + dnormdri[2][2][i] * delz;
fp1[0] = prodnorm1 * normal[i][0] * fpair1;
fp1[1] = prodnorm1 * normal[i][1] * fpair1;
fp1[2] = prodnorm1 * normal[i][2] * fpair1;
fprod1[0] = prodnorm1 * dprodnorm1[0] * fpair1;
fprod1[1] = prodnorm1 * dprodnorm1[1] * fpair1;
fprod1[2] = prodnorm1 * dprodnorm1[2] * fpair1;
fkcx = (delx * fsum - fp1[0]) * Tap - Vilp * dTap * delx / r;
fkcy = (dely * fsum - fp1[1]) * Tap - Vilp * dTap * dely / r;
fkcz = (delz * fsum - fp1[2]) * Tap - Vilp * dTap * delz / r;
f[i][0] += fkcx - fprod1[0] * Tap;
f[i][1] += fkcy - fprod1[1] * Tap;
f[i][2] += fkcz - fprod1[2] * Tap;
f[j][0] -= fkcx;
f[j][1] -= fkcy;
f[j][2] -= fkcz;
// calculate the forces acted on the neighbors of atom i from atom j
ILP_neighs_i = ILP_firstneigh[i];
for (kk = 0; kk < ILP_numneigh[i]; kk++) {
k = ILP_neighs_i[kk];
if (k == i) continue;
// derivatives of the product of rij and ni respect to rk, k=0,1,2, where atom k is the neighbors of atom i
dprodnorm1[0] = dnormal[0][0][kk][i] * delx + dnormal[1][0][kk][i] * dely +
dnormal[2][0][kk][i] * delz;
dprodnorm1[1] = dnormal[0][1][kk][i] * delx + dnormal[1][1][kk][i] * dely +
dnormal[2][1][kk][i] * delz;
dprodnorm1[2] = dnormal[0][2][kk][i] * delx + dnormal[1][2][kk][i] * dely +
dnormal[2][2][kk][i] * delz;
fk[0] = (-prodnorm1 * dprodnorm1[0] * fpair1) * Tap;
fk[1] = (-prodnorm1 * dprodnorm1[1] * fpair1) * Tap;
fk[2] = (-prodnorm1 * dprodnorm1[2] * fpair1) * Tap;
f[k][0] += fk[0];
f[k][1] += fk[1];
f[k][2] += fk[2];
delki[0] = x[k][0] - x[i][0];
delki[1] = x[k][1] - x[i][1];
delki[2] = x[k][2] - x[i][2];
if (evflag)
ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0],
delki[1], delki[2]);
}
if (eflag) pvector[1] += evdwl = Tap * Vilp;
if (evflag)
ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fkcx, fkcy, fkcz, delx, dely, delz);
} // end of speration atoms for bulk materials
} // end of rsq < cutoff
} // loop over jj
} // loop over ii
}

58
src/INTERLAYER/pair_saip_metal.h Normal file
View File

@ -0,0 +1,58 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(saip/metal,PairSAIPMETAL);
// clang-format on
#else
#ifndef LMP_PAIR_SAIP_METAL_H
#define LMP_PAIR_SAIP_METAL_H
#include "pair_ilp_graphene_hbn.h"
namespace LAMMPS_NS {
class PairSAIPMETAL : public PairILPGrapheneHBN {
public:
PairSAIPMETAL(class LAMMPS *);
protected:
void settings(int, char **) override;
void calc_FRep(int, int) override;
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: All pair coeffs are not set
All pair coefficients must be set in the data file or by the
pair_coeff command before running a simulation.
*/