diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 9f2bdbce79..9bbe216dec 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -256,6 +256,7 @@ OPT. * :doc:`rann ` * :doc:`reaxff (ko) ` * :doc:`rebo (io) ` + * :doc:`rebomos (o) ` * :doc:`resquared (go) ` * :doc:`saip/metal (t) ` * :doc:`sdpd/taitwater/isothermal ` diff --git a/doc/src/pair_rebomos.rst b/doc/src/pair_rebomos.rst new file mode 100644 index 0000000000..9f4b8006c1 --- /dev/null +++ b/doc/src/pair_rebomos.rst @@ -0,0 +1,150 @@ +.. index:: pair_style rebomos +.. index:: pair_style rebomos/omp + +pair_style rebomos command +========================== + +Accelerator Variants: *rebomos/omp* + +Syntax +"""""" + +.. code-block:: LAMMPS + + pair_style rebomos + +* rebomos = name of this pair style + +Examples +"""""""" + +.. code-block:: LAMMPS + + pair_style rebomos + pair_coeff * * ../potentials/MoS.rebomos Mo S + +Example input scripts available: examples/threebody/ + +Description +""""""""""" + +.. versionadded:: TBD + +The *rebomos* pair style computes the interactions between molybdenum +and sulfur atoms :ref:`(Stewart) ` utilizing an adaptive +interatomic reactive empirical bond order potential that is similar in +form to the AIREBO potential :ref:`(Stuart) `. The potential +is based on an earlier parameterizations for :math:`\text{MoS}_2` +developed by :ref:`(Liang) `. + +The REBOMoS potential consists of two terms: + +.. math:: + + E & = \frac{1}{2} \sum_i \sum_{j \neq i} + \left[ E^{\text{REBO}}_{ij} + E^{\text{LJ}}_{ij} \right] \\ + +The :math:`E^{\text{REBO}}` term describes the covalently bonded +interactions between Mo and S atoms while the :math:`E^{\text{LJ}}` term +describes longer range dispersion forces between layers. A cubic spline +function is applied to smoothly switch between covalent bonding at short +distances to dispersion interactions at longer distances. This allows +the model to capture bond formation and breaking events which may occur +between adjacent MoS2 layers, edges, defects, and more. + +---------- + +Only a single pair_coeff command is used with the *rebomos* pair style +which specifies an REBOMoS potential file with parameters for Mo and S. +These are mapped to LAMMPS atom types by specifying N additional +arguments after the filename in the pair_coeff command, where N is the +number of LAMMPS atom types: + +* filename +* :math:`N` element names = mapping of REBOMoS elements to atom types + +See the :doc:`pair_coeff ` page for alternate ways +to specify the path for the potential file. + +As an example, if your LAMMPS simulation has three atom types and you want +the first two to be Mo, and the third to be S, you would use the following +pair_coeff command: + +.. code-block:: LAMMPS + + pair_coeff * * MoS.rebomos Mo Mo S + +The first 2 arguments must be \* \* so as to span all LAMMPS atom types. +The first two Mo arguments map LAMMPS atom types 1 and 2 to the Mo +element in the REBOMoS file. The final S argument maps LAMMPS atom type +3 to the S element in the REBOMoS file. If a mapping value is specified +as NULL, the mapping is not performed. This can be used when a +*rebomos* potential is used as part of the *hybrid* pair style. The +NULL values are placeholders for atom types that will be used with other +potentials. + +---------- + +.. include:: accel_styles.rst + +---------- + +Mixing, shift, table, tail correction, restart, rRESPA info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +This pair style does not support the :doc:`pair_modify ` +mix, shift, table, and tail options. + +This pair style does not write their information to :doc:`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. + +This pair styles can only be used via the *pair* keyword of the +:doc:`run_style respa ` command. It does not support the +*inner*, *middle*, *outer* keywords. + +Restrictions +"""""""""""" + +This pair style is part of the MANYBODY package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. + +These pair potentials require the :doc:`newton ` setting to be +"on" for pair interactions. + +The MoS.rebomos potential file provided with LAMMPS (see the potentials +directory) is parameterized for metal :doc:`units `. You can use +the *rebomos* pair style with any LAMMPS units setting, but you would +need to create your own REBOMoS potential file with coefficients listed +in the appropriate units. + +The pair style provided here **only** supports potential files parameterized +for the elements molybdenum and sulfur (designated with "Mo" and "S" in the +*pair_coeff* command. Using potential files for other elements will trigger +an error. + +Related commands +"""""""""""""""" + +:doc:`pair_coeff `, :doc:`pair style rebo ` + +Default +""""""" + +none + +---------- + +.. _Stewart: + +**(Steward)** Stewart, Spearot, Modelling Simul. Mater. Sci. Eng. 21, 045003, (2013). + +.. _Stuart2: + +**(Stuart)** Stuart, Tutein, Harrison, J Chem Phys, 112, 6472-6486, (2000). + +.. _Liang: + +**(Liang)** Liang, Phillpot, Sinnott Phys. Rev. B79 245110, (2009), Erratum: Phys. Rev. B85 199903(E), (2012) diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index a2467bff2b..53bf269e1c 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -333,6 +333,7 @@ accelerated styles exist. * :doc:`rann ` - * :doc:`reaxff ` - ReaxFF potential * :doc:`rebo ` - second generation REBO potential of Brenner +* :doc:`rebomos ` - REBOMoS potential for MoS2 * :doc:`resquared ` - Everaers RE-Squared ellipsoidal potential * :doc:`saip/metal ` - interlayer potential for hetero-junctions formed with hexagonal 2D materials and metal surfaces * :doc:`sdpd/taitwater/isothermal ` - smoothed dissipative particle dynamics for water at isothermal conditions diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 4f5fe6fdaf..e46fb6ca97 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -2266,6 +2266,7 @@ morris Morriss morse Morteza +MoS Mosayebi Moseler Moskalev @@ -3070,6 +3071,7 @@ reaxff ReaxFF REAXFF rebo +rebomos recurse recursing Ree @@ -3597,6 +3599,7 @@ tesselation tesselations Tetot tex +textrm tfac tfmc tfMC diff --git a/examples/threebody/MoS.rebomos b/examples/threebody/MoS.rebomos new file mode 120000 index 0000000000..6146c74c24 --- /dev/null +++ b/examples/threebody/MoS.rebomos @@ -0,0 +1 @@ +../../potentials/MoS.rebomos \ No newline at end of file diff --git a/examples/threebody/in.mos2-bulk b/examples/threebody/in.mos2-bulk new file mode 100644 index 0000000000..032e71fce8 --- /dev/null +++ b/examples/threebody/in.mos2-bulk @@ -0,0 +1,35 @@ +units metal + +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 & + a2 -1.5964590311 2.7651481541 0.0000000000 & + a3 0.0000000000 0.0000000000 13.9827680588 & + basis 0.0000000000 0.000000000 $(3.0/4.0) & + basis 0.0000000000 0.000000000 $(1.0/4.0) & + basis $(2.0/3.0) $(1.0/3.0) 0.862008989 & + basis $(1.0/3.0) $(2.0/3.0) 0.137990996 & + basis $(1.0/3.0) $(2.0/3.0) 0.362008989 & + basis $(2.0/3.0) $(1.0/3.0) 0.637991011 & + origin 0.1 0.1 0.1 + +region box prism 0 4 0 8 0 1 -2.0 0.0 0.0 +create_box 2 box +create_atoms 2 box & + basis 1 1 & + basis 2 1 & + basis 3 2 & + basis 4 2 & + basis 5 2 & + basis 6 2 + +mass 1 95.95 #Mo +mass 2 32.065 #S + +pair_style rebomos +pair_coeff * * MoS.rebomos Mo S + +thermo_style custom step temp press pe ke cellgamma vol +thermo 10 +#dump 1 all atom 10 MoS.lammpstrj +fix 1 all nve +run 20 + diff --git a/examples/threebody/in.mos2.rebomos b/examples/threebody/in.mos2.rebomos new file mode 100644 index 0000000000..ca91f67003 --- /dev/null +++ b/examples/threebody/in.mos2.rebomos @@ -0,0 +1,31 @@ +# monolayer MoS2 +units metal +boundary p p f +processors * * 1 +atom_modify map array + +atom_style atomic +read_data single_layer_MoS2.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 + +########################## Define potentials ################################ +pair_style rebomos +pair_coeff * * MoS.rebomos Mo S S +######################################################################### + +### Simulation settings #### +timestep 0.001 +velocity all create 300.0 12345 loop geom + +############################ + +# Output +thermo 500 +thermo_style custom step etotal pe ke temp +thermo_modify lost warn + +###### Run molecular dynamics ###### +fix thermostat all nve +run 5000 diff --git a/examples/threebody/log.22Feb24.mos2-bulk.g++.1 b/examples/threebody/log.22Feb24.mos2-bulk.g++.1 new file mode 100644 index 0000000000..8218026f3d --- /dev/null +++ b/examples/threebody/log.22Feb24.mos2-bulk.g++.1 @@ -0,0 +1,85 @@ +LAMMPS (7 Feb 2024 - Development - patch_7Feb2024_update1-73-g36fa601fe0) + using 1 OpenMP thread(s) per MPI task +units metal + +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 $(3.0/4.0) basis 0.0000000000 0.000000000 $(1.0/4.0) basis $(2.0/3.0) $(1.0/3.0) 0.862008989 basis $(1.0/3.0) $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 $(1.0/4.0) basis $(2.0/3.0) $(1.0/3.0) 0.862008989 basis $(1.0/3.0) $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis $(2.0/3.0) $(1.0/3.0) 0.862008989 basis $(1.0/3.0) $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 $(1.0/3.0) 0.862008989 basis $(1.0/3.0) $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis $(1.0/3.0) $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 0.66666666666666662966 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 0.66666666666666662966 0.137990996 basis 0.33333333333333331483 $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 0.66666666666666662966 0.137990996 basis 0.33333333333333331483 0.66666666666666662966 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 0.66666666666666662966 0.137990996 basis 0.33333333333333331483 0.66666666666666662966 0.362008989 basis 0.66666666666666662966 $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 0.66666666666666662966 0.137990996 basis 0.33333333333333331483 0.66666666666666662966 0.362008989 basis 0.66666666666666662966 0.33333333333333331483 0.637991011 origin 0.1 0.1 0.1 +Lattice spacing in x,y,z = 4.7867748 2.7651482 13.982768 + +region box prism 0 4 0 8 0 1 -2.0 0.0 0.0 +create_box 2 box +Created triclinic box = (0 0 0) to (19.147099 22.121185 13.982768) with tilt (-9.5735495 0 0) + 1 by 1 by 1 MPI processor grid +create_atoms 2 box basis 1 1 basis 2 1 basis 3 2 basis 4 2 basis 5 2 basis 6 2 +Created 288 atoms + using lattice units in triclinic box = (0 0 0) to (19.147099 22.121185 13.982768) with tilt (-9.5735495 0 0) + create_atoms CPU = 0.000 seconds + +mass 1 95.95 #Mo +mass 2 32.065 #S + +pair_style rebomos +pair_coeff * * MoS.rebomos Mo S +Reading rebomos potential file MoS.rebomos with DATE: 2013-11-04 + +thermo_style custom step temp press pe ke cellgamma vol +thermo 10 +#dump 1 all atom 10 MoS.lammpstrj +fix 1 all nve +run 20 +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 13.4 + ghost atom cutoff = 13.4 + binsize = 6.7, bins = 5 4 3 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair rebomos, perpetual + attributes: full, newton on, ghost + pair build: full/bin/ghost + stencil: full/ghost/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 4.996 | 4.996 | 4.996 Mbytes + Step Temp Press PotEng KinEng CellGamma Volume + 0 0 28799.53 -2061.6112 0 113.40187 5922.4926 + 10 80.776057 13540.088 -2064.6132 2.9966028 113.40187 5922.4926 + 20 146.17503 -20669.371 -2067.0428 5.4227518 113.40187 5922.4926 +Loop time of 0.058071 on 1 procs for 20 steps with 288 atoms + +Performance: 29.757 ns/day, 0.807 hours/ns, 344.406 timesteps/s, 99.189 katom-step/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 | 0.057666 | 0.057666 | 0.057666 | 0.0 | 99.30 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 0.00024654 | 0.00024654 | 0.00024654 | 0.0 | 0.42 +Output | 2.3975e-05 | 2.3975e-05 | 2.3975e-05 | 0.0 | 0.04 +Modify | 3.8394e-05 | 3.8394e-05 | 3.8394e-05 | 0.0 | 0.07 +Other | | 9.596e-05 | | | 0.17 + +Nlocal: 288 ave 288 max 288 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 4285 ave 4285 max 4285 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: 142848 ave 142848 max 142848 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 142848 +Ave neighs/atom = 496 +Neighbor list builds = 0 +Dangerous builds = 0 + +Total wall time: 0:00:00 diff --git a/examples/threebody/log.22Feb24.mos2-bulk.g++.4 b/examples/threebody/log.22Feb24.mos2-bulk.g++.4 new file mode 100644 index 0000000000..0b9cd3ed8a --- /dev/null +++ b/examples/threebody/log.22Feb24.mos2-bulk.g++.4 @@ -0,0 +1,85 @@ +LAMMPS (7 Feb 2024 - Development - patch_7Feb2024_update1-73-g36fa601fe0) + using 1 OpenMP thread(s) per MPI task +units metal + +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 $(3.0/4.0) basis 0.0000000000 0.000000000 $(1.0/4.0) basis $(2.0/3.0) $(1.0/3.0) 0.862008989 basis $(1.0/3.0) $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 $(1.0/4.0) basis $(2.0/3.0) $(1.0/3.0) 0.862008989 basis $(1.0/3.0) $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis $(2.0/3.0) $(1.0/3.0) 0.862008989 basis $(1.0/3.0) $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 $(1.0/3.0) 0.862008989 basis $(1.0/3.0) $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis $(1.0/3.0) $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 $(2.0/3.0) 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 0.66666666666666662966 0.137990996 basis $(1.0/3.0) $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 0.66666666666666662966 0.137990996 basis 0.33333333333333331483 $(2.0/3.0) 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 0.66666666666666662966 0.137990996 basis 0.33333333333333331483 0.66666666666666662966 0.362008989 basis $(2.0/3.0) $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 0.66666666666666662966 0.137990996 basis 0.33333333333333331483 0.66666666666666662966 0.362008989 basis 0.66666666666666662966 $(1.0/3.0) 0.637991011 origin 0.1 0.1 0.1 +lattice custom 1.0 a1 3.1903157234 0.0000000000 0.0000000000 a2 -1.5964590311 2.7651481541 0.0000000000 a3 0.0000000000 0.0000000000 13.9827680588 basis 0.0000000000 0.000000000 0.75 basis 0.0000000000 0.000000000 0.25 basis 0.66666666666666662966 0.33333333333333331483 0.862008989 basis 0.33333333333333331483 0.66666666666666662966 0.137990996 basis 0.33333333333333331483 0.66666666666666662966 0.362008989 basis 0.66666666666666662966 0.33333333333333331483 0.637991011 origin 0.1 0.1 0.1 +Lattice spacing in x,y,z = 4.7867748 2.7651482 13.982768 + +region box prism 0 4 0 8 0 1 -2.0 0.0 0.0 +create_box 2 box +Created triclinic box = (0 0 0) to (19.147099 22.121185 13.982768) with tilt (-9.5735495 0 0) + 2 by 2 by 1 MPI processor grid +create_atoms 2 box basis 1 1 basis 2 1 basis 3 2 basis 4 2 basis 5 2 basis 6 2 +Created 288 atoms + using lattice units in triclinic box = (0 0 0) to (19.147099 22.121185 13.982768) with tilt (-9.5735495 0 0) + create_atoms CPU = 0.000 seconds + +mass 1 95.95 #Mo +mass 2 32.065 #S + +pair_style rebomos +pair_coeff * * MoS.rebomos Mo S +Reading rebomos potential file MoS.rebomos with DATE: 2013-11-04 + +thermo_style custom step temp press pe ke cellgamma vol +thermo 10 +#dump 1 all atom 10 MoS.lammpstrj +fix 1 all nve +run 20 +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 13.4 + ghost atom cutoff = 13.4 + binsize = 6.7, bins = 5 4 3 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair rebomos, perpetual + attributes: full, newton on, ghost + pair build: full/bin/ghost + stencil: full/ghost/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 4.15 | 4.151 | 4.151 Mbytes + Step Temp Press PotEng KinEng CellGamma Volume + 0 0 28799.53 -2061.6112 0 113.40187 5922.4926 + 10 80.776057 13540.088 -2064.6132 2.9966028 113.40187 5922.4926 + 20 146.17503 -20669.371 -2067.0428 5.4227518 113.40187 5922.4926 +Loop time of 0.0219485 on 4 procs for 20 steps with 288 atoms + +Performance: 78.730 ns/day, 0.305 hours/ns, 911.225 timesteps/s, 262.433 katom-step/s +96.3% CPU use with 4 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0.018118 | 0.019372 | 0.020087 | 0.5 | 88.26 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 0.0015635 | 0.0023195 | 0.0035967 | 1.6 | 10.57 +Output | 2.5017e-05 | 4.6834e-05 | 0.00010543 | 0.0 | 0.21 +Modify | 1.3954e-05 | 1.423e-05 | 1.4594e-05 | 0.0 | 0.06 +Other | | 0.0001957 | | | 0.89 + +Nlocal: 72 ave 72 max 72 min +Histogram: 4 0 0 0 0 0 0 0 0 0 +Nghost: 2771.5 ave 2775 max 2768 min +Histogram: 2 0 0 0 0 0 0 0 0 2 +Neighs: 0 ave 0 max 0 min +Histogram: 4 0 0 0 0 0 0 0 0 0 +FullNghs: 35712 ave 35712 max 35712 min +Histogram: 4 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 142848 +Ave neighs/atom = 496 +Neighbor list builds = 0 +Dangerous builds = 0 + +Total wall time: 0:00:00 diff --git a/examples/threebody/log.22Feb24.mos2.rebomos.g++.1 b/examples/threebody/log.22Feb24.mos2.rebomos.g++.1 new file mode 100644 index 0000000000..f7c5b3c74d --- /dev/null +++ b/examples/threebody/log.22Feb24.mos2.rebomos.g++.1 @@ -0,0 +1,95 @@ +LAMMPS (7 Feb 2024 - Development - patch_7Feb2024_update1-73-g36fa601fe0) + using 1 OpenMP thread(s) per MPI task +# monolayer MoS2 +units metal +boundary p p f +processors * * 1 +atom_modify map array + +atom_style atomic +read_data single_layer_MoS2.data +Reading data file ... + triclinic box = (0 0 -100) to (51.15232 44.299209 100) with tilt (25.57616 0 0) +WARNING: Triclinic box skew is large. LAMMPS will run inefficiently. (src/domain.cpp:219) + 1 by 1 by 1 MPI processor grid + reading atoms ... + 768 atoms + read_data CPU = 0.002 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 + +########################## Define potentials ################################ +pair_style rebomos +pair_coeff * * MoS.rebomos Mo S S +Reading rebomos potential file MoS.rebomos with DATE: 2013-11-04 +######################################################################### + +### Simulation settings #### +timestep 0.001 +velocity all create 300.0 12345 loop geom + +############################ + +# Output +thermo 500 +thermo_style custom step etotal pe ke temp +thermo_modify lost warn + +###### Run molecular dynamics ###### +fix thermostat all nve +run 5000 +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 13.4 + ghost atom cutoff = 13.4 + binsize = 6.7, bins = 12 7 30 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair rebomos, perpetual + attributes: full, newton on, ghost + pair build: full/bin/ghost + stencil: full/ghost/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 4.473 | 4.473 | 4.473 Mbytes + Step TotEng PotEng KinEng Temp + 0 -5466.9785 -5496.7212 29.742759 300 + 500 -5466.964 -5482.6985 15.734505 158.7059 + 1000 -5466.9615 -5480.9492 13.98763 141.08607 + 1500 -5466.964 -5482.6912 15.727258 158.63281 + 2000 -5466.9657 -5483.3606 16.394878 165.36675 + 2500 -5466.9624 -5481.6253 14.662948 147.89765 + 3000 -5466.9642 -5482.7515 15.7873 159.23842 + 3500 -5466.9654 -5483.3789 16.413502 165.5546 + 4000 -5466.9628 -5481.848 14.885236 150.13977 + 4500 -5466.9648 -5483.5045 16.539775 166.82825 + 5000 -5466.9649 -5483.4932 16.528298 166.71249 +Loop time of 19.1009 on 1 procs for 5000 steps with 768 atoms + +Performance: 22.617 ns/day, 1.061 hours/ns, 261.768 timesteps/s, 201.038 katom-step/s +99.9% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 19.042 | 19.042 | 19.042 | 0.0 | 99.69 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 0.018451 | 0.018451 | 0.018451 | 0.0 | 0.10 +Output | 0.00015575 | 0.00015575 | 0.00015575 | 0.0 | 0.00 +Modify | 0.023931 | 0.023931 | 0.023931 | 0.0 | 0.13 +Other | | 0.01658 | | | 0.09 + +Nlocal: 768 ave 768 max 768 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 1158 ave 1158 max 1158 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: 141824 ave 141824 max 141824 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 141824 +Ave neighs/atom = 184.66667 +Neighbor list builds = 0 +Dangerous builds = 0 +Total wall time: 0:00:19 diff --git a/examples/threebody/log.22Feb24.mos2.rebomos.g++.4 b/examples/threebody/log.22Feb24.mos2.rebomos.g++.4 new file mode 100644 index 0000000000..dc1cfa84d4 --- /dev/null +++ b/examples/threebody/log.22Feb24.mos2.rebomos.g++.4 @@ -0,0 +1,95 @@ +LAMMPS (7 Feb 2024 - Development - patch_7Feb2024_update1-73-g36fa601fe0) + using 1 OpenMP thread(s) per MPI task +# monolayer MoS2 +units metal +boundary p p f +processors * * 1 +atom_modify map array + +atom_style atomic +read_data single_layer_MoS2.data +Reading data file ... + triclinic box = (0 0 -100) to (51.15232 44.299209 100) with tilt (25.57616 0 0) +WARNING: Triclinic box skew is large. LAMMPS will run inefficiently. (src/domain.cpp:219) + 2 by 2 by 1 MPI processor grid + reading atoms ... + 768 atoms + read_data CPU = 0.002 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 + +########################## Define potentials ################################ +pair_style rebomos +pair_coeff * * MoS.rebomos Mo S S +Reading rebomos potential file MoS.rebomos with DATE: 2013-11-04 +######################################################################### + +### Simulation settings #### +timestep 0.001 +velocity all create 300.0 12345 loop geom + +############################ + +# Output +thermo 500 +thermo_style custom step etotal pe ke temp +thermo_modify lost warn + +###### Run molecular dynamics ###### +fix thermostat all nve +run 5000 +Neighbor list info ... + update: every = 1 steps, delay = 0 steps, check = yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 13.4 + ghost atom cutoff = 13.4 + binsize = 6.7, bins = 12 7 30 + 1 neighbor lists, perpetual/occasional/extra = 1 0 0 + (1) pair rebomos, perpetual + attributes: full, newton on, ghost + pair build: full/bin/ghost + stencil: full/ghost/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 4.045 | 4.045 | 4.045 Mbytes + Step TotEng PotEng KinEng Temp + 0 -5466.9785 -5496.7212 29.742759 300 + 500 -5466.964 -5482.6985 15.734505 158.7059 + 1000 -5466.9615 -5480.9492 13.98763 141.08607 + 1500 -5466.964 -5482.6912 15.727258 158.63281 + 2000 -5466.9657 -5483.3606 16.394878 165.36675 + 2500 -5466.9624 -5481.6253 14.662948 147.89765 + 3000 -5466.9642 -5482.7515 15.7873 159.23842 + 3500 -5466.9654 -5483.3789 16.413502 165.5546 + 4000 -5466.9628 -5481.848 14.885236 150.13977 + 4500 -5466.9648 -5483.5045 16.539775 166.82825 + 5000 -5466.9649 -5483.4932 16.528298 166.71249 +Loop time of 5.69326 on 4 procs for 5000 steps with 768 atoms + +Performance: 75.879 ns/day, 0.316 hours/ns, 878.231 timesteps/s, 674.482 katom-step/s +98.6% CPU use with 4 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 5.2611 | 5.3666 | 5.4358 | 3.0 | 94.26 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 0.23476 | 0.30106 | 0.40642 | 12.8 | 5.29 +Output | 0.00014996 | 0.0004478 | 0.0013353 | 0.0 | 0.01 +Modify | 0.0068861 | 0.0069917 | 0.0072247 | 0.2 | 0.12 +Other | | 0.01814 | | | 0.32 + +Nlocal: 192 ave 194 max 190 min +Histogram: 1 0 0 0 0 2 0 0 0 1 +Nghost: 710 ave 712 max 708 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: 35456 ave 35824 max 35088 min +Histogram: 1 0 0 0 0 2 0 0 0 1 + +Total # of neighbors = 141824 +Ave neighs/atom = 184.66667 +Neighbor list builds = 0 +Dangerous builds = 0 +Total wall time: 0:00:05 diff --git a/potentials/MoS.rebomos b/potentials/MoS.rebomos new file mode 100644 index 0000000000..ea96981a1e --- /dev/null +++ b/potentials/MoS.rebomos @@ -0,0 +1,65 @@ +# DATE: 2013-11-04 UNITS: metal CONTRIBUTOR: J Stewart, K Dang, D Spearot (UArk) CITATION: Stewart J A and Spearot D E, Modelling Simul. Mater. Sci. Eng. 21. +# MoS-S REBO Brenner/Sinnot Potential as published in +# Liang T, Phillpot S R and Sinnott S B (2009) Phys. Rev. B79 245110, Erratum: Phys. Rev. B85 199903(E). + +3.50 rcmin_MM +2.75 rcmin_MS +2.30 rcmin_SS +3.80 rcmax_MM +3.05 rcmax_MS +3.00 rcmax_SS +3.419129390005910 Q_MM +1.505537839153790 Q_MS +0.254959104053671 Q_SS +1.07500712999340 alpha_MM +1.19267902218820 alpha_MS +1.10775022439715 alpha_SS +179.008013654688 A_MM +575.509677721866 A_MS +1228.43233679426 A_SS +706.247903589221 BIJc_MM1 +1344.46820036159 BIJc_MS1 +1498.64815404145 BIJc_SS1 +1.16100322369589 Beta_MM1 +1.26973752204290 Beta_MS1 +1.12673623610320 Beta_SS1 +0.1326842550663270 M_b0 +-0.007642788338017 M_b1 +0.0341395775059370 M_b2 +0.2523050971380870 M_b3 +0.1227287372225670 M_b4 +-0.361387798398897 M_b5 +-0.282577591351457 M_b6 +0.120194301035280 M_bg0 +0.045238287358190 M_bg1 +0.067922807244030 M_bg2 +-0.03672511378682 M_bg3 +0.107516477513860 M_bg4 +0.004964711984940 M_bg5 +-0.12997598358652 M_bg6 +0.006848761596750 S_b0 +-0.02389964401024 S_b1 +0.137457353311170 S_b2 +0.033016467497740 S_b3 +-0.31064291544850 S_b4 +-0.08550273135791 S_b5 +0.149252790306880 S_b6 +-0.2850852 S_bg0 +1.67102480 S_bg1 +-3.5678516 S_bg2 +3.45054990 S_bg3 +-1.2186289 S_bg4 +0.0 S_bg5 +0.0 S_bg6 +0.138040769883614 M_a0 +0.803625443023934 M_a1 +0.292412960851064 M_a2 +0.640588078946224 M_a3 +0.062978539843324 S_a0 +2.478617619878250 S_a1 +0.036666243238154 S_a2 +2.386431372486710 S_a3 +0.00058595 epsilon_MM +0.01386 epsilon_SS +4.200 sigma_MM +3.130 sigma_SS diff --git a/potentials/README b/potentials/README index c234f5f48b..2cb4a383c5 100644 --- a/potentials/README +++ b/potentials/README @@ -84,7 +84,7 @@ Au_u3 = Gold universal 3 The suffix of each file indicates the pair style it is used with: adp ADP angular dependent potential -airebo AI-REBO and REBO potentials +airebo AI-REBO potentials bop.table BOP potential, tabulated form cdeam concentration-dependent EAM comb COMB potential @@ -107,6 +107,8 @@ nb3b.harmonic nonbonded 3-body harmonic potential pod ML potential with proper orthogonal descriptors (POD) poly polymorphic 3-body potential reax ReaxFF potential (see README.reax for more info) +rebo REBO potentials +rebomos REBOMoS potential smtbq second moment tight binding QEq (SMTBQ) potential snap SNAP potential snapcoeff SNAP potential diff --git a/src/.gitignore b/src/.gitignore index 9e2c36bdd0..ba4c4c05b0 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1371,6 +1371,8 @@ /pair_reaxff.h /pair_rebo.cpp /pair_rebo.h +/pair_rebomos.cpp +/pair_rebomos.h /pair_resquared.cpp /pair_resquared.h /pair_saip_metal.cpp diff --git a/src/MANYBODY/pair_rebomos.cpp b/src/MANYBODY/pair_rebomos.cpp new file mode 100644 index 0000000000..d6cea004a7 --- /dev/null +++ b/src/MANYBODY/pair_rebomos.cpp @@ -0,0 +1,1131 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + References: + + This code: + Stewart J A and Spearot D E (2013) Atomistic simulations of nanoindentation on the basal plane of crystalline molybdenum disulfide. Modelling Simul. Mater. Sci. Eng. 21. + + Based on: + Liang T, Phillpot S R and Sinnott S B (2009) Parameterization of a reactive many-body potential for Mo2S systems. Phys. Rev. B79 245110. + Liang T, Phillpot S R and Sinnott S B (2012) Erratum: Parameterization of a reactive many-body potential for Mo-S systems. (Phys. Rev. B79 245110 (2009)) Phys. Rev. B85 199903(E). + + LAMMPS file contributing authors: James Stewart, Khanh Dang and Douglas Spearot (University of Arkansas) +------------------------------------------------------------------------- */ + +// clang-format on + +#include "pair_rebomos.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "math_special.h" +#include "memory.h" +#include "my_page.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "potential_file_reader.h" +#include "text_file_reader.h" + +#include +#include + +using namespace LAMMPS_NS; +using MathSpecial::cube; +using MathSpecial::powint; +using MathSpecial::square; + +static constexpr double TOL = 1.0e-9; +static constexpr int PGDELTA = 1; + +/* ---------------------------------------------------------------------- */ + +PairREBOMoS::PairREBOMoS(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 0; + restartinfo = 0; + one_coeff = 1; + ghostneigh = 1; + manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; + + maxlocal = 0; + REBO_numneigh = nullptr; + REBO_firstneigh = nullptr; + ipage = nullptr; + pgsize = oneatom = 0; + nM = nS = nullptr; +} + +// clang-format off + +/* ---------------------------------------------------------------------- + Check if allocated, since class can be destructed when incomplete +------------------------------------------------------------------------- */ + +PairREBOMoS::~PairREBOMoS() +{ + memory->destroy(REBO_numneigh); + memory->sfree(REBO_firstneigh); + delete[] ipage; + memory->destroy(nM); + memory->destroy(nS); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(cutghost); + + memory->destroy(lj1); + memory->destroy(lj2); + memory->destroy(lj3); + memory->destroy(lj4); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairREBOMoS::compute(int eflag, int vflag) +{ + ev_init(eflag,vflag); + + REBO_neigh(); + FREBO(eflag); + FLJ(eflag); + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairREBOMoS::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutghost,n+1,n+1,"pair:cutghost"); + + // only sized by M,S = 2 types + + memory->create(lj1,2,2,"pair:lj1"); + memory->create(lj2,2,2,"pair:lj2"); + memory->create(lj3,2,2,"pair:lj3"); + memory->create(lj4,2,2,"pair:lj4"); + + map = new int[n+1]; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairREBOMoS::settings(int narg, char ** /* arg */) +{ + if (narg != 0) error->all(FLERR,"Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairREBOMoS::coeff(int narg, char **arg) +{ + if (!allocated) allocate(); + + if (narg != 3 + atom->ntypes) + error->all(FLERR,"Incorrect args for pair coefficients"); + + // insure I,J args are * * + + if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) + error->all(FLERR,"Incorrect args for pair coefficients"); + + // read args that map atom types to Mo and S + // map[i] = which element (0,1) the Ith atom type is, -1 if NULL + + for (int i = 3; i < narg; i++) { + if (strcmp(arg[i],"NULL") == 0) { + map[i-2] = -1; + continue; + } else if (strcmp(arg[i],"Mo") == 0) { + map[i-2] = 0; + } else if (strcmp(arg[i],"M") == 0) { // backward compatibility + map[i-2] = 0; + } else if (strcmp(arg[i],"S") == 0) { + map[i-2] = 1; + } else error->all(FLERR,"Incorrect args for pair coefficients"); + } + + // read potential file and initialize fitting splines + + read_file(arg[2]); + + // clear setflag since coeff() called once with I,J = * * + + int n = atom->ntypes; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + // set setflag i,j for type pairs where both are mapped to elements + + int count = 0; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + if (map[i] >= 0 && map[j] >= 0) { + setflag[i][j] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairREBOMoS::init_style() +{ + if (atom->tag_enable == 0) + error->all(FLERR,"Pair style REBOMoS requires atom IDs"); + if (force->newton_pair == 0) + error->all(FLERR,"Pair style REBOMoS requires newton pair on"); + + // need a full neighbor list, including neighbors of ghosts + + neighbor->add_request(this,NeighConst::REQ_FULL|NeighConst::REQ_GHOST); + + // local REBO neighbor list + // create pages if first time or if neighbor pgsize/oneatom has changed + + int create = 0; + if (ipage == nullptr) create = 1; + if (pgsize != neighbor->pgsize) create = 1; + if (oneatom != neighbor->oneatom) create = 1; + + if (create) { + delete[] ipage; + pgsize = neighbor->pgsize; + oneatom = neighbor->oneatom; + + int nmypage= comm->nthreads; + ipage = new MyPage[nmypage]; + for (int i = 0; i < nmypage; i++) + ipage[i].init(oneatom,pgsize,PGDELTA); + } +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairREBOMoS::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + + // convert to Mo,S types + + int ii = map[i]; + int jj = map[j]; + + // use Mo-Mo values for these cutoffs since M atoms are biggest + + // cut3rebo = 3 REBO distances + + cut3rebo = 3.0 * rcmax[0][0]; + + // cutghost = REBO cutoff used in REBO_neigh() for neighbors of ghosts + + cutghost[i][j] = rcmax[ii][jj]; + lj1[ii][jj] = 48.0 * epsilon[ii][jj] * powint(sigma[ii][jj],12); + lj2[ii][jj] = 24.0 * epsilon[ii][jj] * powint(sigma[ii][jj],6); + lj3[ii][jj] = 4.0 * epsilon[ii][jj] * powint(sigma[ii][jj],12); + lj4[ii][jj] = 4.0 * epsilon[ii][jj] * powint(sigma[ii][jj],6); + + cutghost[j][i] = cutghost[i][j]; + lj1[jj][ii] = lj1[ii][jj]; + lj2[jj][ii] = lj2[ii][jj]; + lj3[jj][ii] = lj3[ii][jj]; + lj4[jj][ii] = lj4[ii][jj]; + + return cut3rebo; +} + +/* ---------------------------------------------------------------------- + create REBO neighbor list from main neighbor list + REBO neighbor list stores neighbors of ghost atoms +------------------------------------------------------------------------- */ + +void PairREBOMoS::REBO_neigh() +{ + int i,j,ii,jj,n,allnum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq,dS; + int *ilist,*jlist,*numneigh,**firstneigh; + int *neighptr; + + double **x = atom->x; + int *type = atom->type; + + if (atom->nmax > maxlocal) { + maxlocal = atom->nmax; + memory->destroy(REBO_numneigh); + memory->sfree(REBO_firstneigh); + memory->destroy(nM); + memory->destroy(nS); + memory->create(REBO_numneigh,maxlocal,"REBOMoS:numneigh"); + REBO_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *), + "REBOMoS:firstneigh"); + memory->create(nM,maxlocal,"REBOMoS:nM"); + memory->create(nS,maxlocal,"REBOMoS:nS"); + } + + allnum = list->inum + list->gnum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // store all REBO neighs of owned and ghost atoms + // scan full neighbor list of I + + ipage->reset(); + + for (ii = 0; ii < allnum; ii++) { + i = ilist[ii]; + + n = 0; + neighptr = ipage->vget(); + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = map[type[i]]; + nM[i] = nS[i] = 0.0; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtype = map[type[j]]; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < rcmaxsq[itype][jtype]) { + neighptr[n++] = j; + if (jtype == 0) + nM[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS); + else + nS[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS); + } + } + + REBO_firstneigh[i] = neighptr; + REBO_numneigh[i] = n; + ipage->vgot(n); + if (ipage->status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } +} + +/* ---------------------------------------------------------------------- + REBO forces and energy +------------------------------------------------------------------------- */ + +void PairREBOMoS::FREBO(int eflag) +{ + int i,j,k,ii,inum,itype,jtype; + tagint itag, jtag; + double delx,dely,delz,evdwl,fpair,xtmp,ytmp,ztmp; + double rsq,rij,wij; + double Qij,Aij,alphaij,VR,pre,dVRdi,VA,bij,dVAdi,dVA; + double dwij,del[3]; + int *ilist,*REBO_neighs; + + evdwl = 0.0; + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + tagint *tag = atom->tag; + int nlocal = atom->nlocal; + + inum = list->inum; + ilist = list->ilist; + + // two-body interactions from REBO neighbor list, skip half of them + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itag = tag[i]; + itype = map[type[i]]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + REBO_neighs = REBO_firstneigh[i]; + + for (k = 0; k < REBO_numneigh[i]; k++) { + j = REBO_neighs[k]; + jtag = tag[j]; + + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + + jtype = map[type[j]]; + + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + rij = sqrt(rsq); + wij = Sp(rij,rcmin[itype][jtype],rcmax[itype][jtype],dwij); + if (wij <= TOL) continue; + + Qij = Q[itype][jtype]; + Aij = A[itype][jtype]; + alphaij = alpha[itype][jtype]; + + VR = wij*(1.0+(Qij/rij)) * Aij*exp(-alphaij*rij); + pre = wij*Aij * exp(-alphaij*rij); + dVRdi = pre * ((-alphaij)-(Qij/rsq)-(Qij*alphaij/rij)); + dVRdi += VR/wij * dwij; + + VA = dVA = 0.0; + VA = -wij * BIJc[itype][jtype] * exp(-Beta[itype][jtype]*rij); + + dVA = -Beta[itype][jtype] * VA; + dVA += VA/wij * dwij; + + del[0] = delx; + del[1] = dely; + del[2] = delz; + bij = bondorder(i,j,del,rij,VA,f); + dVAdi = bij*dVA; + + fpair = -(dVRdi+dVAdi) / rij; + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + + if (eflag) evdwl = VR + bij*VA; + if (evflag) ev_tally(i,j,nlocal,/*newton_pair*/1,evdwl,0.0,fpair,delx,dely,delz); + } + } +} + +/* ---------------------------------------------------------------------- + compute LJ forces and energy +------------------------------------------------------------------------- */ + +void PairREBOMoS::FLJ(int eflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + tagint itag,jtag; + double evdwl,fpair,xtmp,ytmp,ztmp; + double rij,delij[3],rijsq; + double VLJ,dVLJ; + double vdw,dvdw; + double r2inv,r6inv; + int *ilist,*jlist,*numneigh,**firstneigh; + double c2,c3,dr,drp,r6; + + // I-J interaction from full neighbor list + // skip 1/2 of interactions since only consider each pair once + + evdwl = 0.0; + + double **x = atom->x; + double **f = atom->f; + tagint *tag = atom->tag; + int *type = atom->type; + int nlocal = atom->nlocal; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itag = tag[i]; + itype = map[type[i]]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtag = tag[j]; + + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + jtype = map[type[j]]; + + delij[0] = xtmp - x[j][0]; + delij[1] = ytmp - x[j][1]; + delij[2] = ztmp - x[j][2]; + rijsq = delij[0]*delij[0] + delij[1]*delij[1] + delij[2]*delij[2]; + rij = sqrt(rijsq); + + // compute LJ forces and energy + + // Outside Rmax + if (rij > rcLJmax[itype][jtype] || rij < rcLJmin[itype][jtype]){ + VLJ = 0; + dVLJ = 0; + } + + // Inside Rmax and above 0.95*sigma + else if (rij <= rcLJmax[itype][jtype] && rij >= 0.95*sigma[itype][jtype]){ + r2inv = 1.0/rijsq; + r6inv = r2inv*r2inv*r2inv; + VLJ = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + dVLJ = -r6inv*(lj1[itype][jtype]*r6inv - lj2[itype][jtype])/rij; + } + + // Below 0.95*sigma + else if (rij < 0.95*sigma[itype][jtype] && rij >= rcLJmin[itype][jtype]){ + dr = 0.95*sigma[itype][jtype] - rcLJmin[itype][jtype]; + r6 = powint((sigma[itype][jtype]/(0.95*sigma[itype][jtype])),6); + vdw = 4*epsilon[itype][jtype]*r6*(r6 - 1.0); + dvdw = (-4*epsilon[itype][jtype]/(0.95*sigma[itype][jtype]))*r6*(12.0*r6 - 6.0); + c2 = ((3.0/dr)*vdw - dvdw)/dr; + c3 = (vdw/(dr*dr) - c2)/dr; + + drp = rij - rcLJmin[itype][jtype]; + VLJ = drp*drp*(drp*c3 + c2); + dVLJ = drp*(3.0*drp*c3 + 2.0*c2); + } + + fpair = -dVLJ/rij; + f[i][0] += delij[0]*fpair; + f[i][1] += delij[1]*fpair; + f[i][2] += delij[2]*fpair; + f[j][0] -= delij[0]*fpair; + f[j][1] -= delij[1]*fpair; + f[j][2] -= delij[2]*fpair; + + if (eflag) evdwl = VLJ; + if (evflag) ev_tally(i,j,nlocal,/*newton_pair*/1,evdwl,0.0,fpair,delij[0],delij[1],delij[2]); + + } + } +} + +/* ---------------------------------------------------------------------- + Bij function + + The bond order term modified the attractive portion of the REBO + potential based on the number of atoms around a specific pair + and the bond angle between sets of three atoms. + + The functions G(cos(theta)) and P(N) are evaluated and their + derivatives are also computed for use in the force calculation. +------------------------------------------------------------------------- */ + +double PairREBOMoS::bondorder(int i, int j, double rij[3], double rijmag, double VA, double **f) +{ + int atomi,atomj,atomk,atoml; + int k,l; + int itype, jtype, ktype, ltype; + double rik[3], rjl[3], rji[3], rki[3],rlj[3], dwjl, bij; + double NijM,NijS,NjiM,NjiS,wik,dwik,wjl; + double rikmag,rjlmag,cosjik,cosijl,g,tmp2; + double Etmp,pij,tmp,dwij,dS; + double dgdc,pji; + double dcosjikdri[3],dcosijldri[3],dcosjikdrk[3]; + double dp; + double dcosjikdrj[3],dcosijldrj[3],dcosijldrl[3]; + double fi[3],fj[3],fk[3],fl[3]; + double PijS, PjiS; + int *REBO_neighs; + + double **x = atom->x; + int *type = atom->type; + + atomi = i; + atomj = j; + itype = map[type[i]]; + jtype = map[type[j]]; + Sp(rijmag,rcmin[itype][jtype],rcmax[itype][jtype],dwij); + NijM = nM[i]; + NijS = nS[i]; + NjiM = nM[j]; + NjiS = nS[j]; + bij = 0.0; + tmp = 0.0; + tmp2 = 0.0; + dgdc = 0.0; + Etmp = 0.0; + + REBO_neighs = REBO_firstneigh[i]; + for (k = 0; k < REBO_numneigh[i]; k++) { + atomk = REBO_neighs[k]; + if (atomk != atomj) { + ktype = map[type[atomk]]; + rik[0] = x[atomi][0]-x[atomk][0]; + rik[1] = x[atomi][1]-x[atomk][1]; + rik[2] = x[atomi][2]-x[atomk][2]; + rikmag = sqrt((rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2])); + wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dS); + cosjik = ((rij[0]*rik[0])+(rij[1]*rik[1])+(rij[2]*rik[2])) / (rijmag*rikmag); + cosjik = MIN(cosjik,1.0); + cosjik = MAX(cosjik,-1.0); + + // evaluate g and derivative dg + + g = gSpline(cosjik,itype,dgdc); + Etmp = Etmp+(wik*g); + } + } + + dp = 0.0; + PijS = PijSpline(NijM,NijS,itype,dp); + pij = 1.0/sqrt(1.0+Etmp+PijS); + tmp = -0.5*cube(pij); + + // derivative calculations + + REBO_neighs = REBO_firstneigh[i]; + for (k = 0; k < REBO_numneigh[i]; k++) { + atomk = REBO_neighs[k]; + if (atomk != atomj) { + ktype = map[type[atomk]]; + rik[0] = x[atomi][0]-x[atomk][0]; + rik[1] = x[atomi][1]-x[atomk][1]; + rik[2] = x[atomi][2]-x[atomk][2]; + rikmag = sqrt((rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2])); + wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dwik); + cosjik = (rij[0]*rik[0] + rij[1]*rik[1] + rij[2]*rik[2]) / (rijmag*rikmag); + cosjik = MIN(cosjik,1.0); + cosjik = MAX(cosjik,-1.0); + + dcosjikdri[0] = ((rij[0]+rik[0])/(rijmag*rikmag)) - + (cosjik*((rij[0]/(rijmag*rijmag))+(rik[0]/(rikmag*rikmag)))); + dcosjikdri[1] = ((rij[1]+rik[1])/(rijmag*rikmag)) - + (cosjik*((rij[1]/(rijmag*rijmag))+(rik[1]/(rikmag*rikmag)))); + dcosjikdri[2] = ((rij[2]+rik[2])/(rijmag*rikmag)) - + (cosjik*((rij[2]/(rijmag*rijmag))+(rik[2]/(rikmag*rikmag)))); + dcosjikdrk[0] = (-rij[0]/(rijmag*rikmag)) + + (cosjik*(rik[0]/(rikmag*rikmag))); + dcosjikdrk[1] = (-rij[1]/(rijmag*rikmag)) + + (cosjik*(rik[1]/(rikmag*rikmag))); + dcosjikdrk[2] = (-rij[2]/(rijmag*rikmag)) + + (cosjik*(rik[2]/(rikmag*rikmag))); + dcosjikdrj[0] = (-rik[0]/(rijmag*rikmag)) + + (cosjik*(rij[0]/(rijmag*rijmag))); + dcosjikdrj[1] = (-rik[1]/(rijmag*rikmag)) + + (cosjik*(rij[1]/(rijmag*rijmag))); + dcosjikdrj[2] = (-rik[2]/(rijmag*rikmag)) + + (cosjik*(rij[2]/(rijmag*rijmag))); + + g = gSpline(cosjik,itype,dgdc); + tmp2 = VA*0.5*(tmp*wik*dgdc); + fj[0] = -tmp2*dcosjikdrj[0]; + fj[1] = -tmp2*dcosjikdrj[1]; + fj[2] = -tmp2*dcosjikdrj[2]; + fi[0] = -tmp2*dcosjikdri[0]; + fi[1] = -tmp2*dcosjikdri[1]; + fi[2] = -tmp2*dcosjikdri[2]; + fk[0] = -tmp2*dcosjikdrk[0]; + fk[1] = -tmp2*dcosjikdrk[1]; + fk[2] = -tmp2*dcosjikdrk[2]; + + // coordination forces + + // dwik forces (from partial derivative) + + tmp2 = VA*0.5*(tmp*dwik*g)/rikmag; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; + + // PIJ forces (from coordination P(N) term) + + tmp2 = VA*0.5*(tmp*dp*dwik)/rikmag; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; + + // dgdN forces are removed + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2]; + + if (vflag_either) { + rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; + rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2]; + v_tally3(atomi,atomj,atomk,fj,fk,rji,rki); + } + } + } + + // PIJ force contribution additional term + tmp2 = -VA*0.5*(tmp*dp*dwij)/rijmag; + + f[atomi][0] += rij[0]*tmp2; + f[atomi][1] += rij[1]*tmp2; + f[atomi][2] += rij[2]*tmp2; + f[atomj][0] -= rij[0]*tmp2; + f[atomj][1] -= rij[1]*tmp2; + f[atomj][2] -= rij[2]*tmp2; + + if (vflag_either) v_tally2(atomi,atomj,tmp2,rij); + + tmp = 0.0; + tmp2 = 0.0; + Etmp = 0.0; + + REBO_neighs = REBO_firstneigh[j]; + for (l = 0; l < REBO_numneigh[j]; l++) { + atoml = REBO_neighs[l]; + if (atoml != atomi) { + ltype = map[type[atoml]]; + rjl[0] = x[atomj][0]-x[atoml][0]; + rjl[1] = x[atomj][1]-x[atoml][1]; + rjl[2] = x[atomj][2]-x[atoml][2]; + rjlmag = sqrt((rjl[0]*rjl[0])+(rjl[1]*rjl[1])+(rjl[2]*rjl[2])); + wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dS); + cosijl = -1.0*((rij[0]*rjl[0])+(rij[1]*rjl[1])+(rij[2]*rjl[2])) / (rijmag*rjlmag); + cosijl = MIN(cosijl,1.0); + cosijl = MAX(cosijl,-1.0); + + // evaluate g and derivative dg + + g = gSpline(cosijl,jtype,dgdc); + Etmp = Etmp+(wjl*g); + } + } + + dp = 0.0; + PjiS = PijSpline(NjiM,NjiS,jtype,dp); + pji = 1.0/sqrt(1.0+Etmp+PjiS); + tmp = -0.5*cube(pji); + + REBO_neighs = REBO_firstneigh[j]; + for (l = 0; l < REBO_numneigh[j]; l++) { + atoml = REBO_neighs[l]; + if (atoml != atomi) { + ltype = map[type[atoml]]; + rjl[0] = x[atomj][0]-x[atoml][0]; + rjl[1] = x[atomj][1]-x[atoml][1]; + rjl[2] = x[atomj][2]-x[atoml][2]; + rjlmag = sqrt((rjl[0]*rjl[0])+(rjl[1]*rjl[1])+(rjl[2]*rjl[2])); + wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dwjl); + cosijl = (-1.0*((rij[0]*rjl[0])+(rij[1]*rjl[1])+(rij[2]*rjl[2]))) / (rijmag*rjlmag); + cosijl = MIN(cosijl,1.0); + cosijl = MAX(cosijl,-1.0); + + dcosijldri[0] = (-rjl[0]/(rijmag*rjlmag)) - (cosijl*rij[0]/(rijmag*rijmag)); + dcosijldri[1] = (-rjl[1]/(rijmag*rjlmag)) - (cosijl*rij[1]/(rijmag*rijmag)); + dcosijldri[2] = (-rjl[2]/(rijmag*rjlmag)) - (cosijl*rij[2]/(rijmag*rijmag)); + dcosijldrj[0] = ((-rij[0]+rjl[0])/(rijmag*rjlmag)) + + (cosijl*((rij[0]/square(rijmag))-(rjl[0]/(rjlmag*rjlmag)))); + dcosijldrj[1] = ((-rij[1]+rjl[1])/(rijmag*rjlmag)) + + (cosijl*((rij[1]/square(rijmag))-(rjl[1]/(rjlmag*rjlmag)))); + dcosijldrj[2] = ((-rij[2]+rjl[2])/(rijmag*rjlmag)) + + (cosijl*((rij[2]/square(rijmag))-(rjl[2]/(rjlmag*rjlmag)))); + dcosijldrl[0] = (rij[0]/(rijmag*rjlmag))+(cosijl*rjl[0]/(rjlmag*rjlmag)); + dcosijldrl[1] = (rij[1]/(rijmag*rjlmag))+(cosijl*rjl[1]/(rjlmag*rjlmag)); + dcosijldrl[2] = (rij[2]/(rijmag*rjlmag))+(cosijl*rjl[2]/(rjlmag*rjlmag)); + + // evaluate g and derivatives dg + + g = gSpline(cosijl,jtype,dgdc); + tmp2 = VA*0.5*(tmp*wjl*dgdc); + fi[0] = -tmp2*dcosijldri[0]; + fi[1] = -tmp2*dcosijldri[1]; + fi[2] = -tmp2*dcosijldri[2]; + fj[0] = -tmp2*dcosijldrj[0]; + fj[1] = -tmp2*dcosijldrj[1]; + fj[2] = -tmp2*dcosijldrj[2]; + fl[0] = -tmp2*dcosijldrl[0]; + fl[1] = -tmp2*dcosijldrl[1]; + fl[2] = -tmp2*dcosijldrl[2]; + + // coordination forces + + // dwik forces (from partial derivative) + + tmp2 = VA*0.5*(tmp*dwjl*g)/rjlmag; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; + + // PIJ forces (coordination) + + tmp2 = VA*0.5*(tmp*dp*dwjl)/rjlmag; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; + + // dgdN forces are removed + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2]; + + if (vflag_either) { + rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2]; + v_tally3(atomi,atomj,atoml,fi,fl,rij,rlj); + } + } + } + + // PIJ force contribution additional term + + tmp2 = -VA*0.5*(tmp*dp*dwij)/rijmag; + f[atomi][0] += rij[0]*tmp2; + f[atomi][1] += rij[1]*tmp2; + f[atomi][2] += rij[2]*tmp2; + f[atomj][0] -= rij[0]*tmp2; + f[atomj][1] -= rij[1]*tmp2; + f[atomj][2] -= rij[2]*tmp2; + + if (vflag_either) v_tally2(atomi,atomj,tmp2,rij); + + bij = (0.5*(pij+pji)); + return bij; +} + +/* ---------------------------------------------------------------------- + G calculation +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + read REBO potential file +------------------------------------------------------------------------- */ + +void PairREBOMoS::read_file(char *filename) +{ + // REBO Parameters (Mo-S REBO) + + double rcmin_MM,rcmin_MS,rcmin_SS,rcmax_MM,rcmax_MS,rcmax_SS; + double Q_MM,Q_MS,Q_SS,alpha_MM,alpha_MS,alpha_SS,A_MM,A_MS,A_SS; + double BIJc_MM1,BIJc_MS1,BIJc_SS1; + double Beta_MM1,Beta_MS1,Beta_SS1; + double M_bg0,M_bg1,M_bg2,M_bg3,M_bg4,M_bg5,M_bg6; + double S_bg0,S_bg1,S_bg2,S_bg3,S_bg4,S_bg5,S_bg6; + double M_b0,M_b1,M_b2,M_b3,M_b4,M_b5,M_b6; + double S_b0,S_b1,S_b2,S_b3,S_b4,S_b5,S_b6; + double M_a0,M_a1,M_a2,M_a3; + double S_a0,S_a1,S_a2,S_a3; + + // LJ Parameters (Mo-S REBO) + + double epsilon_MM,epsilon_SS; + double sigma_MM,sigma_SS; + + // read file on proc 0 + + if (comm->me == 0) { + PotentialFileReader reader(lmp, filename, "rebomos"); + + // read parameters + + std::vector params { + &rcmin_MM, + &rcmin_MS, + &rcmin_SS, + &rcmax_MM, + &rcmax_MS, + &rcmax_SS, + &Q_MM, + &Q_MS, + &Q_SS, + &alpha_MM, + &alpha_MS, + &alpha_SS, + &A_MM, + &A_MS, + &A_SS, + &BIJc_MM1, + &BIJc_MS1, + &BIJc_SS1, + &Beta_MM1, + &Beta_MS1, + &Beta_SS1, + &M_b0, + &M_b1, + &M_b2, + &M_b3, + &M_b4, + &M_b5, + &M_b6, + &M_bg0, + &M_bg1, + &M_bg2, + &M_bg3, + &M_bg4, + &M_bg5, + &M_bg6, + &S_b0, + &S_b1, + &S_b2, + &S_b3, + &S_b4, + &S_b5, + &S_b6, + &S_bg0, + &S_bg1, + &S_bg2, + &S_bg3, + &S_bg4, + &S_bg5, + &S_bg6, + &M_a0, + &M_a1, + &M_a2, + &M_a3, + &S_a0, + &S_a1, + &S_a2, + &S_a3, + + // LJ parameters + &epsilon_MM, + &epsilon_SS, + &sigma_MM, + &sigma_SS, + }; + + try { + for (auto ¶m : params) { + *param = reader.next_double(); + } + } catch (TokenizerException &e) { + error->one(FLERR, "reading rebomos potential file {}\nREASON: {}\n", filename, e.what()); + } catch (FileReaderException &fre) { + error->one(FLERR, "reading rebomos potential file {}\nREASON: {}\n", filename, fre.what()); + } + + // store read-in values in arrays + + // REBO + + rcmin[0][0] = rcmin_MM; + rcmin[0][1] = rcmin_MS; + rcmin[1][0] = rcmin[0][1]; + rcmin[1][1] = rcmin_SS; + + rcmax[0][0] = rcmax_MM; + rcmax[0][1] = rcmax_MS; + rcmax[1][0] = rcmax[0][1]; + rcmax[1][1] = rcmax_SS; + + rcmaxsq[0][0] = rcmax[0][0]*rcmax[0][0]; + rcmaxsq[1][0] = rcmax[1][0]*rcmax[1][0]; + rcmaxsq[0][1] = rcmax[0][1]*rcmax[0][1]; + rcmaxsq[1][1] = rcmax[1][1]*rcmax[1][1]; + + Q[0][0] = Q_MM; + Q[0][1] = Q_MS; + Q[1][0] = Q[0][1]; + Q[1][1] = Q_SS; + + alpha[0][0] = alpha_MM; + alpha[0][1] = alpha_MS; + alpha[1][0] = alpha[0][1]; + alpha[1][1] = alpha_SS; + + A[0][0] = A_MM; + A[0][1] = A_MS; + A[1][0] = A[0][1]; + A[1][1] = A_SS; + + BIJc[0][0] = BIJc_MM1; + BIJc[0][1] = BIJc_MS1; + BIJc[1][0] = BIJc_MS1; + BIJc[1][1] = BIJc_SS1; + + Beta[0][0] = Beta_MM1; + Beta[0][1] = Beta_MS1; + Beta[1][0] = Beta_MS1; + Beta[1][1] = Beta_SS1; + + b0[0] = M_b0; + b1[0] = M_b1; + b2[0] = M_b2; + b3[0] = M_b3; + b4[0] = M_b4; + b5[0] = M_b5; + b6[0] = M_b6; + + bg0[0] = M_bg0; + bg1[0] = M_bg1; + bg2[0] = M_bg2; + bg3[0] = M_bg3; + bg4[0] = M_bg4; + bg5[0] = M_bg5; + bg6[0] = M_bg6; + + b0[1] = S_b0; + b1[1] = S_b1; + b2[1] = S_b2; + b3[1] = S_b3; + b4[1] = S_b4; + b5[1] = S_b5; + b6[1] = S_b6; + + bg0[1] = S_bg0; + bg1[1] = S_bg1; + bg2[1] = S_bg2; + bg3[1] = S_bg3; + bg4[1] = S_bg4; + bg5[1] = S_bg5; + bg6[1] = S_bg6; + + a0[0] = M_a0; + a1[0] = M_a1; + a2[0] = M_a2; + a3[0] = M_a3; + + a0[1] = S_a0; + a1[1] = S_a1; + a2[1] = S_a2; + a3[1] = S_a3; + + // LJ + + sigma[0][0] = sigma_MM; + sigma[0][1] = (sigma_MM + sigma_SS)/2; + sigma[1][0] = sigma[0][1]; + sigma[1][1] = sigma_SS; + + epsilon[0][0] = epsilon_MM; + epsilon[0][1] = sqrt(epsilon_MM*epsilon_SS); + epsilon[1][0] = epsilon[0][1]; + epsilon[1][1] = epsilon_SS; + + rcLJmin[0][0] = rcmin_MM; + rcLJmin[0][1] = rcmin_MS; + rcLJmin[1][0] = rcmin[0][1]; + rcLJmin[1][1] = rcmin_SS; + + rcLJmax[0][0] = 2.5*sigma[0][0]; + rcLJmax[0][1] = 2.5*sigma[0][1]; + rcLJmax[1][0] = rcLJmax[0][1]; + rcLJmax[1][1] = 2.5*sigma[1][1]; + + rcLJmaxsq[0][0] = rcLJmax[0][0]*rcLJmax[0][0]; + rcLJmaxsq[1][0] = rcLJmax[1][0]*rcLJmax[1][0]; + rcLJmaxsq[0][1] = rcLJmax[0][1]*rcLJmax[0][1]; + rcLJmaxsq[1][1] = rcLJmax[1][1]*rcLJmax[1][1]; + + } + + // broadcast read-in and setup values + + MPI_Bcast(&rcmin[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&rcmax[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&rcmaxsq[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&rcmaxp[0][0],4,MPI_DOUBLE,0,world); + + MPI_Bcast(&Q[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&alpha[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&A[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&BIJc[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&Beta[0][0],4,MPI_DOUBLE,0,world); + + MPI_Bcast(&b0[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&b1[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&b2[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&b3[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&b4[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&b5[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&b6[0],2,MPI_DOUBLE,0,world); + + MPI_Bcast(&a0[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&a1[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&a2[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&a3[0],2,MPI_DOUBLE,0,world); + + MPI_Bcast(&bg0[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&bg1[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&bg2[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&bg3[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&bg4[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&bg5[0],2,MPI_DOUBLE,0,world); + MPI_Bcast(&bg6[0],2,MPI_DOUBLE,0,world); + + MPI_Bcast(&rcLJmin[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&rcLJmax[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&rcLJmaxsq[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&epsilon[0][0],4,MPI_DOUBLE,0,world); + MPI_Bcast(&sigma[0][0],4,MPI_DOUBLE,0,world); +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double PairREBOMoS::memory_usage() +{ + double bytes = 0.0; + bytes += (double)maxlocal * sizeof(int); + bytes += (double)maxlocal * sizeof(int *); + + for (int i = 0; i < comm->nthreads; i++) + bytes += ipage[i].size(); + + bytes += 3.0 * maxlocal * sizeof(double); + return bytes; +} diff --git a/src/MANYBODY/pair_rebomos.h b/src/MANYBODY/pair_rebomos.h new file mode 100644 index 0000000000..4997c65b87 --- /dev/null +++ b/src/MANYBODY/pair_rebomos.h @@ -0,0 +1,220 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + 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(rebomos,PairREBOMoS); +// clang-format on +#else + +#ifndef LMP_PAIR_REBOMOS_H +#define LMP_PAIR_REBOMOS_H + +#include "pair.h" +#include "math_const.h" + +#include + +namespace LAMMPS_NS { + +class PairREBOMoS : public Pair { + public: + PairREBOMoS(class LAMMPS *); + ~PairREBOMoS() override; + void compute(int, int) override; + void settings(int, char **) override; + void coeff(int, char **) override; + void init_style() override; + double init_one(int, int) override; + double memory_usage() override; + + protected: + double cutljrebosq; // cut for when to compute + // REBO neighs of ghost atoms + + double **lj1, **lj2, **lj3, **lj4; // pre-computed LJ coeffs for M,S types + double cut3rebo; // maximum distance for 3rd REBO neigh + + int maxlocal; // size of numneigh, firstneigh arrays + int pgsize; // size of neighbor page + int oneatom; // max # of neighbors for one atom + MyPage *ipage; // neighbor list pages + int *REBO_numneigh; // # of pair neighbors for each atom + int **REBO_firstneigh; // ptr to 1st neighbor of each atom + + double *closestdistsq; // closest owned atom dist to each ghost + double *nM, *nS; // sum of weighting fns with REBO neighs + + double rcmin[2][2], rcmax[2][2], rcmaxsq[2][2], rcmaxp[2][2]; + double Q[2][2], alpha[2][2], A[2][2], BIJc[2][2], Beta[2][2]; + double b0[2], b1[2], b2[2], b3[2], b4[2], b5[2], b6[2]; + double bg0[2], bg1[2], bg2[2], bg3[2], bg4[2], bg5[2], bg6[2]; + double a0[2], a1[2], a2[2], a3[2]; + double rcLJmin[2][2], rcLJmax[2][2], rcLJmaxsq[2][2]; + double epsilon[2][2], sigma[2][2]; + + void REBO_neigh(); + void FREBO(int); + void FLJ(int); + + double bondorder(int, int, double *, double, double, double **); + + inline double gSpline(const double costh, const int typei, double &dgdc) const + { + const double b0i = b0[typei]; + const double b1i = b1[typei]; + const double b2i = b2[typei]; + const double b3i = b3[typei]; + const double b4i = b4[typei]; + const double b5i = b5[typei]; + const double b6i = b6[typei]; + double g = 0.0; + + if (costh >= -1.0 && costh < 0.5) { + g = b6i * costh; + double dg = 6.0 * b6i * costh; + g += b5i; + dg += 5.0 * b5i; + g *= costh; + dg *= costh; + g += b4i; + dg += 4.0 * b4i; + g *= costh; + dg *= costh; + g += b3i; + dg += 3.0 * b3i; + g *= costh; + dg *= costh; + g += b2i; + dg += 2.0 * b2i; + g *= costh; + dg *= costh; + g += b1i; + dg += b1i; + g *= costh; + g += b0i; + dgdc = dg; + + } else if (costh >= 0.5 && costh <= 1.0) { + double gcos = b6i * costh; + double dgcos = 6.0 * b6i * costh; + gcos += b5i; + dgcos += 5.0 * b5i; + gcos *= costh; + dgcos *= costh; + gcos += b4i; + dgcos += 4.0 * b4i; + gcos *= costh; + dgcos *= costh; + gcos += b3i; + dgcos += 3.0 * b3i; + gcos *= costh; + dgcos *= costh; + gcos += b2i; + dgcos += 2.0 * b2i; + gcos *= costh; + dgcos *= costh; + gcos += b1i; + dgcos += b1i; + gcos *= costh; + gcos += b0i; + + const double bg0i = bg0[typei]; + const double bg1i = bg1[typei]; + const double bg2i = bg2[typei]; + const double bg3i = bg3[typei]; + const double bg4i = bg4[typei]; + const double bg5i = bg5[typei]; + const double bg6i = bg6[typei]; + double gamma = bg6i * costh; + double dgamma = 6.0 * bg6i * costh; + gamma += bg5i; + dgamma += 5.0 * bg5i; + gamma *= costh; + dgamma *= costh; + gamma += bg4i; + dgamma += 4.0 * bg4i; + gamma *= costh; + dgamma *= costh; + gamma += bg3i; + dgamma += 3.0 * bg3i; + gamma *= costh; + dgamma *= costh; + gamma += bg2i; + dgamma += 2.0 * bg2i; + gamma *= costh; + dgamma *= costh; + gamma += bg1i; + dgamma += bg1i; + gamma *= costh; + gamma += bg0i; + + const double tmp = MathConst::MY_2PI * (costh - 0.5); + const double psi = 0.5 * (1 - cos(tmp)); + const double dpsi = MathConst::MY_PI * sin(tmp); + g = gcos + psi * (gamma - gcos); + dgdc = dgcos + dpsi * (gamma - gcos) + psi * (dgamma - dgcos); + } else { + dgdc = 0.0; + } + return g; + } + + /* ---------------------------------------------------------------------- + Pij calculation + ------------------------------------------------------------------------- */ + + inline double PijSpline(const double NM, const double NS, const int typei, double &dp) const + { + const double N = NM + NS; + + dp = -a0[typei] + a1[typei] * a2[typei] * exp(-a2[typei] * N); + return -a0[typei] * (N - 1) - a1[typei] * exp(-a2[typei] * N) + a3[typei]; + } + + void read_file(char *); + void allocate(); + + // ---------------------------------------------------------------------- + // S'(t) and S(t) cutoff functions + // added to header for inlining + // ---------------------------------------------------------------------- + + /* ---------------------------------------------------------------------- + cutoff function Sprime + return cutoff and dX = derivative + no side effects + ------------------------------------------------------------------------- */ + + inline double Sp(double Xij, double Xmin, double Xmax, double &dX) const + { + double cutoff; + + const double t = (Xij - Xmin) / (Xmax - Xmin); + if (t <= 0.0) { + cutoff = 1.0; + dX = 0.0; + } else if (t >= 1.0) { + cutoff = 0.0; + dX = 0.0; + } else { + cutoff = 0.5 * (1.0 + cos(t * MathConst::MY_PI)); + dX = (-0.5 * MathConst::MY_PI * sin(t * MathConst::MY_PI)) / (Xmax - Xmin); + } + return cutoff; + }; +}; +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/OPENMP/pair_rebomos_omp.cpp b/src/OPENMP/pair_rebomos_omp.cpp new file mode 100644 index 0000000000..5143fd0f63 --- /dev/null +++ b/src/OPENMP/pair_rebomos_omp.cpp @@ -0,0 +1,702 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + References: + + This code: + Stewart J A and Spearot D E (2013) Atomistic simulations of nanoindentation on the basal plane of crystalline molybdenum disulfide. Modelling Simul. Mater. Sci. Eng. 21. + + Based on: + Liang T, Phillpot S R and Sinnott S B (2009) Parameterization of a reactive many-body potential for Mo2S systems. Phys. Rev. B79 245110. + Liang T, Phillpot S R and Sinnott S B (2012) Erratum: Parameterization of a reactive many-body potential for Mo-S systems. (Phys. Rev. B79 245110 (2009)) Phys. Rev. B85 199903(E). + + LAMMPS file contributing authors: James Stewart, Khanh Dang and Douglas Spearot (University of Arkansas) +------------------------------------------------------------------------- */ + +// clang-format on + +#include "pair_rebomos_omp.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "math_special.h" +#include "memory.h" +#include "my_page.h" +#include "neigh_list.h" + +#include "suffix.h" + +#include + +#include "omp_compat.h" +#if defined(_OPENMP) +#include +#endif + +using namespace LAMMPS_NS; +using namespace MathConst; +using MathSpecial::cube; +using MathSpecial::powint; +using MathSpecial::square; + +static constexpr double TOL = 1.0e-9; + +/* ---------------------------------------------------------------------- */ + +PairREBOMoSOMP::PairREBOMoSOMP(LAMMPS *lmp) : PairREBOMoS(lmp), ThrOMP(lmp, THR_PAIR) +{ + suffix_flag |= Suffix::OMP; + respa_enable = 0; +} + +// clang-format off + +/* ---------------------------------------------------------------------- */ + +void PairREBOMoSOMP::compute(int eflag, int vflag) +{ + ev_init(eflag,vflag); + + REBO_neigh_thr(); + + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int inum = list->inum; + +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) +#endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); + + FREBO_thr(ifrom,ito,eflag,thr); + FLJ_thr(ifrom,ito,eflag,thr); + + thr->timer(Timer::PAIR); + reduce_thr(this, eflag, vflag, thr); + } // end of omp parallel region +} + +/* ---------------------------------------------------------------------- + create REBO neighbor list from main neighbor list + REBO neighbor list stores neighbors of ghost atoms +------------------------------------------------------------------------- */ + +void PairREBOMoSOMP::REBO_neigh_thr() +{ + const int nthreads = comm->nthreads; + + if (atom->nmax > maxlocal) { + maxlocal = atom->nmax; + memory->destroy(REBO_numneigh); + memory->sfree(REBO_firstneigh); + memory->destroy(nM); + memory->destroy(nS); + memory->create(REBO_numneigh,maxlocal,"REBOMoS:numneigh"); + REBO_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *), + "REBOMoS:firstneigh"); + memory->create(nM,maxlocal,"REBOMoS:nM"); + memory->create(nS,maxlocal,"REBOMoS:nS"); + } + +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE +#endif + { + int i,j,ii,jj,n,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq,dS; + int *ilist,*jlist,*numneigh,**firstneigh; + int *neighptr; + + double **x = atom->x; + int *type = atom->type; + + const int allnum = list->inum + list->gnum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + +#if defined(_OPENMP) + const int tid = omp_get_thread_num(); +#else + const int tid = 0; +#endif + + const int iidelta = 1 + allnum/nthreads; + const int iifrom = tid*iidelta; + const int iito = ((iifrom+iidelta)>allnum) ? allnum : (iifrom+iidelta); + + // store all REBO neighs of owned and ghost atoms + // scan full neighbor list of I + + // each thread has its own page allocator + MyPage &ipg = ipage[tid]; + ipg.reset(); + + for (ii = iifrom; ii < iito; ii++) { + i = ilist[ii]; + + n = 0; + neighptr = ipg.vget(); + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = map[type[i]]; + nM[i] = nS[i] = 0.0; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtype = map[type[j]]; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < rcmaxsq[itype][jtype]) { + neighptr[n++] = j; + if (jtype == 0) + nM[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS); + else + nS[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS); + } + } + + REBO_firstneigh[i] = neighptr; + REBO_numneigh[i] = n; + ipg.vgot(n); + if (ipg.status()) + error->one(FLERR,"REBO list overflow, boost neigh_modify one"); + } + } +} + +/* ---------------------------------------------------------------------- + REBO forces and energy +------------------------------------------------------------------------- */ + +void PairREBOMoSOMP::FREBO_thr(int ifrom, int ito, int eflag, ThrData * const thr) +{ + int i,j,k,ii,itype,jtype; + tagint itag, jtag; + double delx,dely,delz,evdwl,fpair,xtmp,ytmp,ztmp; + double rsq,rij,wij; + double Qij,Aij,alphaij,VR,pre,dVRdi,VA,bij,dVAdi,dVA; + double dwij,del[3]; + int *ilist,*REBO_neighs; + + evdwl = 0.0; + + const double * const * const x = atom->x; + double * const * const f = thr->get_f(); + const int * const type = atom->type; + const tagint * const tag = atom->tag; + const int nlocal = atom->nlocal; + + ilist = list->ilist; + + // two-body interactions from REBO neighbor list, skip half of them + + for (ii = ifrom; ii < ito; ii++) { + i = ilist[ii]; + itag = tag[i]; + itype = map[type[i]]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + REBO_neighs = REBO_firstneigh[i]; + + for (k = 0; k < REBO_numneigh[i]; k++) { + j = REBO_neighs[k]; + jtag = tag[j]; + + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + + jtype = map[type[j]]; + + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + rij = sqrt(rsq); + wij = Sp(rij,rcmin[itype][jtype],rcmax[itype][jtype],dwij); + if (wij <= TOL) continue; + + Qij = Q[itype][jtype]; + Aij = A[itype][jtype]; + alphaij = alpha[itype][jtype]; + + VR = wij*(1.0+(Qij/rij)) * Aij*exp(-alphaij*rij); + pre = wij*Aij * exp(-alphaij*rij); + dVRdi = pre * ((-alphaij)-(Qij/rsq)-(Qij*alphaij/rij)); + dVRdi += VR/wij * dwij; + + VA = dVA = 0.0; + VA = -wij * BIJc[itype][jtype] * exp(-Beta[itype][jtype]*rij); + + dVA = -Beta[itype][jtype] * VA; + dVA += VA/wij * dwij; + + del[0] = delx; + del[1] = dely; + del[2] = delz; + bij = bondorder_thr(i,j,del,rij,VA,thr); + dVAdi = bij*dVA; + + fpair = -(dVRdi+dVAdi) / rij; + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + + if (eflag) evdwl = VR + bij*VA; + if (evflag) ev_tally_thr(this,i,j,nlocal,/* newton_pair */1,evdwl,0.0,fpair,delx,dely,delz,thr); + } + } +} + +/* ---------------------------------------------------------------------- + compute LJ forces and energy +------------------------------------------------------------------------- */ + +void PairREBOMoSOMP::FLJ_thr(int ifrom, int ito, int eflag, ThrData * const thr) +{ + int i,j,ii,jj,jnum,itype,jtype; + tagint itag,jtag; + double evdwl,fpair,xtmp,ytmp,ztmp; + double rij,delij[3],rijsq; + double VLJ,dVLJ; + double vdw,dvdw; + double r2inv,r6inv; + int *ilist,*jlist,*numneigh,**firstneigh; + double c2,c3,dr,drp,r6; + + // I-J interaction from full neighbor list + // skip 1/2 of interactions since only consider each pair once + + evdwl = 0.0; + + const double * const * const x = atom->x; + double * const * const f = thr->get_f(); + const tagint * const tag = atom->tag; + const int * const type = atom->type; + const int nlocal = atom->nlocal; + + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = ifrom; ii < ito; ii++) { + i = ilist[ii]; + itag = tag[i]; + itype = map[type[i]]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtag = tag[j]; + + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + jtype = map[type[j]]; + + delij[0] = xtmp - x[j][0]; + delij[1] = ytmp - x[j][1]; + delij[2] = ztmp - x[j][2]; + rijsq = delij[0]*delij[0] + delij[1]*delij[1] + delij[2]*delij[2]; + rij = sqrt(rijsq); + + // compute LJ forces and energy + + // Outside Rmax + if (rij > rcLJmax[itype][jtype] || rij < rcLJmin[itype][jtype]){ + VLJ = 0; + dVLJ = 0; + } + + // Inside Rmax and above 0.95*sigma + else if (rij <= rcLJmax[itype][jtype] && rij >= 0.95*sigma[itype][jtype]){ + r2inv = 1.0/rijsq; + r6inv = r2inv*r2inv*r2inv; + VLJ = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + dVLJ = -r6inv*(lj1[itype][jtype]*r6inv - lj2[itype][jtype])/rij; + } + + // Below 0.95*sigma + else if (rij < 0.95*sigma[itype][jtype] && rij >= rcLJmin[itype][jtype]){ + dr = 0.95*sigma[itype][jtype] - rcLJmin[itype][jtype]; + r6 = powint((sigma[itype][jtype]/(0.95*sigma[itype][jtype])),6); + vdw = 4*epsilon[itype][jtype]*r6*(r6 - 1.0); + dvdw = (-4*epsilon[itype][jtype]/(0.95*sigma[itype][jtype]))*r6*(12.0*r6 - 6.0); + c2 = ((3.0/dr)*vdw - dvdw)/dr; + c3 = (vdw/(dr*dr) - c2)/dr; + + drp = rij - rcLJmin[itype][jtype]; + VLJ = drp*drp*(drp*c3 + c2); + dVLJ = drp*(3.0*drp*c3 + 2.0*c2); + } + + fpair = -dVLJ/rij; + f[i][0] += delij[0]*fpair; + f[i][1] += delij[1]*fpair; + f[i][2] += delij[2]*fpair; + f[j][0] -= delij[0]*fpair; + f[j][1] -= delij[1]*fpair; + f[j][2] -= delij[2]*fpair; + + if (eflag) evdwl = VLJ; + if (evflag) ev_tally_thr(this,i,j,nlocal,/*newton_pair*/1,evdwl,0.0,fpair,delij[0],delij[1],delij[2],thr); + + } + } +} + +/* ---------------------------------------------------------------------- + Bij function + + The bond order term modified the attractive portion of the REBO + potential based on the number of atoms around a specific pair + and the bond angle between sets of three atoms. + + The functions G(cos(theta)) and P(N) are evaluated and their + derivatives are also computed for use in the force calculation. +------------------------------------------------------------------------- */ + +double PairREBOMoSOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, double VA, ThrData *thr) +{ + int atomi,atomj,atomk,atoml; + int k,l; + int itype, jtype, ktype, ltype; + double rik[3], rjl[3], rji[3], rki[3],rlj[3], dwjl, bij; + double NijM,NijS,NjiM,NjiS,wik,dwik,wjl; + double rikmag,rjlmag,cosjik,cosijl,g,tmp2; + double Etmp,pij,tmp,dwij,dS; + double dgdc,pji; + double dcosjikdri[3],dcosijldri[3],dcosjikdrk[3]; + double dp; + double dcosjikdrj[3],dcosijldrj[3],dcosijldrl[3]; + double fi[3],fj[3],fk[3],fl[3]; + double PijS, PjiS; + int *REBO_neighs; + + const double * const * const x = atom->x; + double * const * const f = thr->get_f(); + const int * const type = atom->type; + + atomi = i; + atomj = j; + itype = map[type[i]]; + jtype = map[type[j]]; + Sp(rijmag,rcmin[itype][jtype],rcmax[itype][jtype],dwij); + NijM = nM[i]; + NijS = nS[i]; + NjiM = nM[j]; + NjiS = nS[j]; + bij = 0.0; + tmp = 0.0; + tmp2 = 0.0; + dgdc = 0.0; + Etmp = 0.0; + + REBO_neighs = REBO_firstneigh[i]; + for (k = 0; k < REBO_numneigh[i]; k++) { + atomk = REBO_neighs[k]; + if (atomk != atomj) { + ktype = map[type[atomk]]; + rik[0] = x[atomi][0]-x[atomk][0]; + rik[1] = x[atomi][1]-x[atomk][1]; + rik[2] = x[atomi][2]-x[atomk][2]; + rikmag = sqrt((rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2])); + wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dS); + cosjik = ((rij[0]*rik[0])+(rij[1]*rik[1])+(rij[2]*rik[2])) / (rijmag*rikmag); + cosjik = MIN(cosjik,1.0); + cosjik = MAX(cosjik,-1.0); + + // evaluate g and derivative dg + + g = gSpline(cosjik,itype,dgdc); + Etmp = Etmp+(wik*g); + } + } + + dp = 0.0; + PijS = PijSpline(NijM,NijS,itype,dp); + pij = 1.0/sqrt(1.0+Etmp+PijS); + tmp = -0.5*cube(pij); + + // derivative calculations + + REBO_neighs = REBO_firstneigh[i]; + for (k = 0; k < REBO_numneigh[i]; k++) { + atomk = REBO_neighs[k]; + if (atomk != atomj) { + ktype = map[type[atomk]]; + rik[0] = x[atomi][0]-x[atomk][0]; + rik[1] = x[atomi][1]-x[atomk][1]; + rik[2] = x[atomi][2]-x[atomk][2]; + rikmag = sqrt((rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2])); + wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dwik); + cosjik = (rij[0]*rik[0] + rij[1]*rik[1] + rij[2]*rik[2]) / (rijmag*rikmag); + cosjik = MIN(cosjik,1.0); + cosjik = MAX(cosjik,-1.0); + + dcosjikdri[0] = ((rij[0]+rik[0])/(rijmag*rikmag)) - + (cosjik*((rij[0]/(rijmag*rijmag))+(rik[0]/(rikmag*rikmag)))); + dcosjikdri[1] = ((rij[1]+rik[1])/(rijmag*rikmag)) - + (cosjik*((rij[1]/(rijmag*rijmag))+(rik[1]/(rikmag*rikmag)))); + dcosjikdri[2] = ((rij[2]+rik[2])/(rijmag*rikmag)) - + (cosjik*((rij[2]/(rijmag*rijmag))+(rik[2]/(rikmag*rikmag)))); + dcosjikdrk[0] = (-rij[0]/(rijmag*rikmag)) + + (cosjik*(rik[0]/(rikmag*rikmag))); + dcosjikdrk[1] = (-rij[1]/(rijmag*rikmag)) + + (cosjik*(rik[1]/(rikmag*rikmag))); + dcosjikdrk[2] = (-rij[2]/(rijmag*rikmag)) + + (cosjik*(rik[2]/(rikmag*rikmag))); + dcosjikdrj[0] = (-rik[0]/(rijmag*rikmag)) + + (cosjik*(rij[0]/(rijmag*rijmag))); + dcosjikdrj[1] = (-rik[1]/(rijmag*rikmag)) + + (cosjik*(rij[1]/(rijmag*rijmag))); + dcosjikdrj[2] = (-rik[2]/(rijmag*rikmag)) + + (cosjik*(rij[2]/(rijmag*rijmag))); + + g = gSpline(cosjik,itype,dgdc); + tmp2 = VA*0.5*(tmp*wik*dgdc); + fj[0] = -tmp2*dcosjikdrj[0]; + fj[1] = -tmp2*dcosjikdrj[1]; + fj[2] = -tmp2*dcosjikdrj[2]; + fi[0] = -tmp2*dcosjikdri[0]; + fi[1] = -tmp2*dcosjikdri[1]; + fi[2] = -tmp2*dcosjikdri[2]; + fk[0] = -tmp2*dcosjikdrk[0]; + fk[1] = -tmp2*dcosjikdrk[1]; + fk[2] = -tmp2*dcosjikdrk[2]; + + // coordination forces + + // dwik forces (from partial derivative) + + tmp2 = VA*0.5*(tmp*dwik*g)/rikmag; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; + + // PIJ forces (from coordination P(N) term) + + tmp2 = VA*0.5*(tmp*dp*dwik)/rikmag; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; + + // dgdN forces are removed + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2]; + + if (vflag_either) { + rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; + rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2]; + v_tally3_thr(this,atomi,atomj,atomk,fj,fk,rji,rki,thr); + } + } + } + + // PIJ force contribution additional term + tmp2 = -VA*0.5*(tmp*dp*dwij)/rijmag; + + f[atomi][0] += rij[0]*tmp2; + f[atomi][1] += rij[1]*tmp2; + f[atomi][2] += rij[2]*tmp2; + f[atomj][0] -= rij[0]*tmp2; + f[atomj][1] -= rij[1]*tmp2; + f[atomj][2] -= rij[2]*tmp2; + + if (vflag_either) v_tally2_thr(this,atomi,atomj,tmp2,rij,thr); + + tmp = 0.0; + tmp2 = 0.0; + Etmp = 0.0; + + REBO_neighs = REBO_firstneigh[j]; + for (l = 0; l < REBO_numneigh[j]; l++) { + atoml = REBO_neighs[l]; + if (atoml != atomi) { + ltype = map[type[atoml]]; + rjl[0] = x[atomj][0]-x[atoml][0]; + rjl[1] = x[atomj][1]-x[atoml][1]; + rjl[2] = x[atomj][2]-x[atoml][2]; + rjlmag = sqrt((rjl[0]*rjl[0])+(rjl[1]*rjl[1])+(rjl[2]*rjl[2])); + wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dS); + cosijl = -1.0*((rij[0]*rjl[0])+(rij[1]*rjl[1])+(rij[2]*rjl[2])) / (rijmag*rjlmag); + cosijl = MIN(cosijl,1.0); + cosijl = MAX(cosijl,-1.0); + + // evaluate g and derivative dg + + g = gSpline(cosijl,jtype,dgdc); + Etmp = Etmp+(wjl*g); + } + } + + dp = 0.0; + PjiS = PijSpline(NjiM,NjiS,jtype,dp); + pji = 1.0/sqrt(1.0+Etmp+PjiS); + tmp = -0.5*cube(pji); + + REBO_neighs = REBO_firstneigh[j]; + for (l = 0; l < REBO_numneigh[j]; l++) { + atoml = REBO_neighs[l]; + if (atoml != atomi) { + ltype = map[type[atoml]]; + rjl[0] = x[atomj][0]-x[atoml][0]; + rjl[1] = x[atomj][1]-x[atoml][1]; + rjl[2] = x[atomj][2]-x[atoml][2]; + rjlmag = sqrt((rjl[0]*rjl[0])+(rjl[1]*rjl[1])+(rjl[2]*rjl[2])); + wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dwjl); + cosijl = (-1.0*((rij[0]*rjl[0])+(rij[1]*rjl[1])+(rij[2]*rjl[2]))) / (rijmag*rjlmag); + cosijl = MIN(cosijl,1.0); + cosijl = MAX(cosijl,-1.0); + + dcosijldri[0] = (-rjl[0]/(rijmag*rjlmag)) - + (cosijl*rij[0]/(rijmag*rijmag)); + dcosijldri[1] = (-rjl[1]/(rijmag*rjlmag)) - + (cosijl*rij[1]/(rijmag*rijmag)); + dcosijldri[2] = (-rjl[2]/(rijmag*rjlmag)) - + (cosijl*rij[2]/(rijmag*rijmag)); + dcosijldrj[0] = ((-rij[0]+rjl[0])/(rijmag*rjlmag)) + + (cosijl*((rij[0]/square(rijmag))-(rjl[0]/(rjlmag*rjlmag)))); + dcosijldrj[1] = ((-rij[1]+rjl[1])/(rijmag*rjlmag)) + + (cosijl*((rij[1]/square(rijmag))-(rjl[1]/(rjlmag*rjlmag)))); + dcosijldrj[2] = ((-rij[2]+rjl[2])/(rijmag*rjlmag)) + + (cosijl*((rij[2]/square(rijmag))-(rjl[2]/(rjlmag*rjlmag)))); + dcosijldrl[0] = (rij[0]/(rijmag*rjlmag))+(cosijl*rjl[0]/(rjlmag*rjlmag)); + dcosijldrl[1] = (rij[1]/(rijmag*rjlmag))+(cosijl*rjl[1]/(rjlmag*rjlmag)); + dcosijldrl[2] = (rij[2]/(rijmag*rjlmag))+(cosijl*rjl[2]/(rjlmag*rjlmag)); + + // evaluate g and derivatives dg + + g = gSpline(cosijl,jtype,dgdc); + tmp2 = VA*0.5*(tmp*wjl*dgdc); + fi[0] = -tmp2*dcosijldri[0]; + fi[1] = -tmp2*dcosijldri[1]; + fi[2] = -tmp2*dcosijldri[2]; + fj[0] = -tmp2*dcosijldrj[0]; + fj[1] = -tmp2*dcosijldrj[1]; + fj[2] = -tmp2*dcosijldrj[2]; + fl[0] = -tmp2*dcosijldrl[0]; + fl[1] = -tmp2*dcosijldrl[1]; + fl[2] = -tmp2*dcosijldrl[2]; + + // coordination forces + + // dwik forces (from partial derivative) + + tmp2 = VA*0.5*(tmp*dwjl*g)/rjlmag; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; + + // PIJ forces (coordination) + + tmp2 = VA*0.5*(tmp*dp*dwjl)/rjlmag; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; + + // dgdN forces are removed + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2]; + + if (vflag_either) { + rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2]; + v_tally3_thr(this,atomi,atomj,atoml,fi,fl,rij,rlj,thr); + } + } + } + + // PIJ force contribution additional term + + tmp2 = -VA*0.5*(tmp*dp*dwij)/rijmag; + f[atomi][0] += rij[0]*tmp2; + f[atomi][1] += rij[1]*tmp2; + f[atomi][2] += rij[2]*tmp2; + f[atomj][0] -= rij[0]*tmp2; + f[atomj][1] -= rij[1]*tmp2; + f[atomj][2] -= rij[2]*tmp2; + + if (vflag_either) v_tally2_thr(this,atomi,atomj,tmp2,rij,thr); + + bij = (0.5*(pij+pji)); + return bij; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double PairREBOMoSOMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += PairREBOMoS::memory_usage(); + + return bytes; +} diff --git a/src/OPENMP/pair_rebomos_omp.h b/src/OPENMP/pair_rebomos_omp.h new file mode 100644 index 0000000000..ea87f51950 --- /dev/null +++ b/src/OPENMP/pair_rebomos_omp.h @@ -0,0 +1,46 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + 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(rebomos/omp,PairREBOMoSOMP); +// clang-format on +#else + +#ifndef LMP_PAIR_REBOMOS_OMP_H +#define LMP_PAIR_REBOMOS_OMP_H + +#include "pair_rebomos.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairREBOMoSOMP : public PairREBOMoS, public ThrOMP { + public: + PairREBOMoSOMP(class LAMMPS *); + + void compute(int, int) override; + double memory_usage() override; + + protected: + void FREBO_thr(int ifrom, int ito, int eflag, ThrData *const thr); + void FLJ_thr(int ifrom, int ito, int eflag, ThrData *const thr); + + void REBO_neigh_thr(); + + double bondorder_thr(int, int, double *, double, double, ThrData *const thr); +}; +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/unittest/force-styles/tests/manybody-pair-rebomos.yaml b/unittest/force-styles/tests/manybody-pair-rebomos.yaml new file mode 100644 index 0000000000..74fbe2b001 --- /dev/null +++ b/unittest/force-styles/tests/manybody-pair-rebomos.yaml @@ -0,0 +1,125 @@ +--- +lammps_version: 7 Feb 2024 +tags: slow +date_generated: Thu Feb 22 09:08:59 2024 +epsilon: 1e-12 +skip_tests: +prerequisites: ! | + pair rebomos +pre_commands: ! | + variable newton_pair delete + if "$(is_active(package,gpu)) > 0.0" then "variable newton_pair index off" else "variable newton_pair index on" +post_commands: ! "" +input_file: in.airebo +pair_style: rebomos +pair_coeff: ! | + * * MoS.rebomos Mo S +extract: ! "" +natoms: 48 +init_vdwl: 3158.017726833385 +init_coul: 0 +init_stress: ! |2- + 6.8398718310371441e+03 6.7325636075141883e+03 6.1154248388685965e+03 -3.2850057579078185e+02 -6.6397329123828470e+01 -3.4208234997867203e+02 +init_forces: ! |2 + 1 -3.7681565852737293e+00 2.4378308384489483e+02 -2.6969740060923279e+01 + 2 2.6266038011491645e+02 -2.0681295165426090e+02 -1.9964225279104706e+01 + 3 -2.5778010496370018e+02 -2.3678496612327552e+02 4.3829487525746238e+01 + 4 -1.0244959408607365e+01 -2.6794342328905179e+02 -2.6250177085935768e+01 + 5 2.2561151172717553e+02 2.0458598509633319e+02 -4.2183143539914880e+01 + 6 -2.4766208404861507e+02 1.8231884233051196e+02 1.2358374858464346e+01 + 7 -2.8974107341390795e+01 1.1031227139363692e+02 -1.2014739688866064e+02 + 8 2.1125868843851663e+02 -2.8863731327540540e+02 1.6869558243398984e+01 + 9 -3.1045083889038222e+02 -1.7569259168198960e+02 7.9653354234911163e+01 + 10 3.9886149586070935e+01 -1.0347554263069082e+02 1.5430519714983259e+02 + 11 1.8975960895859026e+02 1.4943835159807929e+02 -7.3224202319244824e+01 + 12 -2.1525843617159188e+02 2.2144411990369784e+02 -1.4893230739895195e+01 + 13 -4.5214798787473534e+00 1.1303681465405538e+02 -1.1385660017739841e+02 + 14 2.3203040375527121e+02 -3.4422786145586961e+01 1.8438743650405590e+02 + 15 -9.7170644212253990e+01 -2.9205504121924292e+02 1.9237141908302370e+01 + 16 -5.1044169101209107e+01 -1.4503490444304657e+02 1.5798442448724001e+02 + 17 2.8968061091312188e+02 8.7626477818666558e+01 -5.8802357863256212e+01 + 18 -1.7269858357220591e+02 8.9805161881631591e+01 -1.7847636055101773e+02 + 19 1.0760757939777642e+02 1.2642591396366257e+02 1.5914528700154095e+01 + 20 8.4608047558760049e+01 -2.4114295115066687e+02 -8.5792447633922421e+01 + 21 -8.8196263713344706e+01 -1.2718199182512282e+02 -7.0146571930858386e+01 + 22 -9.8957412653883708e+01 -1.4706747771082254e+02 -4.1887992372999022e+01 + 23 1.3739608451158423e+02 2.5576040689729510e+02 3.0309933932861746e+01 + 24 -1.7919184387909232e+02 1.4014151619269020e+02 1.6188399417577682e+02 + 25 -7.7347072523993347e+01 1.5555894832740958e+02 2.6182670200573220e+01 + 26 7.1360176064432650e+01 -2.9346871053231695e+02 -4.0766232856732509e+00 + 27 -2.6964124483003886e+01 9.8094028485618843e+00 2.4520760165936114e+01 + 28 1.2261226472559841e+02 -2.1507437048459283e+02 -3.0529457478077248e+01 + 29 2.2123474818204119e+01 2.8827028167050270e+02 1.0492685172067040e+02 + 30 -3.0573547162322302e+02 6.5880918489915032e+01 -1.1459287271946665e+02 + 31 6.0299759271951014e+00 1.7876996651274851e+02 -1.2964028093837922e+02 + 32 2.3752911820416773e+02 -1.2903240767402318e+02 1.2355068302179771e+02 + 33 -2.7462402312452838e+02 -5.3957296643601545e+01 9.5908146508348565e+01 + 34 -2.4514118579997412e+01 -1.0385194449627924e+02 1.4149428440602534e+02 + 35 1.2613298864651166e+02 6.2171771831118461e+00 -2.0642312645163636e+02 + 36 -6.9996802234927173e+01 2.9354514776980676e+02 -2.0963604805137543e+01 + 37 5.0241789394721707e+01 2.4495589528252060e+02 -1.1814425338523709e+01 + 38 2.2684897209808264e+02 -1.8673318434123004e+02 1.5220428086050066e+00 + 39 -2.5107271651763259e+02 -1.6443685417603410e+02 -3.8081572318769169e+01 + 40 -5.6754727284178912e+01 -2.5956532443118618e+02 -5.9343761520548126e+00 + 41 2.7426370608616236e+02 1.7742153146823480e+02 1.3451000229991498e+01 + 42 -2.2732796257204490e+02 2.5279459742055184e+02 -2.7020668313372529e+01 + 43 1.6217419052721431e+02 1.3946372126978235e+02 -2.0623397237537300e+00 + 44 3.0490806490775839e+01 2.1072542774667330e+01 8.2426559376766466e+00 + 45 -1.8440339212966003e+02 -2.5521581334939563e+02 2.6487508881411735e+01 + 46 -3.5020523296396959e+01 -2.5119105337038474e+02 8.0096253171487604e+00 + 47 3.2722641189388571e+02 9.7949745935413119e+01 -1.1898751818014821e+01 + 48 -1.3785292104885133e+02 3.2239007811982697e+02 2.4602884867061839e+01 +run_vdwl: 838.005031679853 +run_coul: 0 +run_stress: ! |2- + 2.0454406617512559e+03 2.1642975706136021e+03 3.2875013790709349e+03 2.4615001869678809e+02 9.3532964794023911e+01 -1.6795786485689854e+02 +run_forces: ! |2 + 1 -2.5386905249574866e+01 -9.9909927316940710e+01 1.9622022776949393e+00 + 2 -3.6801215363277215e+01 1.3354608905768052e+02 1.3651128931415545e+02 + 3 2.1845136723115344e+00 2.3527273580981478e+00 2.4005941692097457e+00 + 4 -3.7818914621624749e+01 3.2946537814811997e+01 4.3018847892043375e+01 + 5 9.6010433504280890e+00 4.3327442916171194e+01 -1.0513311611209312e+02 + 6 8.4308572792324540e+01 -2.3503352214725012e+01 5.0539884405443850e+01 + 7 -1.1497740491253731e+01 -6.4938790153209212e+01 3.8948855627091397e+01 + 8 -4.2897157525993697e+00 -3.3530261299383611e+01 9.5574489119024229e+00 + 9 4.8767710902858741e-01 -3.5176606358955858e+00 -1.4806227207057554e-01 + 10 -1.4766722110672411e+01 2.3818988804615220e+01 -5.3134339246512809e+01 + 11 -1.3935431402422819e+02 -1.5665571900752980e+02 -3.7627334252485582e+01 + 12 7.7074185304589484e+00 -1.2031768398608595e+01 1.8274050209872662e+01 + 13 2.3670515505343992e+01 -6.4282829210616583e+01 4.6336363412751325e+01 + 14 -2.5980738378663424e+00 8.4167979467127569e+00 -2.3849686761968631e+00 + 15 2.6909213334137871e+01 9.4711670949167956e+00 1.4728430099051328e+00 + 16 2.6072330920174934e+01 1.5861581881675242e+01 -5.1429144994723586e+01 + 17 -4.5656802478750279e+01 -1.0557685195759221e+02 -1.8486644164722456e+01 + 18 1.3345348571390588e+02 5.3791289185967734e+01 3.1330462313925089e+01 + 19 -3.4744179120880638e+01 -2.9727735009574776e+01 1.3758697899848009e+01 + 20 3.5805745710027459e+00 2.3635583506556830e+00 -9.8641913612550636e-01 + 21 -1.9112967021807076e+01 2.4093036642722907e+01 -1.8643238887150801e+01 + 22 -9.0201617792954476e+00 3.7160721932912736e+01 -5.0873237843214003e+00 + 23 -2.0559670924056718e+01 -1.3037172781831362e+01 -4.1555140102832148e+01 + 24 3.5016464856031341e+01 -1.2488249708112150e+02 -5.5132899737529639e+00 + 25 2.4856779395007460e+00 -6.7533673428345381e+01 -1.7153563726726016e+01 + 26 -8.6470517146734750e+01 -7.9057942497429794e+00 3.8434105758101332e+01 + 27 3.1062776626038996e+01 2.5568916700234454e+00 -2.3812255018358194e+01 + 28 -3.6128021960690276e+01 8.8535946500139033e+01 -3.1718315624050986e-01 + 29 -1.7473599694524488e+01 5.1605700760038529e+00 3.1477016370448561e+00 + 30 3.2722165297173231e+00 -5.9807188688706105e+00 3.0903475914489444e+01 + 31 -5.5852599877204412e-01 -3.7635607965114964e+01 5.5407547347738578e+01 + 32 -4.8419811021707950e+01 1.2083134296973620e+02 -1.4299477119499166e+01 + 33 2.5993167832212944e+01 5.4415694550334088e+01 -2.2817585641730194e+01 + 34 1.3632918027769463e+01 1.0165209000060044e+02 -6.0709294862836366e+01 + 35 -1.7896100047169217e+01 5.7265522876438499e+00 6.7274680861552483e+00 + 36 -1.8125137765726237e+01 2.7564787465936047e+00 -9.5278575701404478e+00 + 37 -3.5898678616111592e+00 -7.2362382910010311e+01 -9.3588358281326709e+00 + 38 -3.6726824963125949e+01 2.0008528283799198e+01 4.2387983046244827e+01 + 39 1.4756349734330427e+00 3.4387817449154888e+01 4.6004570693407864e+01 + 40 1.3004687263347805e+01 6.5948251357408196e+01 -9.3870897441198231e+00 + 41 -9.2360908880622592e+00 3.0113497063365648e+00 -1.1785409504259610e+01 + 42 -1.1795417170666214e+01 -1.1546754933247897e+01 -1.2188857610119777e+01 + 43 -6.5511968637617390e-01 -3.9133448932537512e+01 5.5788437875425423e-01 + 44 -4.2976958767216473e+01 1.8814792063050870e+01 -3.2742042602021932e+01 + 45 1.8589630046648591e+02 5.3367997582431663e+01 3.8324469380182308e+01 + 46 2.2463767362468872e+01 7.7485287617829670e+01 2.6113916756613811e+01 + 47 3.0703470673161828e+00 2.4966324680483410e+01 -5.0826014249556524e+00 + 48 7.6310071304830785e+01 -9.3082908173611301e+01 -1.1280958703044767e+02 +...