From a9481733a016fbb3f4fcdf2b6d7f562b000faa5e Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 28 Jan 2022 18:34:30 -0700 Subject: [PATCH] Finished debugging, testing, documenting --- doc/src/Commands_fix.rst | 1 + doc/src/Errors_messages.rst | 7 + doc/src/compute_pressure.rst | 2 +- doc/src/fix.rst | 3 +- doc/src/fix_numdiff.rst | 23 ++- doc/src/fix_numdiff_virial.rst | 115 ++++++++++++ doc/src/fix_nvt_sllod.rst | 9 +- doc/src/fix_viscosity.rst | 4 +- doc/utils/sphinx-config/false_positives.txt | 2 + examples/README | 2 +- examples/numdiff/in.numdiff | 90 ++++++--- .../numdiff/log.28Jan2022.log.numdiff.g++.1 | 175 ++++++++++++++++++ .../numdiff/log.28Jan2022.log.numdiff.g++.4 | 175 ++++++++++++++++++ src/EXTRA-FIX/fix_numdiff.cpp | 5 + ...diff_stress.cpp => fix_numdiff_virial.cpp} | 79 ++++---- ..._numdiff_stress.h => fix_numdiff_virial.h} | 31 ++-- 16 files changed, 627 insertions(+), 96 deletions(-) create mode 100644 doc/src/fix_numdiff_virial.rst create mode 100644 examples/numdiff/log.28Jan2022.log.numdiff.g++.1 create mode 100644 examples/numdiff/log.28Jan2022.log.numdiff.g++.4 rename src/EXTRA-FIX/{fix_numdiff_stress.cpp => fix_numdiff_virial.cpp} (82%) rename src/EXTRA-FIX/{fix_numdiff_stress.h => fix_numdiff_virial.h} (73%) diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index c3fca16bef..35e4f671ef 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -129,6 +129,7 @@ OPT. * :doc:`npt/sphere (o) ` * :doc:`npt/uef ` * :doc:`numdiff ` + * :doc:`numdiff/virial ` * :doc:`nve (giko) ` * :doc:`nve/asphere (gi) ` * :doc:`nve/asphere/noforce ` diff --git a/doc/src/Errors_messages.rst b/doc/src/Errors_messages.rst index c06f4c86e3..4e216828d3 100644 --- a/doc/src/Errors_messages.rst +++ b/doc/src/Errors_messages.rst @@ -1941,6 +1941,9 @@ Doc page with :doc:`WARNING messages ` *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 ` 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. diff --git a/doc/src/compute_pressure.rst b/doc/src/compute_pressure.rst index 6d76713b41..ab57786c05 100644 --- a/doc/src/compute_pressure.rst +++ b/doc/src/compute_pressure.rst @@ -141,7 +141,7 @@ Related commands """""""""""""""" :doc:`compute temp `, :doc:`compute stress/atom `, -:doc:`thermo_style `, +:doc:`thermo_style `, :doc:`fix numdiff/virial `, Default """"""" diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 6afe0fe325..92a00abe46 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -271,7 +271,8 @@ accelerated styles exist. * :doc:`npt/eff ` - NPT for nuclei and electrons in the electron force field model * :doc:`npt/sphere ` - NPT for spherical particles * :doc:`npt/uef ` - NPT style time integration with diagonal flow -* :doc:`numdiff ` - compute derivatives of per-atom data from finite differences +* :doc:`numdiff ` - numerically approximate atomic forces using finite energy differences +* :doc:`numdiff/virial ` - numerically approximate virial stress tensor using finite energy differences * :doc:`nve ` - constant NVE time integration * :doc:`nve/asphere ` - NVE for aspherical particles * :doc:`nve/asphere/noforce ` - NVE for aspherical particles without forces diff --git a/doc/src/fix_numdiff.rst b/doc/src/fix_numdiff.rst index b0686f471f..b5fdfb12e6 100644 --- a/doc/src/fix_numdiff.rst +++ b/doc/src/fix_numdiff.rst @@ -13,16 +13,15 @@ Syntax * ID, group-ID are documented in :doc:`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 `. No parameter of this fix can be used with the *start/stop* keywords of the :doc:`run ` command. This fix is invoked during :doc:`energy @@ -108,7 +113,7 @@ was built with that package. See the :doc:`Build package ` page Related commands """""""""""""""" -:doc:`dynamical_matrix `, +:doc:`dynamical_matrix `, :doc:`fix numdiff/virial `, Default """"""" diff --git a/doc/src/fix_numdiff_virial.rst b/doc/src/fix_numdiff_virial.rst new file mode 100644 index 0000000000..647a0e592a --- /dev/null +++ b/doc/src/fix_numdiff_virial.rst @@ -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 ` 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 +`. None of the :doc:`fix_modify ` options are +relevant to this fix. + +This fix produces a global vector which can be accessed by various +:doc:`output commands `, 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 ` +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 `. + +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. This fix is invoked during :doc:`energy +minimization `. + +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 ` page for more info. + +Related commands +"""""""""""""""" + +:doc:`fix numdiff `, :doc:`compute pressure ` + +Default +""""""" + +none diff --git a/doc/src/fix_nvt_sllod.rst b/doc/src/fix_nvt_sllod.rst index 4bb4478991..04e60057a1 100644 --- a/doc/src/fix_nvt_sllod.rst +++ b/doc/src/fix_nvt_sllod.rst @@ -66,7 +66,10 @@ equivalent to Newton's equations of motion for shear flow by :ref:`(Evans and Morriss) `. 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) -`. As implemented in LAMMPS, they are coupled to a +`. +The LAMMPS implementation corresponds to the p-SLLOD variant +of SLLOD. :ref:`(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 ` command. @@ -180,6 +183,10 @@ Same as :doc:`fix nvt `, 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), diff --git a/doc/src/fix_viscosity.rst b/doc/src/fix_viscosity.rst index e077af19f4..e92d1eab92 100644 --- a/doc/src/fix_viscosity.rst +++ b/doc/src/fix_viscosity.rst @@ -108,7 +108,9 @@ fluid, in appropriate units. See the :ref:`Muller-Plathe paper An alternative method for calculating a viscosity is to run a NEMD simulation, as described on the :doc:`Howto nemd ` doc page. -NEMD simulations deform the simulation box via the :doc:`fix deform ` command. Thus they cannot be run on a charged +NEMD simulations deform the simulation box via the :doc:`fix deform ` command. + +Thus they cannot be run on a charged system using a :doc:`PPPM solver ` since PPPM does not currently support non-orthogonal boxes. Using fix viscosity keeps the box orthogonal; thus it does not suffer from this limitation. diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index fe1e40e8ba..53e8781a9e 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -194,6 +194,7 @@ Baczewski Bagchi Bagi Bagnold +Baig Bajaj Bkappa Bal @@ -1569,6 +1570,7 @@ ke KE Keblinski Keefe +Keffer keflag Keir Kelchner diff --git a/examples/README b/examples/README index c398579cde..0c09b6d847 100644 --- a/examples/README +++ b/examples/README @@ -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 diff --git a/examples/numdiff/in.numdiff b/examples/numdiff/in.numdiff index 07d39f2ddf..cd3d44a9e7 100644 --- a/examples/numdiff/in.numdiff +++ b/examples/numdiff/in.numdiff @@ -1,38 +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 3 0 3 0 3 -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.0 bin +velocity all create ${temp} 2357 mom yes dist gaussian -timestep 0.001 +pair_style lj/cubic +pair_coeff * * 0.0102701 3.42 -fix numforce all numdiff 100 0.0001 -fix numstress all numdiff/stress 100 1.0e-2 -fix nve all nve +neighbor 0.0 bin +neigh_modify every 1 delay 0 check no -variable errx atom fx-f_numforce[1] -variable erry atom fy-f_numforce[2] -variable errz atom fz-f_numforce[3] +timestep 0.001 +fix nve all nve -dump errors all custom 100 force_error.dump v_errx v_erry v_errz +# define numerical force calculation + +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 -variable ferrsq atom (fx-f_numforce[1])^2+(fy-f_numforce[2])^2+(fz-f_numforce[3])^2 -compute faverrsq all reduce ave v_ferrsq compute myvirial all pressure NULL virial -variable ratio11 equal f_numstress[1]/c_myvirial[1] -thermo_style custom step temp pe press c_faverrsq c_myvirial[*] f_numstress[*] v_ratio11 -thermo 100 -run 500 +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} diff --git a/examples/numdiff/log.28Jan2022.log.numdiff.g++.1 b/examples/numdiff/log.28Jan2022.log.numdiff.g++.1 new file mode 100644 index 0000000000..f32e72834f --- /dev/null +++ b/examples/numdiff/log.28Jan2022.log.numdiff.g++.1 @@ -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 diff --git a/examples/numdiff/log.28Jan2022.log.numdiff.g++.4 b/examples/numdiff/log.28Jan2022.log.numdiff.g++.4 new file mode 100644 index 0000000000..fc56c4aee8 --- /dev/null +++ b/examples/numdiff/log.28Jan2022.log.numdiff.g++.4 @@ -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 diff --git a/src/EXTRA-FIX/fix_numdiff.cpp b/src/EXTRA-FIX/fix_numdiff.cpp index 4c9c08b681..757d6adab4 100644 --- a/src/EXTRA-FIX/fix_numdiff.cpp +++ b/src/EXTRA-FIX/fix_numdiff.cpp @@ -233,12 +233,17 @@ void FixNumDiff::calculate_forces() } } + // 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]; } + } /* ---------------------------------------------------------------------- diff --git a/src/EXTRA-FIX/fix_numdiff_stress.cpp b/src/EXTRA-FIX/fix_numdiff_virial.cpp similarity index 82% rename from src/EXTRA-FIX/fix_numdiff_stress.cpp rename to src/EXTRA-FIX/fix_numdiff_virial.cpp index cc1adc7e2a..4d1cf20028 100644 --- a/src/EXTRA-FIX/fix_numdiff_stress.cpp +++ b/src/EXTRA-FIX/fix_numdiff_virial.cpp @@ -16,7 +16,7 @@ Contributing author: Aidan Thompson (Sandia) ------------------------------------------------------------------------- */ -#include "fix_numdiff_stress.h" +#include "fix_numdiff_virial.h" #include #include "atom.h" @@ -40,16 +40,17 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ -FixNumDiffStress::FixNumDiffStress(LAMMPS *lmp, int narg, char **arg) : +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/stress command"); + 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_STRESS; + size_vector = NDIR_VIRIAL; extvector = 0; maxatom = 0; @@ -96,7 +97,7 @@ FixNumDiffStress::FixNumDiffStress(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -FixNumDiffStress::~FixNumDiffStress() +FixNumDiffVirial::~FixNumDiffVirial() { memory->destroy(temp_x); memory->destroy(temp_f); @@ -107,7 +108,7 @@ FixNumDiffStress::~FixNumDiffStress() /* ---------------------------------------------------------------------- */ -int FixNumDiffStress::setmask() +int FixNumDiffVirial::setmask() { datamask_read = datamask_modify = 0; @@ -120,7 +121,7 @@ int FixNumDiffStress::setmask() /* ---------------------------------------------------------------------- */ -void FixNumDiffStress::init() +void FixNumDiffVirial::init() { // check for PE compute @@ -141,7 +142,7 @@ void FixNumDiffStress::init() /* ---------------------------------------------------------------------- */ -void FixNumDiffStress::setup(int vflag) +void FixNumDiffVirial::setup(int vflag) { if (utils::strmatch(update->integrate_style,"^verlet")) post_force(vflag); @@ -154,39 +155,39 @@ void FixNumDiffStress::setup(int vflag) /* ---------------------------------------------------------------------- */ -void FixNumDiffStress::min_setup(int vflag) +void FixNumDiffVirial::min_setup(int vflag) { post_force(vflag); } /* ---------------------------------------------------------------------- */ -void FixNumDiffStress::post_force(int /* vflag */) +void FixNumDiffVirial::post_force(int /* vflag */) { if (update->ntimestep % nevery) return; - calculate_stress(); + calculate_virial(); } /* ---------------------------------------------------------------------- */ -void FixNumDiffStress::post_force_respa(int vflag, int ilevel, int /*iloop*/) +void FixNumDiffVirial::post_force_respa(int vflag, int ilevel, int /*iloop*/) { if (ilevel == ilevel_respa) post_force(vflag); } /* ---------------------------------------------------------------------- */ -void FixNumDiffStress::min_post_force(int vflag) +void FixNumDiffVirial::min_post_force(int vflag) { post_force(vflag); } /* ---------------------------------------------------------------------- - compute finite difference stress tensor + compute finite difference virial stress tensor ------------------------------------------------------------------------- */ -void FixNumDiffStress::calculate_stress() +void FixNumDiffVirial::calculate_virial() { double energy; @@ -198,8 +199,7 @@ void FixNumDiffStress::calculate_stress() double **x = atom->x; double **f = atom->f; - int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; + int nall = atom->nlocal + atom->nghost; for (int i = 0; i < nall; i++) for (int k = 0; k < 3; k++) { @@ -216,17 +216,22 @@ void FixNumDiffStress::calculate_stress() double denominator = -0.5 / delta * inv_volume * nktv2p; - for (int idir = 0; idir < NDIR_STRESS; idir++) { + for (int idir = 0; idir < NDIR_VIRIAL; idir++) { displace_atoms(nall, idir, 1.0); energy = update_energy(); - stress[idir] = energy; - displace_atoms(nall, idir,-1.0); + virial[idir] = energy; + restore_atoms(nall, idir); + displace_atoms(nall, idir, -1.0); energy = update_energy(); - stress[idir] -= energy; - stress[idir] *= denominator; + 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++) @@ -239,13 +244,13 @@ void FixNumDiffStress::calculate_stress() displace position of all owned and ghost atoms ---------------------------------------------------------------------- */ -void FixNumDiffStress::displace_atoms(int nall, int idir, double magnitude) +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* + x[i][k] = temp_x[i][k] + delta*magnitude* (temp_x[i][l]-fixedpoint[l]); } @@ -253,7 +258,7 @@ void FixNumDiffStress::displace_atoms(int nall, int idir, double magnitude) restore position of all owned and ghost atoms ---------------------------------------------------------------------- */ -void FixNumDiffStress::restore_atoms(int nall, int idir) +void FixNumDiffVirial::restore_atoms(int nall, int idir) { double **x = atom->x; int k = dirlist[idir][0]; @@ -267,7 +272,7 @@ void FixNumDiffStress::restore_atoms(int nall, int idir) same logic as in Verlet ------------------------------------------------------------------------- */ -double FixNumDiffStress::update_energy() +double FixNumDiffVirial::update_energy() { int eflag = 1; @@ -286,43 +291,33 @@ double FixNumDiffStress::update_energy() return energy; } -/* ---------------------------------------------------------------------- - clear forces needed -------------------------------------------------------------------------- */ - -void FixNumDiffStress::stress_clear() -{ - size_t nbytes = sizeof(double) * size_vector; - memset(&stress[0],0,nbytes); -} - /* ---------------------------------------------------------------------- return Ith vector value, assume in range of size_vector ------------------------------------------------------------------------- */ -double FixNumDiffStress::compute_vector(int i) +double FixNumDiffVirial::compute_vector(int i) { - return stress[i]; + return virial[i]; } /* ---------------------------------------------------------------------- reallocated local per-atoms arrays ------------------------------------------------------------------------- */ -void FixNumDiffStress::reallocate() +void FixNumDiffVirial::reallocate() { memory->destroy(temp_x); memory->destroy(temp_f); maxatom = atom->nmax; - memory->create(temp_x,maxatom,3,"numdiff/stress:temp_x"); - memory->create(temp_f,maxatom,3,"numdiff/stress:temp_f"); + 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 FixNumDiffStress::memory_usage() +double FixNumDiffVirial::memory_usage() { double bytes = 0.0; bytes += (double)2 * maxatom*3 * sizeof(double); diff --git a/src/EXTRA-FIX/fix_numdiff_stress.h b/src/EXTRA-FIX/fix_numdiff_virial.h similarity index 73% rename from src/EXTRA-FIX/fix_numdiff_stress.h rename to src/EXTRA-FIX/fix_numdiff_virial.h index 4ca0c533ad..69eae142a9 100644 --- a/src/EXTRA-FIX/fix_numdiff_stress.h +++ b/src/EXTRA-FIX/fix_numdiff_virial.h @@ -13,21 +13,21 @@ #ifdef FIX_CLASS // clang-format off -FixStyle(numdiff/stress,FixNumDiffStress); +FixStyle(numdiff/virial,FixNumDiffVirial); // clang-format on #else -#ifndef LMP_FIX_NUMDIFF_STRESS_H -#define LMP_FIX_NUMDIFF_STRESS_H +#ifndef LMP_FIX_NUMDIFF_VIRIAL_H +#define LMP_FIX_NUMDIFF_VIRIAL_H #include "fix.h" namespace LAMMPS_NS { -class FixNumDiffStress : public Fix { +class FixNumDiffVirial : public Fix { public: - FixNumDiffStress(class LAMMPS *, int, char **); - ~FixNumDiffStress(); + FixNumDiffVirial(class LAMMPS *, int, char **); + ~FixNumDiffVirial(); int setmask() override; void init() override; void setup(int) override; @@ -38,7 +38,7 @@ class FixNumDiffStress : public Fix { double memory_usage() override; private: - static const int NDIR_STRESS = 6; // dimension of stress and strain vectors + static const int NDIR_VIRIAL = 6; // dimension of virial and strain vectors double delta; // strain magnitude int maxatom; // allocated size of atom arrays int ilevel_respa; @@ -49,18 +49,18 @@ class FixNumDiffStress : public Fix { char *id_pe; // name of energy compute class Compute *pe; // pointer to energy compute - double stress[NDIR_STRESS]; // finite diff stress components (Voigt order) + 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_STRESS][2];// strain cartesian indices (Voigt order) + int dirlist[NDIR_VIRIAL][2];// strain cartesian indices (Voigt order) - double compute_vector(int) override; // access function for stress - void calculate_stress(); // stress calculation + 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 stress_clear(); // set stress to zero + void virial_clear(); // set virial to zero void reallocate(); // grow the atom arrays }; @@ -77,8 +77,13 @@ 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/stress does not exist +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. + */