Unverified Commit 947827fe authored by Axel Kohlmeyer's avatar Axel Kohlmeyer Committed by GitHub
Browse files

Merge pull request #1052 from jguenole/adaptglok

Optimized implementation of FIRE
parents 56ff70de b524db88
Loading
Loading
Loading
Loading
+1 −0
Original line number Diff line number Diff line
@@ -26,6 +26,7 @@ vgcore.*
.vagrant
\#*#
.#*
.vscode

.DS_Store
.DS_Store?
+90 −54
Original line number Diff line number Diff line
@@ -15,7 +15,7 @@ Syntax
  
  .. parsed-literal::
  
     keyword = *dmax* or *line* or *norm* or *alpha_damp* or *discrete_factor*
     keyword = *dmax* or *line* or *norm* or *alpha_damp* or *discrete_factor* or *integrator* or *tmax*
       *dmax* value = max
         max = maximum distance for line search to move (distance units)
       *line* value = *backtrack* or *quadratic* or *forcezero* or *spin_cubic* or *spin_none*
@@ -28,6 +28,10 @@ Syntax
         damping = fictitious Gilbert damping for spin minimization (adim)
       *discrete_factor* value = factor
         factor = discretization factor for adaptive spin timestep (adim)
       *integrator* value = *eulerimplicit* or *verlet*
         time integration scheme for fire minimization
       *tmax* value = factor
         factor = maximum adaptive timestep for fire minimization (adim)



@@ -38,6 +42,7 @@ Examples
.. parsed-literal::

   min_modify dmax 0.2
   min_modify integrator verlet tmax 4

Description
"""""""""""
@@ -60,32 +65,30 @@ highly overlapped atoms from being moved long distances (e.g. through
another atom) due to large forces.

The choice of line search algorithm for the *cg* and *sd* minimization
styles can be selected via the *line* keyword.
The default *quadratic* line search algorithm starts out using
the robust backtracking method described below. However, once
the system gets close to a local
minimum and the linesearch steps get small, so that the energy
is approximately quadratic in the step length, it uses the
estimated location of zero gradient as the linesearch step,
provided the energy change is downhill.
This becomes more efficient than backtracking
for highly-converged relaxations. The *forcezero*
line search algorithm is similar to *quadratic*\ .
It may be more efficient than *quadratic* on some systems.

The backtracking search is robust and should always find a local energy
minimum.  However, it will "converge" when it can no longer reduce the
energy of the system.  Individual atom forces may still be larger than
desired at this point, because the energy change is measured as the
difference of two large values (energy before and energy after) and
that difference may be smaller than machine epsilon even if atoms
could move in the gradient direction to reduce forces further.

The choice of a norm can be modified for the min styles *cg*\ , *sd*\ , 
*quickmin*\ , *fire*\ , *spin*\ , *spin/cg* and *spin/lbfgs* using 
the *norm* keyword.
The default *two* norm computes the 2-norm (Euclidean length) of the
global force vector:
styles can be selected via the *line* keyword.  The default
*quadratic* line search algorithm starts out using the robust
backtracking method described below. However, once the system gets
close to a local minimum and the linesearch steps get small, so that
the energy is approximately quadratic in the step length, it uses the
estimated location of zero gradient as the linesearch step, provided
the energy change is downhill.  This becomes more efficient than
backtracking for highly-converged relaxations. The *forcezero* line
search algorithm is similar to *quadratic*\ .  It may be more
efficient than *quadratic* on some systems.

The backtracking search is robust and should always find a local
energy minimum.  However, it will "converge" when it can no longer
reduce the energy of the system.  Individual atom forces may still be
larger than desired at this point, because the energy change is
measured as the difference of two large values (energy before and
energy after) and that difference may be smaller than machine epsilon
even if atoms could move in the gradient direction to reduce forces
further.

The choice of a norm can be modified for the min styles *cg*\ , *sd*\
, *quickmin*\ , *fire*\ , *fire/old*\ , *spin*\ , *spin/cg* and 
*spin/lbfgs* using the *norm* keyword.  The default *two* norm computes 
the 2-norm (Euclidean length) of the global force vector:

.. image:: Eqs/norm_two.jpg
   :align: center
@@ -106,31 +109,64 @@ all atoms in the system:
For the min styles *spin*\ , *spin/cg* and *spin/lbfgs*\ , the force
norm is replaced by the spin-torque norm.

Keywords *alpha\_damp* and *discrete\_factor* only make sense when a
:doc:`min_spin <min_spin>` command is declared.  Keyword *alpha\_damp*
defines an analog of a magnetic Gilbert damping. It defines a relaxation
rate toward an equilibrium for a given magnetic system.  Keyword
*discrete\_factor* defines a discretization factor for the adaptive
timestep used in the *spin* minimization.  See :doc:`min_spin
<min_spin>` for more information about those quantities.

The choice of a line search algorithm for the *spin/cg* and *spin/lbfgs*
styles can be specified via the *line* keyword.  The *spin\_cubic* and
*spin\_none* options only make sense when one of those two minimization
styles is declared.  The *spin\_cubic* option performs the line search
based on a cubic interpolation of the energy along the search
direction. The *spin\_none* option deactivates the line search
procedure.  The *spin\_none* option is a default value for *line*
keyword for both *spin/lbfgs* and *spin/cg*\ . Convergence of
*spin/lbfgs* can be more robust if *spin\_cubic* line search is used.
Keywords *alpha\_damp* and *discrete\_factor* only make sense when
a :doc:`min_spin <min_spin>` command is declared.
Keyword *alpha\_damp* defines an analog of a magnetic Gilbert
damping. It defines a relaxation rate toward an equilibrium for
a given magnetic system.
Keyword *discrete\_factor* defines a discretization factor for the
adaptive timestep used in the *spin* minimization.
See :doc:`min_spin <min_spin>` for more information about those
quantities.

The choice of a line search algorithm for the *spin/cg* and
*spin/lbfgs* styles can be specified via the *line* keyword.  The
*spin\_cubic* and *spin\_none* keywords only make sense when one of those two
minimization styles is declared.  The *spin\_cubic* performs the line
search based on a cubic interpolation of the energy along the search
direction. The *spin\_none* keyword deactivates the line search
procedure.  The *spin\_none* is a default value for *line* keyword for
both *spin/lbfgs* and *spin/cg*\ . Convergence of *spin/lbfgs* can be
more robust if *spin\_cubic* line search is used.

The Newton *integrator* used for *fire* minimization can be selected
to be either the symplectic Euler (\ *eulerimplicit*\ ) or velocity
Verlet (\ *verlet*\ ).  *tmax* defines the maximum value for the
adaptive timestep during a *fire* minimization. It is a multiplication
factor applied to the current :doc:`timestep <timestep>` (not in time
unit). For example, *tmax* = 4.0 with a :doc:`timestep <timestep>` of
2fs, means that the maximum value the timestep can reach during a *fire*
minimization is 4fs. 
Note that parameter defaults has been chosen to be reliable in most cases, 
but one should consider adjusting :doc:`timestep <timestep>` and *tmax* to
optimize the minimization for large or complex systems.  Other
parameters of the *fire* minimization can be tuned (\ *tmin*\ ,
*delaystep*\ , *dtgrow*\ , *dtshrink*\ , *alpha0*\ , and
*alphashrink*\ ). Please refer to the references describing the
:doc:`min_style <min_style>` *fire*.
An additional stopping criteria *vdfmax* is used by *fire* in order to avoid
unnecessary looping when it is reasonable to think the system will not
be relaxed further.  Note that in this case the system will NOT have
reached your minimization criteria. This could happen when the system 
comes to be stuck in a local basin of the phase space.  *vdfmax* is 
the maximum number of consecutive iterations with P(t) < 0.

The :doc:`min_style <min_style>` *fire* is an optimized implementation of
:doc:`min_style <min_style>` *fire/old*. It can however behave similarly 
to the *fire/old* style by using the following set of parameters:

.. parsed-literal::
  min_modify integrator eulerexplicit tmax 10.0 tmin 0.0 delaystep 5 &
             dtgrow 1.1 dtshrink 0.5 alpha0 0.1 alphashrink 0.99 &
             vdfmax 100000 halfstepback no initialdelay no

Restrictions
""""""""""""


For magnetic GNEB calculations, only the *spin\_none* value for *line* keyword can be used
when minimization styles *spin/cg* and *spin/lbfgs* are employed.
See :doc:`neb/spin <neb_spin>` for more explanations.
For magnetic GNEB calculations, only *spin\_none* value for *line*
keyword can be used when minimization styles *spin/cg* and *spin/lbfgs* are
employed.  See :doc:`neb/spin <neb_spin>` for more explanation.

Related commands
""""""""""""""""
@@ -142,11 +178,11 @@ Default

The option defaults are dmax = 0.1, line = quadratic and norm = two.

For the *spin*\ , *spin/cg* and *spin/lbfgs* styles, the
option defaults are alpha\_damp = 1.0, discrete\_factor = 10.0,
line = spin\_none, and norm = euclidean.

For the *spin*\ , *spin/cg* and *spin/lbfgs* styles, the option
defaults are alpha\_damp = 1.0, discrete\_factor = 10.0, line =
spin\_none, and norm = euclidean.

.. _lws: http://lammps.sandia.gov
.. _ld: Manual.html
.. _lc: Commands_all.html
For the *fire* style, the option defaults are integrator =
eulerimplicit, tmax = 10.0, tmin = 0.02, delaystep = 20, dtgrow = 1.1,
dtshrink = 0.5, alpha0 = 0.25, alphashrink = 0.99, vdfmax = 2000,
halfstepback = yes and initialdelay = yes.
+70 −50
Original line number Diff line number Diff line
@@ -11,7 +11,7 @@ Syntax

   min_style style

* style = *cg* or *hftn* or *sd* or *quickmin* or *fire* or *spin* or *spin/cg* or *spin/lbfgs*
* style = *cg* or *hftn* or *sd* or *quickmin* or *fire* or *fire/old* or *spin* or *spin/cg* or *spin/lbfgs*

Examples
""""""""
@@ -26,8 +26,8 @@ Examples
Description
"""""""""""

Choose a minimization algorithm to use when a :doc:`minimize <minimize>`
command is performed.
Choose a minimization algorithm to use when a :doc:`minimize
<minimize>` command is performed.

Style *cg* is the Polak-Ribiere version of the conjugate gradient (CG)
algorithm.  At each iteration the force gradient is combined with the
@@ -55,52 +55,68 @@ descent will not converge as quickly as CG, but may be more robust in
some situations.

Style *quickmin* is a damped dynamics method described in
:ref:`(Sheppard) <Sheppard>`, where the damping parameter is related to the
projection of the velocity vector along the current force vector for
each atom.  The velocity of each atom is initialized to 0.0 by this
style, at the beginning of a minimization.

Style *fire* is a damped dynamics method described in
:ref:`(Bitzek) <Bitzek>`, which is similar to *quickmin* but adds a variable
timestep and alters the projection operation to maintain components of
the velocity non-parallel to the current force vector.  The velocity
of each atom is initialized to 0.0 by this style, at the beginning of
a minimization.

Style *spin* is a damped spin dynamics with an adaptive
timestep.

Style *spin/cg* uses an orthogonal spin optimization (OSO)
combined to a conjugate gradient (CG) approach to minimize spin
configurations.

Style *spin/lbfgs* uses an orthogonal spin optimization (OSO)
combined to a limited-memory Broyden-Fletcher-Goldfarb-Shanno
(LBFGS) approach to minimize spin configurations.

See the :doc:`min/spin <min_spin>` doc page for more information
about the *spin*\ , *spin/cg* and *spin/lbfgs* styles.

Either the *quickmin* and *fire* styles are useful in the context of
nudged elastic band (NEB) calculations via the :doc:`neb <neb>` command.

Either the *spin*\ , *spin/cg* and *spin/lbfgs* styles are useful 
in the context of magnetic geodesic nudged elastic band (GNEB) calculations 
via the :doc:`neb/spin <neb_spin>` command.
:ref:`(Sheppard) <Sheppard>`, where the damping parameter is related
to the projection of the velocity vector along the current force
vector for each atom.  The velocity of each atom is initialized to 0.0
by this style, at the beginning of a minimization.

Style *fire* is a damped dynamics method described in :ref:`(Bitzek)
<Bitzek>`, which is similar to *quickmin* but adds a variable timestep
and alters the projection operation to maintain components of the
velocity non-parallel to the current force vector.  The velocity of
each atom is initialized to 0.0 by this style, at the beginning of a
minimization. This style correspond to an optimized version described 
in :ref:`(Guenole) <Guenole>` that include different time integration
schemes and defaults parameters. The default parameters can be 
modified with the command :doc:`min_modify <min_modify>`.


Style *fire/old* is the original implementation of *fire* in Lammps, 
conserved for backward compatibility. The main differences regarding 
the current version *fire* are: time integration by Explicit Euler
only, different sequence in maintaining velocity components non-parallel
to the current force vector and hard-coded minimization parameters.
A complete description of the differences between *fire/old* and *fire* 
can be found in :ref:`(Guenole) <Guenole>` (where the current *fire* 
in lammps is named *fire2.0*). By using an appropriate set of 
parameters, *fire* can behave similarly than *fire/old*, as described
in the :doc:`min_modify <min_modify>` command.

Style *spin* is a damped spin dynamics with an adaptive timestep.

Style *spin/cg* uses an orthogonal spin optimization (OSO) combined to
a conjugate gradient (CG) approach to minimize spin configurations.

Style *spin/lbfgs* uses an orthogonal spin optimization (OSO) combined
to a limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) approach
to minimize spin configurations.

See the :doc:`min/spin <min_spin>` doc page for more information about
the *spin*\ , *spin/cg* and *spin/lbfgs* styles.

Either the *quickmin*\ , *fire* and *fire/old* styles are useful in the
context of nudged elastic band (NEB) calculations via the :doc:`neb
<neb>` command.

Either the *spin*\ , *spin/cg* and *spin/lbfgs* styles are useful in
the context of magnetic geodesic nudged elastic band (GNEB)
calculations via the :doc:`neb/spin <neb_spin>` command.

.. note::

   The damped dynamic minimizers use whatever timestep you have
   defined via the :doc:`timestep <timestep>` command.  Often they will
   converge more quickly if you use a timestep about 10x larger than you
   would normally use for dynamics simulations.
   defined via the :doc:`timestep <timestep>` command.  Often they
   will converge more quickly if you use a timestep about 10x larger
   than you would normally use for dynamics simulations.
   For *fire*, the default timestep is recommended to be equal to
   the one you would normally use for dynamics simulations.

.. note::

   The *quickmin*\ , *fire*\ , *hftn*\ , and *cg/kk* styles do not yet
   support the use of the :doc:`fix box/relax <fix_box_relax>` command or
   minimizations involving the electron radius in :doc:`eFF <pair_eff>`
   models.
   The *quickmin*\ , *fire*\ , *fire/old*\ , *hftn*\ , and *cg/kk* styles do not yet
   support the use of the :doc:`fix box/relax <fix_box_relax>` command
   or minimizations involving the electron radius in :doc:`eFF
   <pair_eff>` models.


----------
@@ -109,17 +125,19 @@ via the :doc:`neb/spin <neb_spin>` command.
Styles with a *gpu*\ , *intel*\ , *kk*\ , *omp*\ , or *opt* suffix are
functionally the same as the corresponding style without the suffix.
They have been optimized to run faster, depending on your available
hardware, as discussed on the :doc:`Speed packages <Speed_packages>` doc
page.  The accelerated styles take the same arguments and should
hardware, as discussed on the :doc:`Speed packages <Speed_packages>`
doc page.  The accelerated styles take the same arguments and should
produce the same results, except for round-off and precision issues.

These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
USER-OMP and OPT packages, respectively.  They are only enabled if
LAMMPS was built with those packages.  See the :doc:`Build package <Build_package>` doc page for more info.
LAMMPS was built with those packages.  See the :doc:`Build package
<Build_package>` doc page for more info.

You can specify the accelerated styles explicitly in your input script
by including their suffix, or you can use the :doc:`-suffix command-line switch <Run_options>` when you invoke LAMMPS, or you can use the
:doc:`suffix <suffix>` command in your input script.
by including their suffix, or you can use the :doc:`-suffix
command-line switch <Run_options>` when you invoke LAMMPS, or you can
use the :doc:`suffix <suffix>` command in your input script.

See the :doc:`Speed packages <Speed_packages>` doc page for more
instructions on how to use the accelerated styles effectively.
@@ -164,7 +182,9 @@ Jonsson, Mills, Jacobsen.
**(Bitzek)** Bitzek, Koskinen, Gahler, Moseler, Gumbsch, Phys Rev Lett,
97, 170201 (2006).

.. _Guenole:



.. _lws: http://lammps.sandia.gov
.. _ld: Manual.html
.. _lc: Commands_all.html
**(Guenole)** Guenole, Noehring, Vaid, Houlle, Xie, Prakash, Bitzek,
Comput Mater Sci, (2020), in press (arXiv:190802038).
+26 −25
Original line number Diff line number Diff line
@@ -52,10 +52,12 @@ or :doc:`limited timestep <fix_nve_limit>`. Or dynamics can be run
using :doc:`fix viscous <fix_viscous>` to impose a damping force that
slowly drains all kinetic energy from the system.  The :doc:`pair_style soft <pair_soft>` potential can be used to un-overlap atoms while
running dynamics.
un-overlap atoms while running dynamics.

Note that you can minimize some atoms in the system while holding the
coordinates of other atoms fixed by applying :doc:`fix setforce <fix_setforce>` to the other atoms.  See a fuller
discussion of using fixes while minimizing below.
coordinates of other atoms fixed by applying :doc:`fix setforce
<fix_setforce>` to the other atoms.  See a fuller discussion of using
fixes while minimizing below.

The :doc:`minimization styles <min_style>` *cg*\ , *sd*\ , and *hftn*
involves an outer iteration loop which sets the search direction along
@@ -68,16 +70,16 @@ be more robust than previous line searches we've tried. The
backtracking method is described in Nocedal and Wright's Numerical
Optimization (Procedure 3.1 on p 41).

The :doc:`minimization styles <min_style>` *quickmin* and *fire* perform
damped dynamics using an Euler integration step.  Thus they require a
:doc:`timestep <timestep>` be defined.
The :doc:`minimization styles <min_style>` *quickmin*\ , *fire* and
*fire/old* perform damped dynamics using an Euler integration step.  Thus
they require a :doc:`timestep <timestep>` be defined.

.. note::

   The damped dynamic minimizers use whatever timestep you have
   defined via the :doc:`timestep <timestep>` command.  Often they will
   converge more quickly if you use a timestep about 10x larger than you
   would normally use for dynamics simulations.
   defined via the :doc:`timestep <timestep>` command.  Often they
   will converge more quickly if you use a timestep about 10x larger
   than you would normally use for dynamics simulations.


----------
@@ -90,13 +92,15 @@ coordinates:
.. image:: Eqs/min_energy.jpg
   :align: center

where the first term is the sum of all non-bonded :doc:`pairwise interactions <pair_style>` including :doc:`long-range Coulombic interactions <kspace_style>`, the 2nd through 5th terms are
:doc:`bond <bond_style>`, :doc:`angle <angle_style>`,
:doc:`dihedral <dihedral_style>`, and :doc:`improper <improper_style>`
interactions respectively, and the last term is energy due to
:doc:`fixes <fix>` which can act as constraints or apply force to atoms,
such as through interaction with a wall.  See the discussion below about
how fix commands affect minimization.
where the first term is the sum of all non-bonded :doc:`pairwise
interactions <pair_style>` including :doc:`long-range Coulombic
interactions <kspace_style>`, the 2nd through 5th terms are :doc:`bond
<bond_style>`, :doc:`angle <angle_style>`, :doc:`dihedral
<dihedral_style>`, and :doc:`improper <improper_style>` interactions
respectively, and the last term is energy due to :doc:`fixes <fix>`
which can act as constraints or apply force to atoms, such as through
interaction with a wall.  See the discussion below about how fix
commands affect minimization.

The starting point for the minimization is the current configuration
of the atoms.
@@ -126,9 +130,9 @@ The minimization procedure stops if any of several criteria are met:
.. note::

   You can also use the :doc:`fix halt <fix_halt>` command to specify
   a general criterion for exiting a minimization, that is a calculation
   performed on the state of the current system, as defined by an
   :doc:`equal-style variable <variable>`.
   a general criterion for exiting a minimization, that is a
   calculation performed on the state of the current system, as
   defined by an :doc:`equal-style variable <variable>`.

For the first criterion, the specified energy tolerance *etol* is
unitless; it is met when the energy change between successive
@@ -163,8 +167,8 @@ freedom, such as from the :doc:`fix box/relax <fix_box_relax>` command.

Following minimization, a statistical summary is printed that lists
which convergence criterion caused the minimizer to stop, as well as
information about the energy, force, final line search, and
iteration counts.  An example is as follows:
information about the energy, force, final line search, and iteration
counts.  An example is as follows:


.. parsed-literal::
@@ -336,9 +340,6 @@ Related commands
:doc:`min_modify <min_modify>`, :doc:`min_style <min_style>`,
:doc:`run_style <run_style>`

**Default:** none
**Default:**


.. _lws: http://lammps.sandia.gov
.. _ld: Manual.html
.. _lc: Commands_all.html
none
+26 −9
Original line number Diff line number Diff line
@@ -72,6 +72,7 @@ allocaters
allosws
AlO
Alonso
alphashrink
amap
Amatrix
AMD
@@ -273,8 +274,8 @@ Bouguet
Bourne
boxcolor
bp
bpls
bpclermont
bpls
br
Branduardi
Branicio
@@ -557,6 +558,7 @@ defgrad
defn
deformable
del
delaystep
deleteIDs
Dellago
delocalization
@@ -660,7 +662,9 @@ DRUDE
dsf
dsmc
dt
dtgrow
dtheta
dtshrink
du
dU
Dudarev
@@ -816,6 +820,7 @@ ethernet
etol
etot
etotal
eulerimplicit
Europhys
ev
eV
@@ -958,10 +963,10 @@ Gahler
gainsboro
Galindo
gamA
gamR
gammaA
gammaFactor
gammaR
gamR
Gan
GaN
ganzenmueller
@@ -1053,6 +1058,7 @@ gstyle
GTL
Gubbins
Guericke
Guenole
gui
Gumbsch
Gunsteren
@@ -1064,6 +1070,7 @@ gz
gzipped
Haak
Hafskjold
halfstepback
Halperin
Halver
Hamaker
@@ -1124,6 +1131,7 @@ Homebrew
hooke
Hookean
hotpink
Houlle
howto
Howto
Hoyt
@@ -1197,8 +1205,9 @@ infty
inhomogeneities
inhomogeneous
init
initio
initialdelay
initializations
initio
InP
inregion
Institut
@@ -1348,7 +1357,6 @@ kcl
Kd
KDevelop
ke
kepler
KE
Keblinski
keflag
@@ -1356,6 +1364,7 @@ Keir
Kelchner
Kelkar
Kemper
kepler
keV
Keyes
Khersonskii
@@ -1451,8 +1460,8 @@ lcbop
ld
ldfftw
ldg
Lebedeva
lebedeva
Lebedeva
Lebold
Lechman
Lehoucq
@@ -1715,8 +1724,8 @@ Mie
Mikami
Militzer
Minary
Mindlin
mincap
Mindlin
mingw
minima
minimizations
@@ -1957,6 +1966,7 @@ Nocedal
nocite
nocoeff
nodeless
Noehring
noforce
Noid
nolib
@@ -2076,8 +2086,8 @@ ons
OO
opencl
openKIM
OpenMP
openmp
OpenMP
openmpi
opls
Oppelstrup
@@ -2250,6 +2260,7 @@ powderblue
ppn
pppm
prd
Prakash
pre
Pre
prec
@@ -2283,6 +2294,7 @@ pthread
pthreads
ptm
PTM
ptol
ptr
pu
purdue
@@ -2388,6 +2400,7 @@ reflectionstyle
regoin
Reinders
reinit
relaxbox
relink
relTol
remappings
@@ -2687,9 +2700,9 @@ Startstep
statcoul
statcoulombs
statvolt
stdlib
stdin
stdio
stdlib
steelblue
Stegailov
Steinbach
@@ -2853,6 +2866,7 @@ tlsph
tmax
Tmax
tmd
tmin
Tmin
tmp
tN
@@ -2992,6 +3006,7 @@ uvm
uwo
Uzdin
vacf
Vaid
valent
Valeriu
valgrind
@@ -3009,6 +3024,7 @@ Vashishta
vasp
vcm
Vcm
vdfmax
vdim
vdisplace
vdW
@@ -3138,6 +3154,7 @@ xhi
xHost
Xiaohu
Xiaowang
Xie
xk
xlat
xlo
@@ -3163,8 +3180,8 @@ xy
xyz
xz
xzhou
YAFF
yaff
YAFF
Yamada
Yazdani
Ybar
Loading