tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3105 from athomps/numdiff-stress

Browse Source 
Add new numdiff/virial fix style
This commit is contained in:
Axel Kohlmeyer
2022-01-29 18:37:57 -05:00
committed by GitHub
parent 2a3a7c387b dad72a612a
commit a1d186b5fa
16 changed files with 1040 additions and 97 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_fix.rst
View File

@ -129,6 +129,7 @@ OPT.
* :doc:`npt/sphere (o) <fix_npt_sphere>`
* :doc:`npt/uef <fix_nh_uef>`
* :doc:`numdiff <fix_numdiff>`
* :doc:`numdiff/virial <fix_numdiff_virial>`
* :doc:`nve (giko) <fix_nve>`
* :doc:`nve/asphere (gi) <fix_nve_asphere>`
* :doc:`nve/asphere/noforce <fix_nve_asphere_noforce>`

7
doc/src/Errors_messages.rst
View File

@ -1941,6 +1941,9 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
*Compute ID for fix numdiff does not exist*
Self-explanatory.
*Compute ID for fix numdiff/virial does not exist*
Self-explanatory.
*Compute ID for fix store/state does not exist*
Self-explanatory.
@ -3796,6 +3799,10 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
Self-explanatory. Efficient loop over all atoms for numerical
difference requires consecutive atom IDs.
*Fix numdiff/virial must use group all*
Virial contributions computed by this fix are
computed on all atoms.
*Fix nve/asphere requires extended particles*
This fix can only be used for particles with a shape setting.

2
doc/src/compute_pressure.rst
View File

@ -141,7 +141,7 @@ Related commands
""""""""""""""""
:doc:`compute temp <compute_temp>`, :doc:`compute stress/atom <compute_stress_atom>`,
:doc:`thermo_style <thermo_style>`,
:doc:`thermo_style <thermo_style>`, :doc:`fix numdiff/virial <fix_numdiff_virial>`,
Default
"""""""

3
doc/src/fix.rst
View File

@ -271,7 +271,8 @@ accelerated styles exist.
* :doc:`npt/eff <fix_nh_eff>` - NPT for nuclei and electrons in the electron force field model
* :doc:`npt/sphere <fix_npt_sphere>` - NPT for spherical particles
* :doc:`npt/uef <fix_nh_uef>` - NPT style time integration with diagonal flow
* :doc:`numdiff <fix_numdiff>` - compute derivatives of per-atom data from finite differences
* :doc:`numdiff <fix_numdiff>` - numerically approximate atomic forces using finite energy differences
* :doc:`numdiff/virial <fix_numdiff_virial>` - numerically approximate virial stress tensor using finite energy differences
* :doc:`nve <fix_nve>` - constant NVE time integration
* :doc:`nve/asphere <fix_nve_asphere>` - NVE for aspherical particles
* :doc:`nve/asphere/noforce <fix_nve_asphere_noforce>` - NVE for aspherical particles without forces

23
doc/src/fix_numdiff.rst
View File

@ -13,16 +13,15 @@ Syntax
* ID, group-ID are documented in :doc:`fix <fix>` command
* numdiff = style name of this fix command
* Nevery = calculate force by finite difference every this many timesteps
* delta = finite difference displacement length (distance units)
* delta = size of atom displacements (distance units)
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all numdiff 1 0.0001
fix 1 all numdiff 10 1e-6
fix 1 all numdiff 100 0.01
fix 1 movegroup numdiff 100 0.01
Description
"""""""""""
@ -67,16 +66,17 @@ by two times *delta*.
The cost of each energy evaluation is essentially the cost of an MD
timestep. Thus invoking this fix once for a 3d system has a cost
of 6N timesteps, where N is the total number of atoms in the system
(assuming all atoms are included in the group). So this fix can be
very expensive to use for large systems.
of 6N timesteps, where N is the total number of atoms in the system.
So this fix can be very expensive to use for large systems.
One expedient alternative is to define the fix for a group containing
only a few atoms.
----------
The *Nevery* argument specifies on what timesteps the force will
be used calculated by finite difference.
The *delta* argument specifies the positional displacement each
The *delta* argument specifies the size of the displacement each
atom will undergo.
----------
@ -93,7 +93,12 @@ This fix produces a per-atom array which can be accessed by various
the force on each atom as calculated by finite difference. The
per-atom values can only be accessed on timesteps that are multiples
of *Nevery* since that is when the finite difference forces are
calculated.
calculated. See the examples in *examples/numdiff* directory
to see how this fix can be used to directly compare with
the analytic forces computed by LAMMPS.
The array values calculated by this compute
will be in force :doc:`units <units>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is invoked during :doc:`energy
@ -108,7 +113,7 @@ was built with that package. See the :doc:`Build package <Build_package>` page
Related commands
""""""""""""""""
:doc:`dynamical_matrix <dynamical_matrix>`,
:doc:`dynamical_matrix <dynamical_matrix>`, :doc:`fix numdiff/virial <fix_numdiff_virial>`,
Default
"""""""

115
doc/src/fix_numdiff_virial.rst Normal file
View File

@ -0,0 +1,115 @@
.. index:: fix numdiff/virial
fix numdiff/virial command
==========================
Syntax
""""""
.. parsed-literal::
fix ID group-ID numdiff/virial Nevery delta
* ID, group-ID are documented in :doc:`fix <fix>` command
* numdiff/virial = style name of this fix command
* Nevery = calculate virial by finite difference every this many timesteps
* delta = magnitude of strain fields (dimensionless)
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all numdiff/stress 10 1e-6
Description
"""""""""""
Calculate the virial stress tensor through a finite difference calculation of
energy versus strain. These values can be compared to the analytic virial
tensor computed by pair styles, bond styles, etc. This can be useful for
debugging or other purposes. The specified group must be "all".
This fix 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 energy of the entire system
is calculated after each. The difference in these two energies
divided by two times *delta*, approximates the corresponding
component of the virial stress tensor, 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 *Nevery* argument specifies on what timesteps the force will
be used calculated by finite difference.
The *delta* argument specifies the size of the displacement each
atom will undergo.
----------
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files
<restart>`. None of the :doc:`fix_modify <fix_modify>` options are
relevant to this fix.
This fix produces a global vector which can be accessed by various
:doc:`output commands <Howto_output>`, which stores the components of
the virial stress tensor as calculated by finite difference. The
global vector can only be accessed on timesteps that are multiples
of *Nevery* since that is when the finite difference virial is
calculated. See the examples in *examples/numdiff* directory
to see how this fix can be used to directly compare with
the analytic virial stress tensor computed by LAMMPS.
The order of the virial stress tensor components is *xx*, *yy*, *zz*,
*yz*, *xz*, and *xy*, consistent with Voigt notation. Note that
the vector produced by :doc:`compute pressure <compute_pressure>`
uses a different ordering, with *yz* and *xy* swapped.
The vector values calculated by this compute are
"intensive". The vector values will be in pressure
:doc:`units <units>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is invoked during :doc:`energy
minimization <minimize>`.
Restrictions
""""""""""""
This fix is part of the EXTRA-FIX package. It is only enabled if LAMMPS
was built with that package. See the :doc:`Build package <Build_package>` page for more info.
Related commands
""""""""""""""""
:doc:`fix numdiff <fix_numdiff>`, :doc:`compute pressure <compute_pressure>`
Default
"""""""
none

9
doc/src/fix_nvt_sllod.rst
View File

@ -66,7 +66,10 @@ equivalent to Newton's equations of motion for shear flow by
:ref:`(Evans and Morriss) <Evans3>`. They were later shown to generate
the desired velocity gradient and the correct production of work by
stresses for all forms of homogeneous flow by :ref:`(Daivis and Todd)
<Daivis>`. As implemented in LAMMPS, they are coupled to a
<Daivis>`.
The LAMMPS implementation corresponds to the p-SLLOD variant
of SLLOD. :ref:`(Edwards) <Edwards>`.
The equations of motion are coupled to a
Nose/Hoover chain thermostat in a velocity Verlet formulation, closely
following the implementation used for the :doc:`fix nvt <fix_nh>`
command.
@ -180,6 +183,10 @@ Same as :doc:`fix nvt <fix_nh>`, except tchain = 1.
**(Daivis and Todd)** Daivis and Todd, J Chem Phys, 124, 194103 (2006).
.. _Edwards:
**(Edwards)** Edwards, Baig, and Keffer, J Chem Phys 124, 194104 (2006).
.. _Daivis-sllod:
**(Daivis and Todd)** Daivis and Todd, Nonequilibrium Molecular Dynamics (book),

9
doc/src/fix_viscosity.rst
View File

@ -108,10 +108,11 @@ fluid, in appropriate units. See the :ref:`Muller-Plathe paper <Muller-Plathe2>
An alternative method for calculating a viscosity is to run a NEMD
simulation, as described on the :doc:`Howto nemd <Howto_nemd>` doc page.
NEMD simulations deform the simulation box via the :doc:`fix deform <fix_deform>` command. Thus they cannot be run on a charged
system using a :doc:`PPPM solver <kspace_style>` since PPPM does not
currently support non-orthogonal boxes. Using fix viscosity keeps the
box orthogonal; thus it does not suffer from this limitation.
NEMD simulations deform the simulation box via the :doc:`fix deform <fix_deform>` command.
Some features or combination of settings in LAMMPS do not support
non-orthogonal boxes. Using fix viscosity keeps the box orthogonal;
thus it does not suffer from these limitations.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

3
doc/utils/sphinx-config/false_positives.txt
View File

@ -194,6 +194,7 @@ Baczewski
Bagchi
Bagi
Bagnold
Baig
Bajaj
Bkappa
Bal
@ -668,6 +669,7 @@ Derlet
Deserno
Destree
destructor
destructors
detils
Devanathan
devel
@ -1569,6 +1571,7 @@ ke
KE
Keblinski
Keefe
Keffer
keflag
Keir
Kelchner

2
examples/README
View File

@ -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
numdiff: numerical difference computation of forces and virial
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

87
examples/numdiff/in.numdiff
View File

@ -1,33 +1,74 @@
# Numerical difference calculation of error in forces
# Numerical difference calculation
# of error in forces and virial stress
units metal
atom_style atomic
# adjustable parameters
atom_modify map yes
lattice fcc 5.358000
region box block 0 6 0 6 0 6
create_box 1 box
create_atoms 1 box
mass 1 39.903
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
velocity all create 10 2357 mom yes dist gaussian
units metal
atom_style atomic
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
atom_modify map yes
lattice fcc 5.358000
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
create_box 1 box
create_atoms 1 box
mass 1 39.903
neighbor 1 bin
velocity all create ${temp} 2357 mom yes dist gaussian
timestep 0.001
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
fix numdiff all numdiff 200 0.0001
fix nve all nve
neighbor 0.0 bin
neigh_modify every 1 delay 0 check no
variable errx atom fx-f_numdiff[1]
variable erry atom fy-f_numdiff[2]
variable errz atom fz-f_numdiff[3]
timestep 0.001
fix nve all nve
write_dump all custom tmp.error f_numdiff[1] f_numdiff[2] f_numdiff[3]
# define numerical force calculation
dump forces all custom 200 force_error.dump v_errx v_erry v_errz
thermo 200
run 2000
fix numforce all numdiff ${nthermo} ${fdelta}
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
# define numerical virial stress tensor calculation
compute myvirial all pressure NULL virial
fix numvirial all numdiff/virial ${nthermo} ${vdelta}
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}
run ${nsteps}

175
examples/numdiff/log.28Jan2022.log.numdiff.g++.1 Normal file
View File

@ -0,0 +1,175 @@
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

175
examples/numdiff/log.28Jan2022.log.numdiff.g++.4 Normal file
View File

@ -0,0 +1,175 @@
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

117
src/EXTRA-FIX/fix_numdiff.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -18,22 +17,23 @@
#include "fix_numdiff.h"
#include <cstring>
#include "atom.h"
#include "domain.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "respa.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "atom.h"
#include "bond.h"
#include "compute.h"
#include "dihedral.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "improper.h"
#include "kspace.h"
#include "memory.h"
#include "error.h"
#include "modify.h"
#include "pair.h"
#include "respa.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -41,20 +41,18 @@ using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixNumDiff::FixNumDiff(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), id_pe(nullptr), numdiff_forces(nullptr),
temp_x(nullptr), temp_f(nullptr)
Fix(lmp, narg, arg), id_pe(nullptr), numdiff_forces(nullptr), temp_x(nullptr), temp_f(nullptr)
{
if (narg < 5) error->all(FLERR,"Illegal fix numdiff command");
if (narg < 5) error->all(FLERR, "Illegal fix numdiff command");
peratom_flag = 1;
peratom_freq = nevery;
size_peratom_cols = 3;
respa_level_support = 1;
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
delta = utils::numeric(FLERR,arg[4],false,lmp);
if (nevery <= 0 || delta <= 0.0)
error->all(FLERR,"Illegal fix numdiff command");
nevery = utils::inumeric(FLERR, arg[3], false, lmp);
delta = utils::numeric(FLERR, arg[4], false, lmp);
if (nevery <= 0 || delta <= 0.0) error->all(FLERR, "Illegal fix numdiff command");
std::string cmd = id + std::string("_pe");
id_pe = utils::strdup(cmd);
@ -65,7 +63,7 @@ FixNumDiff::FixNumDiff(LAMMPS *lmp, int narg, char **arg) :
maxatom = 0;
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR,"Fix numdiff requires an atom map, see atom_modify");
error->all(FLERR, "Fix numdiff requires an atom map, see atom_modify");
// perform initial allocation of atom-based arrays
// zero numdiff_forces since dump may access it on timestep 0
@ -84,7 +82,7 @@ FixNumDiff::~FixNumDiff()
memory->destroy(temp_f);
modify->delete_compute(id_pe);
delete [] id_pe;
delete[] id_pe;
}
/* ---------------------------------------------------------------------- */
@ -107,22 +105,26 @@ void FixNumDiff::init()
// require consecutive atom IDs
if (!atom->tag_enable || !atom->tag_consecutive())
error->all(FLERR,"Fix numdiff requires consecutive atom IDs");
error->all(FLERR, "Fix numdiff requires consecutive atom IDs");
// check for PE compute
int icompute = modify->find_compute(id_pe);
if (icompute < 0) error->all(FLERR,"Compute ID for fix numdiff does not exist");
if (icompute < 0) error->all(FLERR, "Compute ID for fix numdiff does not exist");
pe = modify->compute[icompute];
if (force->pair && force->pair->compute_flag) pair_compute_flag = 1;
else pair_compute_flag = 0;
if (force->kspace && force->kspace->compute_flag) kspace_compute_flag = 1;
else kspace_compute_flag = 0;
if (force->pair && force->pair->compute_flag)
pair_compute_flag = 1;
else
pair_compute_flag = 0;
if (force->kspace && force->kspace->compute_flag)
kspace_compute_flag = 1;
else
kspace_compute_flag = 0;
if (utils::strmatch(update->integrate_style,"^respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
if (utils::strmatch(update->integrate_style, "^respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels - 1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level, ilevel_respa);
}
}
@ -130,11 +132,11 @@ void FixNumDiff::init()
void FixNumDiff::setup(int vflag)
{
if (utils::strmatch(update->integrate_style,"^verlet"))
if (utils::strmatch(update->integrate_style, "^verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
post_force_respa(vflag, ilevel_respa, 0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
@ -175,7 +177,7 @@ void FixNumDiff::min_post_force(int vflag)
void FixNumDiff::calculate_forces()
{
int i,j,ilocal;
int i, j, ilocal;
double energy;
// grow arrays if necessary
@ -202,43 +204,44 @@ void FixNumDiff::calculate_forces()
// loop over all atoms in system
// compute a finite difference force in each dimension
int flag,allflag;
int flag, allflag;
double denominator = 0.5 / delta;
int *mask = atom->mask;
int ntotal = static_cast<tagint> (atom->natoms);
int ntotal = static_cast<tagint>(atom->natoms);
int dimension = domain->dimension;
for (tagint m = 1; m <= ntotal; m++) {
ilocal = atom->map(m);
flag = 0;
if ((ilocal >= 0 && ilocal < nlocal) && (mask[ilocal] & groupbit)) flag = 1;
MPI_Allreduce(&flag,&allflag,1,MPI_INT,MPI_SUM,world);
MPI_Allreduce(&flag, &allflag, 1, MPI_INT, MPI_SUM, world);
if (!allflag) continue;
for (int idim = 0; idim < dimension; idim++) {
displace_atoms(ilocal,idim,1);
displace_atoms(ilocal, idim, 1);
energy = update_energy();
if (ilocal >= 0 && ilocal < nlocal)
numdiff_forces[ilocal][idim] -= energy;
if (ilocal >= 0 && ilocal < nlocal) numdiff_forces[ilocal][idim] -= energy;
displace_atoms(ilocal,idim,-2);
displace_atoms(ilocal, idim, -2);
energy = update_energy();
if (ilocal >= 0 && ilocal < nlocal) {
numdiff_forces[ilocal][idim] += energy;
numdiff_forces[ilocal][idim] *= denominator;
}
restore_atoms(ilocal,idim);
restore_atoms(ilocal, idim);
}
}
// recompute energy so all contributions are as before
energy = update_energy();
// restore original forces for owned and ghost atoms
for (i = 0; i < nall; i++)
for (j = 0; j < 3; j++) {
f[i][j] = temp_f[i][j];
}
for (j = 0; j < 3; j++) { f[i][j] = temp_f[i][j]; }
}
/* ----------------------------------------------------------------------
@ -252,11 +255,11 @@ void FixNumDiff::displace_atoms(int ilocal, int idim, int magnitude)
double **x = atom->x;
int *sametag = atom->sametag;
int j = ilocal;
x[ilocal][idim] += delta*magnitude;
x[ilocal][idim] += delta * magnitude;
while (sametag[j] >= 0) {
j = sametag[j];
x[j][idim] += delta*magnitude;
x[j][idim] += delta * magnitude;
}
}
@ -290,16 +293,16 @@ double FixNumDiff::update_energy()
int eflag = 1;
if (pair_compute_flag) force->pair->compute(eflag,0);
if (pair_compute_flag) force->pair->compute(eflag, 0);
if (atom->molecular != Atom::ATOMIC) {
if (force->bond) force->bond->compute(eflag,0);
if (force->angle) force->angle->compute(eflag,0);
if (force->dihedral) force->dihedral->compute(eflag,0);
if (force->improper) force->improper->compute(eflag,0);
if (force->bond) force->bond->compute(eflag, 0);
if (force->angle) force->angle->compute(eflag, 0);
if (force->dihedral) force->dihedral->compute(eflag, 0);
if (force->improper) force->improper->compute(eflag, 0);
}
if (kspace_compute_flag) force->kspace->compute(eflag,0);
if (kspace_compute_flag) force->kspace->compute(eflag, 0);
double energy = pe->compute_scalar();
return energy;
@ -313,7 +316,7 @@ void FixNumDiff::force_clear(double **forces)
{
size_t nbytes = sizeof(double) * atom->nlocal;
if (force->newton) nbytes += sizeof(double) * atom->nghost;
if (nbytes) memset(&forces[0][0],0,3*nbytes);
if (nbytes) memset(&forces[0][0], 0, 3 * nbytes);
}
/* ----------------------------------------------------------------------
@ -326,9 +329,9 @@ void FixNumDiff::reallocate()
memory->destroy(temp_x);
memory->destroy(temp_f);
maxatom = atom->nmax;
memory->create(numdiff_forces,maxatom,3,"numdiff:numdiff_force");
memory->create(temp_x,maxatom,3,"numdiff:temp_x");
memory->create(temp_f,maxatom,3,"numdiff:temp_f");
memory->create(numdiff_forces, maxatom, 3, "numdiff:numdiff_force");
memory->create(temp_x, maxatom, 3, "numdiff:temp_x");
memory->create(temp_f, maxatom, 3, "numdiff:temp_f");
array_atom = numdiff_forces;
}
@ -339,6 +342,6 @@ void FixNumDiff::reallocate()
double FixNumDiff::memory_usage()
{
double bytes = 0.0;
bytes += (double)3 * maxatom*3 * sizeof(double);
bytes += (double) 3 * maxatom * 3 * sizeof(double);
return bytes;
}

320
src/EXTRA-FIX/fix_numdiff_virial.cpp Normal file
View File

@ -0,0 +1,320 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://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 author: Aidan Thompson (Sandia)
------------------------------------------------------------------------- */
#include "fix_numdiff_virial.h"
#include "angle.h"
#include "atom.h"
#include "bond.h"
#include "compute.h"
#include "dihedral.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "improper.h"
#include "kspace.h"
#include "memory.h"
#include "modify.h"
#include "pair.h"
#include "respa.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixNumDiffVirial::FixNumDiffVirial(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), id_pe(nullptr), temp_x(nullptr), temp_f(nullptr)
{
if (narg < 5) error->all(FLERR, "Illegal fix numdiff/virial command");
if (igroup) error->all(FLERR, "Compute numdiff/virial must use group all");
peratom_freq = nevery;
respa_level_support = 1;
vector_flag = 1;
size_vector = NDIR_VIRIAL;
extvector = 0;
maxatom = 0;
nevery = utils::inumeric(FLERR, arg[3], false, lmp);
delta = utils::numeric(FLERR, arg[4], false, lmp);
if (nevery <= 0 || delta <= 0.0) error->all(FLERR, "Illegal fix numdiff command");
std::string cmd = id + std::string("_pe");
id_pe = utils::strdup(cmd);
cmd += " all pe";
modify->add_compute(cmd);
// perform initial allocation of atom-based arrays
// zero numdiff_forces since dump may access it on timestep 0
// zero numdiff_forces since a variable may access it before first run
reallocate();
// set fixed-point to default = center of cell
fixedpoint[0] = 0.5 * (domain->boxlo[0] + domain->boxhi[0]);
fixedpoint[1] = 0.5 * (domain->boxlo[1] + domain->boxhi[1]);
fixedpoint[2] = 0.5 * (domain->boxlo[2] + domain->boxhi[2]);
// define the cartesian indices for each strain (Voigt order)
dirlist[0][0] = 0;
dirlist[0][1] = 0;
dirlist[1][0] = 1;
dirlist[1][1] = 1;
dirlist[2][0] = 2;
dirlist[2][1] = 2;
dirlist[3][0] = 1;
dirlist[3][1] = 2;
dirlist[4][0] = 0;
dirlist[4][1] = 2;
dirlist[5][0] = 0;
dirlist[5][1] = 1;
}
/* ---------------------------------------------------------------------- */
FixNumDiffVirial::~FixNumDiffVirial()
{
memory->destroy(temp_x);
memory->destroy(temp_f);
modify->delete_compute(id_pe);
delete[] id_pe;
}
/* ---------------------------------------------------------------------- */
int FixNumDiffVirial::setmask()
{
datamask_read = datamask_modify = 0;
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixNumDiffVirial::init()
{
// check for PE compute
int icompute = modify->find_compute(id_pe);
if (icompute < 0) error->all(FLERR, "Compute ID for fix numdiff does not exist");
pe = modify->compute[icompute];
if (force->pair && force->pair->compute_flag)
pair_compute_flag = 1;
else
pair_compute_flag = 0;
if (force->kspace && force->kspace->compute_flag)
kspace_compute_flag = 1;
else
kspace_compute_flag = 0;
if (utils::strmatch(update->integrate_style, "^respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels - 1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level, ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
void FixNumDiffVirial::setup(int vflag)
{
if (utils::strmatch(update->integrate_style, "^verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag, ilevel_respa, 0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
void FixNumDiffVirial::min_setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixNumDiffVirial::post_force(int /* vflag */)
{
if (update->ntimestep % nevery) return;
calculate_virial();
}
/* ---------------------------------------------------------------------- */
void FixNumDiffVirial::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixNumDiffVirial::min_post_force(int vflag)
{
post_force(vflag);
}
/* ----------------------------------------------------------------------
compute finite difference virial stress tensor
------------------------------------------------------------------------- */
void FixNumDiffVirial::calculate_virial()
{
double energy;
// grow arrays if necessary
if (atom->nlocal + atom->nghost > maxatom) reallocate();
// store copy of current forces for owned and ghost atoms
double **x = atom->x;
double **f = atom->f;
int nall = atom->nlocal + atom->nghost;
for (int i = 0; i < nall; i++)
for (int k = 0; k < 3; k++) {
temp_x[i][k] = x[i][k];
temp_f[i][k] = f[i][k];
}
// loop over 6 strain directions
// compute a finite difference force in each dimension
int flag, allflag;
double nktv2p = force->nktv2p;
double inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
double denominator = -0.5 / delta * inv_volume * nktv2p;
for (int idir = 0; idir < NDIR_VIRIAL; idir++) {
displace_atoms(nall, idir, 1.0);
energy = update_energy();
virial[idir] = energy;
restore_atoms(nall, idir);
displace_atoms(nall, idir, -1.0);
energy = update_energy();
virial[idir] -= energy;
virial[idir] *= denominator;
restore_atoms(nall, idir);
}
// recompute energy so all contributions are as before
energy = update_energy();
// restore original forces for owned and ghost atoms
for (int i = 0; i < nall; i++)
for (int k = 0; k < 3; k++) f[i][k] = temp_f[i][k];
}
/* ----------------------------------------------------------------------
displace position of all owned and ghost atoms
---------------------------------------------------------------------- */
void FixNumDiffVirial::displace_atoms(int nall, int idir, double magnitude)
{
double **x = atom->x;
int k = dirlist[idir][0];
int l = dirlist[idir][1];
for (int i = 0; i < nall; i++)
x[i][k] = temp_x[i][k] + delta * magnitude * (temp_x[i][l] - fixedpoint[l]);
}
/* ----------------------------------------------------------------------
restore position of all owned and ghost atoms
---------------------------------------------------------------------- */
void FixNumDiffVirial::restore_atoms(int nall, int idir)
{
double **x = atom->x;
int k = dirlist[idir][0];
for (int i = 0; i < nall; i++) { x[i][k] = temp_x[i][k]; }
}
/* ----------------------------------------------------------------------
evaluate potential energy and forces
same logic as in Verlet
------------------------------------------------------------------------- */
double FixNumDiffVirial::update_energy()
{
int eflag = 1;
if (pair_compute_flag) force->pair->compute(eflag, 0);
if (atom->molecular != Atom::ATOMIC) {
if (force->bond) force->bond->compute(eflag, 0);
if (force->angle) force->angle->compute(eflag, 0);
if (force->dihedral) force->dihedral->compute(eflag, 0);
if (force->improper) force->improper->compute(eflag, 0);
}
if (kspace_compute_flag) force->kspace->compute(eflag, 0);
double energy = pe->compute_scalar();
return energy;
}
/* ----------------------------------------------------------------------
return Ith vector value, assume in range of size_vector
------------------------------------------------------------------------- */
double FixNumDiffVirial::compute_vector(int i)
{
return virial[i];
}
/* ----------------------------------------------------------------------
reallocated local per-atoms arrays
------------------------------------------------------------------------- */
void FixNumDiffVirial::reallocate()
{
memory->destroy(temp_x);
memory->destroy(temp_f);
maxatom = atom->nmax;
memory->create(temp_x, maxatom, 3, "numdiff/virial:temp_x");
memory->create(temp_f, maxatom, 3, "numdiff/virial:temp_f");
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixNumDiffVirial::memory_usage()
{
double bytes = 0.0;
bytes += (double) 2 * maxatom * 3 * sizeof(double);
return bytes;
}

89
src/EXTRA-FIX/fix_numdiff_virial.h Normal file
View File

@ -0,0 +1,89 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://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.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(numdiff/virial,FixNumDiffVirial);
// clang-format on
#else
#ifndef LMP_FIX_NUMDIFF_VIRIAL_H
#define LMP_FIX_NUMDIFF_VIRIAL_H
#include "fix.h"
namespace LAMMPS_NS {
class FixNumDiffVirial : public Fix {
public:
FixNumDiffVirial(class LAMMPS *, int, char **);
~FixNumDiffVirial() override;
int setmask() override;
void init() override;
void setup(int) override;
void min_setup(int) override;
void post_force(int) override;
void post_force_respa(int, int, int) override;
void min_post_force(int) override;
double memory_usage() override;
private:
static constexpr int NDIR_VIRIAL = 6; // dimension of virial and strain vectors
double delta; // strain magnitude
int maxatom; // allocated size of atom arrays
int ilevel_respa;
int pair_compute_flag; // 0 if pair->compute is skipped
int kspace_compute_flag; // 0 if kspace->compute is skipped
char *id_pe; // name of energy compute
class Compute *pe; // pointer to energy compute
double virial[NDIR_VIRIAL]; // finite diff virial components (Voigt order)
double **temp_x; // original coords
double **temp_f; // original forces
double fixedpoint[3]; // define displacement field origin
int dirlist[NDIR_VIRIAL][2]; // strain cartesian indices (Voigt order)
double compute_vector(int) override; // access function for virial
void calculate_virial(); // virial calculation
void displace_atoms(int, int, double); // apply displacement field
void restore_atoms(int, int); // restore original positions
double update_energy(); // calculate new energy
void virial_clear(); // set virial to zero
void reallocate(); // grow the atom arrays
};
} // 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: Compute ID for fix numdiff/virial does not exist
Self-explanatory.
E: Fix numdiff/virial must use group all
Virial contributions computed by this fix are
computed on all atoms.
*/