forked from lijiext/lammps
Merge pull request #2019 from z-gong/viscosity
Fix and compute styles for calculating viscosity with periodic perturbation method
This commit is contained in:
Axel Kohlmeyer
GitHub
commit
1329603184
|
|
@ -161,5 +161,6 @@ KOKKOS, o = USER-OMP, t = OPT.
|
|||
* :doc:`torque/chunk <compute_torque_chunk>`
|
||||
* :doc:`vacf <compute_vacf>`
|
||||
* :doc:`vcm/chunk <compute_vcm_chunk>`
|
||||
* :doc:`viscosity/cos <compute_viscosity_cos>`
|
||||
* :doc:`voronoi/atom <compute_voronoi_atom>`
|
||||
* :doc:`xrd <compute_xrd>`
|
||||
|
|
|
|||
|
|
@ -22,6 +22,7 @@ OPT.
|
|||
.. table_from_list::
|
||||
:columns: 5
|
||||
|
||||
* :doc:`accelerate/cos <fix_accelerate_cos>`
|
||||
* :doc:`adapt <fix_adapt>`
|
||||
* :doc:`adapt/fep <fix_adapt_fep>`
|
||||
* :doc:`addforce <fix_addforce>`
|
||||
|
|
|
|||
|
|
@ -1,10 +1,11 @@
|
|||
Calculate viscosity
|
||||
===================
|
||||
|
||||
The shear viscosity eta of a fluid can be measured in at least 5 ways
|
||||
The shear viscosity eta of a fluid can be measured in at least 6 ways
|
||||
using various options in LAMMPS. See the examples/VISCOSITY directory
|
||||
for scripts that implement the 5 methods discussed here for a simple
|
||||
Lennard-Jones fluid model. Also, see the :doc:`Howto kappa <Howto_kappa>` doc page for an analogous discussion for
|
||||
Lennard-Jones fluid model and 1 method for SPC/E water model.
|
||||
Also, see the :doc:`Howto kappa <Howto_kappa>` doc page for an analogous discussion for
|
||||
thermal conductivity.
|
||||
|
||||
Eta is a measure of the propensity of a fluid to transmit momentum in
|
||||
|
|
@ -130,9 +131,25 @@ time-integrated momentum fluxes play the role of Cartesian
|
|||
coordinates, whose mean-square displacement increases linearly
|
||||
with time at sufficiently long times.
|
||||
|
||||
The sixth is periodic perturbation method. It is also a non-equilibrium MD method.
|
||||
However, instead of measure the momentum flux in response of applied velocity gradient,
|
||||
it measures the velocity profile in response of applied stress.
|
||||
A cosine-shaped periodic acceleration is added to the system via the
|
||||
:doc:`fix accelerate/cos <fix_accelerate_cos>` command,
|
||||
and the :doc:`compute viscosity/cos<compute_viscosity_cos>` command is used to monitor the
|
||||
generated velocity profile and remove the velocity bias before thermostatting.
|
||||
|
||||
.. note::
|
||||
|
||||
An article by :ref:`(Hess) <Hess3>` discussed the accuracy and efficiency of these methods.
|
||||
|
||||
----------
|
||||
|
||||
.. _Daivis-viscosity:
|
||||
|
||||
**(Daivis and Todd)** Daivis and Todd, Nonequilibrium Molecular Dynamics (book),
|
||||
Cambridge University Press, https://doi.org/10.1017/9781139017848, (2017).
|
||||
|
||||
.. _Hess3:
|
||||
|
||||
**(Hess)** Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209-217.
|
||||
|
|
|
|||
|
|
@ -307,6 +307,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
|
|||
* :doc:`torque/chunk <compute_torque_chunk>` - torque applied on each chunk
|
||||
* :doc:`vacf <compute_vacf>` - velocity auto-correlation function of group of atoms
|
||||
* :doc:`vcm/chunk <compute_vcm_chunk>` - velocity of center-of-mass for each chunk
|
||||
* :doc:`viscosity/cos <compute_viscosity_cos>` - velocity profile under cosine-shaped acceleration
|
||||
* :doc:`voronoi/atom <compute_voronoi_atom>` - Voronoi volume and neighbors for each atom
|
||||
* :doc:`xrd <compute_xrd>` - x-ray diffraction intensity on a mesh of reciprocal lattice nodes
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,156 @@
|
|||
.. index:: compute viscosity/cos
|
||||
|
||||
compute viscosity/cos command
|
||||
=============================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
compute ID group-ID viscosity/cos
|
||||
|
||||
* ID, group-ID are documented in :doc:`compute <compute>` command
|
||||
* viscosity/cos = style name of this compute command
|
||||
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute cos all viscosity/cos
|
||||
variable V equal c_cos[7]
|
||||
variable A equal 0.02E-5
|
||||
variable density equal density
|
||||
variable lz equal lz
|
||||
variable reciprocalViscosity equal v_V/${A}/v_density*39.4784/v_lz/v_lz*100
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
Define a computation that calculates the velocity amplitude of a group of atoms
|
||||
with an cosine-shaped velocity profile and the temperature of them
|
||||
after subtracting out the velocity profile before computing the kinetic energy.
|
||||
A compute of this style can be used by any command that computes a temperature,
|
||||
e.g. :doc:`thermo_modify <thermo_modify>`, :doc:`fix npt <fix_nh>`, etc.
|
||||
|
||||
This command together with :doc:`fix_accelerate/cos<fix_accelerate_cos>`
|
||||
enables viscosity calculation with periodic perturbation method,
|
||||
as described by :ref:`Hess<Hess1>`.
|
||||
An acceleration along the x-direction is applied to the simulation system
|
||||
by using :doc:`fix_accelerate/cos<fix_accelerate_cos>` command.
|
||||
The acceleration is a periodic function along the z-direction:
|
||||
|
||||
.. math::
|
||||
|
||||
a_{x}(z) = A \cos \left(\frac{2 \pi z}{l_{z}}\right)
|
||||
|
||||
where :math:`A` is the acceleration amplitude, :math:`l_z` is the z-length
|
||||
of the simulation box. At steady state, the acceleration generates
|
||||
a velocity profile:
|
||||
|
||||
.. math::
|
||||
|
||||
v_{x}(z) = V \cos \left(\frac{2 \pi z}{l_{z}}\right)
|
||||
|
||||
The generated velocity amplitude :math:`V` is related to the
|
||||
shear viscosity :math:`\eta` by:
|
||||
|
||||
.. math::
|
||||
|
||||
V = \frac{A \rho}{\eta}\left(\frac{l_{z}}{2 \pi}\right)^{2}
|
||||
|
||||
|
||||
and it can be obtained from ensemble average of the velocity profile:
|
||||
|
||||
.. math::
|
||||
|
||||
V = \frac{\sum_i 2 m_{i} v_{i, x} \cos \left(\frac{2 \pi z_i}{l_{z}}\right)}{\sum_i m_{i}}
|
||||
|
||||
|
||||
where :math:`m_i`, :math:`v_{i,x}` and :math:`z_i` are the mass,
|
||||
x-component velocity and z coordinate of a particle.
|
||||
|
||||
After the cosine-shaped collective velocity in :math:`x` direction
|
||||
has been subtracted for each atom, the temperature is calculated by the formula
|
||||
KE = dim/2 N k T, 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.
|
||||
|
||||
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.
|
||||
|
||||
The number of atoms contributing to the temperature is assumed to be
|
||||
constant for the duration of the run; use the *dynamic* option of the
|
||||
:doc:`compute_modify <compute_modify>` command if this is not the case.
|
||||
However, in order to get meaningful result, the group ID of this compute should be all.
|
||||
|
||||
The removal of the cosine-shaped velocity component by this command is
|
||||
essentially computing the temperature after a "bias" has been removed
|
||||
from the velocity of the atoms. If this compute is used with a fix
|
||||
command that performs thermostatting then this bias will be subtracted
|
||||
from each atom, thermostatting of the remaining thermal velocity will
|
||||
be performed, and the bias will be added back in. Thermostatting
|
||||
fixes that work in this way include :doc:`fix nvt <fix_nh>`, :doc:`fix temp/rescale <fix_temp_rescale>`, :doc:`fix temp/berendsen <fix_temp_berendsen>`, and :doc:`fix langevin <fix_langevin>`.
|
||||
|
||||
This compute subtracts out degrees-of-freedom due to fixes that
|
||||
constrain molecular motion, such as :doc:`fix shake <fix_shake>` and
|
||||
:doc:`fix rigid <fix_rigid>`. This means the temperature of groups of
|
||||
atoms that include these constraints will be computed correctly. If
|
||||
needed, the subtracted degrees-of-freedom can be altered using the
|
||||
*extra* option of the :doc:`compute_modify <compute_modify>` command.
|
||||
|
||||
See the :doc:`Howto thermostat <Howto_thermostat>` doc page for a
|
||||
discussion of different ways to compute temperature and perform
|
||||
thermostatting.
|
||||
|
||||
----------
|
||||
|
||||
**Output info:**
|
||||
|
||||
This compute calculates a global scalar (the temperature) and a global
|
||||
vector of length 7, which can be accessed by indices 1-7.
|
||||
The first 6 elements of the vector are the KE tensor,
|
||||
and the 7-th is the cosine-shaped velocity amplitude :math:`V`,
|
||||
which can be used to calculate the reciprocal viscosity, as shown in the example.
|
||||
These values can be used by any command that uses global scalar or
|
||||
vector values from a compute as input.
|
||||
See the :doc:`Howto output <Howto_output>` doc page for an overview of LAMMPS output options.
|
||||
|
||||
The scalar value calculated by this compute is "intensive". The
|
||||
first 6 elements of vector values are "extensive",
|
||||
and the 7-th element of vector values is "intensive".
|
||||
|
||||
The scalar value will be in temperature :doc:`units <units>`. The
|
||||
first 6 elements of vector values will be in energy :doc:`units <units>`.
|
||||
The 7-th element of vector value will be in velocity :doc:`units <units>`.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This command is only available when LAMMPS was built with the USER-MISC package.
|
||||
Since this compute depends on :doc:`fix accelerate/cos <fix_accelerate_cos>` which can
|
||||
only work for 3d systems, it cannot be used for 2d systems.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`fix accelerate/cos <fix_accelerate_cos>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
none
|
||||
|
||||
----------
|
||||
|
||||
.. _Hess1:
|
||||
|
||||
**(Hess)** Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209-217.
|
||||
|
|
@ -165,6 +165,7 @@ The individual style names on the :doc:`Commands fix <Commands_fix>` doc
|
|||
page are followed by one or more of (g,i,k,o,t) to indicate which
|
||||
accelerated styles exist.
|
||||
|
||||
* :doc:`accelerate/cos <fix_accelerate_cos>` - apply cosine-shaped acceleration to atoms
|
||||
* :doc:`adapt <fix_adapt>` - change a simulation parameter over time
|
||||
* :doc:`adapt/fep <fix_adapt_fep>` - enhanced version of fix adapt
|
||||
* :doc:`addforce <fix_addforce>` - add a force to each atom
|
||||
|
|
|
|||
|
|
@ -0,0 +1,104 @@
|
|||
.. index:: fix accelerate/cos
|
||||
|
||||
fix accelerate/cos command
|
||||
==========================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID accelerate value
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* accelerate/cos = style name of this fix command
|
||||
* value = amplitude of acceleration (in unit of force/mass)
|
||||
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix 1 all accelerate/cos 0.02e-5
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
Give each atom a acceleration in x-direction based on its z coordinate.
|
||||
The acceleration is a periodic function along the z-direction:
|
||||
|
||||
.. math::
|
||||
|
||||
a_{x}(z) = A \cos \left(\frac{2 \pi z}{l_{z}}\right)
|
||||
|
||||
where :math:`A` is the acceleration amplitude, :math:`l_z` is the z-length
|
||||
of the simulation box.
|
||||
At steady state, the acceleration generates a velocity profile:
|
||||
|
||||
.. math::
|
||||
|
||||
v_{x}(z) = V \cos \left(\frac{2 \pi z}{l_{z}}\right)
|
||||
|
||||
The generated velocity amplitude :math:`V` is related to the
|
||||
shear viscosity :math:`\eta` by:
|
||||
|
||||
.. math::
|
||||
|
||||
V = \frac{A \rho}{\eta}\left(\frac{l_{z}}{2 \pi}\right)^{2}
|
||||
|
||||
|
||||
and it can be obtained from ensemble average of the velocity profile:
|
||||
|
||||
.. math::
|
||||
|
||||
V = \frac{\sum_i 2 m_{i} v_{i, x} \cos \left(\frac{2 \pi z_i}{l_{z}}\right)}{\sum_i m_{i}}
|
||||
|
||||
where :math:`m_i`, :math:`v_{i,x}` and :math:`z_i` are the mass,
|
||||
x-component velocity and z coordinate of a particle.
|
||||
|
||||
The velocity amplitude :math:`V` can be calculated with :doc:`compute viscosity/cos <compute_viscosity_cos>`,
|
||||
which enables viscosity calculation with periodic perturbation method,
|
||||
as described by :ref:`Hess<Hess2>`.
|
||||
Because the applied acceleration drives the system away from equilibration,
|
||||
the calculated shear viscosity is lower than the intrinsic viscosity
|
||||
due to the shear-thinning effect.
|
||||
Extrapolation to zero acceleration should generally be performed to
|
||||
predict the zero-shear viscosity.
|
||||
As the shear stress decreases, the signal-noise ratio decreases rapidly,
|
||||
the simulation time must be extended accordingly to get converged result.
|
||||
|
||||
In order to get meaningful result, the group ID of this fix should be all.
|
||||
|
||||
----------
|
||||
|
||||
**Restart, fix_modify, output, run start/stop, minimize info:**
|
||||
|
||||
No information about this fix is written to binary restart files.
|
||||
None of the fix_modify options are relevant to this fix.
|
||||
No global or per-atom quantities are stored by this fix for access by various output commands.
|
||||
No parameter of this fix can be used with the start/stop keywords of the run command.
|
||||
This fix is not invoked during energy minimization.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This command is only available when LAMMPS was built with the USER-MISC package.
|
||||
Since this fix depends on the z-coordinate of atoms, it cannot be used in 2d simulations.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`compute viscosity/cos <compute_viscosity_cos>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
none
|
||||
|
||||
----------
|
||||
|
||||
.. _Hess2:
|
||||
|
||||
**(Hess)** Hess, B. The Journal of Chemical Physics 2002, 116 (1), 209-217.
|
||||
|
|
@ -1,5 +1,6 @@
|
|||
This directory has 5 scripts that compute the viscosity (eta) of a
|
||||
Lennard-Jones fluid using 5 different methods. See the discussion in
|
||||
This directory has 6 scripts that compute the viscosity (eta) of fluid
|
||||
using 6 different methods. 5 of them are for a Lennard-Jones fluid
|
||||
and the last one is for SPC/E water model. See the discussion in
|
||||
Section 6.21 of the manual for an overview of the methods and pointers
|
||||
to doc pages for the commands which implement them. Citations for the
|
||||
various methods can also be found in the manual.
|
||||
|
|
@ -10,7 +11,7 @@ enough to generate good statistics and highly accurate results.
|
|||
|
||||
-------------
|
||||
|
||||
These are the 5 methods for computing viscosity. The first 3 are
|
||||
These are the 5 methods for computing viscosity of a LJ fluid. The first 3 are
|
||||
non-equilibrium methods; the last 2 are equilibrium methods.
|
||||
|
||||
in.wall = move a wall to shear the fluid between two walls
|
||||
|
|
@ -89,3 +90,18 @@ heat/flux doc page - the resulting value prints at the end of the run
|
|||
and is in the log file
|
||||
|
||||
eta = 1.07
|
||||
|
||||
-------------
|
||||
|
||||
in.cos.1000SPCE is an example script of using cosine periodic perturbation method
|
||||
to calculate the viscosity of SPC/E water model.
|
||||
|
||||
The reciprocal of eta is computed within the script, and printed out as v_invVis
|
||||
in thermo_style command. The result will converge after hundreds of picoseconds.
|
||||
Then eta is obtained from the reciprocal of time average of v_invVis.
|
||||
|
||||
eta = 0.75 mPa*s
|
||||
|
||||
Note that the calculated viscosity by this method decreases with increased acceleration.
|
||||
It is therefore generally necessary to perform calculation at different accelerations
|
||||
and extrapolate the viscosity to zero shear.
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,53 @@
|
|||
# DFF generated Lammps input file
|
||||
|
||||
units real
|
||||
atom_style full
|
||||
boundary p p p
|
||||
|
||||
pair_style lj/cut/coul/long 10.0
|
||||
pair_modify mix arithmetic
|
||||
pair_modify tail yes
|
||||
kspace_style pppm 1.0e-4
|
||||
dielectric 1.0
|
||||
special_bonds amber
|
||||
bond_style harmonic
|
||||
angle_style harmonic
|
||||
dihedral_style none
|
||||
improper_style none
|
||||
|
||||
read_data data.cos.1000SPCE
|
||||
|
||||
variable T equal 300
|
||||
variable P equal 1.0
|
||||
|
||||
velocity all create ${T} 12345 mom yes rot yes dist gaussian
|
||||
|
||||
timestep 1.0
|
||||
|
||||
# Constraint ##################################
|
||||
fix com all momentum 100 linear 1 1 1
|
||||
fix rigid all shake 1e-4 20 0 b 1 a 1
|
||||
|
||||
# Viscosity ##################################
|
||||
variable A equal 0.05e-5 # angstrom/fs^2
|
||||
|
||||
fix cos all accelerate/cos ${A}
|
||||
compute cos all viscosity/cos
|
||||
|
||||
variable density equal density
|
||||
variable lz equal lz
|
||||
variable vMax equal c_cos[7] # velocity of atoms at z=0
|
||||
variable invVis equal v_vMax/${A}/v_density*39.4784/v_lz/v_lz*100 # reciprocal of viscosity 1/Pa/s
|
||||
|
||||
fix npt all npt temp ${T} ${T} 100 iso ${P} ${P} 1000
|
||||
fix_modify npt temp cos
|
||||
|
||||
thermo_style custom step cpu temp press pe density v_vMax v_invVis
|
||||
thermo_modify temp cos
|
||||
thermo 100
|
||||
################################################
|
||||
|
||||
dump 1 all custom 10000 dump.lammpstrj id mol type element q xu yu zu
|
||||
dump_modify 1 sort id element O H
|
||||
|
||||
run 2000
|
||||
|
|
@ -0,0 +1,160 @@
|
|||
LAMMPS (3 Mar 2020)
|
||||
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:94)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
# DFF generated Lammps input file
|
||||
|
||||
units real
|
||||
atom_style full
|
||||
boundary p p p
|
||||
|
||||
pair_style lj/cut/coul/long 10.0
|
||||
pair_modify mix arithmetic
|
||||
pair_modify tail yes
|
||||
kspace_style pppm 1.0e-4
|
||||
dielectric 1.0
|
||||
special_bonds amber
|
||||
bond_style harmonic
|
||||
angle_style harmonic
|
||||
dihedral_style none
|
||||
improper_style none
|
||||
|
||||
read_data data.1000SPCE.lmp
|
||||
orthogonal box = (0 0 0) to (31.043 31.043 31.043)
|
||||
2 by 2 by 2 MPI processor grid
|
||||
reading atoms ...
|
||||
3000 atoms
|
||||
scanning bonds ...
|
||||
2 = max bonds/atom
|
||||
scanning angles ...
|
||||
1 = max angles/atom
|
||||
reading bonds ...
|
||||
2000 bonds
|
||||
reading angles ...
|
||||
1000 angles
|
||||
2 = max # of 1-2 neighbors
|
||||
1 = max # of 1-3 neighbors
|
||||
1 = max # of 1-4 neighbors
|
||||
2 = max # of special neighbors
|
||||
special bonds CPU = 0.00114917 secs
|
||||
read_data CPU = 0.00953543 secs
|
||||
|
||||
variable T equal 300
|
||||
variable P equal 1.0
|
||||
|
||||
velocity all create ${T} 12345 mom yes rot yes dist gaussian
|
||||
velocity all create 300 12345 mom yes rot yes dist gaussian
|
||||
|
||||
timestep 1.0
|
||||
|
||||
# Constraint ##################################
|
||||
fix com all momentum 100 linear 1 1 1
|
||||
fix rigid all shake 1e-4 20 0 b 1 a 1
|
||||
0 = # of size 2 clusters
|
||||
0 = # of size 3 clusters
|
||||
0 = # of size 4 clusters
|
||||
1000 = # of frozen angles
|
||||
find clusters CPU = 0.000442737 secs
|
||||
|
||||
# Viscosity ##################################
|
||||
variable A equal 0.02e-5 # angstrom/fs^2
|
||||
|
||||
fix cos all accelerate/cos ${A}
|
||||
fix cos all accelerate/cos 2e-07
|
||||
compute cos all viscosity/cos
|
||||
|
||||
variable density equal density
|
||||
variable lz equal lz
|
||||
variable vMax equal c_cos[7] # velocity of atoms at z=0
|
||||
variable invVis equal v_vMax/${A}/v_density*39.4784/v_lz/v_lz*100 # reciprocal of viscosity 1/Pa/s
|
||||
variable invVis equal v_vMax/2e-07/v_density*39.4784/v_lz/v_lz*100
|
||||
|
||||
fix npt all npt temp ${T} ${T} 100 iso ${P} ${P} 1000
|
||||
fix npt all npt temp 300 ${T} 100 iso ${P} ${P} 1000
|
||||
fix npt all npt temp 300 300 100 iso ${P} ${P} 1000
|
||||
fix npt all npt temp 300 300 100 iso 1 ${P} 1000
|
||||
fix npt all npt temp 300 300 100 iso 1 1 1000
|
||||
fix_modify npt temp cos
|
||||
|
||||
thermo_style custom step cpu temp press pe density v_vMax v_invVis
|
||||
thermo_modify temp cos
|
||||
thermo 100
|
||||
################################################
|
||||
|
||||
dump 1 all custom 10000 dump.lammpstrj id mol type element q xu yu zu
|
||||
dump_modify 1 sort id element O H
|
||||
|
||||
run 2000
|
||||
PPPM initialization ...
|
||||
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
|
||||
G vector (1/distance) = 0.263539
|
||||
grid = 16 16 16
|
||||
stencil order = 5
|
||||
estimated absolute RMS force accuracy = 0.0325342
|
||||
estimated relative force accuracy = 9.79757e-05
|
||||
using double precision MKL FFT
|
||||
3d grid and FFT values/proc = 3375 512
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 12
|
||||
ghost atom cutoff = 12
|
||||
binsize = 6, bins = 6 6 6
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair lj/cut/coul/long, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 10.14 | 10.14 | 10.14 Mbytes
|
||||
Step CPU Temp Press PotEng Density v_vMax v_invVis
|
||||
0 0 450.04468 9838.6886 -7651.736 0.99999331 0.0001293705 2649.9663
|
||||
100 0.33736925 497.65155 2024.4827 -8342.2499 0.98595028 0.00019602427 4034.2452
|
||||
200 0.67116638 446.48518 27.075082 -8908.9684 0.9652009 0.00027615455 5723.7907
|
||||
300 0.99760895 401.79875 -776.99871 -9381.8986 0.95205822 0.00019626685 4086.6103
|
||||
400 1.3239019 369.65373 -510.5336 -9803.6463 0.94817309 0.00022998514 4795.2142
|
||||
500 1.6488092 343.35807 -936.31982 -10146.023 0.94840581 0.0001434423 2990.5423
|
||||
600 1.9826063 319.86131 -1381.3302 -10405.812 0.95459202 6.6411532e-05 1381.5767
|
||||
700 2.3172637 307.74606 -98.775733 -10643.5 0.96669652 0.00010446317 2164.0664
|
||||
800 2.6562841 305.14214 -540.57804 -10865.742 0.97808367 5.4381233e-05 1122.1765
|
||||
900 2.9938415 288.01316 639.00486 -10925.39 0.98686357 0.00010878474 2238.1355
|
||||
1000 3.327893 295.07773 -226.06503 -11033.826 0.99128496 0.00011935058 2451.8608
|
||||
1100 3.6618862 299.21578 306.34231 -11049.152 0.99552203 8.9538943e-05 1836.8166
|
||||
1200 3.9984287 301.82462 85.804646 -11013.564 0.99713434 0.00015912276 3262.51
|
||||
1300 4.3320735 308.6009 268.08897 -11009.836 0.99695358 0.00026212596 5374.72
|
||||
1400 4.668875 298.36903 -258.75495 -10962.299 0.99503447 0.00033087355 6788.7027
|
||||
1500 5.0003694 299.96073 99.512082 -10980.551 0.99315631 0.00033996557 6979.6425
|
||||
1600 5.3367337 304.18018 -500.65441 -11002.054 0.9914558 0.00039075642 8026.9849
|
||||
1700 5.6780828 301.63978 -499.07458 -10992.88 0.99234354 0.00038101175 7824.4738
|
||||
1800 6.0140638 303.25858 640.03432 -11053.335 0.99553958 0.00041336203 8479.7267
|
||||
1900 6.3532521 301.40882 208.28331 -11119.481 0.99534534 0.00032474734 6662.3144
|
||||
2000 6.6938104 298.0462 -236.47954 -11162.212 0.99421846 0.00023869721 4898.8129
|
||||
Loop time of 6.69387 on 8 procs for 2000 steps with 3000 atoms
|
||||
|
||||
Performance: 25.815 ns/day, 0.930 hours/ns, 298.781 timesteps/s
|
||||
99.7% CPU use with 8 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 4.36 | 4.4981 | 4.6781 | 5.1 | 67.20
|
||||
Bond | 0.00074545 | 0.00088463 | 0.0012464 | 0.0 | 0.01
|
||||
Kspace | 0.86696 | 1.0476 | 1.1863 | 10.5 | 15.65
|
||||
Neigh | 0.37733 | 0.37785 | 0.3784 | 0.1 | 5.64
|
||||
Comm | 0.19874 | 0.20168 | 0.20729 | 0.6 | 3.01
|
||||
Output | 0.0015529 | 0.0015803 | 0.0017546 | 0.2 | 0.02
|
||||
Modify | 0.54083 | 0.55143 | 0.55445 | 0.6 | 8.24
|
||||
Other | | 0.01483 | | | 0.22
|
||||
|
||||
Nlocal: 375 ave 385 max 361 min
|
||||
Histogram: 1 1 0 0 1 0 2 0 1 2
|
||||
Nghost: 5772.25 ave 5789 max 5757 min
|
||||
Histogram: 1 1 2 0 0 0 2 0 0 2
|
||||
Neighs: 135285 ave 144189 max 127550 min
|
||||
Histogram: 1 2 1 1 0 0 0 0 1 2
|
||||
|
||||
Total # of neighbors = 1082280
|
||||
Ave neighs/atom = 360.76
|
||||
Ave special neighs/atom = 2
|
||||
Neighbor list builds = 101
|
||||
Dangerous builds = 1
|
||||
Total wall time: 0:00:06
|
||||
|
|
@ -38,6 +38,7 @@ compute momentum, Rupert Nash (University of Edinburgh), r.nash at epcc.ed.ac.uk
|
|||
compute stress/mop, Romain Vermorel (U Pau) & Laurent Joly (U Lyon), romain.vermorel at univ-pau.fr & ljoly.ulyon at gmail.com, 5 Sep 18
|
||||
compute stress/mop/profile, Romain Vermorel (U Pau) & Laurent Joly (U Lyon), romain.vermorel at univ-pau.fr & ljoly.ulyon at gmail.com, 5 Sep 18
|
||||
compute temp/rotate, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
|
||||
compute viscosity/cos, Zheng Gong (ENS de Lyon), z.gong@outlook.com, 24 Apr 20
|
||||
compute PRESSURE/GREM, David Stelter, dstelter@bu.edu, 22 Nov 16
|
||||
dihedral_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
|
||||
dihedral_style fourier, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12
|
||||
|
|
@ -46,6 +47,7 @@ dihedral_style quadratic, Loukas Peristeras, loukas.peristeras at scienomics.com
|
|||
dihedral_style spherical, Andrew Jewett, jewett.aij@gmail.com, 15 Jul 16
|
||||
dihedral_style table, Andrew Jewett, jewett.aij@gmail.com, 10 Jan 12
|
||||
dihedral_style table/cut, Mike Salerno, ksalerno@pha.jhu.edu, 11 May 18
|
||||
fix accelerate/cos, Zheng Gong (ENS de Lyon), z.gong@outlook.com, 24 Apr 20
|
||||
fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
|
||||
fix ave/correlate/long, Jorge Ramirez (UPM Madrid), jorge.ramirez at upm.es, 21 Oct 2015
|
||||
fix electron/stopping, Konstantin Avchaciov, k.avchachov at gmail.com, 26 Feb 2019
|
||||
|
|
|
|||
|
|
@ -0,0 +1,293 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Zheng GONG (ENS de Lyon, z.gong@outlook.com)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <mpi.h>
|
||||
#include "compute_viscosity_cos.h"
|
||||
#include "atom.h"
|
||||
#include "update.h"
|
||||
#include "force.h"
|
||||
#include "domain.h"
|
||||
#include "comm.h"
|
||||
#include "group.h"
|
||||
#include "error.h"
|
||||
#include "math_const.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeViscosityCos::ComputeViscosityCos(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg) {
|
||||
if (narg != 3) error->all(FLERR, "Illegal compute viscosity/cos command");
|
||||
|
||||
scalar_flag = vector_flag = 1;
|
||||
size_vector = 7;
|
||||
extscalar = 0;
|
||||
extvector = -1;
|
||||
extlist = new int[7]{1,1,1,1,1,1,0};
|
||||
tempflag = 1;
|
||||
tempbias = 1;
|
||||
|
||||
maxbias = 0;
|
||||
vbiasall = NULL;
|
||||
|
||||
vector = new double[7];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeViscosityCos::~ComputeViscosityCos() {
|
||||
if (!copymode)
|
||||
delete[] vector;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeViscosityCos::setup() {
|
||||
dynamic = 0;
|
||||
if (dynamic_user || group->dynamic[igroup]) dynamic = 1;
|
||||
dof_compute();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeViscosityCos::dof_compute() {
|
||||
adjust_dof_fix();
|
||||
natoms_temp = group->count(igroup);
|
||||
dof = domain->dimension * natoms_temp;
|
||||
dof -= extra_dof + fix_dof;
|
||||
if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz);
|
||||
else tfactor = 0.0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void ComputeViscosityCos::calc_V() {
|
||||
double zlo = domain->boxlo[2];
|
||||
double zhi = domain->boxhi[2];
|
||||
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double *mass = atom->mass;
|
||||
double *rmass = atom->rmass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double massone;
|
||||
|
||||
double V_m[2];
|
||||
double V_m_local[2] = {0, 0};
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[type[i]];
|
||||
|
||||
V_m_local[0] +=
|
||||
2 * massone * v[i][0] * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
|
||||
V_m_local[1] += massone;
|
||||
}
|
||||
|
||||
MPI_Allreduce(V_m_local, V_m, 2, MPI_DOUBLE, MPI_SUM, world);
|
||||
V = V_m[0] / V_m[1];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double ComputeViscosityCos::compute_scalar() {
|
||||
invoked_scalar = update->ntimestep;
|
||||
|
||||
double zlo = domain->boxlo[2];
|
||||
double zhi = domain->boxhi[2];
|
||||
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double *mass = atom->mass;
|
||||
double *rmass = atom->rmass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double t = 0.0;
|
||||
double vx_acc;
|
||||
double massone;
|
||||
|
||||
calc_V();
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[type[i]];
|
||||
|
||||
vx_acc = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
|
||||
t += ((v[i][0] - vx_acc) * (v[i][0] - vx_acc) + v[i][1] * v[i][1] +
|
||||
v[i][2] * v[i][2]) * massone;
|
||||
}
|
||||
|
||||
MPI_Allreduce(&t, &scalar, 1, MPI_DOUBLE, MPI_SUM, world);
|
||||
if (dynamic) dof_compute();
|
||||
if (dof < 0.0 && natoms_temp > 0.0)
|
||||
error->all(FLERR, "Temperature compute degrees of freedom < 0");
|
||||
scalar *= tfactor;
|
||||
return scalar;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeViscosityCos::compute_vector() {
|
||||
int i;
|
||||
|
||||
invoked_vector = update->ntimestep;
|
||||
|
||||
double zlo = domain->boxlo[2];
|
||||
double zhi = domain->boxhi[2];
|
||||
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double *mass = atom->mass;
|
||||
double *rmass = atom->rmass;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double vx_acc;
|
||||
|
||||
double massone, t[6];
|
||||
for (i = 0; i < 6; i++) t[i] = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[type[i]];
|
||||
|
||||
vx_acc = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
|
||||
t[0] += massone * (v[i][0] - vx_acc) * (v[i][0] - vx_acc);
|
||||
t[1] += massone * v[i][1] * v[i][1];
|
||||
t[2] += massone * v[i][2] * v[i][2];
|
||||
t[3] += massone * (v[i][0] - vx_acc) * v[i][1];
|
||||
t[4] += massone * (v[i][0] - vx_acc) * v[i][2];
|
||||
t[5] += massone * v[i][1] * v[i][2];
|
||||
}
|
||||
|
||||
MPI_Allreduce(t, vector, 6, MPI_DOUBLE, MPI_SUM, world);
|
||||
for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
|
||||
vector[6] = V;
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
remove velocity bias from atom I to leave thermal velocity
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeViscosityCos::remove_bias(int i, double *v) {
|
||||
|
||||
double zlo = domain->boxlo[2];
|
||||
double zhi = domain->boxhi[2];
|
||||
|
||||
double **x = atom->x;
|
||||
|
||||
vbias[0] = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
|
||||
vbias[1] = 0;
|
||||
vbias[2] = 0;
|
||||
v[0] -= vbias[0];
|
||||
// v[1] -= vbias[1];
|
||||
// v[2] -= vbias[2];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
remove velocity bias from atom I to leave thermal velocity
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeViscosityCos::remove_bias_thr(int i, double *v, double *b) {
|
||||
double zlo = domain->boxlo[2];
|
||||
double zhi = domain->boxhi[2];
|
||||
|
||||
double **x = atom->x;
|
||||
|
||||
b[0] = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
|
||||
b[1] = 0;
|
||||
b[2] = 0;
|
||||
v[0] -= b[0];
|
||||
// v[1] -= b[1];
|
||||
// v[2] -= b[2];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
remove velocity bias from all atoms to leave thermal velocity
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeViscosityCos::remove_bias_all() {
|
||||
double zlo = domain->boxlo[2];
|
||||
double zhi = domain->boxhi[2];
|
||||
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
vbiasall[i][0] = V * cos(MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
|
||||
vbiasall[i][1] = 0;
|
||||
vbiasall[i][2] = 0;
|
||||
v[i][0] -= vbiasall[i][0];
|
||||
// v[i][1] -= vbiasall[i][1];
|
||||
// v[i][2] -= vbiasall[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
add back in velocity bias to atom I removed by remove_bias()
|
||||
assume remove_bias() was previously called
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeViscosityCos::restore_bias(int i, double *v) {
|
||||
v[0] += vbias[0];
|
||||
v[1] += vbias[1];
|
||||
v[2] += vbias[2];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
add back in velocity bias to atom I removed by remove_bias_thr()
|
||||
assume remove_bias_thr() was previously called with the same buffer b
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeViscosityCos::restore_bias_thr(int i, double *v, double *b) {
|
||||
v[0] += b[0];
|
||||
v[1] += b[1];
|
||||
v[2] += b[2];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
add back in velocity bias to all atoms removed by remove_bias_all()
|
||||
assume remove_bias_all() was previously called
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeViscosityCos::restore_bias_all() {
|
||||
double **v = atom->v;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
v[i][0] += vbiasall[i][0];
|
||||
v[i][1] += vbiasall[i][1];
|
||||
v[i][2] += vbiasall[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -0,0 +1,73 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Zheng GONG (ENS de Lyon, z.gong@outlook.com)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef COMPUTE_CLASS
|
||||
|
||||
ComputeStyle(viscosity/cos,ComputeViscosityCos)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_COMPUTE_VISCOSITY_COS_H
|
||||
#define LMP_COMPUTE_VISCOSITY_COS_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputeViscosityCos : public Compute {
|
||||
public:
|
||||
ComputeViscosityCos(class LAMMPS *, int, char **);
|
||||
virtual ~ComputeViscosityCos();
|
||||
void init() {}
|
||||
void setup();
|
||||
virtual double compute_scalar();
|
||||
virtual void compute_vector();
|
||||
|
||||
void remove_bias(int, double *);
|
||||
void remove_bias_thr(int, double *, double *);
|
||||
void remove_bias_all();
|
||||
void restore_bias(int, double *);
|
||||
void restore_bias_thr(int, double *, double *);
|
||||
void restore_bias_all();
|
||||
|
||||
protected:
|
||||
double tfactor;
|
||||
double V;
|
||||
|
||||
void dof_compute();
|
||||
void calc_V();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#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: Temperature compute degrees of freedom < 0
|
||||
|
||||
This should not happen if you are calculating the temperature
|
||||
on a valid set of atoms.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,85 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Zheng GONG (ENS de Lyon, z.gong@outlook.com)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "fix_accelerate_cos.h"
|
||||
#include "atom.h"
|
||||
#include "update.h"
|
||||
#include "modify.h"
|
||||
#include "domain.h"
|
||||
#include "region.h"
|
||||
#include "respa.h"
|
||||
#include "input.h"
|
||||
#include "variable.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixAccelerateCos::FixAccelerateCos(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg) {
|
||||
if (narg < 4) error->all(FLERR, "Illegal fix accelerate/cos command");
|
||||
acceleration = force->numeric(FLERR, arg[3]);
|
||||
if (domain->dimension == 2)
|
||||
error->all(FLERR,"Fix accelerate/cos cannot be used with 2d systems");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixAccelerateCos::setmask() {
|
||||
int mask = 0;
|
||||
mask |= POST_FORCE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixAccelerateCos::setup(int vflag) {
|
||||
post_force(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixAccelerateCos::post_force(int vflag) {
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *type = atom->type;
|
||||
double *mass = atom->mass;
|
||||
double *rmass = atom->rmass;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double massone, force_x, acc_x;
|
||||
double zlo = domain->boxlo[2];
|
||||
double zhi = domain->boxhi[2];
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[type[i]];
|
||||
|
||||
acc_x = acceleration *
|
||||
cos(MathConst::MY_2PI * (x[i][2] - zlo) / (zhi - zlo));
|
||||
force_x = acc_x * massone * force->mvv2e;
|
||||
|
||||
f[i][0] += force_x;
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,61 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Zheng GONG (ENS de Lyon, z.gong@outlook.com)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(accelerate/cos,FixAccelerateCos)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_ACCELERATE_COS_H
|
||||
#define LMP_FIX_ACCELERATE_COS_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixAccelerateCos: public Fix {
|
||||
public:
|
||||
FixAccelerateCos(class LAMMPS *, int, char **);
|
||||
virtual ~FixAccelerateCos() {};
|
||||
int setmask();
|
||||
virtual void init() {};
|
||||
void setup(int);
|
||||
virtual void post_force(int);
|
||||
|
||||
protected:
|
||||
double acceleration;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#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: Fix accelerate/cos cannot be used with 2d systems
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
*/
|
||||
Loading…
Reference in New Issue