forked from lijiext/lammps
Merge pull request #987 from sergeylishchuk/master
Added Axilrod-Teller manybody potential
This commit is contained in:
Steve Plimpton
GitHub
commit
c61f9248f4
|
|
@ -33,6 +33,7 @@ OPT.
|
|||
"agni (o)"_pair_agni.html,
|
||||
"airebo (oi)"_pair_airebo.html,
|
||||
"airebo/morse (oi)"_pair_airebo.html,
|
||||
"atm"_pair_atm.html,
|
||||
"awpmd/cut"_pair_awpmd.html,
|
||||
"beck (go)"_pair_beck.html,
|
||||
"body/nparticle"_pair_body_nparticle.html,
|
||||
|
|
|
|||
Binary file not shown.
|
After Width: | Height: | Size: 5.3 KiB |
|
|
@ -0,0 +1,9 @@
|
|||
\documentclass[12pt]{article}
|
||||
|
||||
\begin{document}
|
||||
|
||||
\begin{equation}
|
||||
E=\nu\frac{1+3\cos\gamma_1\cos\gamma_2\cos\gamma_3}{r_{12}^3r_{23}^3r_{31}^3}
|
||||
\end{equation}
|
||||
|
||||
\end{document}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 31 KiB |
|
|
@ -524,6 +524,7 @@ pair_write.html
|
|||
pair_adp.html
|
||||
pair_agni.html
|
||||
pair_airebo.html
|
||||
pair_atm.html
|
||||
pair_awpmd.html
|
||||
pair_beck.html
|
||||
pair_body_nparticle.html
|
||||
|
|
|
|||
|
|
@ -0,0 +1,164 @@
|
|||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Commands_all.html)
|
||||
|
||||
:line
|
||||
|
||||
pair_style atm command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
pair_style atm cutoff cutoff_triple :pre
|
||||
|
||||
cutoff = cutoff for each pair in 3-body interaction (distance units)
|
||||
cutoff_triple = additional cutoff applied to product of 3 pairwise distances (distance units) :ul
|
||||
|
||||
[Examples:]
|
||||
|
||||
pair_style atm 4.5 2.5
|
||||
pair_coeff * * * 0.072 :pre
|
||||
|
||||
pair_style hybrid/overlay lj/cut 6.5 atm 4.5 2.5
|
||||
pair_coeff * * lj/cut 1.0 1.0
|
||||
pair_coeff 1 1 atm 1 0.064
|
||||
pair_coeff 1 1 atm 2 0.080
|
||||
pair_coeff 1 2 atm 2 0.100
|
||||
pair_coeff 2 2 atm 2 0.125 :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
The {atm} style computes a 3-body "Axilrod-Teller-Muto"_#Axilrod
|
||||
potential for the energy E of a system of atoms as
|
||||
|
||||
:c,image(Eqs/pair_atm.jpg)
|
||||
|
||||
where nu is the three-body interaction strength. The distances
|
||||
between pairs of atoms r12, r23, r31 and the angles gamma1, gamma2,
|
||||
gamma3 are as shown in this diagram:
|
||||
|
||||
:c,image(JPG/pair_atm_dia.jpg)
|
||||
|
||||
Note that for the interaction between a triplet of atoms I,J,K, there
|
||||
is no "central" atom. The interaction is symmetric with respect to
|
||||
permutation of the three atoms. Thus the nu value is
|
||||
the same for all those permutations of the atom types of I,J,K
|
||||
and needs to be specified only once, as discussed below.
|
||||
|
||||
The {atm} potential is typically used in combination with a two-body
|
||||
potential using the "pair_style hybrid/overlay"_pair_hybrid.html
|
||||
command as in the example above.
|
||||
|
||||
The potential for a triplet of atom is calculated only if all 3
|
||||
distances r12, r23, r31 between the 3 atoms satisfy rIJ < cutoff.
|
||||
In addition, the product of the 3 distances r12*r23*r31 <
|
||||
cutoff_triple^3 is required, which excludes from calculation the
|
||||
triplets with small contribution to the interaction.
|
||||
|
||||
The following coefficients must be defined for each pair of atoms
|
||||
types via the "pair_coeff"_pair_coeff.html command as in the examples
|
||||
above, or in the restart files read by the
|
||||
"read_restart"_read_restart.html commands:
|
||||
|
||||
K = atom type of the third atom (1 to Ntypes)
|
||||
nu = prefactor (energy/distance^9 units) :ul
|
||||
|
||||
K can be specified in one of two ways. An explicit numeric value can
|
||||
be used, as in the 2nd example above. J <= K is required. LAMMPS
|
||||
sets the coefficients for the other 5 symmetric interactions to the
|
||||
same values. E.g. if I = 1, J = 2, K = 3, then these 6 values are set
|
||||
to the specified nu: nu123, nu132, nu213, nu231, nu312, nu321. This
|
||||
enforces the symmetry discussed above.
|
||||
|
||||
A wildcard asterisk can be used for K to set the coefficients for
|
||||
multiple triplets of atom types. This takes the form "*" or "*n" or
|
||||
"n*" or "m*n". If N = the number of atom types, then an asterisk with
|
||||
no numeric values means all types from 1 to N. A leading asterisk
|
||||
means all types from 1 to n (inclusive). A trailing asterisk means
|
||||
all types from n to N (inclusive). A middle asterisk means all types
|
||||
from m to n (inclusive). Note that only type triplets with J <= K are
|
||||
considered; if asterisks imply type triplets where K < J, they are
|
||||
ignored.
|
||||
|
||||
Note that a pair_coeff command can override a previous setting for the
|
||||
same I,J,K triplet. For example, these commands set nu for all I,J.K
|
||||
triplets, then overwrite nu for just the I,J,K = 2,3,4 triplet:
|
||||
|
||||
pair_coeff * * * 0.25
|
||||
pair_coeff 2 3 4 0.1 :pre
|
||||
|
||||
Note that for a simulation with a single atom type, only a single
|
||||
entry is required, e.g.
|
||||
|
||||
pair_coeff 1 1 1 0.25 :pre
|
||||
|
||||
For a simulation with two atom types, four pair_coeff commands will
|
||||
specify all possible nu values:
|
||||
|
||||
pair_coeff 1 1 1 nu1
|
||||
pair_coeff 1 1 2 nu2
|
||||
pair_coeff 1 2 2 nu3
|
||||
pair_coeff 2 2 2 nu4 :pre
|
||||
|
||||
For a simulation with three atom types, ten pair_coeff commands will
|
||||
specify all possible nu values:
|
||||
|
||||
pair_coeff 1 1 1 nu1
|
||||
pair_coeff 1 1 2 nu2
|
||||
pair_coeff 1 1 3 nu3
|
||||
pair_coeff 1 2 2 nu4
|
||||
pair_coeff 1 2 3 nu5
|
||||
pair_coeff 1 3 3 nu6
|
||||
pair_coeff 2 2 2 nu7
|
||||
pair_coeff 2 2 3 nu8
|
||||
pair_coeff 2 3 3 nu9
|
||||
pair_coeff 3 3 3 nu10 :pre
|
||||
|
||||
By default the nu value for all triplets is set to 0.0. Thus it is
|
||||
not required to provide pair_coeff commands that enumerate triplet
|
||||
interactions for all K types. If some I,J,K combination is not
|
||||
speficied, then there will be no 3-body ATM interactions for that
|
||||
combination and all its permutations. However, as with all pair
|
||||
styles, it is required to specify a pair_coeff command for all I,J
|
||||
combinations, else an error will result.
|
||||
|
||||
:line
|
||||
|
||||
[Mixing, shift, table, tail correction, restart, rRESPA info]:
|
||||
|
||||
This pair styles do not support the "pair_modify"_pair_modify.html
|
||||
mix, shift, table, and tail options.
|
||||
|
||||
This pair style writes its information to "binary restart
|
||||
files"_restart.html, so pair_style and pair_coeff commands do not need
|
||||
to be specified in an input script that reads a restart file.
|
||||
However, if the {atm} potential is used in combination with other
|
||||
potentials using the "pair_style hybrid/overlay"_pair_hybrid.html
|
||||
command then pair_coeff commands need to be re-specified
|
||||
in the restart input script.
|
||||
|
||||
This pair style can only be used via the {pair} keyword of the
|
||||
"run_style respa"_run_style.html command. It does not support the
|
||||
{inner}, {middle}, {outer} keywords.
|
||||
|
||||
:line
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
This pair style is part of the MANYBODY package. It is only enabled
|
||||
if LAMMPS was built with that package. See the "Build
|
||||
package"_Build_package.html doc page for more info.
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"pair_coeff"_pair_coeff.html
|
||||
|
||||
[Default:] none
|
||||
|
||||
:line
|
||||
|
||||
:link(Axilrod)
|
||||
[(Axilrod)]
|
||||
Axilrod and Teller, J Chem Phys, 11, 299 (1943);
|
||||
Muto, Nippon Sugaku-Buturigakkwaishi 17, 629 (1943).
|
||||
|
|
@ -103,6 +103,7 @@ pair"_Commands_pair.html doc page are followed by one or more of
|
|||
"pair_style adp"_pair_adp.html - angular dependent potential (ADP) of Mishin
|
||||
"pair_style airebo"_pair_airebo.html - AIREBO potential of Stuart
|
||||
"pair_style airebo/morse"_pair_airebo.html - AIREBO with Morse instead of LJ
|
||||
"pair_style atm"_pair_atm.html - Axilrod-Teller-Muto potential
|
||||
"pair_style beck"_pair_beck.html - Beck potential
|
||||
"pair_style body/nparticle"_pair_body_nparticle.html - interactions between body particles
|
||||
"pair_style bop"_pair_bop.html - BOP potential of Pettifor
|
||||
|
|
|
|||
|
|
@ -8,6 +8,7 @@ Pair Styles :h1
|
|||
pair_adp
|
||||
pair_agni
|
||||
pair_airebo
|
||||
pair_atm
|
||||
pair_awpmd
|
||||
pair_beck
|
||||
pair_body_nparticle
|
||||
|
|
|
|||
|
|
@ -59,6 +59,7 @@ sub-directories:
|
|||
|
||||
accelerate: use of all the various accelerator packages
|
||||
airebo: polyethylene with AIREBO potential
|
||||
atm: Axilrod-Teller-Muto potential
|
||||
balance: dynamic load balancing, 2d system
|
||||
body: body particles, 2d system
|
||||
cmap: CMAP 5-body contributions to CHARMM force field
|
||||
|
|
|
|||
|
|
@ -0,0 +1,31 @@
|
|||
# Axilrod-Teller-Muto potential example
|
||||
|
||||
variable x index 1
|
||||
variable y index 1
|
||||
variable z index 1
|
||||
|
||||
variable xx equal 10*$x
|
||||
variable yy equal 10*$y
|
||||
variable zz equal 10*$z
|
||||
|
||||
units lj
|
||||
atom_style atomic
|
||||
|
||||
lattice fcc 0.65
|
||||
region box block 0 ${xx} 0 ${yy} 0 ${zz}
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
|
||||
pair_style hybrid/overlay lj/cut 4.5 atm 4.5 2.5
|
||||
pair_coeff * * lj/cut 1.0 1.0
|
||||
pair_coeff * * atm * 0.072
|
||||
|
||||
mass * 1.0
|
||||
velocity all create 1.033 12345678 loop geom
|
||||
|
||||
fix 1 all nvt temp 1.033 1.033 0.05
|
||||
|
||||
timestep 0.002
|
||||
thermo 5
|
||||
|
||||
run 25
|
||||
|
|
@ -0,0 +1,100 @@
|
|||
LAMMPS (22 Aug 2018)
|
||||
# Axilrod-Teller-Muto potential example
|
||||
|
||||
variable x index 1
|
||||
variable y index 1
|
||||
variable z index 1
|
||||
|
||||
variable xx equal 10*$x
|
||||
variable xx equal 10*1
|
||||
variable yy equal 10*$y
|
||||
variable yy equal 10*1
|
||||
variable zz equal 10*$z
|
||||
variable zz equal 10*1
|
||||
|
||||
units lj
|
||||
atom_style atomic
|
||||
|
||||
lattice fcc 0.65
|
||||
Lattice spacing in x,y,z = 1.83252 1.83252 1.83252
|
||||
region box block 0 ${xx} 0 ${yy} 0 ${zz}
|
||||
region box block 0 10 0 ${yy} 0 ${zz}
|
||||
region box block 0 10 0 10 0 ${zz}
|
||||
region box block 0 10 0 10 0 10
|
||||
create_box 1 box
|
||||
Created orthogonal box = (0 0 0) to (18.3252 18.3252 18.3252)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 4000 atoms
|
||||
Time spent = 0.00139618 secs
|
||||
|
||||
pair_style hybrid/overlay lj/cut 4.5 atm 4.5 2.5
|
||||
pair_coeff * * lj/cut 1.0 1.0
|
||||
pair_coeff * * atm * 0.072
|
||||
|
||||
mass * 1.0
|
||||
velocity all create 1.033 12345678 loop geom
|
||||
|
||||
fix 1 all nvt temp 1.033 1.033 0.05
|
||||
|
||||
timestep 0.002
|
||||
thermo 5
|
||||
|
||||
run 25
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 4.8
|
||||
ghost atom cutoff = 4.8
|
||||
binsize = 2.4, bins = 8 8 8
|
||||
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||
(1) pair lj/cut, perpetual, half/full from (2)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) pair atm, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 11.47 | 11.47 | 11.47 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 1.033 -4.8404387 0 -3.291326 -4.1332095
|
||||
5 1.0337247 -4.8402263 0 -3.290027 -4.1207962
|
||||
10 1.0355935 -4.8425889 0 -3.2895869 -4.0870158
|
||||
15 1.0376519 -4.84599 0 -3.2899013 -4.0278711
|
||||
20 1.0382257 -4.8478854 0 -3.2909361 -3.9368052
|
||||
25 1.0347886 -4.84473 0 -3.2929351 -3.8044469
|
||||
Loop time of 15.95 on 1 procs for 25 steps with 4000 atoms
|
||||
|
||||
Performance: 270.846 tau/day, 1.567 timesteps/s
|
||||
100.0% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 15.946 | 15.946 | 15.946 | 0.0 | 99.97
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.0015042 | 0.0015042 | 0.0015042 | 0.0 | 0.01
|
||||
Output | 0.00013781 | 0.00013781 | 0.00013781 | 0.0 | 0.00
|
||||
Modify | 0.0017776 | 0.0017776 | 0.0017776 | 0.0 | 0.01
|
||||
Other | | 0.0006771 | | | 0.00
|
||||
|
||||
Nlocal: 4000 ave 4000 max 4000 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 10895 ave 10895 max 10895 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 640000 ave 640000 max 640000 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 1.28e+06 ave 1.28e+06 max 1.28e+06 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 1280000
|
||||
Ave neighs/atom = 320
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
|
||||
Please see the log.cite file for references relevant to this simulation
|
||||
|
||||
Total wall time: 0:00:16
|
||||
|
|
@ -0,0 +1,100 @@
|
|||
LAMMPS (22 Aug 2018)
|
||||
# Axilrod-Teller-Muto potential example
|
||||
|
||||
variable x index 1
|
||||
variable y index 1
|
||||
variable z index 1
|
||||
|
||||
variable xx equal 10*$x
|
||||
variable xx equal 10*1
|
||||
variable yy equal 10*$y
|
||||
variable yy equal 10*1
|
||||
variable zz equal 10*$z
|
||||
variable zz equal 10*1
|
||||
|
||||
units lj
|
||||
atom_style atomic
|
||||
|
||||
lattice fcc 0.65
|
||||
Lattice spacing in x,y,z = 1.83252 1.83252 1.83252
|
||||
region box block 0 ${xx} 0 ${yy} 0 ${zz}
|
||||
region box block 0 10 0 ${yy} 0 ${zz}
|
||||
region box block 0 10 0 10 0 ${zz}
|
||||
region box block 0 10 0 10 0 10
|
||||
create_box 1 box
|
||||
Created orthogonal box = (0 0 0) to (18.3252 18.3252 18.3252)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 4000 atoms
|
||||
Time spent = 0.000900984 secs
|
||||
|
||||
pair_style hybrid/overlay lj/cut 4.5 atm 4.5 2.5
|
||||
pair_coeff * * lj/cut 1.0 1.0
|
||||
pair_coeff * * atm * 0.072
|
||||
|
||||
mass * 1.0
|
||||
velocity all create 1.033 12345678 loop geom
|
||||
|
||||
fix 1 all nvt temp 1.033 1.033 0.05
|
||||
|
||||
timestep 0.002
|
||||
thermo 5
|
||||
|
||||
run 25
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 4.8
|
||||
ghost atom cutoff = 4.8
|
||||
binsize = 2.4, bins = 8 8 8
|
||||
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||
(1) pair lj/cut, perpetual, half/full from (2)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) pair atm, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 5.532 | 5.532 | 5.532 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 1.033 -4.8404387 0 -3.291326 -4.1332095
|
||||
5 1.0337247 -4.8402263 0 -3.290027 -4.1207962
|
||||
10 1.0355935 -4.8425889 0 -3.2895869 -4.0870158
|
||||
15 1.0376519 -4.84599 0 -3.2899013 -4.0278711
|
||||
20 1.0382257 -4.8478854 0 -3.2909361 -3.9368052
|
||||
25 1.0347886 -4.84473 0 -3.2929351 -3.8044469
|
||||
Loop time of 4.34636 on 4 procs for 25 steps with 4000 atoms
|
||||
|
||||
Performance: 993.935 tau/day, 5.752 timesteps/s
|
||||
99.6% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 3.9977 | 4.1036 | 4.209 | 4.9 | 94.41
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.13588 | 0.24134 | 0.34722 | 20.4 | 5.55
|
||||
Output | 0.00013757 | 0.00015104 | 0.00016761 | 0.0 | 0.00
|
||||
Modify | 0.00087953 | 0.00091547 | 0.00095582 | 0.0 | 0.02
|
||||
Other | | 0.0003656 | | | 0.01
|
||||
|
||||
Nlocal: 1000 ave 1000 max 1000 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 5835 ave 5835 max 5835 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 160000 ave 160000 max 160000 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 320000 ave 320000 max 320000 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 1280000
|
||||
Ave neighs/atom = 320
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
|
||||
Please see the log.cite file for references relevant to this simulation
|
||||
|
||||
Total wall time: 0:00:04
|
||||
|
|
@ -0,0 +1,374 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Sergey Lishchuk
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <cmath>
|
||||
#include "pair_atm.h"
|
||||
#include "atom.h"
|
||||
#include "citeme.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
static const char cite_atm_package[] =
|
||||
"ATM package:\n\n"
|
||||
"@Article{Lishchuk:2012:164501,\n"
|
||||
" author = {S. V. Lishchuk},\n"
|
||||
" title = {Role of three-body interactions in formation of bulk viscosity in liquid argon},\n"
|
||||
" journal = {J.~Chem.~Phys.},\n"
|
||||
" year = 2012,\n"
|
||||
" volume = 136,\n"
|
||||
" pages = {164501}\n"
|
||||
"}\n\n";
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairATM::PairATM(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
if (lmp->citeme) lmp->citeme->add(cite_atm_package);
|
||||
|
||||
single_enable = 0;
|
||||
restartinfo = 1;
|
||||
one_coeff = 0;
|
||||
manybody_flag = 1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
check if allocated, since class can be destructed when incomplete
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
PairATM::~PairATM()
|
||||
{
|
||||
if (copymode) return;
|
||||
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
|
||||
memory->destroy(nu);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
workhorse routine that computes pairwise interactions
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairATM::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,k,ii,jj,kk,inum,jnum,jnumm1;
|
||||
double xi,yi,zi,evdwl;
|
||||
double rij2,rik2,rjk2;
|
||||
double rij[3],rik[3],rjk[3],fj[3],fk[3];
|
||||
double nu_local;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *type = atom->type;
|
||||
|
||||
double cutoff_squared = cut_global*cut_global;
|
||||
double triple = cut_triple*cut_triple*cut_triple;
|
||||
double cutoff_triple_sixth = triple*triple;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// triple loop over local atoms and neighbors twice
|
||||
// must compute each IJK triplet interaction exactly once
|
||||
// by proc that owns the triplet atom with smallest x coord
|
||||
// special logic to break ties if multiple atoms have same x or y coords
|
||||
// inner two loops for jj=1,Jnum and kk=jj+1,Jnum insure
|
||||
// the pair of other 2 non-minimum-x atoms is only considered once
|
||||
// triplet geometry criteria for calculation:
|
||||
// each pair distance <= cutoff
|
||||
// produce of 3 pair distances <= cutoff_triple^3
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xi = x[i][0];
|
||||
yi = x[i][1];
|
||||
zi = x[i][2];
|
||||
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
jnumm1 = jnum - 1;
|
||||
|
||||
for (jj = 0; jj < jnumm1; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
rij[0] = x[j][0] - xi;
|
||||
if (rij[0] < 0.0) continue;
|
||||
rij[1] = x[j][1] - yi;
|
||||
if (rij[0] == 0.0 and rij[1] < 0.0) continue;
|
||||
rij[2] = x[j][2] - zi;
|
||||
if (rij[0] == 0.0 and rij[1] == 0.0 and rij[2] < 0.0) continue;
|
||||
rij2 = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
if (rij2 > cutoff_squared) continue;
|
||||
|
||||
for (kk = jj+1; kk < jnum; kk++) {
|
||||
k = jlist[kk];
|
||||
k &= NEIGHMASK;
|
||||
|
||||
rik[0] = x[k][0] - xi;
|
||||
if (rik[0] < 0.0) continue;
|
||||
rik[1] = x[k][1] - yi;
|
||||
if (rik[0] == 0.0 and rik[1] < 0.0) continue;
|
||||
rik[2] = x[k][2] - zi;
|
||||
if (rik[0] == 0.0 and rik[1] == 0.0 and rik[2] < 0.0) continue;
|
||||
rik2 = rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2];
|
||||
if (rik2 > cutoff_squared) continue;
|
||||
|
||||
rjk[0] = x[k][0] - x[j][0];
|
||||
rjk[1] = x[k][1] - x[j][1];
|
||||
rjk[2] = x[k][2] - x[j][2];
|
||||
rjk2 = rjk[0]*rjk[0] + rjk[1]*rjk[1] + rjk[2]*rjk[2];
|
||||
if (rjk2 > cutoff_squared) continue;
|
||||
|
||||
double r6 = rij2*rjk2*rik2;
|
||||
if (r6 > cutoff_triple_sixth) continue;
|
||||
|
||||
nu_local = nu[type[i]][type[j]][type[k]];
|
||||
if (nu_local == 0.0) continue;
|
||||
|
||||
interaction_ddd(nu_local,
|
||||
r6,rij2,rik2,rjk2,rij,rik,rjk,fj,fk,eflag,evdwl);
|
||||
|
||||
f[i][0] -= fj[0] + fk[0];
|
||||
f[i][1] -= fj[1] + fk[1];
|
||||
f[i][2] -= fj[2] + fk[2];
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
f[j][2] += fj[2];
|
||||
f[k][0] += fk[0];
|
||||
f[k][1] += fk[1];
|
||||
f[k][2] += fk[2];
|
||||
|
||||
if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,rij,rik);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairATM::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(nu,n+1,n+1,n+1,"pair:nu");
|
||||
|
||||
// initialize all nu values to 0.0
|
||||
|
||||
for (int i = 1; i <= n; i++)
|
||||
for (int j = 1; j <= n; j++)
|
||||
for (int k = 1; k <= n; k++)
|
||||
nu[i][j][k] = 0.0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairATM::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 2) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
cut_global = force->numeric(FLERR,arg[0]);
|
||||
cut_triple = force->numeric(FLERR,arg[1]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coefficients for one I,J,K type triplet
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairATM::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 4) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi,jlo,jhi,klo,khi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
force->bounds(FLERR,arg[2],atom->ntypes,klo,khi);
|
||||
|
||||
double nu_one = force->numeric(FLERR,arg[3]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j<=jhi; j++) {
|
||||
for (int k = MAX(klo,j); k<=khi; k++) {
|
||||
nu[i][j][k] = nu_one;
|
||||
count++;
|
||||
}
|
||||
setflag[i][j] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairATM::init_style()
|
||||
{
|
||||
if (force->newton_pair == 0)
|
||||
error->all(FLERR,"Pair style ATM requires newton pair on");
|
||||
|
||||
// need a full neighbor list
|
||||
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one i,j type pair and corresponding j,i
|
||||
also for all k type permutations
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairATM::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
|
||||
// set all 6 symmetric permutations of I,J,K types to same nu value
|
||||
|
||||
int ntypes = atom->ntypes;
|
||||
for (int k = j; k <= ntypes; k++)
|
||||
nu[i][k][j] = nu[j][i][k] = nu[j][k][i] = nu[k][i][j] = nu[k][j][i] =
|
||||
nu[i][j][k];
|
||||
|
||||
return cut_global;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairATM::write_restart(FILE *fp)
|
||||
{
|
||||
write_restart_settings(fp);
|
||||
|
||||
int i,j,k;
|
||||
for (i = 1; i <= atom->ntypes; i++) {
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
fwrite(&setflag[i][j],sizeof(int),1,fp);
|
||||
if (setflag[i][j])
|
||||
for (k = j; k <= atom->ntypes; k++)
|
||||
fwrite(&nu[i][j][k],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairATM::read_restart(FILE *fp)
|
||||
{
|
||||
read_restart_settings(fp);
|
||||
allocate();
|
||||
|
||||
int i,j,k;
|
||||
int me = comm->me;
|
||||
for (i = 1; i <= atom->ntypes; i++) {
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
|
||||
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
|
||||
if (setflag[i][j]) for (k = j; k <= atom->ntypes; k++) {
|
||||
if (me == 0) fread(&nu[i][j][k],sizeof(double),1,fp);
|
||||
MPI_Bcast(&nu[i][j][k],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairATM::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_global,sizeof(double),1,fp);
|
||||
fwrite(&cut_triple,sizeof(double),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairATM::read_restart_settings(FILE *fp)
|
||||
{
|
||||
int me = comm->me;
|
||||
if (me == 0) {
|
||||
fread(&cut_global,sizeof(double),1,fp);
|
||||
fread(&cut_triple,sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_triple,1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Axilrod-Teller-Muto (dipole-dipole-dipole) potential
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairATM::interaction_ddd(double nu, double r6,
|
||||
double rij2, double rik2, double rjk2,
|
||||
double *rij, double *rik, double *rjk,
|
||||
double *fj, double *fk, int eflag, double &eng)
|
||||
{
|
||||
double r5inv,rri,rrj,rrk,rrr;
|
||||
r5inv = nu / (r6*r6*sqrt(r6));
|
||||
rri = rik[0]*rij[0] + rik[1]*rij[1] + rik[2]*rij[2];
|
||||
rrj = rij[0]*rjk[0] + rij[1]*rjk[1] + rij[2]*rjk[2];
|
||||
rrk = rjk[0]*rik[0] + rjk[1]*rik[1] + rjk[2]*rik[2];
|
||||
rrr = 5.0*rri*rrj*rrk;
|
||||
for (int i = 0; i < 3; i++) {
|
||||
fj[i] = rrj*(rrk - rri)*rik[i] -
|
||||
(rrk*rri - rjk2*rik2 + rrr/rij2) * rij[i] +
|
||||
(rrk*rri - rik2*rij2 + rrr/rjk2) * rjk[i];
|
||||
fj[i] *= 3.0*r5inv;
|
||||
fk[i] = rrk*(rri + rrj)*rij[i] +
|
||||
(rri*rrj + rik2*rij2 - rrr/rjk2) * rjk[i] +
|
||||
(rri*rrj + rij2*rjk2 - rrr/rik2) * rik[i];
|
||||
fk[i] *= 3.0*r5inv;
|
||||
}
|
||||
if (eflag) eng = (r6 - 0.6*rrr)*r5inv;
|
||||
}
|
||||
|
|
@ -0,0 +1,77 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(atm,PairATM)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_ATM_H
|
||||
#define LMP_PAIR_ATM_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairATM : public Pair {
|
||||
public:
|
||||
PairATM(class LAMMPS *);
|
||||
virtual ~PairATM();
|
||||
virtual void compute(int, int);
|
||||
void settings(int, char **);
|
||||
virtual void coeff(int, char **);
|
||||
virtual void init_style();
|
||||
virtual double init_one(int, int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
protected:
|
||||
double cut_global,cut_triple;
|
||||
double ***nu;
|
||||
|
||||
void allocate();
|
||||
void interaction_ddd(double, double, double, double, double, double *,
|
||||
double *, double *, double *, double *, int, double &);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal pair_style command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair style ATM requires newton pair on
|
||||
|
||||
See the newton command. This is a restriction to use the ATM
|
||||
potential.
|
||||
|
||||
E: All pair coeffs are not set
|
||||
|
||||
All pair coefficients must be set in the data file or by the
|
||||
pair_coeff command before running a simulation.
|
||||
|
||||
*/
|
||||
Loading…
Reference in New Issue