Merge branch 'develop' into remove-error-docs-in-header
This commit is contained in:
Axel Kohlmeyer
@ -33,6 +33,7 @@ KOKKOS, o = OPENMP, t = OPT.
|
||||
* :doc:`body/local <compute_body_local>`
|
||||
* :doc:`bond <compute_bond>`
|
||||
* :doc:`bond/local <compute_bond_local>`
|
||||
* :doc:`born/matrix <compute_born_matrix>`
|
||||
* :doc:`centro/atom <compute_centro_atom>`
|
||||
* :doc:`centroid/stress/atom <compute_stress_atom>`
|
||||
* :doc:`chunk/atom <compute_chunk_atom>`
|
||||
|
||||
@ -124,7 +124,7 @@ OPT.
|
||||
* :doc:`hbond/dreiding/lj (o) <pair_hbond_dreiding>`
|
||||
* :doc:`hbond/dreiding/morse (o) <pair_hbond_dreiding>`
|
||||
* :doc:`hdnnp <pair_hdnnp>`
|
||||
* :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>`
|
||||
* :doc:`ilp/graphene/hbn (t) <pair_ilp_graphene_hbn>`
|
||||
* :doc:`ilp/tmd <pair_ilp_tmd>`
|
||||
* :doc:`kolmogorov/crespi/full <pair_kolmogorov_crespi_full>`
|
||||
* :doc:`kolmogorov/crespi/z <pair_kolmogorov_crespi_z>`
|
||||
|
||||
@ -214,6 +214,9 @@ Convenience functions
|
||||
.. doxygenfunction:: errorurl
|
||||
:project: progguide
|
||||
|
||||
.. doxygenfunction:: missing_cmd_args
|
||||
:project: progguide
|
||||
|
||||
.. doxygenfunction:: flush_buffers(LAMMPS *lmp)
|
||||
:project: progguide
|
||||
|
||||
|
||||
@ -3,7 +3,6 @@ Bonded particle models
|
||||
|
||||
The BPM package implements bonded particle models which can be used to
|
||||
simulate mesoscale solids. Solids are constructed as a collection of
|
||||
|
||||
particles which each represent a coarse-grained region of space much
|
||||
larger than the atomistic scale. Particles within a solid region are
|
||||
then connected by a network of bonds to provide solid elasticity.
|
||||
@ -47,33 +46,29 @@ this, LAMMPS requires :doc:`newton <newton>` bond off such that all
|
||||
processors containing an atom know when a bond breaks. Additionally,
|
||||
one must do either (A) or (B).
|
||||
|
||||
(A)
|
||||
A) Use the following special bond settings
|
||||
|
||||
Use the following special bond settings
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
special_bonds lj 0 1 1 coul 1 1 1
|
||||
|
||||
special_bonds lj 0 1 1 coul 1 1 1
|
||||
These settings accomplish two goals. First, they turn off 1-3 and 1-4
|
||||
special bond lists, which are not currently supported for BPMs. As
|
||||
BPMs often have dense bond networks, generating 1-3 and 1-4 special
|
||||
bond lists is expensive. By setting the lj weight for 1-2 bonds to
|
||||
zero, this turns off pairwise interactions. Even though there are no
|
||||
charges in BPM models, setting a nonzero coul weight for 1-2 bonds
|
||||
ensures all bonded neighbors are still included in the neighbor list
|
||||
in case bonds break between neighbor list builds.
|
||||
|
||||
These settings accomplish two goals. First, they turn off 1-3 and 1-4
|
||||
special bond lists, which are not currently supported for BPMs. As
|
||||
BPMs often have dense bond networks, generating 1-3 and 1-4 special
|
||||
bond lists is expensive. By setting the lj weight for 1-2 bonds to
|
||||
zero, this turns off pairwise interactions. Even though there are no
|
||||
charges in BPM models, setting a nonzero coul weight for 1-2 bonds
|
||||
ensures all bonded neighbors are still included in the neighbor list
|
||||
in case bonds break between neighbor list builds.
|
||||
B) Alternatively, one can simply overlay pair interactions such that all
|
||||
bonded particles also feel pair interactions. This can be
|
||||
accomplished by using the *overlay/pair* keyword present in all bpm
|
||||
bond styles and by using the following special bond settings
|
||||
|
||||
(B)
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
Alternatively, one can simply overlay pair interactions such that all
|
||||
bonded particles also feel pair interactions. This can be accomplished
|
||||
by using the *overlay/pair* keyword present in all bpm bond styles and
|
||||
by using the following special bond settings
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
special_bonds lj/coul 1 1 1
|
||||
special_bonds lj/coul 1 1 1
|
||||
|
||||
See the :doc:`Howto <Howto_broken_bonds>` page on broken bonds for
|
||||
more information.
|
||||
|
||||
@ -18,23 +18,52 @@ At zero temperature, it is easy to estimate these derivatives by
|
||||
deforming the simulation box in one of the six directions using the
|
||||
:doc:`change_box <change_box>` command and measuring the change in the
|
||||
stress tensor. A general-purpose script that does this is given in the
|
||||
examples/elastic directory described on the :doc:`Examples <Examples>`
|
||||
examples/ELASTIC directory described on the :doc:`Examples <Examples>`
|
||||
doc page.
|
||||
|
||||
Calculating elastic constants at finite temperature is more
|
||||
challenging, because it is necessary to run a simulation that performs
|
||||
time averages of differential properties. One way to do this is to
|
||||
measure the change in average stress tensor in an NVT simulations when
|
||||
time averages of differential properties. There are at least
|
||||
3 ways to do this in LAMMPS. The most reliable way to do this is
|
||||
by exploiting the relationship between elastic constants, stress
|
||||
fluctuations, and the Born matrix, the second derivatives of energy
|
||||
w.r.t. strain :ref:`(Ray) <Ray>`.
|
||||
The Born matrix calculation has been enabled by
|
||||
the :doc:`compute born/matrix <compute_born_matrix>` command,
|
||||
which works for any bonded or non-bonded potential in LAMMPS.
|
||||
The most expensive part of the calculation is the sampling of
|
||||
the stress fluctuations. Several examples of this method are
|
||||
provided in the examples/ELASTIC_T/BORN_MATRIX directory
|
||||
described on the :doc:`Examples <Examples>` doc page.
|
||||
|
||||
A second way is to measure
|
||||
the change in average stress tensor in an NVT simulations when
|
||||
the cell volume undergoes a finite deformation. In order to balance
|
||||
the systematic and statistical errors in this method, the magnitude of
|
||||
the deformation must be chosen judiciously, and care must be taken to
|
||||
fully equilibrate the deformed cell before sampling the stress
|
||||
tensor. Another approach is to sample the triclinic cell fluctuations
|
||||
tensor. An example of this method is provided in the
|
||||
examples/ELASTIC_T/DEFORMATION directory
|
||||
described on the :doc:`Examples <Examples>` doc page.
|
||||
|
||||
Another approach is to sample the triclinic cell fluctuations
|
||||
that occur in an NPT simulation. This method can also be slow to
|
||||
converge and requires careful post-processing :ref:`(Shinoda) <Shinoda1>`
|
||||
converge and requires careful post-processing :ref:`(Shinoda) <Shinoda1>`.
|
||||
We do not provide an example of this method.
|
||||
|
||||
A nice review of the advantages and disadvantages of all of these methods
|
||||
is provided in the paper by Clavier et al. :ref:`(Clavier) <Clavier>`.
|
||||
|
||||
----------
|
||||
|
||||
.. _Ray:
|
||||
|
||||
**(Ray)** J. R. Ray and A. Rahman, J Chem Phys, 80, 4423 (1984).
|
||||
|
||||
.. _Shinoda1:
|
||||
|
||||
**(Shinoda)** Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
|
||||
|
||||
.. _Clavier:
|
||||
|
||||
**(Clavier)** G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend and B. Rousseau, Mol Sim, 43, 1413 (2017).
|
||||
|
||||
@ -179,6 +179,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
|
||||
* :doc:`body/local <compute_body_local>` - attributes of body sub-particles
|
||||
* :doc:`bond <compute_bond>` - energy of each bond sub-style
|
||||
* :doc:`bond/local <compute_bond_local>` - distance and energy of each bond
|
||||
* :doc:`born/matrix <compute_born_matrix>` - second derivative or potential with respect to strain
|
||||
* :doc:`centro/atom <compute_centro_atom>` - centro-symmetry parameter for each atom
|
||||
* :doc:`centroid/stress/atom <compute_stress_atom>` - centroid based stress tensor for each atom
|
||||
* :doc:`chunk/atom <compute_chunk_atom>` - assign chunk IDs to each atom
|
||||
|
||||
213
doc/src/compute_born_matrix.rst
Normal file
213
doc/src/compute_born_matrix.rst
Normal file
@ -0,0 +1,213 @@
|
||||
.. index:: compute born/matrix
|
||||
|
||||
compute born/matrix command
|
||||
===========================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
compute ID group-ID born/matrix keyword value ...
|
||||
|
||||
* ID, group-ID are documented in :doc:`compute <compute>` command
|
||||
* born/matrix = style name of this compute command
|
||||
* zero or more keyword/value pairs may be appended
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
keyword = *numdiff*
|
||||
*numdiff* values = delta virial-ID
|
||||
delta = magnitude of strain (dimensionless)
|
||||
virial-ID = ID of pressure compute for virial (string)
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute 1 all born/matrix
|
||||
compute 1 all born/matrix bond angle
|
||||
compute 1 all born/matrix numdiff 1.0e-4 myvirial
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
Define a compute that calculates
|
||||
:math:`\frac{\partial{}^2U}{\partial\varepsilon_{i}\partial\varepsilon_{j}}` the
|
||||
second derivatives of the potential energy :math:`U` w.r.t. strain
|
||||
tensor :math:`\varepsilon` elements. These values are related to:
|
||||
|
||||
.. math::
|
||||
|
||||
C^{B}_{i,j}=\frac{1}{V}\frac{\partial{}^2U}{\partial{}\varepsilon_{i}\partial\varepsilon_{j}}
|
||||
|
||||
also called the Born term of elastic constants in the stress-stress fluctuation
|
||||
formalism. This quantity can be used to compute the elastic constant tensor.
|
||||
Using the symmetric Voigt notation, the elastic constant tensor can be written
|
||||
as a 6x6 symmetric matrix:
|
||||
|
||||
.. math::
|
||||
|
||||
C_{i,j} = \langle{}C^{B}_{i,j}\rangle
|
||||
+ \frac{V}{k_{B}T}\left(\langle\sigma_{i}\sigma_{j}\rangle\right.
|
||||
\left.- \langle\sigma_{i}\rangle\langle\sigma_{j}\rangle\right)
|
||||
+ \frac{Nk_{B}T}{V}
|
||||
\left(\delta_{i,j}+(\delta_{1,i}+\delta_{2,i}+\delta_{3,i})\right.
|
||||
\left.*(\delta_{1,j}+\delta_{2,j}+\delta_{3,j})\right)
|
||||
|
||||
In the above expression, :math:`\sigma` stands for the virial stress
|
||||
tensor, :math:`\delta` is the Kronecker delta and the usual notation apply for
|
||||
the number of particle, the temperature and volume respectively :math:`N`,
|
||||
:math:`T` and :math:`V`. :math:`k_{B}` is the Boltzmann constant.
|
||||
|
||||
The Born term is a symmetric 6x6 matrix, as is the matrix of second derivatives
|
||||
of potential energy w.r.t strain,
|
||||
whose 21 independent elements are output in this order:
|
||||
|
||||
.. math::
|
||||
|
||||
\begin{matrix}
|
||||
C_{1} & C_{7} & C_{8} & C_{9} & C_{10} & C_{11} \\
|
||||
C_{7} & C_{2} & C_{12} & C_{13} & C_{14} & C_{15} \\
|
||||
\vdots & C_{12} & C_{3} & C_{16} & C_{17} & C_{18} \\
|
||||
\vdots & C_{13} & C_{16} & C_{4} & C_{19} & C_{20} \\
|
||||
\vdots & \vdots & \vdots & C_{19} & C_{5} & C_{21} \\
|
||||
\vdots & \vdots & \vdots & \vdots & C_{21} & C_{6}
|
||||
\end{matrix}
|
||||
|
||||
in this matrix the indices of :math:`C_{k}` value are the corresponding element
|
||||
:math:`k` in the global vector output by this compute. Each term comes from the sum
|
||||
of the derivatives of every contribution to the potential energy
|
||||
in the system as explained in :ref:`(VanWorkum)
|
||||
<VanWorkum>`.
|
||||
|
||||
The output can be accessed using usual Lammps routines:
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute 1 all born/matrix
|
||||
compute 2 all pressure NULL virial
|
||||
variable S1 equal -c_2[1]
|
||||
variable S2 equal -c_2[2]
|
||||
variable S3 equal -c_2[3]
|
||||
variable S4 equal -c_2[4]
|
||||
variable S5 equal -c_2[5]
|
||||
variable S6 equal -c_2[6]
|
||||
fix 1 all ave/time 1 1 1 v_S1 v_S2 v_S3 v_S4 v_S5 v_S6 c_1[*] file born.out
|
||||
|
||||
In this example, the file *born.out* will contain the information needed to
|
||||
compute the first and second terms of the elastic constant matrix in a post
|
||||
processing procedure. The other required quantities can be accessed using any
|
||||
other *LAMMPS* usual method. Several examples of this method are
|
||||
provided in the examples/ELASTIC_T/BORN_MATRIX directory
|
||||
described on the :doc:`Examples <Examples>` doc page.
|
||||
|
||||
NOTE: In the above :math:`C_{i,j}` computation, the fluctuation
|
||||
term involving the virial stress tensor :math:`\sigma` is the
|
||||
covariance between each elements. In a
|
||||
solid the stress fluctuations can vary rapidly, while average
|
||||
fluctuations can be slow to converge.
|
||||
A detailed analysis of the convergence rate of all the terms in
|
||||
the elastic tensor
|
||||
is provided in the paper by Clavier et al. :ref:`(Clavier) <Clavier2>`.
|
||||
|
||||
Two different computation methods for the Born matrix are implemented in this
|
||||
compute and are mutually exclusive.
|
||||
|
||||
The first one is a direct computation from the analytical formula from the
|
||||
different terms of the potential used for the simulations :ref:`(VanWorkum)
|
||||
<VanWorkum>`. However, the implementation of such derivations must be done
|
||||
for every potential form. This has not been done yet and can be very
|
||||
complicated for complex potentials. At the moment a warning message is
|
||||
displayed for every term that is not supporting the compute at the moment.
|
||||
This method is the default for now.
|
||||
|
||||
The second method uses finite differences of energy to numerically approximate
|
||||
the second derivatives :ref:`(Zhen) <Zhen>`. This is useful when using
|
||||
interaction styles for which the analytical second derivatives have not been
|
||||
implemented. In this cases, the compute applies linear strain fields of
|
||||
magnitude *delta* to all the atoms relative to a point at the center of the
|
||||
box. The strain fields are in six different directions, corresponding to the
|
||||
six Cartesian components of the stress tensor defined by LAMMPS. For each
|
||||
direction it applies the strain field in both the positive and negative senses,
|
||||
and the new stress virial tensor of the entire system is calculated after each.
|
||||
The difference in these two virials divided by two times *delta*, approximates
|
||||
the corresponding components of the second derivative, after applying a
|
||||
suitable unit conversion.
|
||||
|
||||
.. note::
|
||||
|
||||
It is important to choose a suitable value for delta, the magnitude of
|
||||
strains that are used to generate finite difference
|
||||
approximations to the exact virial stress. For typical systems, a value in
|
||||
the range of 1 part in 1e5 to 1e6 will be sufficient.
|
||||
However, the best value will depend on a multitude of factors
|
||||
including the stiffness of the interatomic potential, the thermodynamic
|
||||
state of the material being probed, and so on. The only way to be sure
|
||||
that you have made a good choice is to do a sensitivity study on a
|
||||
representative atomic configuration, sweeping over a wide range of
|
||||
values of delta. If delta is too small, the output values will vary
|
||||
erratically due to truncation effects. If delta is increased beyond a
|
||||
certain point, the output values will start to vary smoothly with
|
||||
delta, due to growing contributions from higher order derivatives. In
|
||||
between these two limits, the numerical virial values should be largely
|
||||
independent of delta.
|
||||
|
||||
The keyword requires the additional arguments *delta* and *virial-ID*.
|
||||
*delta* gives the size of the applied strains. *virial-ID* gives
|
||||
the ID string of the pressure compute that provides the virial stress tensor,
|
||||
requiring that it use the virial keyword e.g.
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute myvirial all pressure NULL virial
|
||||
compute 1 all born/matrix numdiff 1.0e-4 myvirial
|
||||
|
||||
**Output info:**
|
||||
|
||||
This compute calculates a global vector with 21 values that are
|
||||
the second derivatives of the potential energy w.r.t. strain.
|
||||
The values are in energy units.
|
||||
The values are ordered as explained above. These values can be used
|
||||
by any command that uses global values from a compute as input. See
|
||||
the :doc:`Howto output <Howto_output>` doc page for an overview of
|
||||
LAMMPS output options.
|
||||
|
||||
The array values calculated by this compute are all "extensive".
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This compute is part of the EXTRA-COMPUTE package. It is only enabled if
|
||||
LAMMPS was built with that package. See the :doc:`Build package
|
||||
<Build_package>` page for more info. LAMMPS was built with that package. See
|
||||
the :doc:`Build package <Build_package>` page for more info.
|
||||
|
||||
The Born term can be decomposed as a product of two terms. The first one is a
|
||||
general term which depends on the configuration. The second one is specific to
|
||||
every interaction composing your force field (non-bonded, bonds, angle...).
|
||||
Currently not all LAMMPS interaction styles implement the *born_matrix* method
|
||||
giving first and second order derivatives and LAMMPS will exit with an error if
|
||||
this compute is used with such interactions unless the *numdiff* option is
|
||||
also used. The *numdiff* option cannot be used with any other keyword. In this
|
||||
situation, LAMMPS will also exit with an error.
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
none
|
||||
|
||||
----------
|
||||
|
||||
.. _VanWorkum:
|
||||
|
||||
**(Van Workum)** K. Van Workum et al., J. Chem. Phys. 125 144506 (2006)
|
||||
|
||||
.. _Clavier2:
|
||||
|
||||
**(Clavier)** G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend and B. Rousseau, Mol Sim, 43, 1413 (2017).
|
||||
|
||||
.. _Zhen:
|
||||
|
||||
**(Zhen)** Y. Zhen, C. Chu, Computer Physics Communications 183(2012)261-265
|
||||
@ -76,21 +76,28 @@ velocity for each atom. Note that if there is only one atom in the
|
||||
bin, its thermal velocity will thus be 0.0.
|
||||
|
||||
After the spatially-averaged velocity field has been subtracted from
|
||||
each atom, the temperature is calculated by the formula KE = (dim\*N
|
||||
- dim\*Nx\*Ny\*Nz) k T/2, where KE = total kinetic energy of the group of
|
||||
atoms (sum of 1/2 m v\^2), dim = 2 or 3 = dimensionality of the
|
||||
simulation, N = number of atoms in the group, k = Boltzmann constant,
|
||||
and T = temperature. The dim\*Nx\*Ny\*Nz term are degrees of freedom
|
||||
subtracted to adjust for the removal of the center-of-mass velocity in
|
||||
each of Nx\*Ny\*Nz bins, as discussed in the :ref:`(Evans) <Evans1>` paper.
|
||||
each atom, the temperature is calculated by the formula
|
||||
*KE* = (*dim\*N* - *Ns\*Nx\*Ny\*Nz* - *extra* ) *k* *T*/2, where *KE* = total
|
||||
kinetic energy of the group of atoms (sum of 1/2 *m* *v*\^2), *dim* = 2
|
||||
or 3 = dimensionality of the simulation, *Ns* = 0, 1, 2 or 3 for
|
||||
streaming velocity subtracted in 0, 1, 2 or 3 dimensions, *extra* = extra
|
||||
degrees-of-freedom, *N* = number of atoms in the group, *k* = Boltzmann
|
||||
constant, and *T* = temperature. The *Ns\*Nx\*Ny\*Nz* term is degrees
|
||||
of freedom subtracted to adjust for the removal of the center-of-mass
|
||||
velocity in each direction of the *Nx\*Ny\*Nz* bins, as discussed in the
|
||||
:ref:`(Evans) <Evans1>` paper. The extra term defaults to (*dim* - *Ns*)
|
||||
and accounts for overall conservation of center-of-mass velocity across
|
||||
the group in directions where streaming velocity is *not* subtracted. This
|
||||
can be altered using the *extra* option of the
|
||||
:doc:`compute_modify <compute_modify>` command.
|
||||
|
||||
If the *out* keyword is used with a *tensor* value, which is the
|
||||
default, a kinetic energy tensor, stored as a 6-element vector, is
|
||||
also calculated by this compute for use in the computation of a
|
||||
pressure tensor. The formula for the components of the tensor is the
|
||||
same as the above formula, except that v\^2 is replaced by vx\*vy for
|
||||
the xy component, etc. The 6 components of the vector are ordered xx,
|
||||
yy, zz, xy, xz, yz.
|
||||
same as the above formula, except that *v*\^2 is replaced by *vx\*vy* for
|
||||
the xy component, etc. The 6 components of the vector are ordered *xx,
|
||||
yy, zz, xy, xz, yz.*
|
||||
|
||||
If the *out* keyword is used with a *bin* value, the count of atoms
|
||||
and computed temperature for each bin are stored for output, as an
|
||||
@ -123,10 +130,20 @@ needed, the subtracted degrees-of-freedom can be altered using the
|
||||
.. note::
|
||||
|
||||
When using the *out* keyword with a value of *bin*, the
|
||||
calculated temperature for each bin does not include the
|
||||
degrees-of-freedom adjustment described in the preceding paragraph,
|
||||
for fixes that constrain molecular motion. It does include the
|
||||
adjustment due to the *extra* option, which is applied to each bin.
|
||||
calculated temperature for each bin includes the degrees-of-freedom
|
||||
adjustment described in the preceding paragraph for fixes that
|
||||
constrain molecular motion, as well as the adjustment due to
|
||||
the *extra* option (which defaults to *dim* - *Ns* as described above),
|
||||
by fractionally applying them based on the fraction of atoms in each
|
||||
bin. As a result, the bin degrees-of-freedom summed over all bins exactly
|
||||
equals the degrees-of-freedom used in the scalar temperature calculation,
|
||||
:math:`\Sigma N_{DOF_i} = N_{DOF}` and the corresponding relation for temperature
|
||||
is also satisfied :math:`\Sigma N_{DOF_i} T_i = N_{DOF} T`.
|
||||
These relations will breakdown in cases where the adjustment
|
||||
exceeds the actual number of degrees-of-freedom in a bin. This could happen
|
||||
if a bin is empty or in situations where rigid molecules
|
||||
are non-uniformly distributed, in which case the reported
|
||||
temperature within a bin may not be accurate.
|
||||
|
||||
See the :doc:`Howto thermostat <Howto_thermostat>` page for a
|
||||
discussion of different ways to compute temperature and perform
|
||||
|
||||
@ -305,7 +305,7 @@ with fix_adapt are
|
||||
+------------------------------------+-------+-----------------+
|
||||
| :doc:`fene <bond_fene>` | k,r0 | type bonds |
|
||||
+------------------------------------+-------+-----------------+
|
||||
| :doc:`fene/nm <bond_fene_nm>` | k,r0 | type bonds |
|
||||
| :doc:`fene/nm <bond_fene>` | k,r0 | type bonds |
|
||||
+------------------------------------+-------+-----------------+
|
||||
| :doc:`gromos <bond_gromos>` | k,r0 | type bonds |
|
||||
+------------------------------------+-------+-----------------+
|
||||
|
||||
@ -36,7 +36,7 @@ are (full) periodic boundary conditions and no other "manipulations"
|
||||
of the system (e.g. fixes that modify forces or velocities).
|
||||
|
||||
This fix invokes the velocity form of the
|
||||
Störmer-Verlet time integration algorithm (velocity-Verlet). Other
|
||||
Stoermer-Verlet time integration algorithm (velocity-Verlet). Other
|
||||
time integration options can be invoked using the :doc:`run_style <run_style>` command.
|
||||
|
||||
----------
|
||||
|
||||
@ -1,8 +1,11 @@
|
||||
.. index:: pair_style ilp/graphene/hbn
|
||||
.. index:: pair_style ilp/graphene/hbn/opt
|
||||
|
||||
pair_style ilp/graphene/hbn command
|
||||
===================================
|
||||
|
||||
Accelerator Variant: *ilp/graphene/hbn/opt*
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
@ -125,6 +128,10 @@ headings) the following commands could be included in an input script:
|
||||
|
||||
----------
|
||||
|
||||
.. include:: accel_styles.rst
|
||||
|
||||
----------
|
||||
|
||||
Mixing, shift, table, tail correction, restart, rRESPA info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
|
||||
|
||||
@ -68,7 +68,7 @@ Choose the style of time integrator used for molecular dynamics
|
||||
simulations performed by LAMMPS.
|
||||
|
||||
The *verlet* style is the velocity form of the
|
||||
Störmer-Verlet time integration algorithm (velocity-Verlet)
|
||||
Stoermer-Verlet time integration algorithm (velocity-Verlet)
|
||||
|
||||
----------
|
||||
|
||||
|
||||
@ -315,6 +315,7 @@ borophene
|
||||
Botero
|
||||
Botu
|
||||
Bouguet
|
||||
Bourasseau
|
||||
Bourne
|
||||
boxcolor
|
||||
boxhi
|
||||
@ -342,6 +343,7 @@ Broglie
|
||||
brownian
|
||||
brownw
|
||||
Broyden
|
||||
Brusselle
|
||||
Bryantsev
|
||||
Btarget
|
||||
btype
|
||||
@ -674,6 +676,7 @@ Deresiewicz
|
||||
Derjagin
|
||||
Derjaguin
|
||||
Derlet
|
||||
Desbiens
|
||||
Deserno
|
||||
Destree
|
||||
destructor
|
||||
@ -787,6 +790,7 @@ Dullweber
|
||||
dumpfile
|
||||
Dunbrack
|
||||
Dunweg
|
||||
Dupend
|
||||
Dupont
|
||||
dUs
|
||||
dV
|
||||
@ -1669,6 +1673,7 @@ Kusters
|
||||
Kutta
|
||||
Kuznetsov
|
||||
kx
|
||||
Lachet
|
||||
Lackmann
|
||||
Ladd
|
||||
lagrangian
|
||||
@ -3188,6 +3193,7 @@ stochastically
|
||||
stochasticity
|
||||
Stockmayer
|
||||
Stoddard
|
||||
Stoermer
|
||||
stoichiometric
|
||||
stoichiometry
|
||||
Stokesian
|
||||
@ -3599,6 +3605,7 @@ Voronoi
|
||||
VORONOI
|
||||
Vorselaars
|
||||
Voth
|
||||
Voyiatzis
|
||||
vpz
|
||||
vratio
|
||||
Vries
|
||||
@ -3665,6 +3672,7 @@ wn
|
||||
Wolde
|
||||
workflow
|
||||
workflows
|
||||
Workum
|
||||
Worley
|
||||
Wriggers
|
||||
Wuppertal
|
||||
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
3.36316 1.87373 1.87607 -0.00346 -0.00172 -0.00104
|
||||
1.87373 3.36170 1.87425 0.00443 0.00033 0.00014
|
||||
1.87607 1.87425 3.36573 0.00143 0.00155 0.00127
|
||||
-0.00346 0.00443 0.00143 1.87425 0.00127 0.00033
|
||||
-0.00172 0.00033 0.00155 0.00127 1.87607 -0.00346
|
||||
-0.00104 0.00014 0.00127 0.00033 -0.00346 1.87373
|
||||
118
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py
Normal file
118
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py
Normal file
@ -0,0 +1,118 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
|
||||
import sys
|
||||
import numpy as np
|
||||
|
||||
def reduce_Born(Cf):
|
||||
C = np.zeros((6,6), dtype=np.float64)
|
||||
C[0,0] = Cf[0]
|
||||
C[1,1] = Cf[1]
|
||||
C[2,2] = Cf[2]
|
||||
C[3,3] = Cf[3]
|
||||
C[4,4] = Cf[4]
|
||||
C[5,5] = Cf[5]
|
||||
C[0,1] = Cf[6]
|
||||
C[0,2] = Cf[7]
|
||||
C[0,3] = Cf[8]
|
||||
C[0,4] = Cf[9]
|
||||
C[0,5] = Cf[10]
|
||||
C[1,2] = Cf[11]
|
||||
C[1,3] = Cf[12]
|
||||
C[1,4] = Cf[13]
|
||||
C[1,5] = Cf[14]
|
||||
C[2,3] = Cf[15]
|
||||
C[2,4] = Cf[16]
|
||||
C[2,5] = Cf[17]
|
||||
C[3,4] = Cf[18]
|
||||
C[3,5] = Cf[19]
|
||||
C[4,5] = Cf[20]
|
||||
C = np.where(C,C,C.T)
|
||||
return C
|
||||
|
||||
def compute_delta():
|
||||
D = np.zeros((3,3,3,3))
|
||||
for a in range(3):
|
||||
for b in range(3):
|
||||
for m in range(3):
|
||||
for n in range(3):
|
||||
D[a,b,m,n] = (a==m)*(b==n) + (a==n)*(b==m)
|
||||
d = np.zeros((6,6))
|
||||
d[0,0] = D[0,0,0,0]
|
||||
d[1,1] = D[1,1,1,1]
|
||||
d[2,2] = D[2,2,2,2]
|
||||
d[3,3] = D[1,2,1,2]
|
||||
d[4,4] = D[0,2,0,2]
|
||||
d[5,5] = D[0,1,0,1]
|
||||
d[0,1] = D[0,0,1,1]
|
||||
d[0,2] = D[0,0,2,2]
|
||||
d[0,3] = D[0,0,1,2]
|
||||
d[0,4] = D[0,0,0,2]
|
||||
d[0,5] = D[0,0,0,1]
|
||||
d[1,2] = D[1,1,2,2]
|
||||
d[1,3] = D[1,1,1,2]
|
||||
d[1,4] = D[1,1,0,2]
|
||||
d[1,5] = D[1,1,0,1]
|
||||
d[2,3] = D[2,2,1,2]
|
||||
d[2,4] = D[2,2,0,2]
|
||||
d[2,5] = D[2,2,0,1]
|
||||
d[3,4] = D[1,2,0,2]
|
||||
d[3,5] = D[1,2,0,1]
|
||||
d[4,5] = D[0,2,0,1]
|
||||
d = np.where(d,d,d.T)
|
||||
return d
|
||||
|
||||
|
||||
def write_matrix(C, filename):
|
||||
with open(filename, 'w') as f:
|
||||
f.write("# Cij Matrix from post process computation\n")
|
||||
for i in C:
|
||||
f.write("{:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f}\n".format(
|
||||
i[0]*10**-9, i[1]*10**-9, i[2]*10**-9, i[3]*10**-9, i[4]*10**-9, i[5]*10**-9,
|
||||
)
|
||||
)
|
||||
return
|
||||
|
||||
def main():
|
||||
|
||||
N = 500
|
||||
vol = 27.047271**3 * 10**-30 # m^3
|
||||
T = 60 # K
|
||||
kb = 1.380649 * 10**-23 # J/K
|
||||
kbT = T*kb # J
|
||||
kcalmol2J = 4183.9954/(6.022*10**23)
|
||||
|
||||
born = np.loadtxt('born.out')
|
||||
stre = np.loadtxt('vir.out')
|
||||
stre[:, 1:] = -stre[:, 1:]*101325 # -> Pa
|
||||
try:
|
||||
mean_born = np.mean(born[:, 1:], axis=0)
|
||||
except IndexError:
|
||||
mean_born = born[1:]
|
||||
|
||||
CB = kcalmol2J/vol*reduce_Born(mean_born) # -> J/m^3=Pa
|
||||
Cs = vol/kbT*np.cov(stre[:,1:].T)
|
||||
Ct = N*kbT/vol * compute_delta()
|
||||
C = CB - Cs + Ct
|
||||
write_matrix(CB, 'born_matrix.out')
|
||||
write_matrix(Cs, 'stre_matrix.out')
|
||||
write_matrix(Ct, 'temp_matrix.out')
|
||||
write_matrix(C, 'full_matrix.out')
|
||||
C11 = np.mean([C[0,0], C[1,1], C[2,2]]) * 10**-9
|
||||
C12 = np.mean([C[0,1], C[0,2], C[1,2]]) * 10**-9
|
||||
C44 = np.mean([C[3,3], C[4,4], C[5,5]]) * 10**-9
|
||||
eC11 = np.std([C[0,0], C[1,1], C[2,2]]) * 10**-9
|
||||
eC12 = np.std([C[0,1], C[0,2], C[1,2]]) * 10**-9
|
||||
eC44 = np.std([C[3,3], C[4,4], C[5,5]]) * 10**-9
|
||||
print(C*10**-9)
|
||||
print("C11 = {:f} ± {:f}; C12 = {:f} ± {:f}; C44 = {:f} ± {:f}".format(C11, eC11, C12, eC12, C44, eC44))
|
||||
|
||||
return
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
try:
|
||||
main()
|
||||
except KeyboardInterrupt:
|
||||
raise SystemExit("User interruption.")
|
||||
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
2.18161 1.13726 1.16596 -0.01607 -0.02637 0.00291
|
||||
1.13726 2.20242 1.16714 0.00386 -0.05820 0.02644
|
||||
1.16596 1.16714 2.24704 -0.00354 -0.00368 0.02714
|
||||
-0.01607 0.00386 -0.00354 1.43706 0.00210 0.01003
|
||||
-0.02637 -0.05820 -0.00368 0.00210 1.37530 0.01401
|
||||
0.00291 0.02644 0.02714 0.01003 0.01401 1.42403
|
||||
154
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/in.ljcov
Normal file
154
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/in.ljcov
Normal file
@ -0,0 +1,154 @@
|
||||
# Analytical calculation
|
||||
# of Born matrix
|
||||
|
||||
# Note that because of cubic symmetry and central forces, we have:
|
||||
# C11, pure axial == positive mean value: 1,2,3
|
||||
# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7)
|
||||
# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19)
|
||||
# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17
|
||||
|
||||
# adjustable parameters
|
||||
|
||||
units real
|
||||
variable nsteps index 100000 # length of run
|
||||
variable nthermo index 1000 # thermo output interval
|
||||
variable nlat equal 5 # size of box
|
||||
variable T equal 60. # Temperature in K
|
||||
variable rho equal 5.405 # Lattice spacing in A
|
||||
|
||||
atom_style atomic
|
||||
|
||||
lattice fcc ${rho}
|
||||
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
|
||||
mass * 39.948
|
||||
|
||||
velocity all create ${T} 87287 loop geom
|
||||
velocity all zero linear
|
||||
|
||||
pair_style lj/cut 12.0
|
||||
pair_coeff 1 1 0.238067 3.405
|
||||
|
||||
neighbor 0.0 bin
|
||||
neigh_modify every 1 delay 0 check no
|
||||
|
||||
variable vol equal vol
|
||||
thermo 100
|
||||
fix aL all ave/time 1 1 1 v_vol ave running
|
||||
fix NPT all npt temp $T $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0.
|
||||
|
||||
run 20000
|
||||
|
||||
unfix NPT
|
||||
|
||||
variable newL equal "f_aL^(1./3.)"
|
||||
change_box all x final 0 ${newL} y final 0. ${newL} z final 0. ${newL} remap units box
|
||||
|
||||
unfix aL
|
||||
|
||||
reset_timestep 0
|
||||
|
||||
# Conversion variables
|
||||
variable kb equal 1.38065e-23 # J/K
|
||||
variable Myvol equal "vol*10^-30" # Volume in m^3
|
||||
variable kbt equal "v_kb*v_T"
|
||||
variable Nat equal atoms
|
||||
variable Rhokbt equal "v_kbt*v_Nat/v_Myvol"
|
||||
variable at2Pa equal 101325
|
||||
variable kcalmol2J equal "4183.9954/(6.022e23)"
|
||||
variable C1 equal "v_kcalmol2J/v_Myvol" # Convert Cb from energy to pressure units
|
||||
variable C2 equal "v_Myvol/v_kbt" # Factor for Cfl terms
|
||||
variable Pa2GPa equal 1e-9
|
||||
|
||||
# Born compute giving <C^b> terms
|
||||
compute born all born/matrix
|
||||
# The six virial stress component to compute <C^fl>
|
||||
compute VIR all pressure NULL virial
|
||||
variable s1 equal "-c_VIR[1]*v_at2Pa"
|
||||
variable s2 equal "-c_VIR[2]*v_at2Pa"
|
||||
variable s3 equal "-c_VIR[3]*v_at2Pa"
|
||||
variable s6 equal "-c_VIR[4]*v_at2Pa"
|
||||
variable s5 equal "-c_VIR[5]*v_at2Pa"
|
||||
variable s4 equal "-c_VIR[6]*v_at2Pa"
|
||||
variable press equal press
|
||||
|
||||
|
||||
# Average of Born term and vector to store stress
|
||||
# for post processing
|
||||
fix CB all ave/time 1 ${nthermo} ${nthermo} c_born[*] ave running file born.out overwrite
|
||||
fix CPR all ave/time 1 1 1 c_VIR[*] file vir.out
|
||||
fix APR all ave/time 1 1 1 v_press ave running
|
||||
fix VEC all vector 1 v_s1 v_s2 v_s3 v_s4 v_s5 v_s6
|
||||
|
||||
thermo ${nthermo}
|
||||
thermo_style custom step temp press f_APR c_born[1] f_CB[1] c_born[12] f_CB[12] c_born[4] f_CB[4]
|
||||
thermo_modify line multi
|
||||
|
||||
fix 1 all nvt temp $T $T 100
|
||||
|
||||
run ${nsteps}
|
||||
|
||||
# Compute vector averages
|
||||
# Note the indice switch.
|
||||
# LAMMPS convention is NOT the Voigt notation.
|
||||
variable aves1 equal "ave(f_VEC[1])"
|
||||
variable aves2 equal "ave(f_VEC[2])"
|
||||
variable aves3 equal "ave(f_VEC[3])"
|
||||
variable aves4 equal "ave(f_VEC[6])"
|
||||
variable aves5 equal "ave(f_VEC[5])"
|
||||
variable aves6 equal "ave(f_VEC[4])"
|
||||
|
||||
# Computing the covariance through the <s_{i}s_{j}>-<s_i><s_j>
|
||||
# is numerically instable. Here we go through the <(s-<s>)^2>
|
||||
# definition.
|
||||
|
||||
# Computing difference relative to average values
|
||||
variable ds1 vector "f_VEC[1]-v_aves1"
|
||||
variable ds2 vector "f_VEC[2]-v_aves2"
|
||||
variable ds3 vector "f_VEC[3]-v_aves3"
|
||||
variable ds4 vector "f_VEC[4]-v_aves4"
|
||||
variable ds5 vector "f_VEC[5]-v_aves5"
|
||||
variable ds6 vector "f_VEC[6]-v_aves6"
|
||||
|
||||
# Squaring and averaging
|
||||
variable dds1 vector "v_ds1*v_ds1"
|
||||
variable dds2 vector "v_ds2*v_ds2"
|
||||
variable dds3 vector "v_ds3*v_ds3"
|
||||
variable vars1 equal "ave(v_dds1)"
|
||||
variable vars2 equal "ave(v_dds2)"
|
||||
variable vars3 equal "ave(v_dds3)"
|
||||
variable C11 equal "v_Pa2GPa*(v_C1*f_CB[1] - v_C2*v_vars1 + 2*v_Rhokbt)"
|
||||
variable C22 equal "v_Pa2GPa*(v_C1*f_CB[2] - v_C2*v_vars2 + 2*v_Rhokbt)"
|
||||
variable C33 equal "v_Pa2GPa*(v_C1*f_CB[3] - v_C2*v_vars3 + 2*v_Rhokbt)"
|
||||
|
||||
variable dds12 vector "v_ds1*v_ds2"
|
||||
variable dds13 vector "v_ds1*v_ds3"
|
||||
variable dds23 vector "v_ds2*v_ds3"
|
||||
variable vars12 equal "ave(v_dds12)"
|
||||
variable vars13 equal "ave(v_dds13)"
|
||||
variable vars23 equal "ave(v_dds23)"
|
||||
variable C12 equal "v_Pa2GPa*(v_C1*f_CB[7] - v_C2*v_vars12)"
|
||||
variable C13 equal "v_Pa2GPa*(v_C1*f_CB[8] - v_C2*v_vars13)"
|
||||
variable C23 equal "v_Pa2GPa*(v_C1*f_CB[12] - v_C2*v_vars23)"
|
||||
|
||||
variable dds4 vector "v_ds4*v_ds4"
|
||||
variable dds5 vector "v_ds5*v_ds5"
|
||||
variable dds6 vector "v_ds6*v_ds6"
|
||||
variable vars4 equal "ave(v_dds4)"
|
||||
variable vars5 equal "ave(v_dds5)"
|
||||
variable vars6 equal "ave(v_dds6)"
|
||||
variable C44 equal "v_Pa2GPa*(v_C1*f_CB[4] - v_C2*v_vars4 + v_Rhokbt)"
|
||||
variable C55 equal "v_Pa2GPa*(v_C1*f_CB[5] - v_C2*v_vars5 + v_Rhokbt)"
|
||||
variable C66 equal "v_Pa2GPa*(v_C1*f_CB[6] - v_C2*v_vars6 + v_Rhokbt)"
|
||||
|
||||
variable aC11 equal "(v_C11 + v_C22 + v_C33)/3."
|
||||
variable aC12 equal "(v_C12 + v_C13 + v_C23)/3."
|
||||
variable aC44 equal "(v_C44 + v_C55 + v_C66)/3."
|
||||
|
||||
print """
|
||||
C11 = ${aC11}
|
||||
C12 = ${aC12}
|
||||
C44 = ${aC44}
|
||||
"""
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
1.22342 0.73647 0.71011 0.01261 0.02465 -0.00395
|
||||
0.73647 1.20115 0.70711 0.00057 0.05854 -0.02630
|
||||
0.71011 0.70711 1.16055 0.00497 0.00524 -0.02587
|
||||
0.01261 0.00057 0.00497 0.45813 -0.00083 -0.00969
|
||||
0.02465 0.05854 0.00524 -0.00083 0.52170 -0.01747
|
||||
-0.00395 -0.02630 -0.02587 -0.00969 -0.01747 0.47064
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
0.04187 0.00000 0.00000 0.00000 0.00000 0.00000
|
||||
0.00000 0.04187 0.00000 0.00000 0.00000 0.00000
|
||||
0.00000 0.00000 0.04187 0.00000 0.00000 0.00000
|
||||
0.00000 0.00000 0.00000 0.02093 0.00000 0.00000
|
||||
0.00000 0.00000 0.00000 0.00000 0.02093 0.00000
|
||||
0.00000 0.00000 0.00000 0.00000 0.00000 0.02093
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
3.35855 1.86892 1.87139 0.00233 0.00218 -0.00179
|
||||
1.86892 3.37104 1.87285 0.00112 0.00085 -0.00007
|
||||
1.87139 1.87285 3.37707 -0.00058 0.00038 -0.00057
|
||||
0.00233 0.00112 -0.00058 1.88326 -0.00039 0.00065
|
||||
0.00218 0.00085 0.00038 -0.00039 1.88229 0.00242
|
||||
-0.00179 -0.00007 -0.00057 0.00065 0.00242 1.87968
|
||||
118
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py
Normal file
118
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py
Normal file
@ -0,0 +1,118 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
|
||||
import sys
|
||||
import numpy as np
|
||||
|
||||
def reduce_Born(Cf):
|
||||
C = np.zeros((6,6), dtype=np.float64)
|
||||
C[0,0] = Cf[0]
|
||||
C[1,1] = Cf[1]
|
||||
C[2,2] = Cf[2]
|
||||
C[3,3] = Cf[3]
|
||||
C[4,4] = Cf[4]
|
||||
C[5,5] = Cf[5]
|
||||
C[0,1] = Cf[6]
|
||||
C[0,2] = Cf[7]
|
||||
C[0,3] = Cf[8]
|
||||
C[0,4] = Cf[9]
|
||||
C[0,5] = Cf[10]
|
||||
C[1,2] = Cf[11]
|
||||
C[1,3] = Cf[12]
|
||||
C[1,4] = Cf[13]
|
||||
C[1,5] = Cf[14]
|
||||
C[2,3] = Cf[15]
|
||||
C[2,4] = Cf[16]
|
||||
C[2,5] = Cf[17]
|
||||
C[3,4] = Cf[18]
|
||||
C[3,5] = Cf[19]
|
||||
C[4,5] = Cf[20]
|
||||
C = np.where(C,C,C.T)
|
||||
return C
|
||||
|
||||
def compute_delta():
|
||||
D = np.zeros((3,3,3,3))
|
||||
for a in range(3):
|
||||
for b in range(3):
|
||||
for m in range(3):
|
||||
for n in range(3):
|
||||
D[a,b,m,n] = (a==m)*(b==n) + (a==n)*(b==m)
|
||||
d = np.zeros((6,6))
|
||||
d[0,0] = D[0,0,0,0]
|
||||
d[1,1] = D[1,1,1,1]
|
||||
d[2,2] = D[2,2,2,2]
|
||||
d[3,3] = D[1,2,1,2]
|
||||
d[4,4] = D[0,2,0,2]
|
||||
d[5,5] = D[0,1,0,1]
|
||||
d[0,1] = D[0,0,1,1]
|
||||
d[0,2] = D[0,0,2,2]
|
||||
d[0,3] = D[0,0,1,2]
|
||||
d[0,4] = D[0,0,0,2]
|
||||
d[0,5] = D[0,0,0,1]
|
||||
d[1,2] = D[1,1,2,2]
|
||||
d[1,3] = D[1,1,1,2]
|
||||
d[1,4] = D[1,1,0,2]
|
||||
d[1,5] = D[1,1,0,1]
|
||||
d[2,3] = D[2,2,1,2]
|
||||
d[2,4] = D[2,2,0,2]
|
||||
d[2,5] = D[2,2,0,1]
|
||||
d[3,4] = D[1,2,0,2]
|
||||
d[3,5] = D[1,2,0,1]
|
||||
d[4,5] = D[0,2,0,1]
|
||||
d = np.where(d,d,d.T)
|
||||
return d
|
||||
|
||||
|
||||
def write_matrix(C, filename):
|
||||
with open(filename, 'w') as f:
|
||||
f.write("# Cij Matrix from post process computation\n")
|
||||
for i in C:
|
||||
f.write("{:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f}\n".format(
|
||||
i[0]*10**-9, i[1]*10**-9, i[2]*10**-9, i[3]*10**-9, i[4]*10**-9, i[5]*10**-9,
|
||||
)
|
||||
)
|
||||
return
|
||||
|
||||
def main():
|
||||
|
||||
N = 500
|
||||
vol = 27.047271**3 * 10**-30 # m^3
|
||||
T = 60 # K
|
||||
kb = 1.380649 * 10**-23 # J/K
|
||||
kbT = T*kb # J
|
||||
kcalmol2J = 4183.9954/(6.022*10**23)
|
||||
|
||||
born = np.loadtxt('born.out')
|
||||
stre = np.loadtxt('vir.out')
|
||||
stre[:, 1:] = -stre[:, 1:]*101325 # -> Pa
|
||||
try:
|
||||
mean_born = np.mean(born[:, 1:], axis=0)
|
||||
except IndexError:
|
||||
mean_born = born[1:]
|
||||
|
||||
CB = kcalmol2J/vol*reduce_Born(mean_born) # -> J/m^3=Pa
|
||||
Cs = vol/kbT*np.cov(stre[:,1:].T)
|
||||
Ct = N*kbT/vol * compute_delta()
|
||||
C = CB - Cs + Ct
|
||||
write_matrix(CB, 'born_matrix.out')
|
||||
write_matrix(Cs, 'stre_matrix.out')
|
||||
write_matrix(Ct, 'temp_matrix.out')
|
||||
write_matrix(C, 'full_matrix.out')
|
||||
C11 = np.mean([C[0,0], C[1,1], C[2,2]]) * 10**-9
|
||||
C12 = np.mean([C[0,1], C[0,2], C[1,2]]) * 10**-9
|
||||
C44 = np.mean([C[3,3], C[4,4], C[5,5]]) * 10**-9
|
||||
eC11 = np.std([C[0,0], C[1,1], C[2,2]]) * 10**-9
|
||||
eC12 = np.std([C[0,1], C[0,2], C[1,2]]) * 10**-9
|
||||
eC44 = np.std([C[3,3], C[4,4], C[5,5]]) * 10**-9
|
||||
print(C*10**-9)
|
||||
print("C11 = {:f} ± {:f}; C12 = {:f} ± {:f}; C44 = {:f} ± {:f}".format(C11, eC11, C12, eC12, C44, eC44))
|
||||
|
||||
return
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
try:
|
||||
main()
|
||||
except KeyboardInterrupt:
|
||||
raise SystemExit("User interruption.")
|
||||
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
2.29922 1.17439 1.20854 -0.01653 -0.01684 0.01188
|
||||
1.17439 2.20673 1.21718 -0.00781 -0.00753 0.00867
|
||||
1.20854 1.21718 2.30804 0.01535 -0.01596 0.00426
|
||||
-0.01653 -0.00781 0.01535 1.47647 -0.01355 -0.01601
|
||||
-0.01684 -0.00753 -0.01596 -0.01355 1.37905 0.01975
|
||||
0.01188 0.00867 0.00426 -0.01601 0.01975 1.40170
|
||||
154
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/in.ljcov
Normal file
154
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/in.ljcov
Normal file
@ -0,0 +1,154 @@
|
||||
# Numerical difference calculation
|
||||
# of Born matrix
|
||||
|
||||
# Note that because of cubic symmetry and central forces, we have:
|
||||
# C11, pure axial == positive mean value: 1,2,3
|
||||
# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7)
|
||||
# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19)
|
||||
# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17
|
||||
|
||||
# adjustable parameters
|
||||
|
||||
units real
|
||||
variable nsteps index 100000 # length of run
|
||||
variable nthermo index 1000 # thermo output interval
|
||||
variable nlat equal 5 # size of box
|
||||
variable T equal 60. # Temperature in K
|
||||
variable rho equal 5.405 # Lattice spacing in A
|
||||
|
||||
atom_style atomic
|
||||
|
||||
lattice fcc ${rho}
|
||||
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
|
||||
mass * 39.948
|
||||
|
||||
velocity all create ${T} 87287 loop geom
|
||||
velocity all zero linear
|
||||
|
||||
pair_style lj/cut 12.0
|
||||
pair_coeff 1 1 0.238067 3.405
|
||||
|
||||
neighbor 0.0 bin
|
||||
neigh_modify every 1 delay 0 check no
|
||||
|
||||
variable vol equal vol
|
||||
thermo 100
|
||||
fix aL all ave/time 1 1 1 v_vol ave running
|
||||
fix NPT all npt temp $T $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0.
|
||||
|
||||
run 20000
|
||||
|
||||
unfix NPT
|
||||
|
||||
variable newL equal "f_aL^(1./3.)"
|
||||
change_box all x final 0 ${newL} y final 0. ${newL} z final 0. ${newL} remap units box
|
||||
|
||||
unfix aL
|
||||
|
||||
reset_timestep 0
|
||||
|
||||
# Conversion variables
|
||||
variable kb equal 1.38065e-23 # J/K
|
||||
variable Myvol equal "vol*10^-30" # Volume in m^3
|
||||
variable kbt equal "v_kb*v_T"
|
||||
variable Nat equal atoms
|
||||
variable Rhokbt equal "v_kbt*v_Nat/v_Myvol"
|
||||
variable at2Pa equal 101325
|
||||
variable kcalmol2J equal "4183.9954/(6.022e23)"
|
||||
variable C1 equal "v_kcalmol2J/v_Myvol" # Convert Cb from energy to pressure units
|
||||
variable C2 equal "v_Myvol/v_kbt" # Factor for Cfl terms
|
||||
variable Pa2GPa equal 1e-9
|
||||
|
||||
# Born compute giving <C^b> terms
|
||||
# The six virial stress component to compute <C^fl>
|
||||
compute VIR all pressure NULL virial
|
||||
compute born all born/matrix numdiff 1e-6 VIR
|
||||
variable s1 equal "-c_VIR[1]*v_at2Pa"
|
||||
variable s2 equal "-c_VIR[2]*v_at2Pa"
|
||||
variable s3 equal "-c_VIR[3]*v_at2Pa"
|
||||
variable s6 equal "-c_VIR[4]*v_at2Pa"
|
||||
variable s5 equal "-c_VIR[5]*v_at2Pa"
|
||||
variable s4 equal "-c_VIR[6]*v_at2Pa"
|
||||
variable press equal press
|
||||
|
||||
|
||||
# Average of Born term and vector to store stress
|
||||
# for post processing
|
||||
fix CB all ave/time 1 ${nthermo} ${nthermo} c_born[*] ave running file born.out overwrite
|
||||
fix CPR all ave/time 1 1 1 c_VIR[*] file vir.out
|
||||
fix APR all ave/time 1 1 1 v_press ave running
|
||||
fix VEC all vector 1 v_s1 v_s2 v_s3 v_s4 v_s5 v_s6
|
||||
|
||||
thermo ${nthermo}
|
||||
thermo_style custom step temp press f_APR c_born[1] f_CB[1] c_born[12] f_CB[12] c_born[4] f_CB[4]
|
||||
thermo_modify line multi
|
||||
|
||||
fix 1 all nvt temp $T $T 100
|
||||
|
||||
run ${nsteps}
|
||||
|
||||
# Compute vector averages
|
||||
# Note the indice switch.
|
||||
# LAMMPS convention is NOT the Voigt notation.
|
||||
variable aves1 equal "ave(f_VEC[1])"
|
||||
variable aves2 equal "ave(f_VEC[2])"
|
||||
variable aves3 equal "ave(f_VEC[3])"
|
||||
variable aves4 equal "ave(f_VEC[6])"
|
||||
variable aves5 equal "ave(f_VEC[5])"
|
||||
variable aves6 equal "ave(f_VEC[4])"
|
||||
|
||||
# Computing the covariance through the <s_{i}s_{j}>-<s_i><s_j>
|
||||
# is numerically instable. Here we go through the <(s-<s>)^2>
|
||||
# definition.
|
||||
|
||||
# Computing difference relative to average values
|
||||
variable ds1 vector "f_VEC[1]-v_aves1"
|
||||
variable ds2 vector "f_VEC[2]-v_aves2"
|
||||
variable ds3 vector "f_VEC[3]-v_aves3"
|
||||
variable ds4 vector "f_VEC[4]-v_aves4"
|
||||
variable ds5 vector "f_VEC[5]-v_aves5"
|
||||
variable ds6 vector "f_VEC[6]-v_aves6"
|
||||
|
||||
# Squaring and averaging
|
||||
variable dds1 vector "v_ds1*v_ds1"
|
||||
variable dds2 vector "v_ds2*v_ds2"
|
||||
variable dds3 vector "v_ds3*v_ds3"
|
||||
variable vars1 equal "ave(v_dds1)"
|
||||
variable vars2 equal "ave(v_dds2)"
|
||||
variable vars3 equal "ave(v_dds3)"
|
||||
variable C11 equal "v_Pa2GPa*(v_C1*f_CB[1] - v_C2*v_vars1 + 2*v_Rhokbt)"
|
||||
variable C22 equal "v_Pa2GPa*(v_C1*f_CB[2] - v_C2*v_vars2 + 2*v_Rhokbt)"
|
||||
variable C33 equal "v_Pa2GPa*(v_C1*f_CB[3] - v_C2*v_vars3 + 2*v_Rhokbt)"
|
||||
|
||||
variable dds12 vector "v_ds1*v_ds2"
|
||||
variable dds13 vector "v_ds1*v_ds3"
|
||||
variable dds23 vector "v_ds2*v_ds3"
|
||||
variable vars12 equal "ave(v_dds12)"
|
||||
variable vars13 equal "ave(v_dds13)"
|
||||
variable vars23 equal "ave(v_dds23)"
|
||||
variable C12 equal "v_Pa2GPa*(v_C1*f_CB[7] - v_C2*v_vars12)"
|
||||
variable C13 equal "v_Pa2GPa*(v_C1*f_CB[8] - v_C2*v_vars13)"
|
||||
variable C23 equal "v_Pa2GPa*(v_C1*f_CB[12] - v_C2*v_vars23)"
|
||||
|
||||
variable dds4 vector "v_ds4*v_ds4"
|
||||
variable dds5 vector "v_ds5*v_ds5"
|
||||
variable dds6 vector "v_ds6*v_ds6"
|
||||
variable vars4 equal "ave(v_dds4)"
|
||||
variable vars5 equal "ave(v_dds5)"
|
||||
variable vars6 equal "ave(v_dds6)"
|
||||
variable C44 equal "v_Pa2GPa*(v_C1*f_CB[4] - v_C2*v_vars4 + v_Rhokbt)"
|
||||
variable C55 equal "v_Pa2GPa*(v_C1*f_CB[5] - v_C2*v_vars5 + v_Rhokbt)"
|
||||
variable C66 equal "v_Pa2GPa*(v_C1*f_CB[6] - v_C2*v_vars6 + v_Rhokbt)"
|
||||
|
||||
variable aC11 equal "(v_C11 + v_C22 + v_C33)/3."
|
||||
variable aC12 equal "(v_C12 + v_C13 + v_C23)/3."
|
||||
variable aC44 equal "(v_C44 + v_C55 + v_C66)/3."
|
||||
|
||||
print """
|
||||
C11 = ${aC11}
|
||||
C12 = ${aC12}
|
||||
C44 = ${aC44}
|
||||
"""
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
1.10120 0.69454 0.66285 0.01885 0.01902 -0.01367
|
||||
0.69454 1.20617 0.65567 0.00893 0.00839 -0.00873
|
||||
0.66285 0.65567 1.11090 -0.01593 0.01634 -0.00483
|
||||
0.01885 0.00893 -0.01593 0.42772 0.01316 0.01666
|
||||
0.01902 0.00839 0.01634 0.01316 0.52416 -0.01733
|
||||
-0.01367 -0.00873 -0.00483 0.01666 -0.01733 0.49891
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
0.04187 0.00000 0.00000 0.00000 0.00000 0.00000
|
||||
0.00000 0.04187 0.00000 0.00000 0.00000 0.00000
|
||||
0.00000 0.00000 0.04187 0.00000 0.00000 0.00000
|
||||
0.00000 0.00000 0.00000 0.02093 0.00000 0.00000
|
||||
0.00000 0.00000 0.00000 0.00000 0.02093 0.00000
|
||||
0.00000 0.00000 0.00000 0.00000 0.00000 0.02093
|
||||
13
examples/ELASTIC_T/BORN_MATRIX/Argon/README.md
Normal file
13
examples/ELASTIC_T/BORN_MATRIX/Argon/README.md
Normal file
@ -0,0 +1,13 @@
|
||||
This repository is a test case for the compute born/matrix. It provides short
|
||||
scripts creating argon fcc crystal and computing the Born term using the two
|
||||
methods described in the documentation.
|
||||
|
||||
In the __Analytical__ directory the terms are computed using the analytical
|
||||
derivation of the Born term for the lj/cut pair style.
|
||||
|
||||
In the __Numdiff__ directory, the Born term is evaluated through small
|
||||
numerical differences of the stress tensor. This method can be used with any
|
||||
interaction potential.
|
||||
|
||||
Both script show examples on how to compute the full Cij elastic stiffness
|
||||
tensor in LAMMPS.
|
||||
1
examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw
Symbolic link
1
examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw
Symbolic link
@ -0,0 +1 @@
|
||||
../../../../potentials/Si.sw
|
||||
68
examples/ELASTIC_T/BORN_MATRIX/Silicon/final_output.in
Normal file
68
examples/ELASTIC_T/BORN_MATRIX/Silicon/final_output.in
Normal file
@ -0,0 +1,68 @@
|
||||
# Average moduli for cubic crystals
|
||||
|
||||
variable C11cubic equal (${C11}+${C22}+${C33})/3.0
|
||||
variable C12cubic equal (${C12}+${C13}+${C23})/3.0
|
||||
variable C44cubic equal (${C44}+${C55}+${C66})/3.0
|
||||
|
||||
variable bulkmodulus equal (${C11cubic}+2*${C12cubic})/3.0
|
||||
variable shearmodulus1 equal ${C44cubic}
|
||||
variable shearmodulus2 equal (${C11cubic}-${C12cubic})/2.0
|
||||
variable poissonratio equal 1.0/(1.0+${C11cubic}/${C12cubic})
|
||||
|
||||
# For Stillinger-Weber silicon, the analytical results
|
||||
# are known to be (E. R. Cowley, 1988):
|
||||
# C11 = 151.4 GPa
|
||||
# C12 = 76.4 GPa
|
||||
# C44 = 56.4 GPa
|
||||
|
||||
#print "========================================="
|
||||
#print "Components of the Elastic Constant Tensor"
|
||||
#print "========================================="
|
||||
|
||||
print "Elastic Constant C11 = ${C11} ${cunits}"
|
||||
print "Elastic Constant C22 = ${C22} ${cunits}"
|
||||
print "Elastic Constant C33 = ${C33} ${cunits}"
|
||||
|
||||
print "Elastic Constant C12 = ${C12} ${cunits}"
|
||||
print "Elastic Constant C13 = ${C13} ${cunits}"
|
||||
print "Elastic Constant C23 = ${C23} ${cunits}"
|
||||
|
||||
print "Elastic Constant C44 = ${C44} ${cunits}"
|
||||
print "Elastic Constant C55 = ${C55} ${cunits}"
|
||||
print "Elastic Constant C66 = ${C66} ${cunits}"
|
||||
|
||||
print "Elastic Constant C14 = ${C14} ${cunits}"
|
||||
print "Elastic Constant C15 = ${C15} ${cunits}"
|
||||
print "Elastic Constant C16 = ${C16} ${cunits}"
|
||||
|
||||
print "Elastic Constant C24 = ${C24} ${cunits}"
|
||||
print "Elastic Constant C25 = ${C25} ${cunits}"
|
||||
print "Elastic Constant C26 = ${C26} ${cunits}"
|
||||
|
||||
print "Elastic Constant C34 = ${C34} ${cunits}"
|
||||
print "Elastic Constant C35 = ${C35} ${cunits}"
|
||||
print "Elastic Constant C36 = ${C36} ${cunits}"
|
||||
|
||||
print "Elastic Constant C45 = ${C45} ${cunits}"
|
||||
print "Elastic Constant C46 = ${C46} ${cunits}"
|
||||
print "Elastic Constant C56 = ${C56} ${cunits}"
|
||||
|
||||
print "========================================="
|
||||
print "Average properties for a cubic crystal"
|
||||
print "========================================="
|
||||
|
||||
print "Bulk Modulus = ${bulkmodulus} ${cunits}"
|
||||
print "Shear Modulus 1 = ${shearmodulus1} ${cunits}"
|
||||
print "Shear Modulus 2 = ${shearmodulus2} ${cunits}"
|
||||
print "Poisson Ratio = ${poissonratio}"
|
||||
|
||||
# summarize sampling protocol
|
||||
|
||||
variable tmp equal atoms
|
||||
print "Number of atoms = ${tmp}"
|
||||
print "Stress sampling interval = ${nevery}"
|
||||
variable tmp equal ${nrun}/${nevery}
|
||||
print "Stress sample count = ${tmp}"
|
||||
print "Born sampling interval = ${neveryborn}"
|
||||
variable tmp equal ${nrun}/${neveryborn}
|
||||
print "Born sample count = ${tmp}"
|
||||
25
examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic
Normal file
25
examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic
Normal file
@ -0,0 +1,25 @@
|
||||
# Compute elastic constant tensor for a crystal at finite temperature
|
||||
# These settings replicate the 1477~K benchmark of
|
||||
# Kluge, Ray, and Rahman (1986) that is Ref.[15] in:
|
||||
# Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
|
||||
|
||||
# here: Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
|
||||
|
||||
include init.in
|
||||
|
||||
# Compute initial state
|
||||
|
||||
include potential.in
|
||||
thermo_style custom step temp pe press density
|
||||
run ${nequil}
|
||||
|
||||
# Run dynamics
|
||||
|
||||
include potential.in
|
||||
include output.in
|
||||
|
||||
run ${nrun}
|
||||
|
||||
# Output final values
|
||||
|
||||
include final_output.in
|
||||
59
examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in
Normal file
59
examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in
Normal file
@ -0,0 +1,59 @@
|
||||
# NOTE: This script can be modified for different atomic structures,
|
||||
# units, etc. See in.elastic for more info.
|
||||
#
|
||||
|
||||
# Define MD parameters
|
||||
# These can be modified by the user
|
||||
# These settings replicate the 1477~K benchmark of
|
||||
# Kluge, Ray, and Rahman (1986) that is Ref.[15] in:
|
||||
# Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
|
||||
|
||||
# select temperature and pressure (lattice constant)
|
||||
|
||||
variable temp index 1477.0 # temperature of initial sample
|
||||
variable a index 5.457 # lattice constant
|
||||
|
||||
# select sampling parameters, important for speed/convergence
|
||||
|
||||
variable nthermo index 1500 # interval for thermo output
|
||||
variable nevery index 10 # stress sampling interval
|
||||
variable neveryborn index 100 # Born sampling interval
|
||||
variable timestep index 0.000766 # timestep
|
||||
variable nlat index 3 # number of lattice unit cells
|
||||
|
||||
# other settings
|
||||
|
||||
variable mass1 index 28.06 # mass
|
||||
variable tdamp index 0.01 # time constant for thermostat
|
||||
variable seed index 123457 # seed for thermostat
|
||||
variable thermostat index 1 # 0 if NVE, 1 if NVT
|
||||
variable delta index 1.0e-6 # Born numdiff strain magnitude
|
||||
|
||||
# hard-coded rules-of-thumb for run length, etc.
|
||||
|
||||
variable nfreq equal ${nthermo} # interval for averaging output
|
||||
variable nrepeat equal floor(${nfreq}/${nevery}) # number of samples
|
||||
variable nrepeatborn equal floor(${nfreq}/${neveryborn}) # number of samples
|
||||
variable nequil equal 10*${nthermo} # length of equilibration run
|
||||
variable nrun equal 100*${nthermo} # length of equilibrated run
|
||||
|
||||
# generate the box and atom positions using a diamond lattice
|
||||
|
||||
units metal
|
||||
|
||||
boundary p p p
|
||||
|
||||
# this generates a standard 8-atom cubic cell
|
||||
|
||||
lattice diamond $a
|
||||
region box prism 0 1 0 1 0 1 0 0 0
|
||||
|
||||
# this generates a 2-atom triclinic cell
|
||||
#include tri.in
|
||||
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
mass 1 ${mass1}
|
||||
replicate ${nlat} ${nlat} ${nlat}
|
||||
velocity all create ${temp} 87287
|
||||
|
||||
140
examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in
Normal file
140
examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in
Normal file
@ -0,0 +1,140 @@
|
||||
# Setup output
|
||||
|
||||
# Stress fluctuation term F
|
||||
|
||||
compute stress all pressure thermo_temp
|
||||
variable s1 equal c_stress[1]
|
||||
variable s2 equal c_stress[2]
|
||||
variable s3 equal c_stress[3]
|
||||
variable s4 equal c_stress[6]
|
||||
variable s5 equal c_stress[5]
|
||||
variable s6 equal c_stress[4]
|
||||
|
||||
variable s11 equal v_s1*v_s1
|
||||
variable s22 equal v_s2*v_s2
|
||||
variable s33 equal v_s3*v_s3
|
||||
variable s44 equal v_s4*v_s4
|
||||
variable s55 equal v_s5*v_s5
|
||||
variable s66 equal v_s6*v_s6
|
||||
variable s33 equal v_s3*v_s3
|
||||
variable s12 equal v_s1*v_s2
|
||||
variable s13 equal v_s1*v_s3
|
||||
variable s14 equal v_s1*v_s4
|
||||
variable s15 equal v_s1*v_s5
|
||||
variable s16 equal v_s1*v_s6
|
||||
variable s23 equal v_s2*v_s3
|
||||
variable s24 equal v_s2*v_s4
|
||||
variable s25 equal v_s2*v_s5
|
||||
variable s26 equal v_s2*v_s6
|
||||
variable s34 equal v_s3*v_s4
|
||||
variable s35 equal v_s3*v_s5
|
||||
variable s36 equal v_s3*v_s6
|
||||
variable s45 equal v_s4*v_s5
|
||||
variable s46 equal v_s4*v_s6
|
||||
variable s56 equal v_s5*v_s6
|
||||
|
||||
variable mytemp equal temp
|
||||
variable mypress equal press
|
||||
variable mype equal pe/atoms
|
||||
fix avt all ave/time ${nevery} ${nrepeat} ${nfreq} v_mytemp ave running
|
||||
fix avp all ave/time ${nevery} ${nrepeat} ${nfreq} v_mypress ave running
|
||||
fix avpe all ave/time ${nevery} ${nrepeat} ${nfreq} v_mype ave running
|
||||
fix avs all ave/time ${nevery} ${nrepeat} ${nfreq} v_s1 v_s2 v_s3 v_s4 v_s5 v_s6 ave running
|
||||
fix avssq all ave/time ${nevery} ${nrepeat} ${nfreq} &
|
||||
v_s11 v_s22 v_s33 v_s44 v_s55 v_s66 &
|
||||
v_s12 v_s13 v_s14 v_s15 v_s16 &
|
||||
v_s23 v_s24 v_s25 v_s26 &
|
||||
v_s34 v_s35 v_s36 &
|
||||
v_s45 v_s46 &
|
||||
v_s56 &
|
||||
ave running
|
||||
|
||||
# bar to GPa
|
||||
variable pconv equal 1.0e5/1.0e9
|
||||
variable cunits index GPa
|
||||
# metal unit constants from LAMMPS
|
||||
# force->nktv2p = 1.6021765e6;
|
||||
# force->boltz = 8.617343e-5;
|
||||
variable boltz equal 8.617343e-5
|
||||
variable nktv2p equal 1.6021765e6
|
||||
variable vkt equal vol/(${boltz}*${temp})/${nktv2p}
|
||||
variable ffac equal ${pconv}*${vkt}
|
||||
|
||||
variable F11 equal -(f_avssq[1]-f_avs[1]*f_avs[1])*${ffac}
|
||||
variable F22 equal -(f_avssq[2]-f_avs[2]*f_avs[2])*${ffac}
|
||||
variable F33 equal -(f_avssq[3]-f_avs[3]*f_avs[3])*${ffac}
|
||||
variable F44 equal -(f_avssq[4]-f_avs[4]*f_avs[4])*${ffac}
|
||||
variable F55 equal -(f_avssq[5]-f_avs[5]*f_avs[5])*${ffac}
|
||||
variable F66 equal -(f_avssq[6]-f_avs[6]*f_avs[6])*${ffac}
|
||||
|
||||
variable F12 equal -(f_avssq[7]-f_avs[1]*f_avs[2])*${ffac}
|
||||
variable F13 equal -(f_avssq[8]-f_avs[1]*f_avs[3])*${ffac}
|
||||
variable F14 equal -(f_avssq[9]-f_avs[1]*f_avs[4])*${ffac}
|
||||
variable F15 equal -(f_avssq[10]-f_avs[1]*f_avs[5])*${ffac}
|
||||
variable F16 equal -(f_avssq[11]-f_avs[1]*f_avs[6])*${ffac}
|
||||
|
||||
variable F23 equal -(f_avssq[12]-f_avs[2]*f_avs[3])*${ffac}
|
||||
variable F24 equal -(f_avssq[13]-f_avs[2]*f_avs[4])*${ffac}
|
||||
variable F25 equal -(f_avssq[14]-f_avs[2]*f_avs[5])*${ffac}
|
||||
variable F26 equal -(f_avssq[15]-f_avs[2]*f_avs[6])*${ffac}
|
||||
|
||||
variable F34 equal -(f_avssq[16]-f_avs[3]*f_avs[4])*${ffac}
|
||||
variable F35 equal -(f_avssq[17]-f_avs[3]*f_avs[5])*${ffac}
|
||||
variable F36 equal -(f_avssq[18]-f_avs[3]*f_avs[6])*${ffac}
|
||||
|
||||
variable F45 equal -(f_avssq[19]-f_avs[4]*f_avs[5])*${ffac}
|
||||
variable F46 equal -(f_avssq[20]-f_avs[4]*f_avs[6])*${ffac}
|
||||
|
||||
variable F56 equal -(f_avssq[21]-f_avs[5]*f_avs[6])*${ffac}
|
||||
|
||||
# Born term
|
||||
|
||||
compute virial all pressure NULL virial
|
||||
compute born all born/matrix numdiff ${delta} virial
|
||||
fix avborn all ave/time ${neveryborn} ${nrepeatborn} ${nfreq} c_born[*] ave running
|
||||
|
||||
variable bfac equal ${pconv}*${nktv2p}/vol
|
||||
variable B vector f_avborn*${bfac}
|
||||
|
||||
# Kinetic term
|
||||
|
||||
variable kfac equal ${pconv}*${nktv2p}*atoms*${boltz}*${temp}/vol
|
||||
variable K11 equal 4.0*${kfac}
|
||||
variable K22 equal 4.0*${kfac}
|
||||
variable K33 equal 4.0*${kfac}
|
||||
variable K44 equal 2.0*${kfac}
|
||||
variable K55 equal 2.0*${kfac}
|
||||
variable K66 equal 2.0*${kfac}
|
||||
|
||||
# Add F, K, and B together
|
||||
|
||||
variable C11 equal v_F11+v_B[1]+v_K11
|
||||
variable C22 equal v_F22+v_B[2]+v_K22
|
||||
variable C33 equal v_F33+v_B[3]+v_K33
|
||||
variable C44 equal v_F44+v_B[4]+v_K44
|
||||
variable C55 equal v_F55+v_B[5]+v_K55
|
||||
variable C66 equal v_F66+v_B[6]+v_K66
|
||||
|
||||
variable C12 equal v_F12+v_B[7]
|
||||
variable C13 equal v_F13+v_B[8]
|
||||
variable C14 equal v_F14+v_B[9]
|
||||
variable C15 equal v_F15+v_B[10]
|
||||
variable C16 equal v_F16+v_B[11]
|
||||
|
||||
variable C23 equal v_F23+v_B[12]
|
||||
variable C24 equal v_F24+v_B[13]
|
||||
variable C25 equal v_F25+v_B[14]
|
||||
variable C26 equal v_F26+v_B[15]
|
||||
|
||||
variable C34 equal v_F34+v_B[16]
|
||||
variable C35 equal v_F35+v_B[17]
|
||||
variable C36 equal v_F36+v_B[18]
|
||||
|
||||
variable C45 equal v_F45+v_B[19]
|
||||
variable C46 equal v_F46+v_B[20]
|
||||
|
||||
variable C56 equal v_F56+v_B[21]
|
||||
|
||||
thermo ${nthermo}
|
||||
thermo_style custom step temp pe press density f_avt f_avp f_avpe v_F11 v_F22 v_F33 v_F44 v_F55 v_F66 v_F12 v_F13 v_F23 v_B[1] v_B[2] v_B[3] v_B[4] v_B[5] v_B[6] v_B[7] v_B[8] v_B[12]
|
||||
thermo_modify norm no
|
||||
21
examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in
Normal file
21
examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in
Normal file
@ -0,0 +1,21 @@
|
||||
# NOTE: This script can be modified for different pair styles
|
||||
# See in.elastic for more info.
|
||||
|
||||
reset_timestep 0
|
||||
|
||||
# Choose potential
|
||||
pair_style sw
|
||||
pair_coeff * * Si.sw Si
|
||||
|
||||
# Setup neighbor style
|
||||
neighbor 1.0 nsq
|
||||
neigh_modify once no every 1 delay 0 check yes
|
||||
|
||||
# Setup MD
|
||||
|
||||
timestep ${timestep}
|
||||
fix 4 all nve
|
||||
if "${thermostat} == 1" then &
|
||||
"fix 5 all langevin ${temp} ${temp} ${tdamp} ${seed}"
|
||||
|
||||
|
||||
26
examples/ELASTIC_T/BORN_MATRIX/Silicon/tri.in
Normal file
26
examples/ELASTIC_T/BORN_MATRIX/Silicon/tri.in
Normal file
@ -0,0 +1,26 @@
|
||||
# this generates a 2-atom triclinic cell
|
||||
# due to rotation on to x-axis,
|
||||
# elastic constant analysis is not working yet
|
||||
|
||||
# unit lattice vectors are
|
||||
# a1 = (1 0 0)
|
||||
# a2 = (1/2 sqrt3/2 0)
|
||||
# a3 = (1/2 1/(2sqrt3) sqrt2/sqrt3)
|
||||
|
||||
variable a1x equal 1
|
||||
variable a2x equal 1/2
|
||||
variable a2y equal sqrt(3)/2
|
||||
variable a3x equal 1/2
|
||||
variable a3y equal 1/(2*sqrt(3))
|
||||
variable a3z equal sqrt(2/3)
|
||||
variable l equal $a/sqrt(2)
|
||||
|
||||
lattice custom ${l} &
|
||||
a1 ${a1x} 0 0 &
|
||||
a2 ${a2x} ${a2y} 0.0 &
|
||||
a3 ${a3x} ${a3y} ${a3z} &
|
||||
basis 0 0 0 &
|
||||
basis 0.25 0.25 0.25 &
|
||||
spacing 1 1 1
|
||||
|
||||
region box prism 0 ${a1x} 0 ${a2y} 0 ${a3z} ${a2x} ${a3x} ${a3y}
|
||||
1
examples/ELASTIC_T/DEFORMATION/Silicon/Si.sw
Symbolic link
1
examples/ELASTIC_T/DEFORMATION/Silicon/Si.sw
Symbolic link
@ -0,0 +1 @@
|
||||
../../../../potentials/Si.sw
|
||||
@ -2,7 +2,7 @@
|
||||
# See in.elastic for more info.
|
||||
|
||||
# we must undefine any fix ave/* fix before using reset_timestep
|
||||
if "$(is_defined(fix,avp)" then "unfix avp"
|
||||
if "$(is_defined(fix,avp))" then "unfix avp"
|
||||
reset_timestep 0
|
||||
|
||||
# Choose potential
|
||||
@ -1,18 +0,0 @@
|
||||
# DATE: 2007-06-11 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Stillinger and Weber, Phys Rev B, 31, 5262, (1985)
|
||||
# Stillinger-Weber parameters for various elements and mixtures
|
||||
# multiple entries can be added to this file, LAMMPS reads the ones it needs
|
||||
# these entries are in LAMMPS "metal" units:
|
||||
# epsilon = eV; sigma = Angstroms
|
||||
# other quantities are unitless
|
||||
|
||||
# format of a single entry (one or more lines):
|
||||
# element 1, element 2, element 3,
|
||||
# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol
|
||||
|
||||
# Here are the original parameters in metal units, for Silicon from:
|
||||
#
|
||||
# Stillinger and Weber, Phys. Rev. B, v. 31, p. 5262, (1985)
|
||||
#
|
||||
|
||||
Si Si Si 2.1683 2.0951 1.80 21.0 1.20 -0.333333333333
|
||||
7.049556277 0.6022245584 4.0 0.0 0.0
|
||||
@ -94,7 +94,7 @@ msst: MSST shock dynamics
|
||||
nb3b: use of nonbonded 3-body harmonic pair style
|
||||
neb: nudged elastic band (NEB) calculation for barrier finding
|
||||
nemd: non-equilibrium MD of 2d sheared system
|
||||
numdiff: numerical difference computation of forces and virial
|
||||
numdiff: numerical difference computation of forces, virial, and Born matrix
|
||||
obstacle: flow around two voids in a 2d channel
|
||||
peptide: dynamics of a small solvated peptide chain (5-mer)
|
||||
peri: Peridynamic model of cylinder impacted by indenter
|
||||
@ -153,12 +153,19 @@ either by itself or in tandem with another code or library. See the
|
||||
COUPLE/README file to get started.
|
||||
|
||||
The ELASTIC directory has an example script for computing elastic
|
||||
constants at zero temperature, using an Si example. See the
|
||||
stiffness tensor (elastic constants)
|
||||
at zero temperature, using an Si example. See the
|
||||
ELASTIC/in.elastic file for more info.
|
||||
|
||||
The ELASTIC_T directory has an example script for computing elastic
|
||||
constants at finite temperature, using an Si example. See the
|
||||
ELASTIC_T/in.elastic file for more info.
|
||||
The ELASTIC_T directory has example scripts for the computing elastic
|
||||
stiffness tensor at finite temperature. Two different methods are
|
||||
demonstrated. DEFORMATION estimates the change in the average
|
||||
stress tensor between multiple simulations
|
||||
in which small finite deformations are made to the simulation cell.
|
||||
BORN_MATRIX runs a single simulation in which the Born matrix and stress
|
||||
fluctuations are averaged. The second method
|
||||
is newer in LAMMPS and is generally more efficient and
|
||||
more reliable.
|
||||
|
||||
The HEAT directory has example scripts for heat exchange algorithms
|
||||
(e.g. used for establishing a thermal gradient), using two different
|
||||
|
||||
@ -1,8 +1,11 @@
|
||||
# LAMMPS FIX NUMDIFF EXAMPLE
|
||||
# LAMMPS NUMDIFF EXAMPLES FOR FORCES, VIRIAL, and BORN MATRIX
|
||||
|
||||
## Numerical Difference Fix
|
||||
## Numerical Difference Fixes and Computes
|
||||
|
||||
This directory contains the ingredients to run an NVE simulation using the numerical difference fix and calculate error in forces.
|
||||
This directory contains the input script for an NVE simulation with
|
||||
fix numdiff, fix numdiff/virial, and compute born/matrix numdiff.
|
||||
In each cases, results are compared to exact analytic expressions
|
||||
and the small relative differences are reported.
|
||||
|
||||
Example:
|
||||
```
|
||||
@ -10,4 +13,4 @@ NP=4 #number of processors
|
||||
mpirun -np $NP lmp_mpi -in.numdiff
|
||||
```
|
||||
|
||||
## Required LAMMPS packages: MOLECULE package
|
||||
## Required LAMMPS packages: EXTRA-FIX, EXTRA-COMPUTE
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
# Numerical difference calculation
|
||||
# of error in forces and virial stress
|
||||
# of error in forces, virial stress, and Born matrix
|
||||
|
||||
# adjustable parameters
|
||||
|
||||
@ -8,8 +8,10 @@ variable nthermo index 10 # thermo output interval
|
||||
variable ndump index 500 # dump output interval
|
||||
variable nlat index 3 # size of box
|
||||
variable fdelta index 1.0e-4 # displacement size
|
||||
variable vdelta index 1.0e-6 # strain size
|
||||
variable vdelta index 1.0e-6 # strain size for numdiff/virial
|
||||
variable bdelta index 1.0e-8 # strain size for numdiff Born matrix
|
||||
variable temp index 10.0 # temperature
|
||||
variable nugget equal 1.0e-6 # regularization for relerr
|
||||
|
||||
units metal
|
||||
atom_style atomic
|
||||
@ -23,8 +25,8 @@ mass 1 39.903
|
||||
|
||||
velocity all create ${temp} 2357 mom yes dist gaussian
|
||||
|
||||
pair_style lj/cubic
|
||||
pair_coeff * * 0.0102701 3.42
|
||||
pair_style lj/cut 5.0
|
||||
pair_coeff 1 1 0.0102701 3.42
|
||||
|
||||
neighbor 0.0 bin
|
||||
neigh_modify every 1 delay 0 check no
|
||||
@ -42,7 +44,7 @@ variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
|
||||
compute faverrsq all reduce ave v_ferrsq
|
||||
variable fsq atom fx^2+fy^2+fz^2
|
||||
compute favsq all reduce ave v_fsq
|
||||
variable frelerr equal sqrt(c_faverrsq/c_favsq)
|
||||
variable frelerr equal sqrt(c_faverrsq/(c_favsq+${nugget}))
|
||||
dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz
|
||||
|
||||
# define numerical virial stress tensor calculation
|
||||
@ -67,8 +69,59 @@ variable vsq equal "c_myvirial[1]^2 + &
|
||||
c_myvirial[4]^2 + &
|
||||
c_myvirial[5]^2 + &
|
||||
c_myvirial[6]^2"
|
||||
variable vrelerr equal sqrt(v_verrsq/v_vsq)
|
||||
variable vrelerr equal sqrt(v_verrsq/(v_vsq+${nugget}))
|
||||
|
||||
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr
|
||||
# define numerical Born matrix calculation
|
||||
|
||||
compute bornnum all born/matrix numdiff ${bdelta} myvirial
|
||||
compute born all born/matrix
|
||||
variable berr vector c_bornnum-c_born
|
||||
variable berrsq equal "v_berr[1]^2 + &
|
||||
v_berr[2]^2 + &
|
||||
v_berr[3]^2 + &
|
||||
v_berr[4]^2 + &
|
||||
v_berr[5]^2 + &
|
||||
v_berr[6]^2 + &
|
||||
v_berr[7]^2 + &
|
||||
v_berr[8]^2 + &
|
||||
v_berr[9]^2 + &
|
||||
v_berr[10]^2 + &
|
||||
v_berr[11]^2 + &
|
||||
v_berr[12]^2 + &
|
||||
v_berr[13]^2 + &
|
||||
v_berr[14]^2 + &
|
||||
v_berr[15]^2 + &
|
||||
v_berr[16]^2 + &
|
||||
v_berr[17]^2 + &
|
||||
v_berr[18]^2 + &
|
||||
v_berr[19]^2 + &
|
||||
v_berr[20]^2 + &
|
||||
v_berr[21]^2"
|
||||
|
||||
variable bsq equal "c_born[1]^2 + &
|
||||
c_born[2]^2 + &
|
||||
c_born[3]^2 + &
|
||||
c_born[4]^2 + &
|
||||
c_born[5]^2 + &
|
||||
c_born[6]^2 + &
|
||||
c_born[7]^2 + &
|
||||
c_born[8]^2 + &
|
||||
c_born[9]^2 + &
|
||||
c_born[10]^2 + &
|
||||
c_born[11]^2 + &
|
||||
c_born[12]^2 + &
|
||||
c_born[13]^2 + &
|
||||
c_born[14]^2 + &
|
||||
c_born[15]^2 + &
|
||||
c_born[16]^2 + &
|
||||
c_born[17]^2 + &
|
||||
c_born[18]^2 + &
|
||||
c_born[19]^2 + &
|
||||
c_born[20]^2 + &
|
||||
c_born[21]^2"
|
||||
|
||||
variable brelerr equal sqrt(v_berrsq/(v_bsq+${nugget}))
|
||||
|
||||
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr v_brelerr
|
||||
thermo ${nthermo}
|
||||
run ${nsteps}
|
||||
|
||||
197
examples/numdiff/log.19Feb2022.log.numdiff.g++.1
Normal file
197
examples/numdiff/log.19Feb2022.log.numdiff.g++.1
Normal file
@ -0,0 +1,197 @@
|
||||
LAMMPS (17 Feb 2022)
|
||||
# Numerical difference calculation
|
||||
# of error in forces, virial stress, and Born matrix
|
||||
|
||||
# adjustable parameters
|
||||
|
||||
variable nsteps index 500 # length of run
|
||||
variable nthermo index 10 # thermo output interval
|
||||
variable ndump index 500 # dump output interval
|
||||
variable nlat index 3 # size of box
|
||||
variable fdelta index 1.0e-4 # displacement size
|
||||
variable vdelta index 1.0e-6 # strain size for numdiff/virial
|
||||
variable bdelta index 1.0e-8 # strain size for numdiff Born matrix
|
||||
variable temp index 10.0 # temperature
|
||||
variable nugget equal 1.0e-6 # regularization for relerr
|
||||
|
||||
units metal
|
||||
atom_style atomic
|
||||
|
||||
atom_modify map yes
|
||||
lattice fcc 5.358000
|
||||
Lattice spacing in x,y,z = 5.358 5.358 5.358
|
||||
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
|
||||
region box block 0 3 0 ${nlat} 0 ${nlat}
|
||||
region box block 0 3 0 3 0 ${nlat}
|
||||
region box block 0 3 0 3 0 3
|
||||
create_box 1 box
|
||||
Created orthogonal box = (0 0 0) to (16.074 16.074 16.074)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 108 atoms
|
||||
using lattice units in orthogonal box = (0 0 0) to (16.074 16.074 16.074)
|
||||
create_atoms CPU = 0.000 seconds
|
||||
mass 1 39.903
|
||||
|
||||
velocity all create ${temp} 2357 mom yes dist gaussian
|
||||
velocity all create 10.0 2357 mom yes dist gaussian
|
||||
|
||||
pair_style lj/cut 5.0
|
||||
pair_coeff 1 1 0.0102701 3.42
|
||||
|
||||
neighbor 0.0 bin
|
||||
neigh_modify every 1 delay 0 check no
|
||||
|
||||
timestep 0.001
|
||||
fix nve all nve
|
||||
|
||||
# define numerical force calculation
|
||||
|
||||
fix numforce all numdiff ${nthermo} ${fdelta}
|
||||
fix numforce all numdiff 10 ${fdelta}
|
||||
fix numforce all numdiff 10 1.0e-4
|
||||
variable ferrx atom f_numforce[1]-fx
|
||||
variable ferry atom f_numforce[2]-fy
|
||||
variable ferrz atom f_numforce[3]-fz
|
||||
variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
|
||||
compute faverrsq all reduce ave v_ferrsq
|
||||
variable fsq atom fx^2+fy^2+fz^2
|
||||
compute favsq all reduce ave v_fsq
|
||||
variable frelerr equal sqrt(c_faverrsq/(c_favsq+${nugget}))
|
||||
variable frelerr equal sqrt(c_faverrsq/(c_favsq+1e-06))
|
||||
dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz
|
||||
dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz
|
||||
|
||||
# define numerical virial stress tensor calculation
|
||||
|
||||
compute myvirial all pressure NULL virial
|
||||
fix numvirial all numdiff/virial ${nthermo} ${vdelta}
|
||||
fix numvirial all numdiff/virial 10 ${vdelta}
|
||||
fix numvirial all numdiff/virial 10 1.0e-6
|
||||
variable errxx equal f_numvirial[1]-c_myvirial[1]
|
||||
variable erryy equal f_numvirial[2]-c_myvirial[2]
|
||||
variable errzz equal f_numvirial[3]-c_myvirial[3]
|
||||
variable erryz equal f_numvirial[4]-c_myvirial[6]
|
||||
variable errxz equal f_numvirial[5]-c_myvirial[5]
|
||||
variable errxy equal f_numvirial[6]-c_myvirial[4]
|
||||
variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2"
|
||||
variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2"
|
||||
variable vrelerr equal sqrt(v_verrsq/(v_vsq+${nugget}))
|
||||
variable vrelerr equal sqrt(v_verrsq/(v_vsq+1e-06))
|
||||
|
||||
# define numerical Born matrix calculation
|
||||
|
||||
compute bornnum all born/matrix numdiff ${bdelta} myvirial
|
||||
compute bornnum all born/matrix numdiff 1.0e-8 myvirial
|
||||
compute born all born/matrix
|
||||
variable berr vector c_bornnum-c_born
|
||||
variable berrsq equal "v_berr[1]^2 + v_berr[2]^2 + v_berr[3]^2 + v_berr[4]^2 + v_berr[5]^2 + v_berr[6]^2 + v_berr[7]^2 + v_berr[8]^2 + v_berr[9]^2 + v_berr[10]^2 + v_berr[11]^2 + v_berr[12]^2 + v_berr[13]^2 + v_berr[14]^2 + v_berr[15]^2 + v_berr[16]^2 + v_berr[17]^2 + v_berr[18]^2 + v_berr[19]^2 + v_berr[20]^2 + v_berr[21]^2"
|
||||
|
||||
variable bsq equal "c_born[1]^2 + c_born[2]^2 + c_born[3]^2 + c_born[4]^2 + c_born[5]^2 + c_born[6]^2 + c_born[7]^2 + c_born[8]^2 + c_born[9]^2 + c_born[10]^2 + c_born[11]^2 + c_born[12]^2 + c_born[13]^2 + c_born[14]^2 + c_born[15]^2 + c_born[16]^2 + c_born[17]^2 + c_born[18]^2 + c_born[19]^2 + c_born[20]^2 + c_born[21]^2"
|
||||
|
||||
variable brelerr equal sqrt(v_berrsq/(v_bsq+${nugget}))
|
||||
variable brelerr equal sqrt(v_berrsq/(v_bsq+1e-06))
|
||||
|
||||
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr v_brelerr
|
||||
thermo ${nthermo}
|
||||
thermo 10
|
||||
run ${nsteps}
|
||||
run 500
|
||||
generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 0 steps, check no
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 5
|
||||
ghost atom cutoff = 5
|
||||
binsize = 2.5, bins = 7 7 7
|
||||
2 neighbor lists, perpetual/occasional/extra = 1 1 0
|
||||
(1) pair lj/cut, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d
|
||||
bin: standard
|
||||
(2) compute born/matrix, occasional, copy from (1)
|
||||
attributes: half, newton on
|
||||
pair build: copy
|
||||
stencil: none
|
||||
bin: none
|
||||
Per MPI rank memory allocation (min/avg/max) = 6.823 | 6.823 | 6.823 Mbytes
|
||||
Step Temp PotEng TotEng Press v_frelerr v_vrelerr v_brelerr
|
||||
0 10 -6.6101864 -6.471878 947.70558 5.7012262e-09 1.4849734e-08 9.036398e-09
|
||||
10 9.9357441 -6.6092976 -6.471878 949.3486 1.3828998e-08 1.9248385e-09 4.0233493e-09
|
||||
20 9.7444561 -6.6066519 -6.4718779 954.23637 1.385204e-08 1.7476399e-09 4.0061081e-09
|
||||
30 9.4311148 -6.6023181 -6.4718779 962.23331 1.4147226e-08 1.7647816e-09 3.3866543e-09
|
||||
40 9.0043293 -6.5964152 -6.4718778 973.10762 1.4128155e-08 1.6390138e-09 3.2652821e-09
|
||||
50 8.4762135 -6.5891108 -6.4718777 986.53572 1.4168048e-08 2.3910821e-09 4.7266627e-09
|
||||
60 7.8621735 -6.5806179 -6.4718775 1002.1092 1.411958e-08 2.0683414e-09 2.6951001e-09
|
||||
70 7.1805874 -6.5711908 -6.4718773 1019.3448 1.4139911e-08 1.6084571e-09 3.1477301e-09
|
||||
80 6.4523557 -6.5611186 -6.4718771 1037.6974 1.4105096e-08 1.9929271e-09 3.4733802e-09
|
||||
90 5.7003071 -6.5507169 -6.4718769 1056.5767 1.4084183e-08 1.750579e-09 4.310104e-09
|
||||
100 4.9484503 -6.5403179 -6.4718767 1075.3674 1.4063796e-08 1.0250271e-09 2.9213594e-09
|
||||
110 4.221081 -6.5302576 -6.4718765 1093.4526 1.400901e-08 1.389277e-09 4.3909721e-09
|
||||
120 3.5417733 -6.520862 -6.4718763 1110.2388 1.4038158e-08 8.6231891e-10 2.5890696e-09
|
||||
130 2.9323072 -6.5124324 -6.4718762 1125.183 1.4048645e-08 7.0840985e-10 3.388192e-09
|
||||
140 2.411607 -6.5052306 -6.471876 1137.8182 1.3968429e-08 1.8508015e-09 3.2976031e-09
|
||||
150 1.9947801 -6.4994654 -6.4718759 1147.7764 1.395965e-08 1.9484728e-09 4.2924605e-09
|
||||
160 1.6923481 -6.4952825 -6.4718759 1154.8063 1.3948606e-08 1.5275137e-09 4.0204309e-09
|
||||
170 1.5097515 -6.492757 -6.4718759 1158.7853 1.3845523e-08 1.5455e-09 4.8781309e-09
|
||||
180 1.4471795 -6.4918916 -6.4718759 1159.7221 1.3788451e-08 1.578099e-09 3.0795316e-09
|
||||
190 1.4997431 -6.4926187 -6.471876 1157.7529 1.374841e-08 2.142073e-09 2.4376961e-09
|
||||
200 1.6579637 -6.4948072 -6.4718761 1153.1286 1.3674788e-08 2.111894e-09 3.7055708e-09
|
||||
210 1.908522 -6.4982727 -6.4718763 1146.1965 1.3639408e-08 1.2386489e-09 3.160881e-09
|
||||
220 2.23518 -6.5027908 -6.4718764 1137.3775 1.3524209e-08 1.7016573e-09 3.6982265e-09
|
||||
230 2.6197892 -6.5081105 -6.4718766 1127.1415 1.3344007e-08 1.5843477e-09 3.7272821e-09
|
||||
240 3.043298 -6.5139681 -6.4718768 1115.9815 1.3245227e-08 1.5502368e-09 3.898015e-09
|
||||
250 3.4866901 -6.5201007 -6.4718769 1104.3906 1.3080142e-08 1.369987e-09 4.9133863e-09
|
||||
260 3.9318061 -6.5262572 -6.4718771 1092.84 1.2885339e-08 1.0743728e-09 5.7271364e-09
|
||||
270 4.3620216 -6.5322076 -6.4718772 1081.7617 1.2705966e-08 1.3618619e-09 2.3225062e-09
|
||||
280 4.7627723 -6.5377504 -6.4718773 1071.5341 1.2480463e-08 1.4346869e-09 3.281167e-09
|
||||
290 5.1219322 -6.542718 -6.4718774 1062.4716 1.2434727e-08 2.1935942e-09 2.8198924e-09
|
||||
300 5.4300557 -6.5469796 -6.4718774 1054.8177 1.2321314e-08 8.2365886e-10 3.2731015e-09
|
||||
310 5.6804997 -6.5504435 -6.4718774 1048.7409 1.2300884e-08 1.4855741e-09 4.1031988e-09
|
||||
320 5.8694423 -6.5530567 -6.4718774 1044.3341 1.2483087e-08 1.8711589e-09 3.9368436e-09
|
||||
330 5.9958115 -6.5548045 -6.4718774 1041.6165 1.2627617e-08 1.9256986e-09 4.3283764e-09
|
||||
340 6.0611353 -6.555708 -6.4718774 1040.5369 1.2935701e-08 1.6609255e-09 3.8728039e-09
|
||||
350 6.0693222 -6.5558211 -6.4718773 1040.9803 1.3218179e-08 1.985355e-09 2.618577e-09
|
||||
360 6.0263776 -6.5552271 -6.4718773 1042.7755 1.3471701e-08 1.5125203e-09 2.936238e-09
|
||||
370 5.9400629 -6.5540332 -6.4718772 1045.7049 1.3676495e-08 1.7364093e-09 2.9097362e-09
|
||||
380 5.8195019 -6.5523657 -6.4718771 1049.515 1.3859995e-08 1.6834835e-09 2.7416302e-09
|
||||
390 5.6747442 -6.5503635 -6.471877 1053.9288 1.3987553e-08 1.7893896e-09 2.8552537e-09
|
||||
400 5.5162948 -6.5481719 -6.4718769 1058.6583 1.4091878e-08 1.4468098e-09 3.2733654e-09
|
||||
410 5.3546269 -6.5459358 -6.4718768 1063.4182 1.4188438e-08 1.7231047e-09 3.3165187e-09
|
||||
420 5.1996958 -6.5437929 -6.4718768 1067.9384 1.4205207e-08 1.3551982e-09 3.8687611e-09
|
||||
430 5.0604771 -6.5418673 -6.4718767 1071.9767 1.4267199e-08 1.361845e-09 3.1210672e-09
|
||||
440 4.9445529 -6.5402639 -6.4718766 1075.3292 1.4253464e-08 1.3945282e-09 2.6483572e-09
|
||||
450 4.8577717 -6.5390637 -6.4718766 1077.8394 1.4240998e-08 1.8767323e-09 3.2040422e-09
|
||||
460 4.8040023 -6.53832 -6.4718766 1079.4048 1.4242259e-08 1.4785379e-09 3.4402279e-09
|
||||
470 4.7849977 -6.5380571 -6.4718766 1079.9795 1.4227939e-08 1.8623848e-09 4.3634918e-09
|
||||
480 4.8003794 -6.5382699 -6.4718766 1079.5756 1.4215836e-08 1.2821795e-09 2.6846581e-09
|
||||
490 4.8477405 -6.538925 -6.4718767 1078.2596 1.4186541e-08 2.47604e-09 3.2044632e-09
|
||||
500 4.9228588 -6.539964 -6.4718767 1076.1469 1.4099819e-08 1.6653302e-09 3.267113e-09
|
||||
Loop time of 0.458483 on 1 procs for 500 steps with 108 atoms
|
||||
|
||||
Performance: 94.224 ns/day, 0.255 hours/ns, 1090.552 timesteps/s
|
||||
99.7% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.0042278 | 0.0042278 | 0.0042278 | 0.0 | 0.92
|
||||
Neigh | 0.02481 | 0.02481 | 0.02481 | 0.0 | 5.41
|
||||
Comm | 0.002944 | 0.002944 | 0.002944 | 0.0 | 0.64
|
||||
Output | 0.014731 | 0.014731 | 0.014731 | 0.0 | 3.21
|
||||
Modify | 0.41122 | 0.41122 | 0.41122 | 0.0 | 89.69
|
||||
Other | | 0.0005545 | | | 0.12
|
||||
|
||||
Nlocal: 108 ave 108 max 108 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 256 ave 256 max 256 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 648 ave 648 max 648 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 648
|
||||
Ave neighs/atom = 6
|
||||
Neighbor list builds = 500
|
||||
Dangerous builds not checked
|
||||
Total wall time: 0:00:00
|
||||
197
examples/numdiff/log.19Feb2022.log.numdiff.g++.4
Normal file
197
examples/numdiff/log.19Feb2022.log.numdiff.g++.4
Normal file
@ -0,0 +1,197 @@
|
||||
LAMMPS (17 Feb 2022)
|
||||
# Numerical difference calculation
|
||||
# of error in forces, virial stress, and Born matrix
|
||||
|
||||
# adjustable parameters
|
||||
|
||||
variable nsteps index 500 # length of run
|
||||
variable nthermo index 10 # thermo output interval
|
||||
variable ndump index 500 # dump output interval
|
||||
variable nlat index 3 # size of box
|
||||
variable fdelta index 1.0e-4 # displacement size
|
||||
variable vdelta index 1.0e-6 # strain size for numdiff/virial
|
||||
variable bdelta index 1.0e-8 # strain size for numdiff Born matrix
|
||||
variable temp index 10.0 # temperature
|
||||
variable nugget equal 1.0e-6 # regularization for relerr
|
||||
|
||||
units metal
|
||||
atom_style atomic
|
||||
|
||||
atom_modify map yes
|
||||
lattice fcc 5.358000
|
||||
Lattice spacing in x,y,z = 5.358 5.358 5.358
|
||||
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
|
||||
region box block 0 3 0 ${nlat} 0 ${nlat}
|
||||
region box block 0 3 0 3 0 ${nlat}
|
||||
region box block 0 3 0 3 0 3
|
||||
create_box 1 box
|
||||
Created orthogonal box = (0 0 0) to (16.074 16.074 16.074)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 108 atoms
|
||||
using lattice units in orthogonal box = (0 0 0) to (16.074 16.074 16.074)
|
||||
create_atoms CPU = 0.000 seconds
|
||||
mass 1 39.903
|
||||
|
||||
velocity all create ${temp} 2357 mom yes dist gaussian
|
||||
velocity all create 10.0 2357 mom yes dist gaussian
|
||||
|
||||
pair_style lj/cut 5.0
|
||||
pair_coeff 1 1 0.0102701 3.42
|
||||
|
||||
neighbor 0.0 bin
|
||||
neigh_modify every 1 delay 0 check no
|
||||
|
||||
timestep 0.001
|
||||
fix nve all nve
|
||||
|
||||
# define numerical force calculation
|
||||
|
||||
fix numforce all numdiff ${nthermo} ${fdelta}
|
||||
fix numforce all numdiff 10 ${fdelta}
|
||||
fix numforce all numdiff 10 1.0e-4
|
||||
variable ferrx atom f_numforce[1]-fx
|
||||
variable ferry atom f_numforce[2]-fy
|
||||
variable ferrz atom f_numforce[3]-fz
|
||||
variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
|
||||
compute faverrsq all reduce ave v_ferrsq
|
||||
variable fsq atom fx^2+fy^2+fz^2
|
||||
compute favsq all reduce ave v_fsq
|
||||
variable frelerr equal sqrt(c_faverrsq/(c_favsq+${nugget}))
|
||||
variable frelerr equal sqrt(c_faverrsq/(c_favsq+1e-06))
|
||||
dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz
|
||||
dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz
|
||||
|
||||
# define numerical virial stress tensor calculation
|
||||
|
||||
compute myvirial all pressure NULL virial
|
||||
fix numvirial all numdiff/virial ${nthermo} ${vdelta}
|
||||
fix numvirial all numdiff/virial 10 ${vdelta}
|
||||
fix numvirial all numdiff/virial 10 1.0e-6
|
||||
variable errxx equal f_numvirial[1]-c_myvirial[1]
|
||||
variable erryy equal f_numvirial[2]-c_myvirial[2]
|
||||
variable errzz equal f_numvirial[3]-c_myvirial[3]
|
||||
variable erryz equal f_numvirial[4]-c_myvirial[6]
|
||||
variable errxz equal f_numvirial[5]-c_myvirial[5]
|
||||
variable errxy equal f_numvirial[6]-c_myvirial[4]
|
||||
variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2"
|
||||
variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2"
|
||||
variable vrelerr equal sqrt(v_verrsq/(v_vsq+${nugget}))
|
||||
variable vrelerr equal sqrt(v_verrsq/(v_vsq+1e-06))
|
||||
|
||||
# define numerical Born matrix calculation
|
||||
|
||||
compute bornnum all born/matrix numdiff ${bdelta} myvirial
|
||||
compute bornnum all born/matrix numdiff 1.0e-8 myvirial
|
||||
compute born all born/matrix
|
||||
variable berr vector c_bornnum-c_born
|
||||
variable berrsq equal "v_berr[1]^2 + v_berr[2]^2 + v_berr[3]^2 + v_berr[4]^2 + v_berr[5]^2 + v_berr[6]^2 + v_berr[7]^2 + v_berr[8]^2 + v_berr[9]^2 + v_berr[10]^2 + v_berr[11]^2 + v_berr[12]^2 + v_berr[13]^2 + v_berr[14]^2 + v_berr[15]^2 + v_berr[16]^2 + v_berr[17]^2 + v_berr[18]^2 + v_berr[19]^2 + v_berr[20]^2 + v_berr[21]^2"
|
||||
|
||||
variable bsq equal "c_born[1]^2 + c_born[2]^2 + c_born[3]^2 + c_born[4]^2 + c_born[5]^2 + c_born[6]^2 + c_born[7]^2 + c_born[8]^2 + c_born[9]^2 + c_born[10]^2 + c_born[11]^2 + c_born[12]^2 + c_born[13]^2 + c_born[14]^2 + c_born[15]^2 + c_born[16]^2 + c_born[17]^2 + c_born[18]^2 + c_born[19]^2 + c_born[20]^2 + c_born[21]^2"
|
||||
|
||||
variable brelerr equal sqrt(v_berrsq/(v_bsq+${nugget}))
|
||||
variable brelerr equal sqrt(v_berrsq/(v_bsq+1e-06))
|
||||
|
||||
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr v_brelerr
|
||||
thermo ${nthermo}
|
||||
thermo 10
|
||||
run ${nsteps}
|
||||
run 500
|
||||
generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 0 steps, check no
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 5
|
||||
ghost atom cutoff = 5
|
||||
binsize = 2.5, bins = 7 7 7
|
||||
2 neighbor lists, perpetual/occasional/extra = 1 1 0
|
||||
(1) pair lj/cut, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d
|
||||
bin: standard
|
||||
(2) compute born/matrix, occasional, copy from (1)
|
||||
attributes: half, newton on
|
||||
pair build: copy
|
||||
stencil: none
|
||||
bin: none
|
||||
Per MPI rank memory allocation (min/avg/max) = 6.816 | 6.816 | 6.816 Mbytes
|
||||
Step Temp PotEng TotEng Press v_frelerr v_vrelerr v_brelerr
|
||||
0 10 -6.6101864 -6.471878 947.70558 1.9110624e-09 9.4407596e-10 3.1867416e-09
|
||||
10 9.9369961 -6.6093149 -6.471878 949.31222 1.3055176e-08 4.996456e-10 2.7421655e-09
|
||||
20 9.7500224 -6.6067289 -6.4718779 954.07898 1.3721178e-08 5.6039795e-10 2.3689718e-09
|
||||
30 9.4448115 -6.6025075 -6.4718779 961.85502 1.3813156e-08 6.8451692e-10 1.9844663e-09
|
||||
40 9.0305392 -6.5967776 -6.4718777 972.39819 1.3961749e-08 3.1134064e-10 1.7915052e-09
|
||||
50 8.5196068 -6.5897109 -6.4718776 985.38158 1.3996941e-08 7.0149406e-10 2.002272e-09
|
||||
60 7.9273388 -6.5815192 -6.4718775 1000.4024 1.4000005e-08 3.5766629e-10 2.4944703e-09
|
||||
70 7.2715879 -6.5724494 -6.4718773 1016.9932 1.3996503e-08 6.2731503e-10 1.7010533e-09
|
||||
80 6.5722375 -6.5627766 -6.4718771 1034.6361 1.3973603e-08 3.1142917e-10 2.808524e-09
|
||||
90 5.8505991 -6.5527956 -6.4718769 1052.7794 1.3983301e-08 3.9931135e-10 2.6118214e-09
|
||||
100 5.128708 -6.542811 -6.4718767 1070.8561 1.395586e-08 2.3152413e-10 2.8742755e-09
|
||||
110 4.4285344 -6.5331269 -6.4718766 1088.305 1.3938374e-08 4.2173005e-10 2.3059886e-09
|
||||
120 3.7711361 -6.5240343 -6.4718764 1104.5919 1.3915264e-08 2.5458038e-10 1.4864012e-09
|
||||
130 3.1757964 -6.5158002 -6.4718762 1119.2319 1.3858843e-08 5.7490448e-10 2.6191823e-09
|
||||
140 2.6591997 -6.5086551 -6.4718761 1131.8095 1.3814891e-08 3.5434633e-10 2.2009364e-09
|
||||
150 2.2347034 -6.5027839 -6.471876 1141.9961 1.3781115e-08 5.0639594e-10 2.9032558e-09
|
||||
160 1.9117661 -6.4983173 -6.471876 1149.564 1.3734288e-08 3.1954962e-10 2.6097446e-09
|
||||
170 1.6955808 -6.4953273 -6.471876 1154.3946 1.3682252e-08 3.5426781e-10 2.9605676e-09
|
||||
180 1.586949 -6.4938249 -6.471876 1156.4812 1.363e-08 4.0804881e-10 2.1707904e-09
|
||||
190 1.5824056 -6.4937621 -6.4718761 1155.925 1.3532637e-08 4.0767685e-10 3.0091462e-09
|
||||
200 1.6745831 -6.4950371 -6.4718762 1152.926 1.3455927e-08 2.953369e-10 2.5029298e-09
|
||||
210 1.8527803 -6.4975018 -6.4718763 1147.7684 1.335224e-08 3.5042319e-10 3.0550064e-09
|
||||
220 2.1036825 -6.5009721 -6.4718764 1140.8026 1.3239176e-08 3.5988448e-10 2.6852683e-09
|
||||
230 2.4121721 -6.5052389 -6.4718766 1132.4243 1.3090019e-08 3.5004036e-10 2.8838602e-09
|
||||
240 2.7621668 -6.5100798 -6.4718767 1123.0538 1.2946525e-08 4.1216361e-10 2.1105916e-09
|
||||
250 3.1374274 -6.5152701 -6.4718768 1113.1152 1.277789e-08 5.9848318e-10 2.3087106e-09
|
||||
260 3.5222906 -6.5205932 -6.471877 1103.0171 1.2591089e-08 2.0080182e-10 1.6969069e-09
|
||||
270 3.9022942 -6.5258491 -6.4718771 1093.1369 1.2432232e-08 4.2494727e-10 1.7375594e-09
|
||||
280 4.2646753 -6.5308612 -6.4718772 1083.8072 1.2268238e-08 6.1239266e-10 1.7005135e-09
|
||||
290 4.598736 -6.5354816 -6.4718772 1075.306 1.2181179e-08 4.9338341e-10 2.1326848e-09
|
||||
300 4.896078 -6.5395941 -6.4718773 1067.85 1.2098274e-08 3.4564838e-10 2.4199891e-09
|
||||
310 5.150715 -6.543116 -6.4718773 1061.5918 1.2184958e-08 4.2383299e-10 2.2243759e-09
|
||||
320 5.3590742 -6.5459978 -6.4718773 1056.6189 1.2312948e-08 3.5194185e-10 1.3856935e-09
|
||||
330 5.5199009 -6.5482222 -6.4718773 1052.9565 1.2573918e-08 4.2401322e-10 2.9882e-09
|
||||
340 5.6340787 -6.5498013 -6.4718773 1050.5719 1.2821551e-08 5.8802825e-10 2.7333289e-09
|
||||
350 5.7043792 -6.5507736 -6.4718772 1049.3813 1.3067314e-08 4.0014945e-10 2.3564728e-09
|
||||
360 5.7351548 -6.5511992 -6.4718772 1049.2581 1.331283e-08 4.1684815e-10 1.735621e-09
|
||||
370 5.7319891 -6.5511553 -6.4718771 1050.042 1.354018e-08 3.8495426e-10 2.4460056e-09
|
||||
380 5.7013193 -6.5507311 -6.4718771 1051.5496 1.3734888e-08 3.5333605e-10 2.5174342e-09
|
||||
390 5.6500487 -6.5500219 -6.471877 1053.5847 1.3892287e-08 3.8154957e-10 1.77358e-09
|
||||
400 5.5851679 -6.5491245 -6.471877 1055.9489 1.3988171e-08 5.8769536e-10 1.9262201e-09
|
||||
410 5.5134009 -6.5481319 -6.4718769 1058.4508 1.4088779e-08 3.6754739e-10 2.7586362e-09
|
||||
420 5.4408957 -6.547129 -6.4718769 1060.9152 1.4139924e-08 4.9030281e-10 3.2871245e-09
|
||||
430 5.3729707 -6.5461895 -6.4718768 1063.1898 1.4173041e-08 5.2345074e-10 3.5995984e-09
|
||||
440 5.3139284 -6.5453729 -6.4718768 1065.1506 1.4191516e-08 5.9481094e-10 2.5005297e-09
|
||||
450 5.2669383 -6.5447229 -6.4718768 1066.7054 1.4168424e-08 3.0799668e-10 2.0864191e-09
|
||||
460 5.2339881 -6.5442672 -6.4718768 1067.7958 1.4163444e-08 6.3927736e-10 2.2872669e-09
|
||||
470 5.2158979 -6.544017 -6.4718768 1068.3968 1.413819e-08 5.5108262e-10 4.4334972e-09
|
||||
480 5.2123873 -6.5439685 -6.4718768 1068.5155 1.4083227e-08 3.9249548e-10 2.5568235e-09
|
||||
490 5.2221849 -6.544104 -6.4718768 1068.188 1.4035287e-08 2.1988631e-10 2.1264034e-09
|
||||
500 5.2431716 -6.5443943 -6.4718768 1067.4759 1.3968666e-08 3.9100701e-10 3.290368e-09
|
||||
Loop time of 0.170182 on 4 procs for 500 steps with 108 atoms
|
||||
|
||||
Performance: 253.846 ns/day, 0.095 hours/ns, 2938.035 timesteps/s
|
||||
99.7% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.0012069 | 0.0012994 | 0.0013512 | 0.2 | 0.76
|
||||
Neigh | 0.0048233 | 0.0050244 | 0.0053881 | 0.3 | 2.95
|
||||
Comm | 0.0072462 | 0.0078013 | 0.008175 | 0.4 | 4.58
|
||||
Output | 0.0080632 | 0.0081244 | 0.0082899 | 0.1 | 4.77
|
||||
Modify | 0.1476 | 0.14764 | 0.14768 | 0.0 | 86.75
|
||||
Other | | 0.0002961 | | | 0.17
|
||||
|
||||
Nlocal: 27 ave 31 max 24 min
|
||||
Histogram: 1 0 1 0 1 0 0 0 0 1
|
||||
Nghost: 135 ave 138 max 131 min
|
||||
Histogram: 1 0 0 0 0 1 0 1 0 1
|
||||
Neighs: 162 ave 191 max 148 min
|
||||
Histogram: 1 2 0 0 0 0 0 0 0 1
|
||||
|
||||
Total # of neighbors = 648
|
||||
Ave neighs/atom = 6
|
||||
Neighbor list builds = 500
|
||||
Dangerous builds not checked
|
||||
Total wall time: 0:00:00
|
||||
@ -1,175 +0,0 @@
|
||||
LAMMPS (7 Jan 2022)
|
||||
# Numerical difference calculation
|
||||
# of error in forces and virial stress
|
||||
|
||||
# adjustable parameters
|
||||
|
||||
variable nsteps index 500 # length of run
|
||||
variable nthermo index 10 # thermo output interval
|
||||
variable ndump index 500 # dump output interval
|
||||
variable nlat index 3 # size of box
|
||||
variable fdelta index 1.0e-4 # displacement size
|
||||
variable vdelta index 1.0e-6 # strain size
|
||||
variable temp index 10.0 # temperature
|
||||
|
||||
units metal
|
||||
atom_style atomic
|
||||
|
||||
atom_modify map yes
|
||||
lattice fcc 5.358000
|
||||
Lattice spacing in x,y,z = 5.358 5.358 5.358
|
||||
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
|
||||
region box block 0 3 0 ${nlat} 0 ${nlat}
|
||||
region box block 0 3 0 3 0 ${nlat}
|
||||
region box block 0 3 0 3 0 3
|
||||
create_box 1 box
|
||||
Created orthogonal box = (0 0 0) to (16.074 16.074 16.074)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 108 atoms
|
||||
using lattice units in orthogonal box = (0 0 0) to (16.074 16.074 16.074)
|
||||
create_atoms CPU = 0.000 seconds
|
||||
mass 1 39.903
|
||||
|
||||
velocity all create ${temp} 2357 mom yes dist gaussian
|
||||
velocity all create 10.0 2357 mom yes dist gaussian
|
||||
|
||||
pair_style lj/cubic
|
||||
pair_coeff * * 0.0102701 3.42
|
||||
|
||||
neighbor 0.0 bin
|
||||
neigh_modify every 1 delay 0 check no
|
||||
|
||||
timestep 0.001
|
||||
fix nve all nve
|
||||
|
||||
# define numerical force calculation
|
||||
|
||||
fix numforce all numdiff ${nthermo} ${fdelta}
|
||||
fix numforce all numdiff 10 ${fdelta}
|
||||
fix numforce all numdiff 10 1.0e-4
|
||||
variable ferrx atom f_numforce[1]-fx
|
||||
variable ferry atom f_numforce[2]-fy
|
||||
variable ferrz atom f_numforce[3]-fz
|
||||
variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
|
||||
compute faverrsq all reduce ave v_ferrsq
|
||||
variable fsq atom fx^2+fy^2+fz^2
|
||||
compute favsq all reduce ave v_fsq
|
||||
variable frelerr equal sqrt(c_faverrsq/c_favsq)
|
||||
dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz
|
||||
dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz
|
||||
|
||||
# define numerical virial stress tensor calculation
|
||||
|
||||
compute myvirial all pressure NULL virial
|
||||
fix numvirial all numdiff/virial ${nthermo} ${vdelta}
|
||||
fix numvirial all numdiff/virial 10 ${vdelta}
|
||||
fix numvirial all numdiff/virial 10 1.0e-6
|
||||
variable errxx equal f_numvirial[1]-c_myvirial[1]
|
||||
variable erryy equal f_numvirial[2]-c_myvirial[2]
|
||||
variable errzz equal f_numvirial[3]-c_myvirial[3]
|
||||
variable erryz equal f_numvirial[4]-c_myvirial[6]
|
||||
variable errxz equal f_numvirial[5]-c_myvirial[5]
|
||||
variable errxy equal f_numvirial[6]-c_myvirial[4]
|
||||
variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2"
|
||||
variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2"
|
||||
variable vrelerr equal sqrt(v_verrsq/v_vsq)
|
||||
|
||||
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr
|
||||
thermo ${nthermo}
|
||||
thermo 10
|
||||
run ${nsteps}
|
||||
run 500
|
||||
generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 0 steps, check no
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 5.9407173
|
||||
ghost atom cutoff = 5.9407173
|
||||
binsize = 2.9703587, bins = 6 6 6
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair lj/cubic, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 6.083 | 6.083 | 6.083 Mbytes
|
||||
Step Temp PotEng TotEng Press v_frelerr v_vrelerr
|
||||
0 10 -7.0259569 -6.8876486 28.564278 19203.344 1.5660292e-06
|
||||
10 9.9376583 -7.0250947 -6.8876486 30.254762 1.5040965e-08 2.1991382e-07
|
||||
20 9.7520139 -7.022527 -6.8876485 35.28505 1.4756358e-08 2.6265315e-06
|
||||
30 9.4477557 -7.0183188 -6.8876485 43.519863 1.4688198e-08 2.6356166e-07
|
||||
40 9.0330215 -7.0125826 -6.8876484 54.727797 1.4637921e-08 5.2292327e-08
|
||||
50 8.5192918 -7.0054772 -6.8876483 68.585553 1.4587854e-08 7.1324716e-08
|
||||
60 7.9212026 -6.997205 -6.8876481 84.684636 1.4525561e-08 3.1108149e-08
|
||||
70 7.2562592 -6.9880081 -6.8876479 102.54088 1.450885e-08 3.2311094e-08
|
||||
80 6.5444294 -6.9781627 -6.8876478 121.60715 1.4444738e-08 2.1776998e-08
|
||||
90 5.8075961 -6.9679715 -6.8876476 141.2895 1.4493562e-08 2.0400898e-08
|
||||
100 5.0688629 -6.957754 -6.8876474 160.9668 1.445455e-08 1.2636688e-08
|
||||
110 4.3517145 -6.947835 -6.8876472 180.0135 1.4460371e-08 1.2528038e-08
|
||||
120 3.6790589 -6.9385314 -6.887647 197.82486 1.4371757e-08 1.4489522e-08
|
||||
130 3.0721984 -6.9301379 -6.8876468 213.84331 1.4364708e-08 1.2461922e-08
|
||||
140 2.5497991 -6.9229125 -6.8876467 227.58429 1.4330926e-08 9.3913926e-09
|
||||
150 2.1269443 -6.917064 -6.8876466 238.6596 1.4287002e-08 4.1510266e-09
|
||||
160 1.8143642 -6.9127407 -6.8876465 246.79599 1.4282669e-08 7.7048281e-09
|
||||
170 1.6179191 -6.9100237 -6.8876465 251.84748 1.42726e-08 1.2719973e-08
|
||||
180 1.5383946 -6.9089239 -6.8876466 253.79991 1.4236534e-08 8.1200831e-09
|
||||
190 1.5716287 -6.9093836 -6.8876467 252.76745 1.41706e-08 6.5670612e-09
|
||||
200 1.7089493 -6.911283 -6.8876468 248.98142 1.4096463e-08 1.1685863e-08
|
||||
210 1.9378716 -6.9144493 -6.8876469 242.77289 1.4008978e-08 1.1226902e-08
|
||||
220 2.2429731 -6.9186692 -6.887647 234.55055 1.3886901e-08 9.9914102e-09
|
||||
230 2.606862 -6.9237023 -6.8876472 224.77626 1.3864576e-08 1.1540228e-08
|
||||
240 3.0111524 -6.9292941 -6.8876474 213.93996 1.3696314e-08 1.1697747e-08
|
||||
250 3.4373794 -6.9351893 -6.8876475 202.53583 1.3626701e-08 1.0398197e-08
|
||||
260 3.8678047 -6.9411426 -6.8876476 191.04084 1.3489489e-08 6.6603364e-09
|
||||
270 4.2860853 -6.9469279 -6.8876478 179.89646 1.3312014e-08 1.1687917e-08
|
||||
280 4.6777954 -6.9523457 -6.8876479 169.49404 1.3081144e-08 1.1336675e-08
|
||||
290 5.030805 -6.9572282 -6.887648 160.16371 1.2947385e-08 1.7342825e-08
|
||||
300 5.3355278 -6.9614428 -6.887648 152.16682 1.2893673e-08 1.7510534e-08
|
||||
310 5.5850532 -6.964894 -6.887648 145.69148 1.2842022e-08 1.2782546e-08
|
||||
320 5.7751794 -6.9675236 -6.8876481 140.85102 1.2903488e-08 1.5319437e-08
|
||||
330 5.9043601 -6.9693103 -6.887648 137.68497 1.3076809e-08 1.1208999e-08
|
||||
340 5.9735784 -6.9702676 -6.887648 136.16232 1.3296904e-08 1.891087e-08
|
||||
350 5.9861549 -6.9704415 -6.887648 136.18679 1.3504051e-08 2.5783601e-08
|
||||
360 5.947496 -6.9699067 -6.8876479 137.60397 1.3731112e-08 2.0556839e-08
|
||||
370 5.8647874 -6.9687627 -6.8876478 140.2101 1.4009878e-08 2.1771736e-08
|
||||
380 5.7466376 -6.9671285 -6.8876477 143.76234 1.4092054e-08 1.1085162e-08
|
||||
390 5.6026773 -6.9651374 -6.8876477 147.99019 1.4282872e-08 2.0221602e-08
|
||||
400 5.4431231 -6.9629305 -6.8876476 152.60787 1.4317739e-08 1.7076065e-08
|
||||
410 5.2783192 -6.960651 -6.8876475 157.32722 1.4415075e-08 2.5031776e-08
|
||||
420 5.1182723 -6.9584374 -6.8876474 161.87063 1.4441435e-08 2.2519289e-08
|
||||
430 4.9722 -6.956417 -6.8876473 165.98344 1.4550624e-08 2.4512613e-08
|
||||
440 4.8481153 -6.9547008 -6.8876473 169.44527 1.4544672e-08 1.4758301e-08
|
||||
450 4.7524707 -6.9533779 -6.8876472 172.07964 1.4546492e-08 1.324687e-08
|
||||
460 4.6898817 -6.9525122 -6.8876472 173.76132 1.4537475e-08 1.351367e-08
|
||||
470 4.6629495 -6.9521397 -6.8876472 174.42109 1.4530458e-08 1.521106e-08
|
||||
480 4.6721922 -6.9522675 -6.8876472 174.04742 1.4543785e-08 1.0905422e-08
|
||||
490 4.7160887 -6.9528747 -6.8876473 172.68525 1.4545591e-08 2.0128525e-08
|
||||
500 4.7912313 -6.953914 -6.8876473 170.43183 1.4438981e-08 1.6062775e-08
|
||||
Loop time of 0.837333 on 1 procs for 500 steps with 108 atoms
|
||||
|
||||
Performance: 51.592 ns/day, 0.465 hours/ns, 597.134 timesteps/s
|
||||
99.8% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.0097726 | 0.0097726 | 0.0097726 | 0.0 | 1.17
|
||||
Neigh | 0.03095 | 0.03095 | 0.03095 | 0.0 | 3.70
|
||||
Comm | 0.005564 | 0.005564 | 0.005564 | 0.0 | 0.66
|
||||
Output | 0.0042451 | 0.0042451 | 0.0042451 | 0.0 | 0.51
|
||||
Modify | 0.78618 | 0.78618 | 0.78618 | 0.0 | 93.89
|
||||
Other | | 0.0006258 | | | 0.07
|
||||
|
||||
Nlocal: 108 ave 108 max 108 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 558 ave 558 max 558 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 972 ave 972 max 972 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 972
|
||||
Ave neighs/atom = 9
|
||||
Neighbor list builds = 500
|
||||
Dangerous builds not checked
|
||||
Total wall time: 0:00:00
|
||||
@ -1,175 +0,0 @@
|
||||
LAMMPS (7 Jan 2022)
|
||||
# Numerical difference calculation
|
||||
# of error in forces and virial stress
|
||||
|
||||
# adjustable parameters
|
||||
|
||||
variable nsteps index 500 # length of run
|
||||
variable nthermo index 10 # thermo output interval
|
||||
variable ndump index 500 # dump output interval
|
||||
variable nlat index 3 # size of box
|
||||
variable fdelta index 1.0e-4 # displacement size
|
||||
variable vdelta index 1.0e-6 # strain size
|
||||
variable temp index 10.0 # temperature
|
||||
|
||||
units metal
|
||||
atom_style atomic
|
||||
|
||||
atom_modify map yes
|
||||
lattice fcc 5.358000
|
||||
Lattice spacing in x,y,z = 5.358 5.358 5.358
|
||||
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
|
||||
region box block 0 3 0 ${nlat} 0 ${nlat}
|
||||
region box block 0 3 0 3 0 ${nlat}
|
||||
region box block 0 3 0 3 0 3
|
||||
create_box 1 box
|
||||
Created orthogonal box = (0 0 0) to (16.074 16.074 16.074)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 108 atoms
|
||||
using lattice units in orthogonal box = (0 0 0) to (16.074 16.074 16.074)
|
||||
create_atoms CPU = 0.000 seconds
|
||||
mass 1 39.903
|
||||
|
||||
velocity all create ${temp} 2357 mom yes dist gaussian
|
||||
velocity all create 10.0 2357 mom yes dist gaussian
|
||||
|
||||
pair_style lj/cubic
|
||||
pair_coeff * * 0.0102701 3.42
|
||||
|
||||
neighbor 0.0 bin
|
||||
neigh_modify every 1 delay 0 check no
|
||||
|
||||
timestep 0.001
|
||||
fix nve all nve
|
||||
|
||||
# define numerical force calculation
|
||||
|
||||
fix numforce all numdiff ${nthermo} ${fdelta}
|
||||
fix numforce all numdiff 10 ${fdelta}
|
||||
fix numforce all numdiff 10 1.0e-4
|
||||
variable ferrx atom f_numforce[1]-fx
|
||||
variable ferry atom f_numforce[2]-fy
|
||||
variable ferrz atom f_numforce[3]-fz
|
||||
variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
|
||||
compute faverrsq all reduce ave v_ferrsq
|
||||
variable fsq atom fx^2+fy^2+fz^2
|
||||
compute favsq all reduce ave v_fsq
|
||||
variable frelerr equal sqrt(c_faverrsq/c_favsq)
|
||||
dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz
|
||||
dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz
|
||||
|
||||
# define numerical virial stress tensor calculation
|
||||
|
||||
compute myvirial all pressure NULL virial
|
||||
fix numvirial all numdiff/virial ${nthermo} ${vdelta}
|
||||
fix numvirial all numdiff/virial 10 ${vdelta}
|
||||
fix numvirial all numdiff/virial 10 1.0e-6
|
||||
variable errxx equal f_numvirial[1]-c_myvirial[1]
|
||||
variable erryy equal f_numvirial[2]-c_myvirial[2]
|
||||
variable errzz equal f_numvirial[3]-c_myvirial[3]
|
||||
variable erryz equal f_numvirial[4]-c_myvirial[6]
|
||||
variable errxz equal f_numvirial[5]-c_myvirial[5]
|
||||
variable errxy equal f_numvirial[6]-c_myvirial[4]
|
||||
variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2"
|
||||
variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2"
|
||||
variable vrelerr equal sqrt(v_verrsq/v_vsq)
|
||||
|
||||
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr
|
||||
thermo ${nthermo}
|
||||
thermo 10
|
||||
run ${nsteps}
|
||||
run 500
|
||||
generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 0 steps, check no
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 5.9407173
|
||||
ghost atom cutoff = 5.9407173
|
||||
binsize = 2.9703587, bins = 6 6 6
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair lj/cubic, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 6.067 | 6.067 | 6.067 Mbytes
|
||||
Step Temp PotEng TotEng Press v_frelerr v_vrelerr
|
||||
0 10 -7.0259569 -6.8876486 28.564278 10664.391 9.1481187e-08
|
||||
10 9.9388179 -7.0251107 -6.8876486 30.21736 1.4771865e-08 1.3452512e-07
|
||||
20 9.7572185 -7.022599 -6.8876485 35.123527 1.437525e-08 8.0966999e-07
|
||||
30 9.4606673 -7.0184974 -6.8876484 43.132052 1.4375468e-08 1.990012e-08
|
||||
40 9.0579092 -7.0129268 -6.8876483 54.000927 1.4350331e-08 1.7239028e-08
|
||||
50 8.5607685 -7.0060508 -6.8876482 67.403151 1.4353284e-08 7.813181e-09
|
||||
60 7.9838726 -6.9980717 -6.8876481 82.935358 1.4398078e-08 2.022316e-08
|
||||
70 7.3442875 -6.9892255 -6.8876479 100.12892 1.434409e-08 7.5938179e-09
|
||||
80 6.6610579 -6.9797757 -6.8876477 118.46358 1.4324787e-08 7.1972571e-09
|
||||
90 5.9546462 -6.9700053 -6.8876476 137.38365 1.4322718e-08 4.3978378e-09
|
||||
100 5.2462727 -6.9602077 -6.8876474 156.31651 1.4273172e-08 4.6728038e-09
|
||||
110 4.5571706 -6.9506767 -6.8876472 174.69294 1.4266163e-08 3.522225e-09
|
||||
120 3.9077807 -6.9416949 -6.887647 191.96859 1.42241e-08 3.5357511e-09
|
||||
130 3.3169241 -6.9335227 -6.8876469 207.64566 1.4203813e-08 2.5182488e-09
|
||||
140 2.8010028 -6.926387 -6.8876468 221.29333 1.4164215e-08 2.3112509e-09
|
||||
150 2.3732854 -6.9204712 -6.8876467 232.5658 1.4134122e-08 1.9368963e-09
|
||||
160 2.0433329 -6.9159076 -6.8876466 241.21647 1.4095473e-08 3.6806452e-09
|
||||
170 1.8166184 -6.912772 -6.8876466 247.10715 1.4049531e-08 2.5319322e-09
|
||||
180 1.6943727 -6.9110813 -6.8876467 250.21143 1.3997912e-08 1.952404e-09
|
||||
190 1.6736731 -6.910795 -6.8876467 250.61203 1.3915487e-08 1.4271767e-09
|
||||
200 1.7477659 -6.9118199 -6.8876468 248.49223 1.3850618e-08 2.4515718e-09
|
||||
210 1.9065921 -6.9140167 -6.8876469 244.12226 1.3747916e-08 1.7957531e-09
|
||||
220 2.1374676 -6.91721 -6.887647 237.84173 1.3674779e-08 2.6613511e-09
|
||||
230 2.4258607 -6.9211989 -6.8876472 230.0395 1.3565712e-08 3.0777067e-09
|
||||
240 2.7562034 -6.9257679 -6.8876473 221.13265 1.3440442e-08 1.7111501e-09
|
||||
250 3.1126827 -6.9306984 -6.8876474 211.54566 1.3270914e-08 1.6690112e-09
|
||||
260 3.4799641 -6.9357784 -6.8876476 201.69126 1.3105092e-08 2.1637558e-09
|
||||
270 3.8438148 -6.9408108 -6.8876477 191.95361 1.2962846e-08 4.4613506e-09
|
||||
280 4.191607 -6.9456212 -6.8876478 182.6745 1.2846383e-08 3.3730203e-09
|
||||
290 4.5126922 -6.9500621 -6.8876478 174.14285 1.2710692e-08 2.2809889e-09
|
||||
300 4.7986487 -6.9540172 -6.8876479 166.58747 1.2657778e-08 6.9880891e-09
|
||||
310 5.0434083 -6.9574025 -6.8876479 160.17316 1.266381e-08 4.2486217e-09
|
||||
320 5.243275 -6.9601668 -6.8876479 154.99974 1.279856e-08 5.1505673e-09
|
||||
330 5.3968455 -6.9622908 -6.8876479 151.1038 1.2981831e-08 3.3339333e-09
|
||||
340 5.5048468 -6.9637845 -6.8876479 148.46296 1.3159021e-08 1.7881393e-09
|
||||
350 5.569902 -6.9646843 -6.8876479 147.00205 1.3439465e-08 5.6876219e-09
|
||||
360 5.5962378 -6.9650485 -6.8876478 146.60113 1.3645943e-08 7.233847e-09
|
||||
370 5.5893468 -6.9649531 -6.8876478 147.10471 1.3829581e-08 4.5514318e-09
|
||||
380 5.5556199 -6.9644866 -6.8876477 148.33195 1.4005893e-08 4.2971846e-09
|
||||
390 5.5019639 -6.9637444 -6.8876476 150.08725 1.4157454e-08 3.3564262e-09
|
||||
400 5.4354239 -6.962824 -6.8876476 152.17073 1.4226422e-08 4.2393923e-09
|
||||
410 5.3628267 -6.9618199 -6.8876475 154.38825 1.4302679e-08 3.8937698e-09
|
||||
420 5.2904639 -6.960819 -6.8876475 156.56034 1.4381099e-08 4.315875e-09
|
||||
430 5.2238282 -6.9598973 -6.8876474 158.52969 1.4426567e-08 4.2658185e-09
|
||||
440 5.1674149 -6.9591171 -6.8876474 160.16704 1.4453381e-08 5.7055268e-09
|
||||
450 5.1245913 -6.9585248 -6.8876474 161.37513 1.4449488e-08 4.4308801e-09
|
||||
460 5.0975361 -6.9581506 -6.8876474 162.09077 1.4445596e-08 5.8269923e-09
|
||||
470 5.0872416 -6.9580082 -6.8876474 162.28517 1.4444348e-08 4.8263194e-09
|
||||
480 5.0935712 -6.9580957 -6.8876474 161.96268 1.4411666e-08 6.222228e-09
|
||||
490 5.115362 -6.9583971 -6.8876474 161.15816 1.4369716e-08 3.3926077e-09
|
||||
500 5.1505605 -6.958884 -6.8876474 159.9333 1.4288515e-08 3.8845251e-09
|
||||
Loop time of 0.252598 on 4 procs for 500 steps with 108 atoms
|
||||
|
||||
Performance: 171.023 ns/day, 0.140 hours/ns, 1979.430 timesteps/s
|
||||
99.8% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.0021545 | 0.0022845 | 0.0023794 | 0.2 | 0.90
|
||||
Neigh | 0.0063887 | 0.0067604 | 0.0069752 | 0.3 | 2.68
|
||||
Comm | 0.01048 | 0.010916 | 0.011408 | 0.3 | 4.32
|
||||
Output | 0.0026603 | 0.0027399 | 0.0029738 | 0.3 | 1.08
|
||||
Modify | 0.2295 | 0.22952 | 0.22954 | 0.0 | 90.86
|
||||
Other | | 0.0003814 | | | 0.15
|
||||
|
||||
Nlocal: 27 ave 29 max 25 min
|
||||
Histogram: 1 0 1 0 0 0 0 1 0 1
|
||||
Nghost: 325 ave 327 max 323 min
|
||||
Histogram: 1 0 1 0 0 0 0 1 0 1
|
||||
Neighs: 243 ave 273 max 228 min
|
||||
Histogram: 1 1 1 0 0 0 0 0 0 1
|
||||
|
||||
Total # of neighbors = 972
|
||||
Ave neighs/atom = 9
|
||||
Neighbor list builds = 500
|
||||
Dangerous builds not checked
|
||||
Total wall time: 0:00:00
|
||||
@ -387,7 +387,7 @@ double LammpsInterface::atom_quantity_conversion(FundamentalAtomQuantity quantit
|
||||
|
||||
int LammpsInterface::dimension() const { return lammps_->domain->dimension; }
|
||||
|
||||
int LammpsInterface::nregion() const { return lammps_->domain->nregion; }
|
||||
int LammpsInterface::nregion() const { return lammps_->domain->get_region_list().size(); }
|
||||
|
||||
void LammpsInterface::box_bounds(double & boxxlo, double & boxxhi,
|
||||
double & boxylo, double & boxyhi,
|
||||
@ -527,14 +527,15 @@ void LammpsInterface::box_periodicity(int & xperiodic,
|
||||
zperiodic = lammps_->domain->zperiodic;
|
||||
}
|
||||
|
||||
int LammpsInterface::region_id(const char * regionName) const {
|
||||
int nregion = this->nregion();
|
||||
for (int iregion = 0; iregion < nregion; iregion++) {
|
||||
if (strcmp(regionName, region_name(iregion)) == 0) {
|
||||
int LammpsInterface::region_id(const char *regionName) const {
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
int iregion = 0;
|
||||
for (auto reg : regions) {
|
||||
if (strcmp(regionName, reg->id) == 0) {
|
||||
return iregion;
|
||||
}
|
||||
++iregion;
|
||||
}
|
||||
throw ATC_Error("Region has not been defined");
|
||||
return -1;
|
||||
}
|
||||
|
||||
@ -1322,61 +1323,73 @@ int** LammpsInterface::bond_list() const { return lammps_->neighbor->bondlist;
|
||||
|
||||
char * LammpsInterface::region_name(int iRegion) const
|
||||
{
|
||||
return lammps_->domain->regions[iRegion]->id;
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->id;
|
||||
}
|
||||
|
||||
char * LammpsInterface::region_style(int iRegion) const
|
||||
{
|
||||
return lammps_->domain->regions[iRegion]->style;
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->style;
|
||||
}
|
||||
|
||||
double LammpsInterface::region_xlo(int iRegion) const
|
||||
{
|
||||
return lammps_->domain->regions[iRegion]->extent_xlo;
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->extent_xlo;
|
||||
}
|
||||
|
||||
double LammpsInterface::region_xhi(int iRegion) const
|
||||
{
|
||||
return lammps_->domain->regions[iRegion]->extent_xhi;
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->extent_xhi;
|
||||
}
|
||||
|
||||
double LammpsInterface::region_ylo(int iRegion) const
|
||||
{
|
||||
return lammps_->domain->regions[iRegion]->extent_ylo;
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->extent_ylo;
|
||||
}
|
||||
|
||||
double LammpsInterface::region_yhi(int iRegion) const
|
||||
{
|
||||
return lammps_->domain->regions[iRegion]->extent_yhi;
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->extent_yhi;
|
||||
}
|
||||
|
||||
double LammpsInterface::region_zlo(int iRegion) const
|
||||
{
|
||||
return lammps_->domain->regions[iRegion]->extent_zlo;
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->extent_zlo;
|
||||
}
|
||||
|
||||
double LammpsInterface::region_zhi(int iRegion) const
|
||||
{
|
||||
return lammps_->domain->regions[iRegion]->extent_zhi;
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->extent_zhi;
|
||||
}
|
||||
|
||||
double LammpsInterface::region_xscale(int iRegion) const
|
||||
{
|
||||
return lammps_->domain->regions[iRegion]->xscale;
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->xscale;
|
||||
}
|
||||
|
||||
double LammpsInterface::region_yscale(int iRegion) const
|
||||
{
|
||||
return lammps_->domain->regions[iRegion]->yscale;
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->yscale;
|
||||
}
|
||||
|
||||
double LammpsInterface::region_zscale(int iRegion) const
|
||||
{
|
||||
return lammps_->domain->regions[iRegion]->zscale;
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->zscale;
|
||||
}
|
||||
|
||||
int LammpsInterface::region_match(int iRegion, double x, double y, double z) const {
|
||||
return lammps_->domain->regions[iRegion]->match(x,y,z);
|
||||
auto regions = lammps_->domain->get_region_list();
|
||||
return regions[iRegion]->match(x,y,z);
|
||||
}
|
||||
|
||||
// -----------------------------------------------------------------
|
||||
|
||||
4
src/.gitignore
vendored
4
src/.gitignore
vendored
@ -453,6 +453,8 @@
|
||||
/compute_basal_atom.h
|
||||
/compute_body_local.cpp
|
||||
/compute_body_local.h
|
||||
/compute_born_matrix.cpp
|
||||
/compute_born_matrix.h
|
||||
/compute_cnp_atom.cpp
|
||||
/compute_cnp_atom.h
|
||||
/compute_damage_atom.cpp
|
||||
@ -1091,6 +1093,8 @@
|
||||
/pair_hdnnp.h
|
||||
/pair_ilp_graphene_hbn.cpp
|
||||
/pair_ilp_graphene_hbn.h
|
||||
/pair_ilp_graphene_hbn_opt.cpp
|
||||
/pair_ilp_graphene_hbn_opt.h
|
||||
/pair_ilp_tmd.cpp
|
||||
/pair_ilp_tmd.h
|
||||
/pair_kolmogorov_crespi_full.cpp
|
||||
|
||||
@ -134,10 +134,10 @@ void DumpAtomADIOS::write()
|
||||
|
||||
// Now we know the global size and the local subset size and offset
|
||||
// of the atoms table
|
||||
auto nAtomsGlobal = static_cast<size_t>(ntotal);
|
||||
auto startRow = static_cast<size_t>(atomOffset);
|
||||
auto nAtomsLocal = static_cast<size_t>(nme);
|
||||
auto nColumns = static_cast<size_t>(size_one);
|
||||
auto nAtomsGlobal = static_cast<size_t>(ntotal);
|
||||
auto startRow = static_cast<size_t>(atomOffset);
|
||||
auto nAtomsLocal = static_cast<size_t>(nme);
|
||||
auto nColumns = static_cast<size_t>(size_one);
|
||||
internal->varAtoms.SetShape({nAtomsGlobal, nColumns});
|
||||
internal->varAtoms.SetSelection({{startRow, 0}, {nAtomsLocal, nColumns}});
|
||||
|
||||
@ -238,7 +238,7 @@ void DumpAtomADIOS::init_style()
|
||||
columnNames = {"id", "type", "xs", "ys", "zs", "ix", "iy", "iz"};
|
||||
}
|
||||
|
||||
for (int icol = 0; icol < (int)columnNames.size(); ++icol)
|
||||
for (int icol = 0; icol < (int) columnNames.size(); ++icol)
|
||||
if (keyword_user[icol].size()) columnNames[icol] = keyword_user[icol];
|
||||
|
||||
// setup function ptrs
|
||||
@ -296,7 +296,7 @@ void DumpAtomADIOS::init_style()
|
||||
int *boundaryptr = reinterpret_cast<int *>(domain->boundary);
|
||||
internal->io.DefineAttribute<int>("boundary", boundaryptr, 6);
|
||||
|
||||
auto nColumns = static_cast<size_t>(size_one);
|
||||
auto nColumns = static_cast<size_t>(size_one);
|
||||
internal->io.DefineAttribute<std::string>("columns", columnNames.data(), nColumns);
|
||||
internal->io.DefineAttribute<std::string>("columnstr", columns);
|
||||
internal->io.DefineAttribute<std::string>("boundarystr", boundstr);
|
||||
|
||||
@ -146,10 +146,10 @@ void DumpCustomADIOS::write()
|
||||
|
||||
// Now we know the global size and the local subset size and offset
|
||||
// of the atoms table
|
||||
auto nAtomsGlobal = static_cast<size_t>(ntotal);
|
||||
auto startRow = static_cast<size_t>(atomOffset);
|
||||
auto nAtomsLocal = static_cast<size_t>(nme);
|
||||
auto nColumns = static_cast<size_t>(size_one);
|
||||
auto nAtomsGlobal = static_cast<size_t>(ntotal);
|
||||
auto startRow = static_cast<size_t>(atomOffset);
|
||||
auto nAtomsLocal = static_cast<size_t>(nme);
|
||||
auto nColumns = static_cast<size_t>(size_one);
|
||||
internal->varAtoms.SetShape({nAtomsGlobal, nColumns});
|
||||
internal->varAtoms.SetSelection({{startRow, 0}, {nAtomsLocal, nColumns}});
|
||||
|
||||
@ -221,8 +221,10 @@ void DumpCustomADIOS::init_style()
|
||||
int icol = 0;
|
||||
for (auto item : utils::split_words(columns_default)) {
|
||||
if (combined.size()) combined += " ";
|
||||
if (keyword_user[icol].size()) combined += keyword_user[icol];
|
||||
else combined += item;
|
||||
if (keyword_user[icol].size())
|
||||
combined += keyword_user[icol];
|
||||
else
|
||||
combined += item;
|
||||
++icol;
|
||||
}
|
||||
columns = utils::strdup(combined);
|
||||
@ -249,34 +251,30 @@ void DumpCustomADIOS::init_style()
|
||||
*/
|
||||
// find current ptr for each compute,fix,variable
|
||||
// check that fix frequency is acceptable
|
||||
int icompute;
|
||||
for (int i = 0; i < ncompute; i++) {
|
||||
icompute = modify->find_compute(id_compute[i]);
|
||||
if (icompute < 0) error->all(FLERR, "Could not find dump custom compute ID");
|
||||
compute[i] = modify->compute[icompute];
|
||||
compute[i] = modify->get_compute_by_id(id_compute[i]);
|
||||
if (!compute[i])
|
||||
error->all(FLERR, "Could not find dump custom/adios compute ID {}", id_compute[i]);
|
||||
}
|
||||
|
||||
int ifix;
|
||||
for (int i = 0; i < nfix; i++) {
|
||||
ifix = modify->find_fix(id_fix[i]);
|
||||
if (ifix < 0) error->all(FLERR, "Could not find dump custom fix ID");
|
||||
fix[i] = modify->fix[ifix];
|
||||
if (nevery % modify->fix[ifix]->peratom_freq)
|
||||
error->all(FLERR, "Dump custom and fix not computed at compatible times");
|
||||
fix[i] = modify->get_fix_by_id(id_fix[i]);
|
||||
if (!fix[i]) error->all(FLERR, "Could not find dump custom/adios fix ID {}", id_fix[i]);
|
||||
if (nevery % fix[i]->peratom_freq)
|
||||
error->all(FLERR, "dump custom/adios and fix {} with ID {} not computed at compatible times",
|
||||
fix[i]->style, id_fix[i]);
|
||||
}
|
||||
|
||||
int ivariable;
|
||||
for (int i = 0; i < nvariable; i++) {
|
||||
ivariable = input->variable->find(id_variable[i]);
|
||||
if (ivariable < 0) error->all(FLERR, "Could not find dump custom variable name");
|
||||
if (ivariable < 0) error->all(FLERR, "Could not find dump custom/adios variable name");
|
||||
variable[i] = ivariable;
|
||||
}
|
||||
|
||||
// set index and check validity of region
|
||||
if (iregion >= 0) {
|
||||
iregion = domain->find_region(idregion);
|
||||
if (iregion == -1) error->all(FLERR, "Region ID for dump custom does not exist");
|
||||
}
|
||||
if (idregion && !domain->get_region_by_id(idregion))
|
||||
error->all(FLERR, "Region {} for dump custom/adios does not exist", idregion);
|
||||
|
||||
/* Define the group of variables for the atom style here since it's a fixed
|
||||
* set */
|
||||
@ -316,7 +314,7 @@ void DumpCustomADIOS::init_style()
|
||||
int *boundaryptr = reinterpret_cast<int *>(domain->boundary);
|
||||
internal->io.DefineAttribute<int>("boundary", boundaryptr, 6);
|
||||
|
||||
auto nColumns = static_cast<size_t>(size_one);
|
||||
auto nColumns = static_cast<size_t>(size_one);
|
||||
internal->io.DefineAttribute<std::string>("columns", internal->columnNames.data(), nColumns);
|
||||
internal->io.DefineAttribute<std::string>("columnstr", columns);
|
||||
internal->io.DefineAttribute<std::string>("boundarystr", boundstr);
|
||||
|
||||
@ -221,7 +221,7 @@ bigint ReaderADIOS::read_header(double box[3][3], int &boxinfo, int &triclinic,
|
||||
uint64_t rem = nAtomsTotal % comm->nprocs;
|
||||
nAtoms = nAtomsTotal / comm->nprocs;
|
||||
atomOffset = comm->me * nAtoms;
|
||||
if (comm->me < (int)rem) {
|
||||
if (comm->me < (int) rem) {
|
||||
++nAtoms;
|
||||
atomOffset += comm->me;
|
||||
} else {
|
||||
@ -421,7 +421,7 @@ void ReaderADIOS::read_atoms(int n, int nfield, double **fields)
|
||||
|
||||
adios2::Variable<double> varAtoms = internal->io.InquireVariable<double>("atoms");
|
||||
|
||||
if ((uint64_t)n != nAtoms)
|
||||
if ((uint64_t) n != nAtoms)
|
||||
error->one(FLERR,
|
||||
"ReaderADIOS::read_atoms() expects 'n={}' equal to the number of "
|
||||
"atoms (={}) for process {} in ADIOS file {}.",
|
||||
|
||||
@ -40,8 +40,8 @@ class ReaderADIOS : public Reader {
|
||||
|
||||
int read_time(bigint &) override;
|
||||
void skip() override;
|
||||
bigint read_header(double[3][3], int &, int &, int, int, int *, char **, int, int, int &,
|
||||
int &, int &, int &) override;
|
||||
bigint read_header(double[3][3], int &, int &, int, int, int *, char **, int, int, int &, int &,
|
||||
int &, int &) override;
|
||||
void read_atoms(int, int, double **) override;
|
||||
|
||||
void open_file(const std::string &) override;
|
||||
|
||||
@ -317,10 +317,10 @@ void PairRESquared::coeff(int narg, char **arg)
|
||||
|
||||
void PairRESquared::init_style()
|
||||
{
|
||||
avec = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
|
||||
avec = dynamic_cast<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
|
||||
if (!avec) error->all(FLERR, "Pair resquared requires atom style ellipsoid");
|
||||
|
||||
neighbor->add_request(this,NeighConst::REQ_DEFAULT);
|
||||
neighbor->add_request(this, NeighConst::REQ_DEFAULT);
|
||||
|
||||
// per-type shape precalculations
|
||||
// require that atom shapes are identical within each type
|
||||
|
||||
@ -30,7 +30,7 @@ namespace LAMMPS_NS {
|
||||
class FixBocs : public Fix {
|
||||
public:
|
||||
FixBocs(class LAMMPS *, int, char **); // MRD NJD
|
||||
~FixBocs() override; // MRD NJD
|
||||
~FixBocs() override; // MRD NJD
|
||||
int setmask() override;
|
||||
void init() override;
|
||||
void setup(int) override;
|
||||
|
||||
@ -29,7 +29,8 @@ using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
BondBPMSpring::BondBPMSpring(LAMMPS *_lmp) : BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr)
|
||||
BondBPMSpring::BondBPMSpring(LAMMPS *_lmp) :
|
||||
BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr)
|
||||
{
|
||||
partial_flag = 1;
|
||||
smooth_flag = 1;
|
||||
|
||||
@ -56,7 +56,7 @@ FixBrownianAsphere::FixBrownianAsphere(LAMMPS *lmp, int narg, char **arg) :
|
||||
|
||||
void FixBrownianAsphere::init()
|
||||
{
|
||||
avec = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
|
||||
avec = dynamic_cast<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
|
||||
if (!avec) error->all(FLERR, "Compute brownian/asphere requires atom style ellipsoid");
|
||||
|
||||
// check that all particles are finite-size ellipsoids
|
||||
|
||||
@ -48,9 +48,9 @@ class FixBrownianBase : public Fix {
|
||||
int noise_flag; // 0/1 for noise off/on
|
||||
int gaussian_noise_flag; // 0/1 for uniform/gaussian noise
|
||||
|
||||
double temp; // temperature
|
||||
double rot_temp; // temperature
|
||||
double g1, g2; // prefactors in time stepping
|
||||
double temp; // temperature
|
||||
double rot_temp; // temperature
|
||||
double g1, g2; // prefactors in time stepping
|
||||
|
||||
class RanMars *rng;
|
||||
};
|
||||
|
||||
@ -100,7 +100,7 @@ void FixPropelSelf::init()
|
||||
error->all(FLERR, "Fix propel/self requires atom attribute mu with option dipole");
|
||||
|
||||
if (mode == QUAT) {
|
||||
avec = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
|
||||
avec = dynamic_cast<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
|
||||
if (!avec) error->all(FLERR, "Fix propel/self requires atom style ellipsoid with option quat");
|
||||
|
||||
// check that all particles are finite-size ellipsoids
|
||||
|
||||
@ -145,16 +145,16 @@ void BondOxdnaFene::ev_tally_xyz(int i, int j, int nlocal, int newton_bond, doub
|
||||
------------------------------------------------------------------------- */
|
||||
void BondOxdnaFene::compute(int eflag, int vflag)
|
||||
{
|
||||
int a,b,in,type;
|
||||
double delf[3],delta[3],deltb[3]; // force, torque increment;;
|
||||
double delr[3],ebond,fbond;
|
||||
double rsq,Deltasq,rlogarg;
|
||||
double r,rr0,rr0sq;
|
||||
int a, b, in, type;
|
||||
double delf[3], delta[3], deltb[3]; // force, torque increment;;
|
||||
double delr[3], ebond, fbond;
|
||||
double rsq, Deltasq, rlogarg;
|
||||
double r, rr0, rr0sq;
|
||||
// vectors COM-backbone site in lab frame
|
||||
double ra_cs[3],rb_cs[3];
|
||||
double ra_cs[3], rb_cs[3];
|
||||
// Cartesian unit vectors in lab frame
|
||||
double ax[3],ay[3],az[3];
|
||||
double bx[3],by[3],bz[3];
|
||||
double ax[3], ay[3], az[3];
|
||||
double bx[3], by[3], bz[3];
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
@ -170,9 +170,9 @@ void BondOxdnaFene::compute(int eflag, int vflag)
|
||||
|
||||
// n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv
|
||||
int dim;
|
||||
nx_xtrct = (double **) force->pair->extract("nx",dim);
|
||||
ny_xtrct = (double **) force->pair->extract("ny",dim);
|
||||
nz_xtrct = (double **) force->pair->extract("nz",dim);
|
||||
nx_xtrct = (double **) force->pair->extract("nx", dim);
|
||||
ny_xtrct = (double **) force->pair->extract("ny", dim);
|
||||
nz_xtrct = (double **) force->pair->extract("nz", dim);
|
||||
|
||||
// loop over FENE bonds
|
||||
|
||||
|
||||
@ -39,8 +39,8 @@ class BondOxdnaFene : public Bond {
|
||||
double single(int, double, int, int, double &) override;
|
||||
|
||||
protected:
|
||||
double *k, *Delta, *r0; // FENE
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
double *k, *Delta, *r0; // FENE
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
|
||||
void allocate();
|
||||
void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double);
|
||||
|
||||
@ -55,7 +55,7 @@ class PairOxdna2Coaxstk : public Pair {
|
||||
double **a_cxst6, **theta_cxst6_0, **dtheta_cxst6_ast;
|
||||
double **b_cxst6, **dtheta_cxst6_c;
|
||||
double **AA_cxst1, **BB_cxst1;
|
||||
double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
|
||||
virtual void allocate();
|
||||
};
|
||||
|
||||
@ -45,7 +45,7 @@ class PairOxdna2Dh : public Pair {
|
||||
protected:
|
||||
double **qeff_dh_pf, **kappa_dh;
|
||||
double **b_dh, **cut_dh_ast, **cutsq_dh_ast, **cut_dh_c, **cutsq_dh_c;
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
|
||||
virtual void allocate();
|
||||
};
|
||||
|
||||
@ -57,7 +57,7 @@ class PairOxdnaCoaxstk : public Pair {
|
||||
double **b_cxst6, **dtheta_cxst6_c;
|
||||
double **a_cxst3p, **cosphi_cxst3p_ast, **b_cxst3p, **cosphi_cxst3p_c;
|
||||
double **a_cxst4p, **cosphi_cxst4p_ast, **b_cxst4p, **cosphi_cxst4p_c;
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
|
||||
virtual void allocate();
|
||||
};
|
||||
|
||||
@ -54,7 +54,7 @@ class PairOxdnaExcv : public Pair {
|
||||
double **lj1_sb, **lj2_sb, **b_sb, **cut_sb_c, **cutsq_sb_c;
|
||||
double **epsilon_bb, **sigma_bb, **cut_bb_ast, **cutsq_bb_ast;
|
||||
double **lj1_bb, **lj2_bb, **b_bb, **cut_bb_c, **cutsq_bb_c;
|
||||
double **nx, **ny, **nz; // per-atom arrays for local unit vectors
|
||||
double **nx, **ny, **nz; // per-atom arrays for local unit vectors
|
||||
|
||||
virtual void allocate();
|
||||
};
|
||||
|
||||
@ -60,7 +60,7 @@ class PairOxdnaHbond : public Pair {
|
||||
double **b_hb7, **dtheta_hb7_c;
|
||||
double **a_hb8, **theta_hb8_0, **dtheta_hb8_ast;
|
||||
double **b_hb8, **dtheta_hb8_c;
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
int seqdepflag;
|
||||
|
||||
virtual void allocate();
|
||||
|
||||
@ -59,7 +59,7 @@ class PairOxdnaStk : public Pair {
|
||||
double **b_st6, **dtheta_st6_c;
|
||||
double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c;
|
||||
double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c;
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
int seqdepflag;
|
||||
|
||||
virtual void allocate();
|
||||
|
||||
@ -59,7 +59,7 @@ class PairOxdnaXstk : public Pair {
|
||||
double **b_xst7, **dtheta_xst7_c;
|
||||
double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast;
|
||||
double **b_xst8, **dtheta_xst8_c;
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
|
||||
virtual void allocate();
|
||||
};
|
||||
|
||||
@ -60,7 +60,7 @@ class PairOxrna2Stk : public Pair {
|
||||
double **b_st10, **dtheta_st10_c;
|
||||
double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c;
|
||||
double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c;
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
int seqdepflag;
|
||||
|
||||
virtual void allocate();
|
||||
|
||||
@ -56,7 +56,7 @@ class PairOxrna2Xstk : public Pair {
|
||||
double **b_xst7, **dtheta_xst7_c;
|
||||
double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast;
|
||||
double **b_xst8, **dtheta_xst8_c;
|
||||
double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
|
||||
|
||||
virtual void allocate();
|
||||
};
|
||||
|
||||
@ -484,7 +484,7 @@ void PairLJClass2::init_style()
|
||||
int list_style = NeighConst::REQ_DEFAULT;
|
||||
|
||||
if (update->whichflag == 1 && utils::strmatch(update->integrate_style, "^respa")) {
|
||||
auto respa = dynamic_cast<Respa *>( update->integrate);
|
||||
auto respa = dynamic_cast<Respa *>(update->integrate);
|
||||
if (respa->level_inner >= 0) list_style = NeighConst::REQ_RESPA_INOUT;
|
||||
if (respa->level_middle >= 0) list_style = NeighConst::REQ_RESPA_ALL;
|
||||
}
|
||||
@ -493,8 +493,8 @@ void PairLJClass2::init_style()
|
||||
// set rRESPA cutoffs
|
||||
|
||||
if (utils::strmatch(update->integrate_style, "^respa") &&
|
||||
(dynamic_cast<Respa *>( update->integrate))->level_inner >= 0)
|
||||
cut_respa = (dynamic_cast<Respa *>( update->integrate))->cutoff;
|
||||
(dynamic_cast<Respa *>(update->integrate))->level_inner >= 0)
|
||||
cut_respa = (dynamic_cast<Respa *>(update->integrate))->cutoff;
|
||||
else
|
||||
cut_respa = nullptr;
|
||||
}
|
||||
|
||||
@ -672,7 +672,7 @@ void PairLJClass2CoulLong::init_style()
|
||||
int list_style = NeighConst::REQ_DEFAULT;
|
||||
|
||||
if (update->whichflag == 1 && utils::strmatch(update->integrate_style, "^respa")) {
|
||||
auto respa = dynamic_cast<Respa *>( update->integrate);
|
||||
auto respa = dynamic_cast<Respa *>(update->integrate);
|
||||
if (respa->level_inner >= 0) list_style = NeighConst::REQ_RESPA_INOUT;
|
||||
if (respa->level_middle >= 0) list_style = NeighConst::REQ_RESPA_ALL;
|
||||
}
|
||||
@ -683,8 +683,8 @@ void PairLJClass2CoulLong::init_style()
|
||||
// set rRESPA cutoffs
|
||||
|
||||
if (utils::strmatch(update->integrate_style, "^respa") &&
|
||||
(dynamic_cast<Respa *>( update->integrate))->level_inner >= 0)
|
||||
cut_respa = (dynamic_cast<Respa *>( update->integrate))->cutoff;
|
||||
(dynamic_cast<Respa *>(update->integrate))->level_inner >= 0)
|
||||
cut_respa = (dynamic_cast<Respa *>(update->integrate))->cutoff;
|
||||
else
|
||||
cut_respa = nullptr;
|
||||
|
||||
|
||||
@ -494,7 +494,7 @@ void PairBrownian::init_style()
|
||||
else if (strstr(modify->fix[i]->style, "wall") != nullptr) {
|
||||
if (flagwall) error->all(FLERR, "Cannot use multiple fix wall commands with pair brownian");
|
||||
flagwall = 1; // Walls exist
|
||||
wallfix = dynamic_cast<FixWall *>( modify->fix[i]);
|
||||
wallfix = dynamic_cast<FixWall *>(modify->fix[i]);
|
||||
if (wallfix->xflag) flagwall = 2; // Moving walls exist
|
||||
}
|
||||
}
|
||||
|
||||
@ -101,7 +101,8 @@ class colvarproxy_lammps : public colvarproxy {
|
||||
void log(std::string const &message) override;
|
||||
void error(std::string const &message) override;
|
||||
|
||||
cvm::rvector position_distance(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) const override;
|
||||
cvm::rvector position_distance(cvm::atom_pos const &pos1,
|
||||
cvm::atom_pos const &pos2) const override;
|
||||
|
||||
int backup_file(char const *filename) override;
|
||||
|
||||
|
||||
@ -91,8 +91,10 @@ void DumpCFGGZ::write_header(bigint n)
|
||||
// so molecules are not split across periodic box boundaries
|
||||
|
||||
double scale = 1.0;
|
||||
if (atom->peri_flag) scale = atom->pdscale;
|
||||
else if (unwrapflag == 1) scale = UNWRAPEXPAND;
|
||||
if (atom->peri_flag)
|
||||
scale = atom->pdscale;
|
||||
else if (unwrapflag == 1)
|
||||
scale = UNWRAPEXPAND;
|
||||
|
||||
std::string header = fmt::format("Number of particles = {}\n", n);
|
||||
header += fmt::format("A = {:g} Angstrom (basic length-scale)\n", scale);
|
||||
|
||||
@ -26,7 +26,7 @@ namespace LAMMPS_NS {
|
||||
|
||||
class PairLJSFDipoleSF : public Pair {
|
||||
public:
|
||||
PairLJSFDipoleSF(class LAMMPS *_lmp) : Pair(_lmp) {};
|
||||
PairLJSFDipoleSF(class LAMMPS *_lmp) : Pair(_lmp){};
|
||||
~PairLJSFDipoleSF() override;
|
||||
void compute(int, int) override;
|
||||
void settings(int, char **) override;
|
||||
|
||||
@ -25,7 +25,7 @@ FixStyle(drude/transform/inverse,FixDrudeTransform<true>);
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
template <bool inverse> class FixDrudeTransform: public Fix {
|
||||
template <bool inverse> class FixDrudeTransform : public Fix {
|
||||
public:
|
||||
FixDrudeTransform(class LAMMPS *, int, char **);
|
||||
~FixDrudeTransform() override;
|
||||
|
||||
@ -97,6 +97,10 @@ if (test $1 = "GRANULAR") then
|
||||
depend OPENMP
|
||||
fi
|
||||
|
||||
if (test $1 = "INTERLAYER") then
|
||||
depend OPT
|
||||
fi
|
||||
|
||||
if (test $1 = "KSPACE") then
|
||||
depend CG-SDK
|
||||
depend CORESHELL
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -16,7 +15,6 @@
|
||||
Contributing author: Andres Jaramillo-Botero (Caltech)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
|
||||
#include "compute_temp_region_eff.h"
|
||||
|
||||
#include "atom.h"
|
||||
@ -33,16 +31,15 @@ using namespace LAMMPS_NS;
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeTempRegionEff::ComputeTempRegionEff(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg)
|
||||
Compute(lmp, narg, arg), region(nullptr), idregion(nullptr)
|
||||
{
|
||||
if (!atom->electron_flag)
|
||||
error->all(FLERR,"Compute temp/region/eff requires atom style electron");
|
||||
error->all(FLERR, "Compute temp/region/eff requires atom style electron");
|
||||
|
||||
if (narg != 4) error->all(FLERR,"Illegal compute temp/region/eff command");
|
||||
if (narg != 4) error->all(FLERR, "Illegal compute temp/region/eff command");
|
||||
|
||||
iregion = domain->find_region(arg[3]);
|
||||
if (iregion == -1)
|
||||
error->all(FLERR,"Region ID for compute temp/region/eff does not exist");
|
||||
region = domain->get_region_by_id(arg[3]);
|
||||
if (!region) error->all(FLERR, "Region {} for compute temp/region/eff does not exist", arg[3]);
|
||||
idregion = utils::strdup(arg[3]);
|
||||
|
||||
scalar_flag = vector_flag = 1;
|
||||
@ -61,9 +58,9 @@ ComputeTempRegionEff::ComputeTempRegionEff(LAMMPS *lmp, int narg, char **arg) :
|
||||
|
||||
ComputeTempRegionEff::~ComputeTempRegionEff()
|
||||
{
|
||||
delete [] idregion;
|
||||
delete[] idregion;
|
||||
memory->destroy(vbiasall);
|
||||
delete [] vector;
|
||||
delete[] vector;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -72,9 +69,8 @@ void ComputeTempRegionEff::init()
|
||||
{
|
||||
// set index and check validity of region
|
||||
|
||||
iregion = domain->find_region(idregion);
|
||||
if (iregion == -1)
|
||||
error->all(FLERR,"Region ID for compute temp/region/eff does not exist");
|
||||
region = domain->get_region_by_id(idregion);
|
||||
if (!region) error->all(FLERR, "Region {} for compute temp/region/eff does not exist", idregion);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -90,7 +86,7 @@ void ComputeTempRegionEff::setup()
|
||||
|
||||
void ComputeTempRegionEff::dof_remove_pre()
|
||||
{
|
||||
domain->regions[iregion]->prematch();
|
||||
region->prematch();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -98,7 +94,7 @@ void ComputeTempRegionEff::dof_remove_pre()
|
||||
int ComputeTempRegionEff::dof_remove(int i)
|
||||
{
|
||||
double *x = atom->x[i];
|
||||
if (domain->regions[iregion]->match(x[0],x[1],x[2])) return 0;
|
||||
if (region->match(x[0], x[1], x[2])) return 0;
|
||||
return 1;
|
||||
}
|
||||
|
||||
@ -116,9 +112,8 @@ double ComputeTempRegionEff::compute_scalar()
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
double mefactor = domain->dimension/4.0;
|
||||
double mefactor = domain->dimension / 4.0;
|
||||
|
||||
Region *region = domain->regions[iregion];
|
||||
region->prematch();
|
||||
|
||||
int count = 0;
|
||||
@ -127,34 +122,35 @@ double ComputeTempRegionEff::compute_scalar()
|
||||
|
||||
if (mass) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
|
||||
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) {
|
||||
count++;
|
||||
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
|
||||
mass[type[i]];
|
||||
if (abs(spin[i])==1) {
|
||||
t += mefactor*mass[type[i]]*ervel[i]*ervel[i];
|
||||
t += (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]) * mass[type[i]];
|
||||
if (abs(spin[i]) == 1) {
|
||||
t += mefactor * mass[type[i]] * ervel[i] * ervel[i];
|
||||
ecount++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double tarray[2],tarray_all[2];
|
||||
double tarray[2], tarray_all[2];
|
||||
// Assume 3/2 k T per nucleus
|
||||
tarray[0] = count-ecount;
|
||||
tarray[0] = count - ecount;
|
||||
tarray[1] = t;
|
||||
MPI_Allreduce(tarray,tarray_all,2,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(tarray, tarray_all, 2, MPI_DOUBLE, MPI_SUM, world);
|
||||
dof = domain->dimension * tarray_all[0] - extra_dof;
|
||||
if (dof < 0.0 && tarray_all[0] > 0.0)
|
||||
error->all(FLERR,"Temperature compute degrees of freedom < 0");
|
||||
error->all(FLERR, "Temperature compute degrees of freedom < 0");
|
||||
|
||||
int one = 0;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
|
||||
if (abs(spin[i])==1) one++;
|
||||
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) {
|
||||
if (abs(spin[i]) == 1) one++;
|
||||
}
|
||||
|
||||
if (dof > 0.0) scalar = force->mvv2e * tarray_all[1] / (dof * force->boltz);
|
||||
else scalar = 0.0;
|
||||
if (dof > 0.0)
|
||||
scalar = force->mvv2e * tarray_all[1] / (dof * force->boltz);
|
||||
else
|
||||
scalar = 0.0;
|
||||
return scalar;
|
||||
}
|
||||
|
||||
@ -174,33 +170,32 @@ void ComputeTempRegionEff::compute_vector()
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
double mefactor = domain->dimension/4.0;
|
||||
double mefactor = domain->dimension / 4.0;
|
||||
|
||||
Region *region = domain->regions[iregion];
|
||||
region->prematch();
|
||||
|
||||
double massone,t[6];
|
||||
double massone, t[6];
|
||||
for (i = 0; i < 6; i++) t[i] = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
|
||||
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) {
|
||||
massone = mass[type[i]];
|
||||
|
||||
t[0] += massone * v[i][0]*v[i][0];
|
||||
t[1] += massone * v[i][1]*v[i][1];
|
||||
t[2] += massone * v[i][2]*v[i][2];
|
||||
t[3] += massone * v[i][0]*v[i][1];
|
||||
t[4] += massone * v[i][0]*v[i][2];
|
||||
t[5] += massone * v[i][1]*v[i][2];
|
||||
t[0] += massone * v[i][0] * v[i][0];
|
||||
t[1] += massone * v[i][1] * v[i][1];
|
||||
t[2] += massone * v[i][2] * v[i][2];
|
||||
t[3] += massone * v[i][0] * v[i][1];
|
||||
t[4] += massone * v[i][0] * v[i][2];
|
||||
t[5] += massone * v[i][1] * v[i][2];
|
||||
|
||||
if (abs(spin[i])==1) {
|
||||
t[0] += mefactor * massone * ervel[i]*ervel[i];
|
||||
t[1] += mefactor * massone * ervel[i]*ervel[i];
|
||||
t[2] += mefactor * massone * ervel[i]*ervel[i];
|
||||
if (abs(spin[i]) == 1) {
|
||||
t[0] += mefactor * massone * ervel[i] * ervel[i];
|
||||
t[1] += mefactor * massone * ervel[i] * ervel[i];
|
||||
t[2] += mefactor * massone * ervel[i] * ervel[i];
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(t, vector, 6, MPI_DOUBLE, MPI_SUM, world);
|
||||
for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
|
||||
}
|
||||
|
||||
@ -212,7 +207,7 @@ void ComputeTempRegionEff::compute_vector()
|
||||
void ComputeTempRegionEff::remove_bias(int i, double *v)
|
||||
{
|
||||
double *x = atom->x[i];
|
||||
if (domain->regions[iregion]->match(x[0],x[1],x[2]))
|
||||
if (region->match(x[0], x[1], x[2]))
|
||||
vbias[0] = vbias[1] = vbias[2] = 0.0;
|
||||
else {
|
||||
vbias[0] = v[0];
|
||||
@ -236,14 +231,12 @@ void ComputeTempRegionEff::remove_bias_all()
|
||||
if (atom->nmax > maxbias) {
|
||||
memory->destroy(vbiasall);
|
||||
maxbias = atom->nmax;
|
||||
memory->create(vbiasall,maxbias,3,"temp/region:vbiasall");
|
||||
memory->create(vbiasall, maxbias, 3, "temp/region:vbiasall");
|
||||
}
|
||||
|
||||
Region *region = domain->regions[iregion];
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (region->match(x[i][0],x[i][1],x[i][2]))
|
||||
if (region->match(x[i][0], x[i][1], x[i][2]))
|
||||
vbiasall[i][0] = vbiasall[i][1] = vbiasall[i][2] = 0.0;
|
||||
else {
|
||||
vbiasall[i][0] = v[i][0];
|
||||
@ -289,6 +282,6 @@ void ComputeTempRegionEff::restore_bias_all()
|
||||
|
||||
double ComputeTempRegionEff::memory_usage()
|
||||
{
|
||||
double bytes = (double)maxbias * sizeof(double);
|
||||
double bytes = (double) maxbias * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
||||
@ -42,7 +42,7 @@ class ComputeTempRegionEff : public Compute {
|
||||
double memory_usage() override;
|
||||
|
||||
protected:
|
||||
int iregion;
|
||||
class Region *region;
|
||||
char *idregion;
|
||||
};
|
||||
|
||||
|
||||
@ -30,14 +30,12 @@
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeAveSphereAtom::ComputeAveSphereAtom(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg),
|
||||
result(nullptr)
|
||||
Compute(lmp, narg, arg), result(nullptr)
|
||||
{
|
||||
if (narg < 3 || narg > 5) error->all(FLERR,"Illegal compute ave/sphere/atom command");
|
||||
if (narg < 3 || narg > 5) error->all(FLERR, "Illegal compute ave/sphere/atom command");
|
||||
|
||||
// process optional args
|
||||
|
||||
@ -45,12 +43,13 @@ ComputeAveSphereAtom::ComputeAveSphereAtom(LAMMPS *lmp, int narg, char **arg) :
|
||||
|
||||
int iarg = 3;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"cutoff") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal compute ave/sphere/atom command");
|
||||
cutoff = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
if (cutoff <= 0.0) error->all(FLERR,"Illegal compute ave/sphere/atom command");
|
||||
if (strcmp(arg[iarg], "cutoff") == 0) {
|
||||
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute ave/sphere/atom command");
|
||||
cutoff = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
|
||||
if (cutoff <= 0.0) error->all(FLERR, "Illegal compute ave/sphere/atom command");
|
||||
iarg += 2;
|
||||
} else error->all(FLERR,"Illegal compute ave/sphere/atom command");
|
||||
} else
|
||||
error->all(FLERR, "Illegal compute ave/sphere/atom command");
|
||||
}
|
||||
|
||||
peratom_flag = 1;
|
||||
@ -74,32 +73,32 @@ ComputeAveSphereAtom::~ComputeAveSphereAtom()
|
||||
void ComputeAveSphereAtom::init()
|
||||
{
|
||||
if (!force->pair && cutoff == 0.0)
|
||||
error->all(FLERR,"Compute ave/sphere/atom requires a cutoff be specified "
|
||||
error->all(FLERR,
|
||||
"Compute ave/sphere/atom requires a cutoff be specified "
|
||||
"or a pair style be defined");
|
||||
|
||||
double skin = neighbor->skin;
|
||||
if (cutoff != 0.0) {
|
||||
double cutghost; // as computed by Neighbor and Comm
|
||||
double cutghost; // as computed by Neighbor and Comm
|
||||
if (force->pair)
|
||||
cutghost = MAX(force->pair->cutforce+skin,comm->cutghostuser);
|
||||
cutghost = MAX(force->pair->cutforce + skin, comm->cutghostuser);
|
||||
else
|
||||
cutghost = comm->cutghostuser;
|
||||
|
||||
if (cutoff > cutghost)
|
||||
error->all(FLERR,"Compute ave/sphere/atom cutoff exceeds ghost atom range - "
|
||||
error->all(FLERR,
|
||||
"Compute ave/sphere/atom cutoff exceeds ghost atom range - "
|
||||
"use comm_modify cutoff command");
|
||||
}
|
||||
|
||||
int cutflag = 1;
|
||||
if (force->pair) {
|
||||
if (cutoff == 0.0) {
|
||||
cutoff = force->pair->cutforce;
|
||||
}
|
||||
if (cutoff <= force->pair->cutforce+skin) cutflag = 0;
|
||||
if (cutoff == 0.0) { cutoff = force->pair->cutforce; }
|
||||
if (cutoff <= force->pair->cutforce + skin) cutflag = 0;
|
||||
}
|
||||
|
||||
cutsq = cutoff*cutoff;
|
||||
sphere_vol = 4.0/3.0*MY_PI*cutsq*cutoff;
|
||||
cutsq = cutoff * cutoff;
|
||||
sphere_vol = 4.0 / 3.0 * MY_PI * cutsq * cutoff;
|
||||
|
||||
// need an occasional full neighbor list
|
||||
|
||||
@ -118,11 +117,11 @@ void ComputeAveSphereAtom::init_list(int /*id*/, NeighList *ptr)
|
||||
|
||||
void ComputeAveSphereAtom::compute_peratom()
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int i, j, ii, jj, inum, jnum;
|
||||
double xtmp, ytmp, ztmp, delx, dely, delz, rsq;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
int count;
|
||||
double vsum[3],vavg[3],vnet[3];
|
||||
double vsum[3], vavg[3], vnet[3];
|
||||
|
||||
invoked_peratom = update->ntimestep;
|
||||
|
||||
@ -131,7 +130,7 @@ void ComputeAveSphereAtom::compute_peratom()
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(result);
|
||||
nmax = atom->nmax;
|
||||
memory->create(result,nmax,2,"ave/sphere/atom:result");
|
||||
memory->create(result, nmax, 2, "ave/sphere/atom:result");
|
||||
array_atom = result;
|
||||
}
|
||||
|
||||
@ -179,7 +178,7 @@ void ComputeAveSphereAtom::compute_peratom()
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
if (rsq < cutsq) {
|
||||
count++;
|
||||
vsum[0] += v[j][0];
|
||||
@ -188,9 +187,9 @@ void ComputeAveSphereAtom::compute_peratom()
|
||||
}
|
||||
}
|
||||
|
||||
vavg[0] = vsum[0]/count;
|
||||
vavg[1] = vsum[1]/count;
|
||||
vavg[2] = vsum[2]/count;
|
||||
vavg[0] = vsum[0] / count;
|
||||
vavg[1] = vsum[1] / count;
|
||||
vavg[2] = vsum[2] / count;
|
||||
|
||||
// i atom contribution
|
||||
|
||||
@ -198,7 +197,7 @@ void ComputeAveSphereAtom::compute_peratom()
|
||||
vnet[0] = v[i][0] - vavg[0];
|
||||
vnet[1] = v[i][1] - vavg[1];
|
||||
vnet[2] = v[i][2] - vavg[2];
|
||||
double ke_sum = vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2];
|
||||
double ke_sum = vnet[0] * vnet[0] + vnet[1] * vnet[1] + vnet[2] * vnet[2];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
@ -207,17 +206,17 @@ void ComputeAveSphereAtom::compute_peratom()
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
if (rsq < cutsq) {
|
||||
count++;
|
||||
vnet[0] = v[j][0] - vavg[0];
|
||||
vnet[1] = v[j][1] - vavg[1];
|
||||
vnet[2] = v[j][2] - vavg[2];
|
||||
ke_sum += vnet[0]*vnet[0] + vnet[1]*vnet[1] + vnet[2]*vnet[2];
|
||||
ke_sum += vnet[0] * vnet[0] + vnet[1] * vnet[1] + vnet[2] * vnet[2];
|
||||
}
|
||||
}
|
||||
double density = count/sphere_vol;
|
||||
double temp = ke_sum/3.0/count;
|
||||
double density = count / sphere_vol;
|
||||
double temp = ke_sum / 3.0 / count;
|
||||
result[i][0] = density;
|
||||
result[i][1] = temp;
|
||||
}
|
||||
@ -226,12 +225,12 @@ void ComputeAveSphereAtom::compute_peratom()
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int ComputeAveSphereAtom::pack_forward_comm(int n, int *list, double *buf,
|
||||
int /*pbc_flag*/, int * /*pbc*/)
|
||||
int ComputeAveSphereAtom::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
|
||||
int * /*pbc*/)
|
||||
{
|
||||
double **v = atom->v;
|
||||
|
||||
int i,m=0;
|
||||
int i, m = 0;
|
||||
for (i = 0; i < n; ++i) {
|
||||
buf[m++] = v[list[i]][0];
|
||||
buf[m++] = v[list[i]][1];
|
||||
@ -247,7 +246,7 @@ void ComputeAveSphereAtom::unpack_forward_comm(int n, int first, double *buf)
|
||||
{
|
||||
double **v = atom->v;
|
||||
|
||||
int i,last,m=0;
|
||||
int i, last, m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; ++i) {
|
||||
v[i][0] = buf[m++];
|
||||
@ -262,6 +261,6 @@ void ComputeAveSphereAtom::unpack_forward_comm(int n, int first, double *buf)
|
||||
|
||||
double ComputeAveSphereAtom::memory_usage()
|
||||
{
|
||||
double bytes = (double)2*nmax * sizeof(double);
|
||||
double bytes = (double) 2 * nmax * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
||||
@ -12,9 +12,9 @@
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef COMPUTE_CLASS
|
||||
|
||||
ComputeStyle(ave/sphere/atom,ComputeAveSphereAtom)
|
||||
|
||||
// clang-format off
|
||||
ComputeStyle(ave/sphere/atom,ComputeAveSphereAtom);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_COMPUTE_AVE_SPHERE_ATOM_H
|
||||
@ -37,13 +37,13 @@ class ComputeAveSphereAtom : public Compute {
|
||||
|
||||
protected:
|
||||
int nmax;
|
||||
double cutoff,cutsq,sphere_vol;
|
||||
double cutoff, cutsq, sphere_vol;
|
||||
class NeighList *list;
|
||||
|
||||
double **result;
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
1188
src/EXTRA-COMPUTE/compute_born_matrix.cpp
Normal file
1188
src/EXTRA-COMPUTE/compute_born_matrix.cpp
Normal file
File diff suppressed because it is too large
Load Diff
97
src/EXTRA-COMPUTE/compute_born_matrix.h
Normal file
97
src/EXTRA-COMPUTE/compute_born_matrix.h
Normal file
@ -0,0 +1,97 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://www.lammps.org/, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/*------------------------------------------------------------------------
|
||||
Contributing Authors : Germain Clavier (TUe), Aidan Thompson (Sandia)
|
||||
--------------------------------------------------------------------------*/
|
||||
|
||||
#ifdef COMPUTE_CLASS
|
||||
// clang-format off
|
||||
ComputeStyle(born/matrix,ComputeBornMatrix);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_COMPUTE_BORN_MATRIX_H
|
||||
#define LMP_COMPUTE_BORN_MATRIX_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputeBornMatrix : public Compute {
|
||||
public:
|
||||
ComputeBornMatrix(class LAMMPS *, int, char **);
|
||||
~ComputeBornMatrix() override;
|
||||
void init() override;
|
||||
void init_list(int, class NeighList *) override;
|
||||
void compute_vector() override;
|
||||
double memory_usage() override;
|
||||
|
||||
private:
|
||||
// Born matrix contributions
|
||||
|
||||
void compute_pairs(); // pair and manybody
|
||||
void compute_bonds(); // bonds
|
||||
void compute_angles(); // angles
|
||||
void compute_dihedrals(); // dihedrals
|
||||
void compute_numdiff(); // stress virial finite differences
|
||||
void displace_atoms(int, int, double); // displace atoms
|
||||
void force_clear(int); // zero out force array
|
||||
void update_virial(); // recalculate the virial
|
||||
void restore_atoms(int, int); // restore atom positions
|
||||
void virial_addon(); // restore atom positions
|
||||
void reallocate(); // grow the atom arrays
|
||||
|
||||
int me; // process rank
|
||||
int nvalues; // length of elastic tensor
|
||||
int numflag; // 1 if using finite differences
|
||||
double numdelta; // size of finite strain
|
||||
int maxatom; // allocated size of atom arrays
|
||||
|
||||
int pairflag, bondflag, angleflag;
|
||||
int dihedflag, impflag, kspaceflag;
|
||||
|
||||
double *values_local, *values_global;
|
||||
double pos, pos1, dt, nktv2p, ftm2v;
|
||||
class NeighList *list;
|
||||
|
||||
char *id_virial; // name of virial compute
|
||||
class Compute *compute_virial; // pointer to virial compute
|
||||
|
||||
static constexpr int NDIR_VIRIAL = 6; // dimension of virial and strain vectors
|
||||
static constexpr int NXYZ_VIRIAL = 3; // number of Cartesian coordinates
|
||||
int revalbe[NDIR_VIRIAL][NDIR_VIRIAL];
|
||||
int virialVtoV[NDIR_VIRIAL];
|
||||
double **temp_x; // original coords
|
||||
double **temp_f; // original forces
|
||||
double fixedpoint[NXYZ_VIRIAL]; // displacement field origin
|
||||
int dirlist[NDIR_VIRIAL][2]; // strain cartesian indices
|
||||
};
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: ... style does not support compute born/matrix
|
||||
|
||||
Some component of the force field (pair, bond, angle...) does not provide
|
||||
a function to return the Born term contribution.
|
||||
*/
|
||||
@ -452,7 +452,8 @@ void ComputeStressCartesian::compute_pressure_1d(double fpair, double xi, double
|
||||
}
|
||||
|
||||
void ComputeStressCartesian::compute_pressure_2d(double fpair, double xi, double yi, double /*xj*/,
|
||||
double /*yj*/, double delx, double dely, double delz)
|
||||
double /*yj*/, double delx, double dely,
|
||||
double delz)
|
||||
{
|
||||
int bin1, bin2, next_bin1, next_bin2;
|
||||
double la = 0.0, lb = 0.0, l_sum = 0.0;
|
||||
|
||||
@ -28,7 +28,7 @@ class DumpYAML : public DumpCustom {
|
||||
public:
|
||||
DumpYAML(class LAMMPS *, int, char **);
|
||||
|
||||
protected:
|
||||
protected:
|
||||
bool thermo;
|
||||
|
||||
void init_style() override;
|
||||
|
||||
@ -59,8 +59,8 @@ extern "C" {
|
||||
|
||||
typedef int bool_t;
|
||||
|
||||
#if defined(_WIN32) || defined(__APPLE__) || defined(__FreeBSD__) || \
|
||||
defined(__DragonFly__) || defined(__OpenBSD__) || defined(__NetBSD__)
|
||||
#if defined(_WIN32) || defined(__APPLE__) || defined(__FreeBSD__) || defined(__DragonFly__) || \
|
||||
defined(__OpenBSD__) || defined(__NetBSD__)
|
||||
typedef char *caddr_t;
|
||||
typedef unsigned int u_int;
|
||||
#endif
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -39,63 +38,53 @@
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
#define MAXLINE 1024
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixElectronStopping::FixElectronStopping(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
Fix(lmp, narg, arg), elstop_ranges(nullptr), idregion(nullptr), region(nullptr), list(nullptr)
|
||||
{
|
||||
scalar_flag = 1; // Has compute_scalar
|
||||
global_freq = 1; // SeLoss computed every step
|
||||
extscalar = 0; // SeLoss compute_scalar is intensive
|
||||
nevery = 1; // Run fix every step
|
||||
|
||||
scalar_flag = 1; // Has compute_scalar
|
||||
global_freq = 1; // SeLoss computed every step
|
||||
extscalar = 0; // SeLoss compute_scalar is intensive
|
||||
nevery = 1; // Run fix every step
|
||||
|
||||
// args: 0 = fix ID, 1 = group ID, 2 = "electron/stopping"
|
||||
// 3 = Ecut, 4 = file path
|
||||
// optional rest: "region" <region name>
|
||||
// "minneigh" <min number of neighbors>
|
||||
|
||||
if (narg < 5) error->all(FLERR,
|
||||
"Illegal fix electron/stopping command: too few arguments");
|
||||
if (narg < 5) error->all(FLERR, "Illegal fix electron/stopping command: too few arguments");
|
||||
|
||||
Ecut = utils::numeric(FLERR, arg[3],false,lmp);
|
||||
if (Ecut <= 0.0) error->all(FLERR,
|
||||
"Illegal fix electron/stopping command: Ecut <= 0");
|
||||
Ecut = utils::numeric(FLERR, arg[3], false, lmp);
|
||||
if (Ecut <= 0.0) error->all(FLERR, "Illegal fix electron/stopping command: Ecut <= 0");
|
||||
|
||||
int iarg = 5;
|
||||
iregion = -1;
|
||||
minneigh = 1;
|
||||
bool minneighflag = false;
|
||||
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg], "region") == 0) {
|
||||
if (iregion >= 0) error->all(FLERR,
|
||||
"Illegal fix electron/stopping command: region given twice");
|
||||
if (iarg+2 > narg) error->all(FLERR,
|
||||
"Illegal fix electron/stopping command: region name missing");
|
||||
iregion = domain->find_region(arg[iarg+1]);
|
||||
if (iregion < 0) error->all(FLERR,
|
||||
"Region ID for fix electron/stopping does not exist");
|
||||
if (region) error->all(FLERR, "Illegal fix electron/stopping command: region given twice");
|
||||
if (iarg + 2 > narg)
|
||||
error->all(FLERR, "Illegal fix electron/stopping command: region name missing");
|
||||
region = domain->get_region_by_id(arg[iarg + 1]);
|
||||
if (!region)
|
||||
error->all(FLERR, "Region {} for fix electron/stopping does not exist", arg[iarg + 1]);
|
||||
idregion = utils::strdup(arg[iarg + 1]);
|
||||
iarg += 2;
|
||||
}
|
||||
else if (strcmp(arg[iarg], "minneigh") == 0) {
|
||||
if (minneighflag) error->all(FLERR,
|
||||
"Illegal fix electron/stopping command: minneigh given twice");
|
||||
} else if (strcmp(arg[iarg], "minneigh") == 0) {
|
||||
if (minneighflag)
|
||||
error->all(FLERR, "Illegal fix electron/stopping command: minneigh given twice");
|
||||
minneighflag = true;
|
||||
if (iarg+2 > narg) error->all(FLERR,
|
||||
"Illegal fix electron/stopping command: minneigh number missing");
|
||||
minneigh = utils::inumeric(FLERR, arg[iarg+1],false,lmp);
|
||||
if (minneigh < 0) error->all(FLERR,
|
||||
"Illegal fix electron/stopping command: minneigh < 0");
|
||||
if (iarg + 2 > narg)
|
||||
error->all(FLERR, "Illegal fix electron/stopping command: minneigh number missing");
|
||||
minneigh = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
|
||||
if (minneigh < 0) error->all(FLERR, "Illegal fix electron/stopping command: minneigh < 0");
|
||||
iarg += 2;
|
||||
}
|
||||
else error->all(FLERR,
|
||||
"Illegal fix electron/stopping command: unknown argument");
|
||||
} else
|
||||
error->all(FLERR, "Illegal fix electron/stopping command: unknown argument");
|
||||
}
|
||||
|
||||
|
||||
// Read the input file for energy ranges and stopping powers.
|
||||
// First proc 0 reads the file, then bcast to others.
|
||||
const int ncol = atom->ntypes + 1;
|
||||
@ -105,13 +94,12 @@ FixElectronStopping::FixElectronStopping(LAMMPS *lmp, int narg, char **arg) :
|
||||
read_table(arg[4]);
|
||||
}
|
||||
|
||||
MPI_Bcast(&maxlines, 1 , MPI_INT, 0, world);
|
||||
MPI_Bcast(&table_entries, 1 , MPI_INT, 0, world);
|
||||
MPI_Bcast(&maxlines, 1, MPI_INT, 0, world);
|
||||
MPI_Bcast(&table_entries, 1, MPI_INT, 0, world);
|
||||
|
||||
if (comm->me != 0)
|
||||
memory->create(elstop_ranges, ncol, maxlines, "electron/stopping:table");
|
||||
if (comm->me != 0) memory->create(elstop_ranges, ncol, maxlines, "electron/stopping:table");
|
||||
|
||||
MPI_Bcast(&elstop_ranges[0][0], ncol*maxlines, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&elstop_ranges[0][0], ncol * maxlines, MPI_DOUBLE, 0, world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -136,6 +124,10 @@ void FixElectronStopping::init()
|
||||
{
|
||||
SeLoss_sync_flag = 0;
|
||||
SeLoss = 0.0;
|
||||
if (idregion) {
|
||||
region = domain->get_region_by_id(idregion);
|
||||
if (!region) error->all(FLERR, "Region {} for fix electron/stopping does not exist", idregion);
|
||||
}
|
||||
|
||||
// need an occasional full neighbor list
|
||||
neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
|
||||
@ -176,18 +168,17 @@ void FixElectronStopping::post_force(int /*vflag*/)
|
||||
|
||||
int itype = type[i];
|
||||
double massone = (atom->rmass) ? atom->rmass[i] : atom->mass[itype];
|
||||
double v2 = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
|
||||
double v2 = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
|
||||
double energy = 0.5 * force->mvv2e * massone * v2;
|
||||
|
||||
if (energy < Ecut) continue;
|
||||
if (energy < elstop_ranges[0][0]) continue;
|
||||
if (energy > elstop_ranges[0][table_entries - 1]) error->one(FLERR,
|
||||
"Atom kinetic energy too high for fix electron/stopping");
|
||||
if (energy > elstop_ranges[0][table_entries - 1])
|
||||
error->one(FLERR, "Atom kinetic energy too high for fix electron/stopping");
|
||||
|
||||
if (iregion >= 0) {
|
||||
if (region) {
|
||||
// Only apply in the given region
|
||||
if (domain->regions[iregion]->match(x[i][0], x[i][1], x[i][2]) != 1)
|
||||
continue;
|
||||
if (region->match(x[i][0], x[i][1], x[i][2]) != 1) continue;
|
||||
}
|
||||
|
||||
// Binary search to find correct energy range
|
||||
@ -196,8 +187,10 @@ void FixElectronStopping::post_force(int /*vflag*/)
|
||||
while (true) {
|
||||
int ihalf = idown + (iup - idown) / 2;
|
||||
if (ihalf == idown) break;
|
||||
if (elstop_ranges[0][ihalf] < energy) idown = ihalf;
|
||||
else iup = ihalf;
|
||||
if (elstop_ranges[0][ihalf] < energy)
|
||||
idown = ihalf;
|
||||
else
|
||||
iup = ihalf;
|
||||
}
|
||||
|
||||
double Se_lo = elstop_ranges[itype][idown];
|
||||
@ -215,7 +208,7 @@ void FixElectronStopping::post_force(int /*vflag*/)
|
||||
f[i][1] += v[i][1] * factor;
|
||||
f[i][2] += v[i][2] * factor;
|
||||
|
||||
SeLoss += Se * vabs * dt; // very rough approx
|
||||
SeLoss += Se * vabs * dt; // very rough approx
|
||||
}
|
||||
}
|
||||
|
||||
@ -254,19 +247,17 @@ void FixElectronStopping::read_table(const char *file)
|
||||
ValueTokenizer values(line);
|
||||
elstop_ranges[0][nlines] = values.next_double();
|
||||
if (elstop_ranges[0][nlines] <= oldvalue)
|
||||
throw TokenizerException("energy values must be positive and in ascending order",line);
|
||||
throw TokenizerException("energy values must be positive and in ascending order", line);
|
||||
|
||||
oldvalue = elstop_ranges[0][nlines];
|
||||
for (int i = 1; i < ncol; ++i)
|
||||
elstop_ranges[i][nlines] = values.next_double();
|
||||
for (int i = 1; i < ncol; ++i) elstop_ranges[i][nlines] = values.next_double();
|
||||
|
||||
++nlines;
|
||||
}
|
||||
} catch (std::exception &e) {
|
||||
error->one(FLERR, "Problem parsing electron stopping data: {}", e.what());
|
||||
}
|
||||
if (nlines == 0)
|
||||
error->one(FLERR, "Did not find any data in electron/stopping table file");
|
||||
if (nlines == 0) error->one(FLERR, "Did not find any data in electron/stopping table file");
|
||||
|
||||
table_entries = nlines;
|
||||
}
|
||||
@ -281,8 +272,7 @@ void FixElectronStopping::grow_table()
|
||||
double **new_array;
|
||||
memory->create(new_array, ncol, new_maxlines, "electron/stopping:table");
|
||||
|
||||
for (int i = 0; i < ncol; i++)
|
||||
memcpy(new_array[i], elstop_ranges[i], maxlines*sizeof(double));
|
||||
for (int i = 0; i < ncol; i++) memcpy(new_array[i], elstop_ranges[i], maxlines * sizeof(double));
|
||||
|
||||
memory->destroy(elstop_ranges);
|
||||
elstop_ranges = new_array;
|
||||
|
||||
@ -52,8 +52,9 @@ class FixElectronStopping : public Fix {
|
||||
double **elstop_ranges; // [ 0][i]: energies
|
||||
// [>0][i]: stopping powers per type
|
||||
|
||||
int iregion; // region index if used, else -1
|
||||
int minneigh; // minimum number of neighbors
|
||||
char *idregion; // region id
|
||||
class Region *region; // region pointer if used, else NULL
|
||||
int minneigh; // minimum number of neighbors
|
||||
|
||||
class NeighList *list;
|
||||
};
|
||||
|
||||
@ -41,7 +41,8 @@ using namespace FixConst;
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixNumDiff::FixNumDiff(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg), id_pe(nullptr), pe(nullptr), numdiff_forces(nullptr), temp_x(nullptr), temp_f(nullptr)
|
||||
Fix(lmp, narg, arg), id_pe(nullptr), pe(nullptr), numdiff_forces(nullptr), temp_x(nullptr),
|
||||
temp_f(nullptr)
|
||||
{
|
||||
if (narg < 5) error->all(FLERR, "Illegal fix numdiff command");
|
||||
|
||||
@ -122,7 +123,7 @@ void FixNumDiff::init()
|
||||
kspace_compute_flag = 0;
|
||||
|
||||
if (utils::strmatch(update->integrate_style, "^respa")) {
|
||||
ilevel_respa = (dynamic_cast<Respa *>( update->integrate))->nlevels - 1;
|
||||
ilevel_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels - 1;
|
||||
if (respa_level >= 0) ilevel_respa = MIN(respa_level, ilevel_respa);
|
||||
}
|
||||
}
|
||||
@ -134,9 +135,9 @@ void FixNumDiff::setup(int vflag)
|
||||
if (utils::strmatch(update->integrate_style, "^verlet"))
|
||||
post_force(vflag);
|
||||
else {
|
||||
(dynamic_cast<Respa *>( update->integrate))->copy_flevel_f(ilevel_respa);
|
||||
(dynamic_cast<Respa *>(update->integrate))->copy_flevel_f(ilevel_respa);
|
||||
post_force_respa(vflag, ilevel_respa, 0);
|
||||
(dynamic_cast<Respa *>( update->integrate))->copy_f_flevel(ilevel_respa);
|
||||
(dynamic_cast<Respa *>(update->integrate))->copy_f_flevel(ilevel_respa);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -132,7 +132,7 @@ void FixNumDiffVirial::init()
|
||||
kspace_compute_flag = 0;
|
||||
|
||||
if (utils::strmatch(update->integrate_style, "^respa")) {
|
||||
ilevel_respa = (dynamic_cast<Respa *>( update->integrate))->nlevels - 1;
|
||||
ilevel_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels - 1;
|
||||
if (respa_level >= 0) ilevel_respa = MIN(respa_level, ilevel_respa);
|
||||
}
|
||||
}
|
||||
@ -144,9 +144,9 @@ void FixNumDiffVirial::setup(int vflag)
|
||||
if (utils::strmatch(update->integrate_style, "^verlet"))
|
||||
post_force(vflag);
|
||||
else {
|
||||
(dynamic_cast<Respa *>( update->integrate))->copy_flevel_f(ilevel_respa);
|
||||
(dynamic_cast<Respa *>(update->integrate))->copy_flevel_f(ilevel_respa);
|
||||
post_force_respa(vflag, ilevel_respa, 0);
|
||||
(dynamic_cast<Respa *>( update->integrate))->copy_f_flevel(ilevel_respa);
|
||||
(dynamic_cast<Respa *>(update->integrate))->copy_f_flevel(ilevel_respa);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -28,35 +27,36 @@
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
enum{NONE=-1,X=0,Y=1,Z=2,XYZMASK=3,MINUS=4,PLUS=0};
|
||||
enum { NONE = -1, X = 0, Y = 1, Z = 2, XYZMASK = 3, MINUS = 4, PLUS = 0 };
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixOneWay::FixOneWay(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
|
||||
FixOneWay::FixOneWay(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg), region(nullptr), idregion(nullptr)
|
||||
{
|
||||
direction = NONE;
|
||||
regionidx = 0;
|
||||
regionstr = nullptr;
|
||||
|
||||
if (narg < 6) error->all(FLERR,"Illegal fix oneway command");
|
||||
if (narg < 6) error->all(FLERR, "Illegal fix oneway command");
|
||||
|
||||
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
|
||||
if (nevery < 1) error->all(FLERR,"Illegal fix oneway command");
|
||||
nevery = utils::inumeric(FLERR, arg[3], false, lmp);
|
||||
if (nevery < 1) error->all(FLERR, "Illegal fix oneway command");
|
||||
|
||||
regionstr = utils::strdup(arg[4]);
|
||||
idregion = utils::strdup(arg[4]);
|
||||
if (!domain->get_region_by_id(idregion))
|
||||
error->all(FLERR, "Region {} for fix oneway does not exist", idregion);
|
||||
|
||||
if (strcmp(arg[5], "x") == 0) direction = X|PLUS;
|
||||
if (strcmp(arg[5], "X") == 0) direction = X|PLUS;
|
||||
if (strcmp(arg[5], "y") == 0) direction = Y|PLUS;
|
||||
if (strcmp(arg[5], "Y") == 0) direction = Y|PLUS;
|
||||
if (strcmp(arg[5], "z") == 0) direction = Z|PLUS;
|
||||
if (strcmp(arg[5], "Z") == 0) direction = Z|PLUS;
|
||||
if (strcmp(arg[5],"-x") == 0) direction = X|MINUS;
|
||||
if (strcmp(arg[5],"-X") == 0) direction = X|MINUS;
|
||||
if (strcmp(arg[5],"-y") == 0) direction = Y|MINUS;
|
||||
if (strcmp(arg[5],"-Y") == 0) direction = Y|MINUS;
|
||||
if (strcmp(arg[5],"-z") == 0) direction = Z|MINUS;
|
||||
if (strcmp(arg[5],"-Z") == 0) direction = Z|MINUS;
|
||||
if (strcmp(arg[5], "x") == 0) direction = X | PLUS;
|
||||
if (strcmp(arg[5], "X") == 0) direction = X | PLUS;
|
||||
if (strcmp(arg[5], "y") == 0) direction = Y | PLUS;
|
||||
if (strcmp(arg[5], "Y") == 0) direction = Y | PLUS;
|
||||
if (strcmp(arg[5], "z") == 0) direction = Z | PLUS;
|
||||
if (strcmp(arg[5], "Z") == 0) direction = Z | PLUS;
|
||||
if (strcmp(arg[5], "-x") == 0) direction = X | MINUS;
|
||||
if (strcmp(arg[5], "-X") == 0) direction = X | MINUS;
|
||||
if (strcmp(arg[5], "-y") == 0) direction = Y | MINUS;
|
||||
if (strcmp(arg[5], "-Y") == 0) direction = Y | MINUS;
|
||||
if (strcmp(arg[5], "-z") == 0) direction = Z | MINUS;
|
||||
if (strcmp(arg[5], "-Z") == 0) direction = Z | MINUS;
|
||||
|
||||
global_freq = nevery;
|
||||
}
|
||||
@ -65,7 +65,7 @@ FixOneWay::FixOneWay(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
|
||||
|
||||
FixOneWay::~FixOneWay()
|
||||
{
|
||||
delete[] regionstr;
|
||||
delete[] idregion;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -79,26 +79,24 @@ int FixOneWay::setmask()
|
||||
|
||||
void FixOneWay::init()
|
||||
{
|
||||
regionidx = domain->find_region(regionstr);
|
||||
if (regionidx < 0)
|
||||
error->all(FLERR,"Region for fix oneway does not exist");
|
||||
region = domain->get_region_by_id(idregion);
|
||||
if (!region) error->all(FLERR, "Region {} for fix oneway does not exist", idregion);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixOneWay::end_of_step()
|
||||
{
|
||||
Region *region = domain->regions[regionidx];
|
||||
region->prematch();
|
||||
|
||||
const int idx = direction & XYZMASK;
|
||||
const double * const * const x = atom->x;
|
||||
double * const * const v = atom->v;
|
||||
const double *const *const x = atom->x;
|
||||
double *const *const v = atom->v;
|
||||
const int *mask = atom->mask;
|
||||
const int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 0; i < nlocal; ++i) {
|
||||
if ((mask[i] & groupbit) && region->match(x[i][0],x[i][1],x[i][2])) {
|
||||
if ((mask[i] & groupbit) && region->match(x[i][0], x[i][1], x[i][2])) {
|
||||
if (direction & MINUS) {
|
||||
if (v[i][idx] > 0.0) v[i][idx] = -v[i][idx];
|
||||
} else {
|
||||
@ -107,4 +105,3 @@ void FixOneWay::end_of_step()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -34,8 +34,8 @@ class FixOneWay : public Fix {
|
||||
|
||||
protected:
|
||||
int direction;
|
||||
int regionidx;
|
||||
char *regionstr;
|
||||
class Region *region;
|
||||
char *idregion;
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
@ -297,7 +297,7 @@ void FixTTMGrid::read_electron_temperatures(const std::string &filename)
|
||||
try {
|
||||
ValueTokenizer values(utils::trim_comment(line));
|
||||
if (values.count() == 0) {
|
||||
; // ignore comment only lines
|
||||
; // ignore comment only lines
|
||||
} else if (values.count() == 4) {
|
||||
++nread;
|
||||
|
||||
@ -306,18 +306,18 @@ void FixTTMGrid::read_electron_temperatures(const std::string &filename)
|
||||
int iz = values.next_int();
|
||||
|
||||
if (ix < 0 || ix >= nxgrid || iy < 0 || iy >= nygrid || iz < 0 || iz >= nzgrid)
|
||||
throw TokenizerException("Fix ttm/grid invalid grid index in input","");
|
||||
throw TokenizerException("Fix ttm/grid invalid grid index in input", "");
|
||||
|
||||
if (ix >= nxlo_in && ix <= nxhi_in && iy >= nylo_in && iy <= nyhi_in
|
||||
&& iz >= nzlo_in && iz <= nzhi_in) {
|
||||
if (ix >= nxlo_in && ix <= nxhi_in && iy >= nylo_in && iy <= nyhi_in && iz >= nzlo_in &&
|
||||
iz <= nzhi_in) {
|
||||
T_electron[iz][iy][ix] = values.next_double();
|
||||
T_initial_set[iz][iy][ix] = 1;
|
||||
}
|
||||
} else {
|
||||
throw TokenizerException("Incorrect format in fix ttm electron grid file","");
|
||||
throw TokenizerException("Incorrect format in fix ttm electron grid file", "");
|
||||
}
|
||||
} catch (std::exception &e) {
|
||||
error->one(FLERR,e.what());
|
||||
error->one(FLERR, e.what());
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -356,9 +356,11 @@ void FixTTMGrid::write_electron_temperatures(const std::string &filename)
|
||||
FPout = fopen(filename.c_str(), "w");
|
||||
if (!FPout) error->one(FLERR, "Fix ttm/grid could not open output file");
|
||||
|
||||
fmt::print(FPout,"# DATE: {} UNITS: {} COMMENT: Electron temperature "
|
||||
"{}x{}x{} grid at step {}. Created by fix {}\n", utils::current_date(),
|
||||
update->unit_style, nxgrid, nygrid, nzgrid, update->ntimestep, style);
|
||||
fmt::print(FPout,
|
||||
"# DATE: {} UNITS: {} COMMENT: Electron temperature "
|
||||
"{}x{}x{} grid at step {}. Created by fix {}\n",
|
||||
utils::current_date(), update->unit_style, nxgrid, nygrid, nzgrid, update->ntimestep,
|
||||
style);
|
||||
}
|
||||
|
||||
gc->gather(GridComm::FIX, this, 1, sizeof(double), 1, nullptr, MPI_DOUBLE);
|
||||
|
||||
@ -48,16 +48,16 @@ class FixTTMGrid : public FixTTM {
|
||||
double memory_usage() override;
|
||||
|
||||
private:
|
||||
int ngridmine,ngridout;
|
||||
int nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in;
|
||||
int nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out;
|
||||
double delxinv,delyinv,delzinv;
|
||||
int ngridmine, ngridout;
|
||||
int nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in;
|
||||
int nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out;
|
||||
double delxinv, delyinv, delzinv;
|
||||
double skin_original;
|
||||
FILE *FPout;
|
||||
|
||||
class GridComm *gc;
|
||||
int ngc_buf1,ngc_buf2;
|
||||
double *gc_buf1,*gc_buf2;
|
||||
int ngc_buf1, ngc_buf2;
|
||||
double *gc_buf1, *gc_buf2;
|
||||
|
||||
void allocate_grid() override;
|
||||
void deallocate_grid() override;
|
||||
|
||||
@ -114,7 +114,7 @@ void FixViscousSphere::init()
|
||||
int max_respa = 0;
|
||||
|
||||
if (utils::strmatch(update->integrate_style, "^respa")) {
|
||||
ilevel_respa = max_respa = (dynamic_cast<Respa *>( update->integrate))->nlevels - 1;
|
||||
ilevel_respa = max_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels - 1;
|
||||
if (respa_level >= 0) ilevel_respa = MIN(respa_level, max_respa);
|
||||
}
|
||||
|
||||
@ -135,9 +135,9 @@ void FixViscousSphere::setup(int vflag)
|
||||
if (utils::strmatch(update->integrate_style, "^verlet"))
|
||||
post_force(vflag);
|
||||
else {
|
||||
(dynamic_cast<Respa *>( update->integrate))->copy_flevel_f(ilevel_respa);
|
||||
(dynamic_cast<Respa *>(update->integrate))->copy_flevel_f(ilevel_respa);
|
||||
post_force_respa(vflag, ilevel_respa, 0);
|
||||
(dynamic_cast<Respa *>( update->integrate))->copy_f_flevel(ilevel_respa);
|
||||
(dynamic_cast<Respa *>(update->integrate))->copy_f_flevel(ilevel_respa);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -23,6 +22,7 @@
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "math_extra.h"
|
||||
#include "math_special.h"
|
||||
#include "region.h"
|
||||
#include "respa.h"
|
||||
#include "update.h"
|
||||
@ -31,13 +31,14 @@
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
using MathSpecial::powint;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixWallRegionEES::FixWallRegionEES(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
Fix(lmp, narg, arg), idregion(nullptr), region(nullptr)
|
||||
{
|
||||
if (narg != 7) error->all(FLERR,"Illegal fix wall/region/ees command");
|
||||
if (narg != 7) error->all(FLERR, "Illegal fix wall/region/ees command");
|
||||
|
||||
scalar_flag = 1;
|
||||
vector_flag = 1;
|
||||
@ -49,15 +50,14 @@ FixWallRegionEES::FixWallRegionEES(LAMMPS *lmp, int narg, char **arg) :
|
||||
|
||||
// parse args
|
||||
|
||||
iregion = domain->find_region(arg[3]);
|
||||
if (iregion == -1)
|
||||
error->all(FLERR,"Region ID for fix wall/region/ees does not exist");
|
||||
region = domain->get_region_by_id(arg[3]);
|
||||
if (!region) error->all(FLERR, "Region {} for fix wall/region/ees does not exist", arg[3]);
|
||||
idregion = utils::strdup(arg[3]);
|
||||
epsilon = utils::numeric(FLERR,arg[4],false,lmp);
|
||||
sigma = utils::numeric(FLERR,arg[5],false,lmp);
|
||||
cutoff = utils::numeric(FLERR,arg[6],false,lmp);
|
||||
epsilon = utils::numeric(FLERR, arg[4], false, lmp);
|
||||
sigma = utils::numeric(FLERR, arg[5], false, lmp);
|
||||
cutoff = utils::numeric(FLERR, arg[6], false, lmp);
|
||||
|
||||
if (cutoff <= 0.0) error->all(FLERR,"Fix wall/region/ees cutoff <= 0.0");
|
||||
if (cutoff <= 0.0) error->all(FLERR, "Fix wall/region/ees cutoff <= 0.0");
|
||||
|
||||
eflag = 0;
|
||||
ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
|
||||
@ -67,7 +67,7 @@ FixWallRegionEES::FixWallRegionEES(LAMMPS *lmp, int narg, char **arg) :
|
||||
|
||||
FixWallRegionEES::~FixWallRegionEES()
|
||||
{
|
||||
delete [] idregion;
|
||||
delete[] idregion;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -87,13 +87,11 @@ void FixWallRegionEES::init()
|
||||
{
|
||||
// set index and check validity of region
|
||||
|
||||
iregion = domain->find_region(idregion);
|
||||
if (iregion == -1)
|
||||
error->all(FLERR,"Region ID for fix wall/region/ees does not exist");
|
||||
region = domain->get_region_by_id(idregion);
|
||||
if (!region) error->all(FLERR, "Region {} for fix wall/region/ees does not exist", idregion);
|
||||
|
||||
avec = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
|
||||
if (!avec)
|
||||
error->all(FLERR,"Fix wall/region/ees requires atom style ellipsoid");
|
||||
avec = dynamic_cast<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
|
||||
if (!avec) error->all(FLERR, "Fix wall/region/ees requires atom style ellipsoid");
|
||||
|
||||
// check that all particles are finite-size ellipsoids
|
||||
// no point particles allowed, spherical is OK
|
||||
@ -105,33 +103,33 @@ void FixWallRegionEES::init()
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
if (ellipsoid[i] < 0)
|
||||
error->one(FLERR,"Fix wall/region/ees requires extended particles");
|
||||
error->one(FLERR, "Fix wall/region/ees requires only extended particles");
|
||||
|
||||
// setup coefficients
|
||||
|
||||
coeff1 = ( 2. / 4725. ) * epsilon * pow(sigma,12.0);
|
||||
coeff2 = ( 1. / 24. ) * epsilon * pow(sigma,6.0);
|
||||
coeff3 = ( 2. / 315. ) * epsilon * pow(sigma,12.0);
|
||||
coeff4 = ( 1. / 3. ) * epsilon * pow(sigma,6.0);
|
||||
coeff5 = ( 4. / 315. ) * epsilon * pow(sigma,12.0);
|
||||
coeff6 = ( 1. / 12. ) * epsilon * pow(sigma,6.0);
|
||||
coeff1 = (2.0 / 4725.0) * epsilon * powint(sigma, 12);
|
||||
coeff2 = (1.0 / 24.0) * epsilon * powint(sigma, 6);
|
||||
coeff3 = (2.0 / 315.0) * epsilon * powint(sigma, 12);
|
||||
coeff4 = (1.0 / 3.0) * epsilon * powint(sigma, 6);
|
||||
coeff5 = (4.0 / 315.0) * epsilon * powint(sigma, 12);
|
||||
coeff6 = (1.0 / 12.0) * epsilon * powint(sigma, 6);
|
||||
offset = 0;
|
||||
|
||||
|
||||
if (utils::strmatch(update->integrate_style,"^respa"))
|
||||
nlevels_respa = (dynamic_cast<Respa *>( update->integrate))->nlevels;
|
||||
if (utils::strmatch(update->integrate_style, "^respa"))
|
||||
nlevels_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixWallRegionEES::setup(int vflag)
|
||||
{
|
||||
if (utils::strmatch(update->integrate_style,"^verlet"))
|
||||
if (utils::strmatch(update->integrate_style, "^respa")) {
|
||||
auto respa = dynamic_cast<Respa *>(update->integrate);
|
||||
respa->copy_flevel_f(nlevels_respa - 1);
|
||||
post_force_respa(vflag, nlevels_respa - 1, 0);
|
||||
respa->copy_f_flevel(nlevels_respa - 1);
|
||||
} else {
|
||||
post_force(vflag);
|
||||
else {
|
||||
(dynamic_cast<Respa *>( update->integrate))->copy_flevel_f(nlevels_respa-1);
|
||||
post_force_respa(vflag,nlevels_respa-1,0);
|
||||
(dynamic_cast<Respa *>( update->integrate))->copy_f_flevel(nlevels_respa-1);
|
||||
}
|
||||
}
|
||||
|
||||
@ -139,7 +137,7 @@ void FixWallRegionEES::setup(int vflag)
|
||||
|
||||
void FixWallRegionEES::min_setup(int vflag)
|
||||
{
|
||||
post_force(vflag);
|
||||
post_force(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -149,8 +147,8 @@ void FixWallRegionEES::post_force(int /*vflag*/)
|
||||
//sth is needed here, but I dont know what
|
||||
//that is calculation of sn
|
||||
|
||||
int i,m,n;
|
||||
double rinv,fx,fy,fz,sn,tooclose[3];
|
||||
int i, m, n;
|
||||
double rinv, fx, fy, fz, sn, tooclose[3];
|
||||
|
||||
eflag = 0;
|
||||
ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
|
||||
@ -165,7 +163,6 @@ void FixWallRegionEES::post_force(int /*vflag*/)
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
Region *region = domain->regions[iregion];
|
||||
region->prematch();
|
||||
|
||||
int onflag = 0;
|
||||
@ -176,33 +173,34 @@ void FixWallRegionEES::post_force(int /*vflag*/)
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (!region->match(x[i][0],x[i][1],x[i][2])) {
|
||||
if (!region->match(x[i][0], x[i][1], x[i][2])) {
|
||||
onflag = 1;
|
||||
continue;
|
||||
}
|
||||
|
||||
double A[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
|
||||
double tempvec[3]= {0,0,0};
|
||||
double A[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
|
||||
double tempvec[3] = {0, 0, 0};
|
||||
double sn2 = 0.0;
|
||||
double nhat[3] = {0,0,0};
|
||||
double* shape = bonus[ellipsoid[i]].shape;;
|
||||
MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat,A);
|
||||
double nhat[3] = {0, 0, 0};
|
||||
double *shape = bonus[ellipsoid[i]].shape;
|
||||
;
|
||||
MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat, A);
|
||||
|
||||
for (int which = 0 ; which < 3; which ++) {//me
|
||||
nhat[which]=1;
|
||||
nhat[(which+1)%3] = 0 ;
|
||||
nhat[(which+2)%3] = 0 ;
|
||||
sn2 = 0 ;
|
||||
MathExtra::transpose_matvec(A,nhat,tempvec);
|
||||
for (int k = 0; k<3; k++) {
|
||||
for (int which = 0; which < 3; which++) { //me
|
||||
nhat[which] = 1;
|
||||
nhat[(which + 1) % 3] = 0;
|
||||
nhat[(which + 2) % 3] = 0;
|
||||
sn2 = 0;
|
||||
MathExtra::transpose_matvec(A, nhat, tempvec);
|
||||
for (int k = 0; k < 3; k++) {
|
||||
tempvec[k] *= shape[k];
|
||||
sn2 += tempvec[k]*tempvec[k];
|
||||
sn2 += tempvec[k] * tempvec[k];
|
||||
}
|
||||
sn = sqrt(sn2);
|
||||
tooclose[which] = sn;
|
||||
}
|
||||
|
||||
n = region->surface(x[i][0],x[i][1],x[i][2],cutoff);
|
||||
n = region->surface(x[i][0], x[i][1], x[i][2], cutoff);
|
||||
|
||||
for (m = 0; m < n; m++) {
|
||||
|
||||
@ -212,12 +210,13 @@ void FixWallRegionEES::post_force(int /*vflag*/)
|
||||
} else if (region->contact[m].dely != 0 && region->contact[m].r <= tooclose[1]) {
|
||||
onflag = 1;
|
||||
continue;
|
||||
} else if (region->contact[m].delz !=0 && region->contact[m].r <= tooclose[2]) {
|
||||
} else if (region->contact[m].delz != 0 && region->contact[m].r <= tooclose[2]) {
|
||||
onflag = 1;
|
||||
continue;
|
||||
} else rinv = 1.0/region->contact[m].r;
|
||||
} else
|
||||
rinv = 1.0 / region->contact[m].r;
|
||||
|
||||
ees(m,i);
|
||||
ees(m, i);
|
||||
|
||||
ewall[0] += eng;
|
||||
fx = fwall * region->contact[m].delx * rinv;
|
||||
@ -237,15 +236,15 @@ void FixWallRegionEES::post_force(int /*vflag*/)
|
||||
}
|
||||
}
|
||||
|
||||
if (onflag) error->one(FLERR,"Particle on or inside surface of region "
|
||||
"used in fix wall/region/ees");
|
||||
if (onflag)
|
||||
error->one(FLERR, "Particle on or inside surface of region used in fix wall/region/ees");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixWallRegionEES::post_force_respa(int vflag, int ilevel, int /*iloop*/)
|
||||
{
|
||||
if (ilevel == nlevels_respa-1) post_force(vflag);
|
||||
if (ilevel == nlevels_respa - 1) post_force(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -264,7 +263,7 @@ double FixWallRegionEES::compute_scalar()
|
||||
// only sum across procs one time
|
||||
|
||||
if (eflag == 0) {
|
||||
MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(ewall, ewall_all, 4, MPI_DOUBLE, MPI_SUM, world);
|
||||
eflag = 1;
|
||||
}
|
||||
return ewall_all[0];
|
||||
@ -279,10 +278,10 @@ double FixWallRegionEES::compute_vector(int n)
|
||||
// only sum across procs one time
|
||||
|
||||
if (eflag == 0) {
|
||||
MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(ewall, ewall_all, 4, MPI_DOUBLE, MPI_SUM, world);
|
||||
eflag = 1;
|
||||
}
|
||||
return ewall_all[n+1];
|
||||
return ewall_all[n + 1];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -292,24 +291,23 @@ double FixWallRegionEES::compute_vector(int n)
|
||||
|
||||
void FixWallRegionEES::ees(int m, int i)
|
||||
{
|
||||
Region *region = domain->regions[iregion];
|
||||
region->prematch();
|
||||
|
||||
double delta, delta2, delta3, delta4, delta5, delta6;
|
||||
double sigman, sigman2 , sigman3, sigman4, sigman5, sigman6;
|
||||
double hhss, hhss2, hhss4, hhss7, hhss8; //h^2 - s_n^2
|
||||
double hps; //h+s_n
|
||||
double hms; //h-s_n
|
||||
double sigman, sigman2, sigman3, sigman4, sigman5, sigman6;
|
||||
double hhss, hhss2, hhss4, hhss7, hhss8; //h^2 - s_n^2
|
||||
double hps; //h+s_n
|
||||
double hms; //h-s_n
|
||||
double twall;
|
||||
|
||||
double A[3][3], nhat[3], SAn[3], that[3];
|
||||
|
||||
double tempvec[3]= {0,0,0};
|
||||
double tempvec2[3]= {0,0,0};
|
||||
double tempvec[3] = {0, 0, 0};
|
||||
double tempvec2[3] = {0, 0, 0};
|
||||
|
||||
double Lx[3][3] = {{0,0,0},{0,0,-1},{0,1,0}};
|
||||
double Ly[3][3] = {{0,0,1},{0,0,0},{-1,0,0}};
|
||||
double Lz[3][3] = {{0,-1,0},{1,0,0},{0,0,0}};
|
||||
double Lx[3][3] = {{0, 0, 0}, {0, 0, -1}, {0, 1, 0}};
|
||||
double Ly[3][3] = {{0, 0, 1}, {0, 0, 0}, {-1, 0, 0}};
|
||||
double Lz[3][3] = {{0, -1, 0}, {1, 0, 0}, {0, 0, 0}};
|
||||
|
||||
nhat[0] = region->contact[m].delx / region->contact[m].r;
|
||||
nhat[1] = region->contact[m].dely / region->contact[m].r;
|
||||
@ -318,14 +316,15 @@ void FixWallRegionEES::ees(int m, int i)
|
||||
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
|
||||
double* shape = bonus[ellipsoid[i]].shape;;
|
||||
MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat,A);
|
||||
double *shape = bonus[ellipsoid[i]].shape;
|
||||
;
|
||||
MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat, A);
|
||||
|
||||
sigman2 = 0.0;
|
||||
MathExtra::transpose_matvec(A,nhat,tempvec);
|
||||
for (int k = 0; k<3; k++) {
|
||||
MathExtra::transpose_matvec(A, nhat, tempvec);
|
||||
for (int k = 0; k < 3; k++) {
|
||||
tempvec[k] *= shape[k];
|
||||
sigman2 += tempvec[k]*tempvec[k];
|
||||
sigman2 += tempvec[k] * tempvec[k];
|
||||
SAn[k] = tempvec[k];
|
||||
}
|
||||
|
||||
@ -337,14 +336,14 @@ void FixWallRegionEES::ees(int m, int i)
|
||||
sigman5 = sigman4 * sigman;
|
||||
sigman6 = sigman3 * sigman3;
|
||||
|
||||
delta2 = delta * delta;
|
||||
delta2 = delta * delta;
|
||||
delta3 = delta2 * delta;
|
||||
delta4 = delta2 * delta2;
|
||||
delta5 = delta3 * delta2;
|
||||
delta6 = delta3 * delta3;
|
||||
|
||||
hhss = delta2 - sigman2;
|
||||
hhss2 = hhss * hhss;
|
||||
hhss2 = hhss * hhss;
|
||||
hhss4 = hhss2 * hhss2;
|
||||
hhss8 = hhss4 * hhss4;
|
||||
hhss7 = hhss4 * hhss2 * hhss;
|
||||
@ -352,31 +351,31 @@ void FixWallRegionEES::ees(int m, int i)
|
||||
hps = delta + sigman;
|
||||
hms = delta - sigman;
|
||||
|
||||
fwall = -1*coeff4/hhss2 + coeff3
|
||||
* (21*delta6 + 63*delta4*sigman2 + 27*delta2*sigman4 + sigman6) / hhss8;
|
||||
fwall = -1 * coeff4 / hhss2 +
|
||||
coeff3 * (21 * delta6 + 63 * delta4 * sigman2 + 27 * delta2 * sigman4 + sigman6) / hhss8;
|
||||
|
||||
eng = -1*coeff2 * (4*delta/sigman2/hhss + 2*log(hms/hps)/sigman3) +
|
||||
coeff1 * (35*delta5 + 70*delta3*sigman2 + 15*delta*sigman4) / hhss7;
|
||||
eng = -1 * coeff2 * (4 * delta / sigman2 / hhss + 2 * log(hms / hps) / sigman3) +
|
||||
coeff1 * (35 * delta5 + 70 * delta3 * sigman2 + 15 * delta * sigman4) / hhss7;
|
||||
|
||||
twall = coeff6 * (6*delta3/sigman4/hhss2 - 10*delta/sigman2/hhss2
|
||||
+ 3*log(hms/hps)/sigman5)
|
||||
+ coeff5 * (21.*delta5 + 30.*delta3*sigman2 + 5.*delta*sigman4) / hhss8;
|
||||
twall = coeff6 *
|
||||
(6 * delta3 / sigman4 / hhss2 - 10 * delta / sigman2 / hhss2 +
|
||||
3 * log(hms / hps) / sigman5) +
|
||||
coeff5 * (21. * delta5 + 30. * delta3 * sigman2 + 5. * delta * sigman4) / hhss8;
|
||||
|
||||
MathExtra::matvec(Lx,nhat,tempvec);
|
||||
MathExtra::transpose_matvec(A,tempvec,tempvec2);
|
||||
for (int k = 0; k<3; k++) tempvec2[k] *= shape[k];
|
||||
that[0] = MathExtra::dot3(SAn,tempvec2);
|
||||
|
||||
MathExtra::matvec(Ly,nhat,tempvec);
|
||||
MathExtra::transpose_matvec(A,tempvec,tempvec2);
|
||||
for (int k = 0; k<3; k++) tempvec2[k] *= shape[k];
|
||||
that[1] = MathExtra::dot3(SAn,tempvec2);
|
||||
|
||||
MathExtra::matvec(Lz,nhat,tempvec);
|
||||
MathExtra::transpose_matvec(A,tempvec,tempvec2);
|
||||
MathExtra::matvec(Lx, nhat, tempvec);
|
||||
MathExtra::transpose_matvec(A, tempvec, tempvec2);
|
||||
for (int k = 0; k < 3; k++) tempvec2[k] *= shape[k];
|
||||
that[2] = MathExtra::dot3(SAn,tempvec2);
|
||||
that[0] = MathExtra::dot3(SAn, tempvec2);
|
||||
|
||||
for (int j = 0; j<3 ; j++)
|
||||
torque[j] = twall * that[j];
|
||||
MathExtra::matvec(Ly, nhat, tempvec);
|
||||
MathExtra::transpose_matvec(A, tempvec, tempvec2);
|
||||
for (int k = 0; k < 3; k++) tempvec2[k] *= shape[k];
|
||||
that[1] = MathExtra::dot3(SAn, tempvec2);
|
||||
|
||||
MathExtra::matvec(Lz, nhat, tempvec);
|
||||
MathExtra::transpose_matvec(A, tempvec, tempvec2);
|
||||
for (int k = 0; k < 3; k++) tempvec2[k] *= shape[k];
|
||||
that[2] = MathExtra::dot3(SAn, tempvec2);
|
||||
|
||||
for (int j = 0; j < 3; j++) torque[j] = twall * that[j];
|
||||
}
|
||||
|
||||
@ -41,12 +41,12 @@ class FixWallRegionEES : public Fix {
|
||||
private:
|
||||
class AtomVecEllipsoid *avec;
|
||||
|
||||
int iregion;
|
||||
double epsilon, sigma, cutoff;
|
||||
int eflag;
|
||||
double ewall[4], ewall_all[4];
|
||||
int nlevels_respa;
|
||||
char *idregion;
|
||||
class Region *region;
|
||||
|
||||
double coeff1, coeff2, coeff3, coeff4, offset;
|
||||
double coeff5, coeff6;
|
||||
|
||||
@ -89,7 +89,7 @@ void BondFENENM::compute(int eflag, int vflag)
|
||||
|
||||
fbond = -k[type] / rlogarg;
|
||||
// force from n-m term
|
||||
if (rsq < sigma[type]*sigma[type]) {
|
||||
if (rsq < sigma[type] * sigma[type]) {
|
||||
r = sqrt(rsq);
|
||||
fbond += epsilon[type] * (nn[type] * mm[type] / (nn[type] - mm[type])) *
|
||||
(pow(sigma[type] / r, nn[type]) - pow(sigma[type] / r, mm[type])) / rsq;
|
||||
@ -99,7 +99,7 @@ void BondFENENM::compute(int eflag, int vflag)
|
||||
|
||||
if (eflag) {
|
||||
ebond = -0.5 * k[type] * r0sq * log(rlogarg);
|
||||
if (rsq < sigma[type]*sigma[type])
|
||||
if (rsq < sigma[type] * sigma[type])
|
||||
ebond += (epsilon[type] / (nn[type] - mm[type])) *
|
||||
(mm[type] * pow(sigma[type] / r, nn[type]) - nn[type] * pow(sigma[type] / r, mm[type]));
|
||||
}
|
||||
@ -257,7 +257,7 @@ double BondFENENM::single(int type, double rsq, int /*i*/, int /*j*/, double &ff
|
||||
double eng = -0.5 * k[type] * r0sq * log(rlogarg);
|
||||
fforce = -k[type] / rlogarg;
|
||||
|
||||
if (rsq < sigma[type]*sigma[type]) {
|
||||
if (rsq < sigma[type] * sigma[type]) {
|
||||
r = sqrt(rsq);
|
||||
fforce += epsilon[type] * (nn[type] * mm[type] / (nn[type] - mm[type])) *
|
||||
(pow(sigma[type] / r, nn[type]) - pow(sigma[type] / r, mm[type])) / rsq;
|
||||
|
||||
@ -34,7 +34,8 @@ using namespace MathConst;
|
||||
BondGaussian::BondGaussian(LAMMPS *lmp) :
|
||||
Bond(lmp), nterms(nullptr), bond_temperature(nullptr), alpha(nullptr), width(nullptr),
|
||||
r0(nullptr)
|
||||
{}
|
||||
{
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
|
||||
@ -25,6 +25,7 @@
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
@ -39,6 +40,7 @@ DihedralNHarmonic::DihedralNHarmonic(LAMMPS *lmp) : Dihedral(lmp)
|
||||
{
|
||||
writedata = 1;
|
||||
a = nullptr;
|
||||
born_matrix_enable = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -336,3 +338,62 @@ void DihedralNHarmonic::write_data(FILE *fp)
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------*/
|
||||
|
||||
void DihedralNHarmonic::born_matrix(int nd, int i1, int i2, int i3, int i4,
|
||||
double &dudih, double &du2dih) {
|
||||
int i,type;
|
||||
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
|
||||
double c,ax,ay,az,bx,by,bz,rasq,rbsq,ra2inv,rb2inv,rabinv;
|
||||
|
||||
int **dihedrallist = neighbor->dihedrallist;
|
||||
double **x = atom->x;
|
||||
|
||||
type = dihedrallist[nd][4];
|
||||
|
||||
vb1x = x[i1][0] - x[i2][0];
|
||||
vb1y = x[i1][1] - x[i2][1];
|
||||
vb1z = x[i1][2] - x[i2][2];
|
||||
|
||||
vb2x = x[i3][0] - x[i2][0];
|
||||
vb2y = x[i3][1] - x[i2][1];
|
||||
vb2z = x[i3][2] - x[i2][2];
|
||||
|
||||
vb2xm = -vb2x;
|
||||
vb2ym = -vb2y;
|
||||
vb2zm = -vb2z;
|
||||
|
||||
vb3x = x[i4][0] - x[i3][0];
|
||||
vb3y = x[i4][1] - x[i3][1];
|
||||
vb3z = x[i4][2] - x[i3][2];
|
||||
|
||||
// c,s calculation
|
||||
|
||||
ax = vb1y*vb2zm - vb1z*vb2ym;
|
||||
ay = vb1z*vb2xm - vb1x*vb2zm;
|
||||
az = vb1x*vb2ym - vb1y*vb2xm;
|
||||
bx = vb3y*vb2zm - vb3z*vb2ym;
|
||||
by = vb3z*vb2xm - vb3x*vb2zm;
|
||||
bz = vb3x*vb2ym - vb3y*vb2xm;
|
||||
|
||||
rasq = ax*ax + ay*ay + az*az;
|
||||
rbsq = bx*bx + by*by + bz*bz;
|
||||
|
||||
ra2inv = rb2inv = 0.0;
|
||||
if (rasq > 0) ra2inv = 1.0/rasq;
|
||||
if (rbsq > 0) rb2inv = 1.0/rbsq;
|
||||
rabinv = sqrt(ra2inv*rb2inv);
|
||||
|
||||
c = (ax*bx + ay*by + az*bz)*rabinv;
|
||||
|
||||
dudih = 0.0;
|
||||
du2dih = 0.0;
|
||||
for (i = 1; i < nterms[type]; i++) {
|
||||
dudih += this->a[type][i]*i*pow(c,i-1);
|
||||
}
|
||||
|
||||
for (i = 2; i < nterms[type]; i++) {
|
||||
du2dih += this->a[type][i]*i*(i-1)*pow(c, i-2);
|
||||
}
|
||||
}
|
||||
|
||||
Some files were not shown because too many files have changed in this diff Show More
Reference in New Issue
Block a user