forked from lijiext/lammps
Merge pull request #2126 from vmohles/add-ECO-DF
Add fix orient/eco to USER-MISC
This commit is contained in:
Axel Kohlmeyer
GitHub
commit
c3831f8971
|
|
@ -147,6 +147,7 @@ OPT.
|
|||
* :doc:`oneway <fix_oneway>`
|
||||
* :doc:`orient/bcc <fix_orient>`
|
||||
* :doc:`orient/fcc <fix_orient>`
|
||||
* :doc:`orient/eco <fix_orient_eco>`
|
||||
* :doc:`phonon <fix_phonon>`
|
||||
* :doc:`pimd <fix_pimd>`
|
||||
* :doc:`planeforce <fix_planeforce>`
|
||||
|
|
|
|||
|
|
@ -1,37 +1,51 @@
|
|||
Submitting new features for inclusion in LAMMPS
|
||||
===============================================
|
||||
|
||||
We encourage users to submit new features or modifications for LAMMPS
|
||||
to `the core developers <https://lammps.sandia.gov/authors.html>`_ so they
|
||||
can be added to the LAMMPS distribution. The preferred way to manage
|
||||
and coordinate this is as of Fall 2016 via the LAMMPS project on
|
||||
`GitHub <https://github.com/lammps/lammps>`_. An alternative is to
|
||||
contact the LAMMPS developers or the indicated developer of a package
|
||||
or feature directly and send in your contribution via e-mail.
|
||||
We encourage users to submit new features or modifications for LAMMPS to
|
||||
`the core developers <https://lammps.sandia.gov/authors.html>`_ so they
|
||||
can be added to the LAMMPS distribution. The preferred way to manage and
|
||||
coordinate this is via the LAMMPS project on `GitHub
|
||||
<https://github.com/lammps/lammps>`_. Please see the :doc:`GitHub
|
||||
Tutorial <Howto_github>` for a demonstration on how to do that. An
|
||||
alternative is to contact the LAMMPS developers or the indicated
|
||||
developer of a package or feature directly and send in your contribution
|
||||
via e-mail, but that can add a significant delay on getting your
|
||||
contribution included, depending on how busy the developer is you
|
||||
contact, how complex a task it would be to integrate that code, and how
|
||||
many - if any - changes are required before the code can be included.
|
||||
|
||||
For any larger modifications or programming project, you are
|
||||
encouraged to contact the LAMMPS developers ahead of time, in order to
|
||||
discuss implementation strategies and coding guidelines, that will
|
||||
make it easier to integrate your contribution and result in less work
|
||||
for everybody involved. You are also encouraged to search through the
|
||||
list of `open issues on GitHub <https://github.com/lammps/lammps/issues>`_ and submit a new issue
|
||||
for a planned feature, so you would not duplicate the work of others
|
||||
(and possibly get scooped by them) or have your work duplicated by
|
||||
others.
|
||||
For any larger modifications or programming project, you are encouraged
|
||||
to contact the LAMMPS developers ahead of time, in order to discuss
|
||||
implementation strategies and coding guidelines, that will make it
|
||||
easier to integrate your contribution and result in less work for
|
||||
everybody involved. You are also encouraged to search through the list
|
||||
of `open issues on GitHub <https://github.com/lammps/lammps/issues>`_
|
||||
and submit a new issue for a planned feature, so you would not duplicate
|
||||
the work of others (and possibly get scooped by them) or have your work
|
||||
duplicated by others.
|
||||
|
||||
How quickly your contribution will be integrated depends largely on
|
||||
how much effort it will cause to integrate and test it, how much it
|
||||
requires changes to the core codebase, and of how much interest it is
|
||||
to the larger LAMMPS community. Please see below for a checklist of
|
||||
typical requirements. Once you have prepared everything, see the
|
||||
:doc:`Using GitHub with LAMMPS Howto <Howto_github>` doc page for instructions on how to
|
||||
submit your changes or new files through a GitHub pull request. If you
|
||||
prefer to submit patches or full files, you should first make certain,
|
||||
that your code works correctly with the latest patch-level version of
|
||||
LAMMPS and contains all bug fixes from it. Then create a gzipped tar
|
||||
file of all changed or added files or a corresponding patch file using
|
||||
'diff -u' or 'diff -c' and compress it with gzip. Please only use gzip
|
||||
compression, as this works well on all platforms.
|
||||
For informal communication with (some of) the LAMMPS developers you may
|
||||
ask to join the `LAMMPS developers on Slack <https://lammps.slack.com>`_.
|
||||
This slack work space is by invitation only. Thus for access, please
|
||||
send an e-mail to ``slack@lammps.org`` explaining what part of LAMMPS
|
||||
you are working on. Only discussions related to LAMMPS development are
|
||||
tolerated, so this is **NOT** for people that look for help with compiling,
|
||||
installing, or using LAMMPS. Please contact the `lammps-users mailing
|
||||
list <https://lammps.sandia.gov>`_ for those purposes instead.
|
||||
|
||||
How quickly your contribution will be integrated depends largely on how
|
||||
much effort it will cause to integrate and test it, how much it requires
|
||||
changes to the core codebase, and of how much interest it is to the
|
||||
larger LAMMPS community. Please see below for a checklist of typical
|
||||
requirements. Once you have prepared everything, see the :doc:`Using
|
||||
GitHub with LAMMPS Howto <Howto_github>` doc page for instructions on
|
||||
how to submit your changes or new files through a GitHub pull
|
||||
request. If you prefer to submit patches or full files, you should first
|
||||
make certain, that your code works correctly with the latest patch-level
|
||||
version of LAMMPS and contains all bug fixes from it. Then create a
|
||||
gzipped tar file of all changed or added files or a corresponding patch
|
||||
file using 'diff -u' or 'diff -c' and compress it with gzip. Please only
|
||||
use gzip compression, as this works well on all platforms.
|
||||
|
||||
If the new features/files are broadly useful we may add them as core
|
||||
files to LAMMPS or as part of a :doc:`standard package <Packages_standard>`. Else we will add them as a
|
||||
|
|
|
|||
|
|
@ -290,6 +290,7 @@ accelerated styles exist.
|
|||
* :doc:`oneway <fix_oneway>` - constrain particles on move in one direction
|
||||
* :doc:`orient/bcc <fix_orient>` - add grain boundary migration force for BCC
|
||||
* :doc:`orient/fcc <fix_orient>` - add grain boundary migration force for FCC
|
||||
* :doc:`orient/eco <fix_orient_eco>` - add generalized grain boundary migration force
|
||||
* :doc:`phonon <fix_phonon>` - calculate dynamical matrix from MD simulations
|
||||
* :doc:`pimd <fix_pimd>` - Feynman path integral molecular dynamics
|
||||
* :doc:`planeforce <fix_planeforce>` - constrain atoms to move in a plane
|
||||
|
|
|
|||
|
|
@ -0,0 +1,154 @@
|
|||
.. index:: fix orient/eco
|
||||
|
||||
fix orient/eco command
|
||||
======================
|
||||
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID orient/eco u0 eta cutoff orientationsFile
|
||||
|
||||
|
||||
* ID, group-ID are documented in fix command
|
||||
* u0 = energy added to each atom (energy units)
|
||||
* eta = cutoff value (usually 0.25)
|
||||
* cutoff = cutoff radius for orientation parameter calculation
|
||||
* orientationsFile = file that specifies orientation of each grain
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix gb all orient/eco 0.08 0.25 3.524 sigma5.ori
|
||||
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
The fix applies a synthetic driving force to a grain boundary which can
|
||||
be used for the investigation of grain boundary motion. The affiliation
|
||||
of atoms to either of the two grains forming the grain boundary is
|
||||
determined from an orientation-dependent order parameter as described
|
||||
in :ref:`(Ulomek) <Ulomek>`. The potential energy of atoms is either increased by an amount
|
||||
of 0.5*\ *u0* or -0.5*\ *u0* according to the orientation of the surrounding
|
||||
crystal. This creates a potential energy gradient which pushes atoms near
|
||||
the grain boundary to orient according to the energetically favorable
|
||||
grain orientation. This fix is designed for applications in bicrystal system
|
||||
with one grain boundary and open ends, or two opposite grain boundaries in
|
||||
a periodic system. In either case, the entire system can experience a
|
||||
displacement during the simulation which needs to be accounted for in the
|
||||
evaluation of the grain boundary velocity. While the basic method is
|
||||
described in :ref:`(Ulomek) <Ulomek>`, the implementation follows the efficient
|
||||
implementation from :ref:`(Schratt & Mohles) <Schratt>`. The synthetic potential energy added to an
|
||||
atom j is given by the following formulas
|
||||
|
||||
.. math::
|
||||
|
||||
w(|\vec{r}_{jk}|) = w_{jk} & = \left\{\begin{array}{lc} \frac{|\vec{r}_{jk}|^{4}}{r_{\mathrm{cut}}^{4}}
|
||||
-2\frac{|\vec{r}_{jk}|^{2}}{r_{\mathrm{cut}}^{2}}+1, & |\vec{r}_{jk}|<r_{\mathrm{cut}} \\
|
||||
0, & |\vec{r}_{jk}|\ge r_{\mathrm{cut}}
|
||||
\end{array}\right. \\
|
||||
\chi_{j} & = \frac{1}{N}\sum_{l=1}^{3}\left\lbrack\left\vert\psi_{l}^{\mathrm{I}}(\vec{r}_{j})\right\vert^{2}-\left\vert\psi_{l}^{\mathrm{II}}(\vec{r}_{j})\right\vert^{2}\right\rbrack \\
|
||||
\psi_{l}^{\mathrm{X}}(\vec{r}_{j}) & = \sum_{k\in\mathit{\Gamma}_{j}}w_{jk}\exp\left(\mathrm{i}\vec{r}_{jk}\cdot\vec{q}_{l}^{\mathrm{X}}\right) \\
|
||||
u(\chi_{j}) & = \frac{u_{0}}{2}\left\{\begin{array}{lc}
|
||||
1, & \chi_{j}\ge\eta\\
|
||||
\sin\left(\frac{\pi\chi_{j}}{2\eta}\right), & -\eta<\chi_{j}<\eta\\
|
||||
-1, & \chi_{j}\le-\eta
|
||||
\end{array}\right.
|
||||
|
||||
which are fully explained in :ref:`(Ulomek) <Ulomek>`
|
||||
and :ref:`(Schratt & Mohles) <Schratt>`.
|
||||
|
||||
The force on each atom is the negative gradient of the synthetic potential energy. It
|
||||
depends on the surrounding of this atom. An atom far from the grain boundary does not
|
||||
experience a synthetic force as its surrounding is that of an oriented single crystal
|
||||
and thermal fluctuations are masked by the parameter *eta*\ . Near the grain boundary
|
||||
however, the gradient is nonzero and synthetic force terms are computed.
|
||||
The orientationsFile specifies the perfect oriented crystal basis vectors for the
|
||||
two adjoining crystals. The first three lines (line=row vector) for the energetically penalized and the
|
||||
last three lines for the energetically favored grain assuming *u0* is positive. For
|
||||
negative *u0*, this is reversed. With the *cutoff* parameter, the size of the region around
|
||||
each atom which is used in the order parameter computation is defined. The cutoff must be
|
||||
smaller than the interaction range of the MD potential. It should at
|
||||
least include the nearest neighbor shell. For high temperatures or low angle
|
||||
grain boundaries, it might be beneficial to increase the cutoff in order to get a more
|
||||
precise identification of the atoms surrounding. However, computation time will
|
||||
increase as more atoms are considered in the order parameter and force computation.
|
||||
It is also worth noting that the cutoff radius must not exceed the communication
|
||||
distance for ghost atoms in LAMMPS. With orientationsFile, the
|
||||
6 oriented crystal basis vectors is specified. Each line of the input file
|
||||
contains the three components of a primitive lattice vector oriented according to
|
||||
the grain orientation in the simulation box. The first (last) three lines correspond
|
||||
to the primitive lattice vectors of the first (second) grain. An example for
|
||||
a :math:`\Sigma\langle001\rangle` mis-orientation is given at the end.
|
||||
|
||||
If no synthetic energy difference between the grains is created, :math:`u0=0`, the
|
||||
force computation is omitted. In this case, still, the order parameter of the
|
||||
driving force is computed and can be used to track the grain boundary motion throughout the
|
||||
simulation.
|
||||
|
||||
|
||||
|
||||
**Restart, fix_modify, output, run start/stop, minimize info:**
|
||||
|
||||
No information about this fix is written to :doc: `binary restart files <restart>`.
|
||||
|
||||
The :doc:`fix_modify <fix_modify>` *energy* option is supported by this fix to
|
||||
add the potential energy of atom interactions with the grain boundary
|
||||
driving force to the system's potential energy as part of thermodynamic output.
|
||||
The total sum of added synthetic potential energy is computed and can be accessed
|
||||
by various output options. The order parameter as well as the thermally masked
|
||||
output parameter are stored in per-atom arrays and can also be accessed by various
|
||||
:doc:`output commands <Howto_output>`.
|
||||
|
||||
No parameter of this fix can be used with the start/stop keywords of the run command. This fix is
|
||||
not invoked during energy minimization.
|
||||
|
||||
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This fix is part of the USER-MISC package. It is only enabled if LAMMPS was
|
||||
built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
|
||||
|
||||
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`fix_modify <fix_modify>`
|
||||
|
||||
:doc:`fix_orient <fix_orient>`
|
||||
|
||||
**Default:** none
|
||||
|
||||
----------
|
||||
|
||||
.. _Ulomek:
|
||||
|
||||
**(Ulomek)** Ulomek, Brien, Foiles, Mohles, Modelling Simul. Mater. Sci. Eng. 23 (2015) 025007
|
||||
|
||||
.. _Schratt:
|
||||
|
||||
**(Schratt & Mohles)** Schratt, Mohles. Comp. Mat. Sci. 182 (2020) 109774
|
||||
|
||||
----------
|
||||
|
||||
|
||||
For illustration purposes, here is an example file that specifies a
|
||||
:math:`\Sigma=5 \langle 001 \rangle` tilt grain boundary. This is for a lattice constant of 3.52 Angstrom:
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
sigma5.ori:
|
||||
|
||||
1.671685 0.557228 1.76212
|
||||
0.557228 -1.671685 1.76212
|
||||
2.228913 -1.114456 0.00000
|
||||
0.557228 1.671685 1.76212
|
||||
1.671685 -0.557228 1.76212
|
||||
2.228913 1.114456 0.00000
|
||||
|
||||
|
||||
|
|
@ -223,6 +223,7 @@ Bfrac
|
|||
bgq
|
||||
Bh
|
||||
Bialke
|
||||
bicrystal
|
||||
Biersack
|
||||
bigbig
|
||||
bigint
|
||||
|
|
@ -284,6 +285,7 @@ br
|
|||
Branduardi
|
||||
Branicio
|
||||
brennan
|
||||
Brien
|
||||
Brilliantov
|
||||
Broadwell
|
||||
Broglie
|
||||
|
|
@ -720,6 +722,7 @@ ebook
|
|||
ebt
|
||||
ec
|
||||
Ec
|
||||
eco
|
||||
ecoul
|
||||
ecp
|
||||
Ecut
|
||||
|
|
@ -1822,6 +1825,7 @@ Modine
|
|||
mofff
|
||||
MOFFF
|
||||
Mohd
|
||||
Mohles
|
||||
mol
|
||||
Mol
|
||||
molfile
|
||||
|
|
@ -2163,6 +2167,7 @@ optimizations
|
|||
orangered
|
||||
organometallic
|
||||
orientational
|
||||
orientationsFile
|
||||
orientorder
|
||||
Orlikowski
|
||||
ornl
|
||||
|
|
@ -2628,6 +2633,7 @@ Schimansky
|
|||
Schiotz
|
||||
Schlitter
|
||||
Schmid
|
||||
Schratt
|
||||
Schoen
|
||||
Schotte
|
||||
Schulten
|
||||
|
|
@ -3040,6 +3046,7 @@ ul
|
|||
ulb
|
||||
Uleft
|
||||
uloop
|
||||
Ulomek
|
||||
ulsph
|
||||
uMech
|
||||
umin
|
||||
|
|
|
|||
|
|
@ -0,0 +1 @@
|
|||
../../../../potentials/Ni_u3.eam
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,20 @@
|
|||
|
||||
units metal
|
||||
atom_style atomic
|
||||
read_data data.sigma5
|
||||
|
||||
pair_style eam
|
||||
pair_coeff * * Ni_u3.eam
|
||||
timestep 0.001
|
||||
|
||||
fix integrator all npt temp 750 750 0.1 iso 0 0 0.1
|
||||
fix eco all orient/eco 0.08 0.25 3.6 sigma5.ori
|
||||
|
||||
thermo 100
|
||||
thermo_style custom step temp etotal press vol f_eco
|
||||
velocity all create 750 18527782
|
||||
|
||||
#dump save all custom 100 orient_eco.dump xs ys zs f_eco[1] f_eco[2]
|
||||
#dump_modify save sort id
|
||||
|
||||
run 1000
|
||||
|
|
@ -0,0 +1,91 @@
|
|||
LAMMPS (2 Jun 2020)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
|
||||
units metal
|
||||
atom_style atomic
|
||||
read_data data.sigma5
|
||||
orthogonal box = (-90.27200752837744 -4.427188724235731e-07 0.0) to (90.27200752837744 88.54377448471462 17.5)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
25600 atoms
|
||||
read_data CPU = 0.020 secs
|
||||
|
||||
pair_style eam
|
||||
pair_coeff * * Ni_u3.eam
|
||||
Reading potential file Ni_u3.eam with DATE: 2007-06-11
|
||||
timestep 0.001
|
||||
|
||||
fix integrator all npt temp 750 750 0.1 iso 0 0 0.1
|
||||
fix eco all orient/eco 0.08 0.25 3.6 sigma5.ori
|
||||
|
||||
thermo 100
|
||||
thermo_style custom step temp etotal press vol f_eco
|
||||
velocity all create 750 18527782
|
||||
|
||||
#dump save all custom 100 orient_eco.dump xs ys zs f_eco[1] f_eco[2]
|
||||
#dump_modify save sort id
|
||||
|
||||
run 1000
|
||||
fix orient/eco: cutoff=3.6 norm_fac=30.012843706295556 neighbors=18
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 6.8
|
||||
ghost atom cutoff = 6.8
|
||||
binsize = 3.4, bins = 54 27 6
|
||||
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||
(1) pair eam, perpetual, half/full from (2)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) fix orient/eco, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 42.15 | 42.15 | 42.15 Mbytes
|
||||
Step Temp TotEng Press Volume f_eco
|
||||
0 750 -111027.95 15148.458 279755.85 8.7507522e-10
|
||||
100 422.2475 -110928.85 214.27204 285332.91 -1.0089186
|
||||
200 388.79514 -110721.02 -481.9574 286443.13 -0.78446905
|
||||
300 441.87285 -110478.56 -174.88288 286640.75 -1.1348357
|
||||
400 486.71094 -110211.6 65.075638 287023.56 -1.495318
|
||||
500 531.12815 -109923.33 97.309245 287552.73 -2.9498072
|
||||
600 569.82126 -109622.61 -85.73157 288229.35 -4.2855812
|
||||
700 615.84724 -109317.16 -52.508824 288799.86 -6.5214427
|
||||
800 666.09015 -109018.62 38.120383 289305.37 -8.4745641
|
||||
900 705.03939 -108738.66 263.39673 289867.87 -12.4514
|
||||
1000 735.59866 -108500.52 242.46405 290353.65 -15.427653
|
||||
Loop time of 109.966 on 1 procs for 1000 steps with 25600 atoms
|
||||
|
||||
Performance: 0.786 ns/day, 30.546 hours/ns, 9.094 timesteps/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 | 27.002 | 27.002 | 27.002 | 0.0 | 24.55
|
||||
Neigh | 0.83228 | 0.83228 | 0.83228 | 0.0 | 0.76
|
||||
Comm | 0.12623 | 0.12623 | 0.12623 | 0.0 | 0.11
|
||||
Output | 0.00080323 | 0.00080323 | 0.00080323 | 0.0 | 0.00
|
||||
Modify | 81.917 | 81.917 | 81.917 | 0.0 | 74.49
|
||||
Other | | 0.08839 | | | 0.08
|
||||
|
||||
Nlocal: 25600 ave 25600 max 25600 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 28259 ave 28259 max 28259 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 1.59326e+06 ave 1.59326e+06 max 1.59326e+06 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 3.18651e+06 ave 3.18651e+06 max 3.18651e+06 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 3186510
|
||||
Ave neighs/atom = 124.473
|
||||
Neighbor list builds = 10
|
||||
Dangerous builds = 0
|
||||
|
||||
Please see the log.cite file for references relevant to this simulation
|
||||
|
||||
Total wall time: 0:01:50
|
||||
|
|
@ -0,0 +1,91 @@
|
|||
LAMMPS (2 Jun 2020)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
|
||||
units metal
|
||||
atom_style atomic
|
||||
read_data data.sigma5
|
||||
orthogonal box = (-90.27200752837744 -4.427188724235731e-07 0.0) to (90.27200752837744 88.54377448471462 17.5)
|
||||
4 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
25600 atoms
|
||||
read_data CPU = 0.025 secs
|
||||
|
||||
pair_style eam
|
||||
pair_coeff * * Ni_u3.eam
|
||||
Reading potential file Ni_u3.eam with DATE: 2007-06-11
|
||||
timestep 0.001
|
||||
|
||||
fix integrator all npt temp 750 750 0.1 iso 0 0 0.1
|
||||
fix eco all orient/eco 0.08 0.25 3.6 sigma5.ori
|
||||
|
||||
thermo 100
|
||||
thermo_style custom step temp etotal press vol f_eco
|
||||
velocity all create 750 18527782
|
||||
|
||||
#dump save all custom 100 orient_eco.dump xs ys zs f_eco[1] f_eco[2]
|
||||
#dump_modify save sort id
|
||||
|
||||
run 1000
|
||||
fix orient/eco: cutoff=3.6 norm_fac=30.012843706295556 neighbors=18
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 6.8
|
||||
ghost atom cutoff = 6.8
|
||||
binsize = 3.4, bins = 54 27 6
|
||||
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||
(1) pair eam, perpetual, half/full from (2)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) fix orient/eco, 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.69 | 11.7 | 11.72 Mbytes
|
||||
Step Temp TotEng Press Volume f_eco
|
||||
0 750 -111027.95 15148.458 279755.85 8.749339e-10
|
||||
100 422.2475 -110928.85 214.27204 285332.91 -1.0089186
|
||||
200 388.79514 -110721.02 -481.9574 286443.13 -0.78446905
|
||||
300 441.87285 -110478.56 -174.88288 286640.75 -1.1348357
|
||||
400 486.71094 -110211.6 65.075638 287023.56 -1.495318
|
||||
500 531.12815 -109923.33 97.309245 287552.73 -2.9498072
|
||||
600 569.82126 -109622.61 -85.73157 288229.35 -4.2855812
|
||||
700 615.84724 -109317.16 -52.508824 288799.86 -6.5214427
|
||||
800 666.09015 -109018.62 38.120383 289305.37 -8.4745641
|
||||
900 705.03939 -108738.66 263.39673 289867.87 -12.4514
|
||||
1000 735.59866 -108500.52 242.46405 290353.65 -15.427653
|
||||
Loop time of 29.6634 on 4 procs for 1000 steps with 25600 atoms
|
||||
|
||||
Performance: 2.913 ns/day, 8.240 hours/ns, 33.712 timesteps/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 | 7.0794 | 7.1053 | 7.1325 | 0.7 | 23.95
|
||||
Neigh | 0.20947 | 0.21088 | 0.21231 | 0.2 | 0.71
|
||||
Comm | 0.15708 | 0.18471 | 0.22338 | 6.7 | 0.62
|
||||
Output | 0.00032616 | 0.00064683 | 0.0015936 | 0.0 | 0.00
|
||||
Modify | 22.105 | 22.118 | 22.139 | 0.3 | 74.56
|
||||
Other | | 0.0437 | | | 0.15
|
||||
|
||||
Nlocal: 6400 ave 6421 max 6384 min
|
||||
Histogram: 1 1 0 0 0 1 0 0 0 1
|
||||
Nghost: 9872.5 ave 9892 max 9855 min
|
||||
Histogram: 1 0 0 1 0 1 0 0 0 1
|
||||
Neighs: 398314 ave 400737 max 395743 min
|
||||
Histogram: 1 0 0 1 0 0 0 1 0 1
|
||||
FullNghs: 796628 ave 801194 max 792566 min
|
||||
Histogram: 1 0 1 0 0 0 1 0 0 1
|
||||
|
||||
Total # of neighbors = 3186510
|
||||
Ave neighs/atom = 124.473
|
||||
Neighbor list builds = 10
|
||||
Dangerous builds = 0
|
||||
|
||||
Please see the log.cite file for references relevant to this simulation
|
||||
|
||||
Total wall time: 0:00:29
|
||||
|
|
@ -0,0 +1,6 @@
|
|||
1.671685 0.557228 1.76212
|
||||
0.557228 -1.671685 1.76212
|
||||
2.228913 -1.114456 0.000000
|
||||
0.557228 1.671685 1.76212
|
||||
1.671685 -0.557228 1.76212
|
||||
2.228913 1.114456 0.000000
|
||||
|
|
@ -637,6 +637,8 @@
|
|||
/fix_oneway.h
|
||||
/fix_orient_bcc.cpp
|
||||
/fix_orient_bcc.h
|
||||
/fix_orient_eco.cpp
|
||||
/fix_orient_eco.h
|
||||
/fix_orient_fcc.cpp
|
||||
/fix_orient_fcc.h
|
||||
/fix_peri_neigh.cpp
|
||||
|
|
|
|||
|
|
@ -7,7 +7,7 @@ SHELL = /bin/sh
|
|||
# specify flags and libraries needed for your compiler
|
||||
|
||||
CC = mpicxx
|
||||
CCFLAGS = -g -O3
|
||||
CCFLAGS = -g -O3
|
||||
SHFLAGS = -fPIC
|
||||
DEPFLAGS = -M
|
||||
|
||||
|
|
|
|||
|
|
@ -60,6 +60,7 @@ fix imd, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 2009
|
|||
fix ipi, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014
|
||||
fix npt/cauchy, R. E. Miller (Carleton University), F. Pavia and S. Pattamatta, 12 Jan 2020
|
||||
fix nvk, Efrem Braun (UC Berkeley), efrem.braun at gmail.com, https://github.com/lammps/lammps/pull/310
|
||||
fix orient/eco Adrian A. Schratt and Volker Mohles (Ruhr-Uni Bochum), volker.mohles at rub.de, 6 Jun 2020
|
||||
fix pimd, Yuxing Peng (U Chicago), yuxing at uchicago.edu, 24 Nov 2014
|
||||
fix propel/self, Stefan Paquay (Brandeis U), stefanpaquay at gmail.com, 20 Jan 2020
|
||||
fix rhok, Ulf Pedersen (Roskilde U), ulf at urp.dk, 25 Sep 2017
|
||||
|
|
|
|||
|
|
@ -0,0 +1,586 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratdir_veces
|
||||
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 authors: Adrian A. Schratt and Volker Mohles (ICAMS)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "fix_orient_eco.h"
|
||||
#include <mpi.h>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include "atom.h"
|
||||
#include "citeme.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
#include "pair.h"
|
||||
#include "respa.h"
|
||||
#include "update.h"
|
||||
#include "utils.h"
|
||||
#include "fmt/format.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
using namespace MathConst;
|
||||
|
||||
static const char cite_fix_orient_eco[] =
|
||||
"fix orient/eco command:\n\n"
|
||||
"@Article{Schratt20,\n"
|
||||
" author = {A. A. Schratt, V. Mohles},\n"
|
||||
" title = {Efficient calculation of the ECO driving force for atomistic simulations of grain boundary motion},\n"
|
||||
" journal = {Computational Materials Science},\n"
|
||||
" volume = {182},\n"
|
||||
" year = {2020},\n"
|
||||
" pages = {109774},\n"
|
||||
" doi = {j.commatsci.2020.109774},\n"
|
||||
" url = {https://doi.org/10.1016/j.commatsci.2020.109774}\n"
|
||||
"}\n\n";
|
||||
|
||||
struct FixOrientECO::Nbr {
|
||||
double duchi; // potential derivative
|
||||
double real_phi[2][3]; // real part of wave function
|
||||
double imag_phi[2][3]; // imaginary part of wave function
|
||||
};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixOrientECO::FixOrientECO(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg),
|
||||
dir_filename(NULL), order(NULL), nbr(NULL), list(NULL)
|
||||
{
|
||||
if (lmp->citeme) lmp->citeme->add(cite_fix_orient_eco);
|
||||
|
||||
// get rank of this processor
|
||||
MPI_Comm_rank(world, &me);
|
||||
|
||||
// check for illegal command
|
||||
if (narg != 7) error->all(FLERR, "Illegal fix orient/eco command");
|
||||
|
||||
// set fix flags
|
||||
scalar_flag = 1; // computes scalar
|
||||
global_freq = 1; // values can be computed at every timestep
|
||||
extscalar = 1; // scalar scales with # of atoms
|
||||
peratom_flag = 1; // quantities are per atom quantities
|
||||
size_peratom_cols = 2; // # of per atom quantities
|
||||
peratom_freq = 1; //
|
||||
|
||||
// parse input parameters
|
||||
u_0 = force->numeric(FLERR, arg[3]);
|
||||
sign = (u_0 >= 0.0 ? 1 : -1);
|
||||
eta = force->numeric(FLERR, arg[4]);
|
||||
r_cut = force->numeric(FLERR, arg[5]);
|
||||
|
||||
// read reference orientations from file
|
||||
// work on rank 0 only
|
||||
int n = strlen(arg[6]) + 1;
|
||||
dir_filename = new char[n];
|
||||
strcpy(dir_filename, arg[6]);
|
||||
if (me == 0) {
|
||||
char line[IMGMAX];
|
||||
char *result;
|
||||
int count;
|
||||
|
||||
FILE *infile = force->open_potential(dir_filename);
|
||||
if (infile == NULL)
|
||||
error->one(FLERR,fmt::format("Cannot open fix orient/eco file {}: {}",
|
||||
dir_filename, utils::getsyserror()));
|
||||
for (int i = 0; i < 6; ++i) {
|
||||
result = fgets(line, IMGMAX, infile);
|
||||
if (!result) error->one(FLERR, "Fix orient/eco file read failed");
|
||||
count = sscanf(line, "%lg %lg %lg", &dir_vec[i][0], &dir_vec[i][1], &dir_vec[i][2]);
|
||||
if (count != 3) error->one(FLERR, "Fix orient/eco file read failed");
|
||||
}
|
||||
fclose(infile);
|
||||
|
||||
// calculate reciprocal lattice vectors
|
||||
get_reciprocal();
|
||||
|
||||
squared_cutoff = r_cut * r_cut;
|
||||
inv_squared_cutoff = 1.0 / squared_cutoff;
|
||||
half_u = 0.5 * u_0;
|
||||
inv_eta = 1.0 / eta;
|
||||
}
|
||||
|
||||
// initializations
|
||||
MPI_Bcast(&dir_vec[0][0], 18, MPI_DOUBLE, 0, world); // communicate direct lattice vectors
|
||||
MPI_Bcast(&reciprocal_vectors[0][0][0], 18, MPI_DOUBLE, 0, world); // communicate reciprocal vectors
|
||||
MPI_Bcast(&squared_cutoff, 1, MPI_DOUBLE, 0, world); // communicate squared cutoff radius
|
||||
MPI_Bcast(&inv_squared_cutoff, 1, MPI_DOUBLE, 0, world); // communicate inverse squared cutoff radius
|
||||
MPI_Bcast(&half_u, 1, MPI_DOUBLE, 0, world); // communicate half potential energy
|
||||
MPI_Bcast(&inv_eta, 1, MPI_DOUBLE, 0, world); // communicate inverse threshold
|
||||
|
||||
// set comm size needed by this Fix
|
||||
if (u_0 != 0) comm_forward = sizeof(Nbr) / sizeof(double);
|
||||
|
||||
added_energy = 0.0; // initial energy
|
||||
|
||||
nmax = atom->nmax;
|
||||
nbr = (Nbr *) memory->smalloc(nmax * sizeof(Nbr), "orient/eco:nbr");
|
||||
memory->create(order, nmax, 2, "orient/eco:order");
|
||||
array_atom = order;
|
||||
|
||||
// zero the array since a variable may access it before first run
|
||||
for (int i = 0; i < atom->nlocal; ++i) order[i][0] = order[i][1] = 0.0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixOrientECO::~FixOrientECO() {
|
||||
memory->destroy(order);
|
||||
memory->sfree(nbr);
|
||||
delete[] dir_filename;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixOrientECO::setmask() {
|
||||
int mask = 0;
|
||||
mask |= POST_FORCE;
|
||||
mask |= THERMO_ENERGY;
|
||||
mask |= POST_FORCE_RESPA;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixOrientECO::init() {
|
||||
// get this processors rank
|
||||
MPI_Comm_rank(world, &me);
|
||||
|
||||
// compute normalization factor
|
||||
int neigh = get_norm();
|
||||
if (me == 0) {
|
||||
utils::logmesg(lmp,fmt::format(" fix orient/eco: cutoff={} norm_fac={} "
|
||||
"neighbors={}\n", r_cut, norm_fac, neigh));
|
||||
}
|
||||
|
||||
inv_norm_fac = 1.0 / norm_fac;
|
||||
|
||||
// check parameters
|
||||
if (r_cut > force->pair->cutforce) {
|
||||
error->all(FLERR, "Cutoff radius used by fix orient/eco must be smaller than force cutoff");
|
||||
}
|
||||
|
||||
// communicate normalization factor
|
||||
MPI_Bcast(&norm_fac, 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&inv_norm_fac, 1, MPI_DOUBLE, 0, world);
|
||||
|
||||
if (strstr(update->integrate_style, "respa")) {
|
||||
ilevel_respa = ((Respa *) update->integrate)->nlevels - 1;
|
||||
if (respa_level >= 0) ilevel_respa = MIN(respa_level, ilevel_respa);
|
||||
}
|
||||
|
||||
// need a full neighbor list
|
||||
// perpetual list, built whenever re-neighboring occurs
|
||||
|
||||
int irequest = neighbor->request(this, instance_me);
|
||||
neighbor->requests[irequest]->pair = 0;
|
||||
neighbor->requests[irequest]->fix = 1;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixOrientECO::init_list(int /* id */, NeighList *ptr) {
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixOrientECO::setup(int vflag) {
|
||||
if (strstr(update->integrate_style, "verlet"))
|
||||
post_force(vflag);
|
||||
else {
|
||||
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
|
||||
post_force_respa(vflag,ilevel_respa, 0);
|
||||
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixOrientECO::post_force(int /* vflag */) {
|
||||
// set local variables
|
||||
int ii, i, jj, j; // loop variables and atom IDs
|
||||
int k; // variable to loop over 3 reciprocal directions
|
||||
int lambda; // variable to loop over 2 crystals
|
||||
int dim; // variable to loop over 3 spatial components
|
||||
double dx, dy, dz; // stores current interatomic vector
|
||||
double squared_distance; // stores current squared distance
|
||||
double chi; // stores current order parameter
|
||||
double weight; // stores current weight function
|
||||
double scalar_product; // stores current scalar product
|
||||
double omega; // phase of sine transition
|
||||
double omega_pre = MY_PI2 * inv_eta; // prefactor for omega
|
||||
double duchi_pre = half_u * MY_PI * inv_eta * inv_norm_fac; // prefactor for duchi
|
||||
double sin_om; // stores the value of the sine transition
|
||||
|
||||
// initialize added energy at this step
|
||||
added_energy = 0.0;
|
||||
|
||||
// set local pointers and neighbor lists
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *mask = atom->mask;
|
||||
int nall = atom->nlocal + atom->nghost;
|
||||
|
||||
int inum = list->inum;
|
||||
int jnum;
|
||||
int *ilist = list->ilist;
|
||||
int *jlist;
|
||||
int *numneigh = list->numneigh;
|
||||
int **firstneigh = list->firstneigh;
|
||||
|
||||
// insure nbr and order data structures are adequate size
|
||||
if (nall > nmax) {
|
||||
nmax = nall;
|
||||
memory->destroy(nbr);
|
||||
memory->destroy(order);
|
||||
nbr = (Nbr *) memory->smalloc(nmax * sizeof(Nbr), "orient/eco:nbr");
|
||||
memory->create(order, nmax, 2, "orient/eco:order");
|
||||
array_atom = order;
|
||||
}
|
||||
|
||||
// loop over owned atoms and compute order parameter
|
||||
for (ii = 0; ii < inum; ++ii) {
|
||||
i = ilist[ii];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
// initializations
|
||||
chi = 0.0;
|
||||
for (k = 0; k < 3; ++k) {
|
||||
nbr[i].real_phi[0][k] = nbr[i].real_phi[1][k] = 0.0;
|
||||
nbr[i].imag_phi[0][k] = nbr[i].imag_phi[1][k] = 0.0;
|
||||
}
|
||||
|
||||
// loop over all neighbors of atom i
|
||||
for (jj = 0; jj < jnum; ++jj) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
// vector difference
|
||||
dx = x[i][0] - x[j][0];
|
||||
dy = x[i][1] - x[j][1];
|
||||
dz = x[i][2] - x[j][2];
|
||||
squared_distance = dx * dx + dy * dy + dz * dz;
|
||||
|
||||
if (squared_distance < squared_cutoff) {
|
||||
squared_distance *= inv_squared_cutoff;
|
||||
weight = squared_distance * (squared_distance - 2.0) + 1.0;
|
||||
|
||||
for (lambda = 0; lambda < 2; ++lambda) {
|
||||
for (k = 0; k < 3; ++k) {
|
||||
scalar_product = reciprocal_vectors[lambda][k][0] * dx + reciprocal_vectors[lambda][k][1] * dy + reciprocal_vectors[lambda][k][2] * dz;
|
||||
nbr[i].real_phi[lambda][k] += weight * cos(scalar_product);
|
||||
nbr[i].imag_phi[lambda][k] += weight * sin(scalar_product);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// collect contributions
|
||||
for (k = 0; k < 3; ++k) {
|
||||
chi += (nbr[i].real_phi[0][k] * nbr[i].real_phi[0][k] + nbr[i].imag_phi[0][k] * nbr[i].imag_phi[0][k] -
|
||||
nbr[i].real_phi[1][k] * nbr[i].real_phi[1][k] - nbr[i].imag_phi[1][k] * nbr[i].imag_phi[1][k]);
|
||||
}
|
||||
chi *= inv_norm_fac;
|
||||
order[i][0] = chi;
|
||||
|
||||
// compute normalized order parameter
|
||||
// and potential energy
|
||||
if (chi > eta) {
|
||||
added_energy += half_u;
|
||||
nbr[i].duchi = 0.0;
|
||||
order[i][1] = sign;
|
||||
} else if (chi < -eta) {
|
||||
added_energy -= half_u;
|
||||
nbr[i].duchi = 0.0;
|
||||
order[i][1] = -sign;
|
||||
} else {
|
||||
omega = omega_pre * chi;
|
||||
sin_om = sin(omega);
|
||||
|
||||
added_energy += half_u * sin_om;
|
||||
nbr[i].duchi = duchi_pre * cos(omega);
|
||||
order[i][1] = sign * sin_om;
|
||||
}
|
||||
|
||||
// compute product with potential derivative
|
||||
for (k = 0; k < 3; ++k) {
|
||||
for (lambda = 0; lambda < 2; ++lambda) {
|
||||
nbr[i].real_phi[lambda][k] *= nbr[i].duchi;
|
||||
nbr[i].imag_phi[lambda][k] *= nbr[i].duchi;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// set additional local variables
|
||||
double gradient_ii_cos[2][3][3]; // gradient ii cosine term
|
||||
double gradient_ii_sin[2][3][3]; // gradient ii sine term
|
||||
double gradient_ij_vec[2][3][3]; // gradient ij vector term
|
||||
double gradient_ij_sca[2][3]; // gradient ij scalar term
|
||||
double weight_gradient_prefactor; // gradient prefactor
|
||||
double weight_gradient[3]; // gradient of weight
|
||||
double cos_scalar_product; // cosine of scalar product
|
||||
double sin_scalar_product; // sine of scalar product
|
||||
double gcos_scalar_product; // gradient weight function * cosine of scalar product
|
||||
double gsin_scalar_product; // gradient weight function * sine of scalar product
|
||||
|
||||
// compute force only if synthetic
|
||||
// potential is not zero
|
||||
if (u_0 != 0.0) {
|
||||
// communicate to acquire nbr data for ghost atoms
|
||||
comm->forward_comm_fix(this);
|
||||
|
||||
// loop over all atoms
|
||||
for (ii = 0; ii < inum; ++ii) {
|
||||
i = ilist[ii];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
const bool no_boundary_atom = (nbr[i].duchi == 0.0);
|
||||
|
||||
// skip atoms not in group
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
|
||||
// initializations
|
||||
for (k = 0; k < 3; ++k) {
|
||||
for (lambda = 0; lambda < 2; ++lambda) {
|
||||
for (dim = 0; dim < 3; ++dim) {
|
||||
gradient_ii_cos[lambda][k][dim] = 0.0;
|
||||
gradient_ii_sin[lambda][k][dim] = 0.0;
|
||||
gradient_ij_vec[lambda][k][dim] = 0.0;
|
||||
}
|
||||
gradient_ij_sca[lambda][k] = 0.0;
|
||||
}
|
||||
}
|
||||
// loop over all neighbors of atom i
|
||||
// for those within squared_cutoff, compute force
|
||||
for (jj = 0; jj < jnum; ++jj) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
// do not compute force on atom i if it is far from boundary
|
||||
if ((nbr[j].duchi == 0.0) && no_boundary_atom) continue;
|
||||
|
||||
// vector difference
|
||||
dx = x[i][0] - x[j][0];
|
||||
dy = x[i][1] - x[j][1];
|
||||
dz = x[i][2] - x[j][2];
|
||||
squared_distance = dx * dx + dy * dy + dz * dz;
|
||||
|
||||
if (squared_distance < squared_cutoff) {
|
||||
// compute force on atom i
|
||||
// need weight and its gradient
|
||||
squared_distance *= inv_squared_cutoff;
|
||||
weight = squared_distance * (squared_distance - 2.0) + 1.0;
|
||||
weight_gradient_prefactor = 4.0 * (squared_distance - 1.0) * inv_squared_cutoff;
|
||||
weight_gradient[0] = weight_gradient_prefactor * dx;
|
||||
weight_gradient[1] = weight_gradient_prefactor * dy;
|
||||
weight_gradient[2] = weight_gradient_prefactor * dz;
|
||||
|
||||
// (1) compute scalar product and sine and cosine of it
|
||||
// (2) compute product of sine and cosine with gradient of weight function
|
||||
// (3) compute gradient_ii_cos and gradient_ii_sin by summing up result of (2)
|
||||
// (4) compute gradient_ij_vec and gradient_ij_sca
|
||||
for (lambda = 0; lambda < 2; ++lambda) {
|
||||
for (k = 0; k < 3; ++k) {
|
||||
scalar_product = reciprocal_vectors[lambda][k][0] * dx + reciprocal_vectors[lambda][k][1] * dy + reciprocal_vectors[lambda][k][2] * dz;
|
||||
cos_scalar_product = cos(scalar_product);
|
||||
sin_scalar_product = sin(scalar_product);
|
||||
for (dim = 0; dim < 3; ++dim) {
|
||||
gradient_ii_cos[lambda][k][dim] += (gcos_scalar_product = weight_gradient[dim] * cos_scalar_product);
|
||||
gradient_ii_sin[lambda][k][dim] += (gsin_scalar_product = weight_gradient[dim] * sin_scalar_product);
|
||||
gradient_ij_vec[lambda][k][dim] += (nbr[j].real_phi[lambda][k] * gcos_scalar_product - nbr[j].imag_phi[lambda][k] * gsin_scalar_product);
|
||||
}
|
||||
gradient_ij_sca[lambda][k] += weight * (nbr[j].real_phi[lambda][k] * sin_scalar_product + nbr[j].imag_phi[lambda][k] * cos_scalar_product);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// sum contributions
|
||||
for (k = 0; k < 3; ++k) {
|
||||
for (dim = 0; dim < 3; ++dim) {
|
||||
f[i][dim] -= (nbr[i].real_phi[0][k] * gradient_ii_cos[0][k][dim] + nbr[i].imag_phi[0][k] * gradient_ii_sin[0][k][dim] + gradient_ij_vec[0][k][dim] + reciprocal_vectors[1][k][dim] * gradient_ij_sca[1][k]);
|
||||
f[i][dim] += (nbr[i].real_phi[1][k] * gradient_ii_cos[1][k][dim] + nbr[i].imag_phi[1][k] * gradient_ii_sin[1][k][dim] + gradient_ij_vec[1][k][dim] + reciprocal_vectors[0][k][dim] * gradient_ij_sca[0][k]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixOrientECO::post_force_respa(int vflag, int ilevel, int /* iloop */) {
|
||||
if (ilevel == ilevel_respa) post_force(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double FixOrientECO::compute_scalar() {
|
||||
double added_energy_total;
|
||||
MPI_Allreduce(&added_energy, &added_energy_total, 1, MPI_DOUBLE, MPI_SUM, world);
|
||||
return added_energy_total;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixOrientECO::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int */*pbc*/) {
|
||||
int i, j;
|
||||
int m = 0;
|
||||
|
||||
for (i = 0; i < n; ++i) {
|
||||
j = list[i];
|
||||
*((Nbr*)&buf[m]) = nbr[j];
|
||||
m += sizeof(Nbr)/sizeof(double);
|
||||
}
|
||||
|
||||
return m;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixOrientECO::unpack_forward_comm(int n, int first, double *buf) {
|
||||
int i;
|
||||
int last = first + n;
|
||||
int m = 0;
|
||||
|
||||
for (i = first; i < last; ++i) {
|
||||
nbr[i] = *((Nbr*) &buf[m]);
|
||||
m += sizeof(Nbr) / sizeof(double);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of local atom-based arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixOrientECO::memory_usage() {
|
||||
double bytes = nmax * sizeof(Nbr);
|
||||
bytes += 2 * nmax * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
reciprocal lattice vectors from file input
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixOrientECO::get_reciprocal() {
|
||||
// volume of unit cell A
|
||||
double vol = 0.5 * (dir_vec[0][0] * (dir_vec[1][1] * dir_vec[2][2] - dir_vec[2][1] * dir_vec[1][2]) + dir_vec[1][0] * (dir_vec[2][1] * dir_vec[0][2] - dir_vec[0][1] * dir_vec[2][2]) + dir_vec[2][0] * (dir_vec[0][1] * dir_vec[1][2] - dir_vec[1][1] * dir_vec[0][2])) / MY_PI;
|
||||
double i_vol = 1.0 / vol;
|
||||
|
||||
// grain A: reciprocal_vectors 0
|
||||
reciprocal_vectors[0][0][0] = (dir_vec[1][1] * dir_vec[2][2] - dir_vec[2][1] * dir_vec[1][2]) * i_vol;
|
||||
reciprocal_vectors[0][0][1] = (dir_vec[1][2] * dir_vec[2][0] - dir_vec[2][2] * dir_vec[1][0]) * i_vol;
|
||||
reciprocal_vectors[0][0][2] = (dir_vec[1][0] * dir_vec[2][1] - dir_vec[2][0] * dir_vec[1][1]) * i_vol;
|
||||
|
||||
// grain A: reciprocal_vectors 1
|
||||
reciprocal_vectors[0][1][0] = (dir_vec[2][1] * dir_vec[0][2] - dir_vec[0][1] * dir_vec[2][2]) * i_vol;
|
||||
reciprocal_vectors[0][1][1] = (dir_vec[2][2] * dir_vec[0][0] - dir_vec[0][2] * dir_vec[2][0]) * i_vol;
|
||||
reciprocal_vectors[0][1][2] = (dir_vec[2][0] * dir_vec[0][1] - dir_vec[0][0] * dir_vec[2][1]) * i_vol;
|
||||
|
||||
// grain A: reciprocal_vectors 2
|
||||
reciprocal_vectors[0][2][0] = (dir_vec[0][1] * dir_vec[1][2] - dir_vec[1][1] * dir_vec[0][2]) * i_vol;
|
||||
reciprocal_vectors[0][2][1] = (dir_vec[0][2] * dir_vec[1][0] - dir_vec[1][2] * dir_vec[0][0]) * i_vol;
|
||||
reciprocal_vectors[0][2][2] = (dir_vec[0][0] * dir_vec[1][1] - dir_vec[1][0] * dir_vec[0][1]) * i_vol;
|
||||
|
||||
// volume of unit cell B
|
||||
vol = 0.5 * (dir_vec[3][0] * (dir_vec[4][1] * dir_vec[5][2] - dir_vec[5][1] * dir_vec[4][2]) + dir_vec[4][0] * (dir_vec[5][1] * dir_vec[3][2] - dir_vec[3][1] * dir_vec[5][2]) + dir_vec[5][0] * (dir_vec[3][1] * dir_vec[4][2] - dir_vec[4][1] * dir_vec[3][2])) / MY_PI;
|
||||
i_vol = 1.0 / vol;
|
||||
|
||||
// grain B: reciprocal_vectors 0
|
||||
reciprocal_vectors[1][0][0] = (dir_vec[4][1] * dir_vec[5][2] - dir_vec[5][1] * dir_vec[4][2]) * i_vol;
|
||||
reciprocal_vectors[1][0][1] = (dir_vec[4][2] * dir_vec[5][0] - dir_vec[5][2] * dir_vec[4][0]) * i_vol;
|
||||
reciprocal_vectors[1][0][2] = (dir_vec[4][0] * dir_vec[5][1] - dir_vec[5][0] * dir_vec[4][1]) * i_vol;
|
||||
|
||||
// grain B: reciprocal_vectors 1
|
||||
reciprocal_vectors[1][1][0] = (dir_vec[5][1] * dir_vec[3][2] - dir_vec[3][1] * dir_vec[5][2]) * i_vol;
|
||||
reciprocal_vectors[1][1][1] = (dir_vec[5][2] * dir_vec[3][0] - dir_vec[3][2] * dir_vec[5][0]) * i_vol;
|
||||
reciprocal_vectors[1][1][2] = (dir_vec[5][0] * dir_vec[3][1] - dir_vec[3][0] * dir_vec[5][1]) * i_vol;
|
||||
|
||||
// grain B: reciprocal_vectors 2
|
||||
reciprocal_vectors[1][2][0] = (dir_vec[3][1] * dir_vec[4][2] - dir_vec[4][1] * dir_vec[3][2]) * i_vol;
|
||||
reciprocal_vectors[1][2][1] = (dir_vec[3][2] * dir_vec[4][0] - dir_vec[4][2] * dir_vec[3][0]) * i_vol;
|
||||
reciprocal_vectors[1][2][2] = (dir_vec[3][0] * dir_vec[4][1] - dir_vec[4][0] * dir_vec[3][1]) * i_vol;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
normalization factor
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int FixOrientECO::get_norm() {
|
||||
// set up local variables
|
||||
double delta[3]; // relative position
|
||||
double squared_distance; // squared distance of atoms i and j
|
||||
double weight; // weight function for atoms i and j
|
||||
double wsum = 0.0; // sum of all weight functions
|
||||
double scalar_product; // scalar product
|
||||
double reesum[3] = {0.0, 0.0, 0.0}; // sum of real part
|
||||
double imesum[3] = {0.0, 0.0, 0.0}; // sum of imaginary part
|
||||
|
||||
int max_co = 4; // will produce wrong results for rcut > 3 * lattice constant
|
||||
|
||||
int neigh = 0; // count number of neighbors used
|
||||
|
||||
// loop over ideal lattice positions
|
||||
int i, k, idx[3];
|
||||
for (idx[0] = -max_co; idx[0] <= max_co; ++idx[0]) {
|
||||
for (idx[1] = -max_co; idx[1] <= max_co; ++idx[1]) {
|
||||
for (idx[2] = -max_co; idx[2] <= max_co; ++idx[2]) {
|
||||
// distance of atoms
|
||||
for (i = 0; i < 3; ++i) {
|
||||
delta[i] = dir_vec[0][i] * idx[0] + dir_vec[1][i] * idx[1] + dir_vec[2][i] * idx[2];
|
||||
}
|
||||
squared_distance = delta[0] * delta[0] + delta[1] * delta[1] + delta[2] * delta[2];
|
||||
|
||||
// check if atom is within cutoff region
|
||||
if ((squared_distance != 0.0) and (squared_distance < squared_cutoff)) {
|
||||
++neigh;
|
||||
squared_distance *= inv_squared_cutoff;
|
||||
|
||||
// weight
|
||||
weight = squared_distance * (squared_distance - 2.0) + 1.0;
|
||||
wsum += weight;
|
||||
|
||||
// three reciprocal directions
|
||||
for (k = 0; k < 3; ++k) {
|
||||
scalar_product = reciprocal_vectors[1][k][0] * delta[0] + reciprocal_vectors[1][k][1] * delta[1] + reciprocal_vectors[1][k][2] * delta[2];
|
||||
reesum[k] += weight * cos(scalar_product);
|
||||
imesum[k] -= weight * sin(scalar_product);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// compute normalization
|
||||
norm_fac = 3.0 * wsum * wsum;
|
||||
for (k = 0; k < 3; ++k) {
|
||||
norm_fac -= reesum[k] * reesum[k] + imesum[k] * imesum[k];
|
||||
}
|
||||
return neigh;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,80 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 FIX_CLASS
|
||||
|
||||
FixStyle(orient/eco,FixOrientECO)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_ORIENT_ECO_H
|
||||
#define LMP_FIX_ORIENT_ECO_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixOrientECO : public Fix {
|
||||
public:
|
||||
FixOrientECO(class LAMMPS *, int, char **);
|
||||
~FixOrientECO();
|
||||
int setmask();
|
||||
void init();
|
||||
void init_list(int, class NeighList *);
|
||||
void setup(int);
|
||||
void post_force(int);
|
||||
void post_force_respa(int, int, int);
|
||||
double compute_scalar();
|
||||
int pack_forward_comm(int, int *, double *, int, int *);
|
||||
void unpack_forward_comm(int, int, double *);
|
||||
double memory_usage();
|
||||
|
||||
private:
|
||||
struct Nbr; // forward declaration. private struct for managing precomputed terms
|
||||
|
||||
int me; // this processors rank
|
||||
int nmax; // maximal # of owned + ghost atoms on this processor
|
||||
int ilevel_respa; // used for RESPA integrator only
|
||||
|
||||
int sign; // from sign of u
|
||||
double u_0; // synthetic potential energy
|
||||
double half_u; // half synthetic potential energy
|
||||
double eta; // threshold for thermal effects
|
||||
double inv_eta; // inverse threshold for thermal effects
|
||||
double r_cut; // cutoff radius
|
||||
double squared_cutoff; // squared cutoff radius
|
||||
double inv_squared_cutoff; // inverse squared cutoff radius
|
||||
char *dir_filename; // filename of reference grain input
|
||||
double dir_vec[6][3]; // direct lattice vectors
|
||||
double reciprocal_vectors[2][3][3]; // reciprocal lattice vectors
|
||||
|
||||
double added_energy; // energy added by fix
|
||||
|
||||
double **order; // order parameter and normalized order
|
||||
// parameter per atom
|
||||
|
||||
double norm_fac; // normalization constant
|
||||
double inv_norm_fac; // inverse normalization constant
|
||||
|
||||
Nbr *nbr; // pointer on array of precomputed terms
|
||||
class NeighList *list; // LAMMPS' neighbor list
|
||||
|
||||
void get_reciprocal(); // calculate reciprocal lattice vectors
|
||||
int get_norm(); // compute normalization factor
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // LMP_FIX_ORIENT_ECO_H
|
||||
#endif // FIX_CLASS
|
||||
|
||||
|
|
@ -150,7 +150,7 @@ TextFileReader * PotentialFileReader::open_potential(const std::string& path) {
|
|||
date = utils::get_potential_date(filepath, filetype);
|
||||
|
||||
if(!date.empty()) {
|
||||
utils::logmesg(lmp, fmt::format("Reading potential file {} with DATE: {}", filename, date));
|
||||
utils::logmesg(lmp, fmt::format("Reading potential file {} with DATE: {}\n", filename, date));
|
||||
}
|
||||
return new TextFileReader(filepath, filetype);
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue