Merge remote-tracking branch 'github/develop' into collected-small-changes

This commit is contained in:
Axel Kohlmeyer
2025-05-06 15:06:43 -04:00
23 changed files with 1839 additions and 532 deletions

View File

@ -77,6 +77,7 @@ OPT.
* :doc:`flow/gauss <fix_flow_gauss>`
* :doc:`freeze (k) <fix_freeze>`
* :doc:`gcmc <fix_gcmc>`
* :doc:`gjf <fix_gjf>`
* :doc:`gld <fix_gld>`
* :doc:`gle <fix_gle>`
* :doc:`gravity (ko) <fix_gravity>`

View File

@ -12,6 +12,17 @@ stop LAMMPS and print a suitable error message in most cases, when a
style/command is used that has been removed or will replace the command
with the direct alternative (if available) and print a warning.
GJF formulation in fix langevin
-------------------------------
.. deprecated:: TBD
The *gjf* keyword in fix langevin is deprecated and will be removed
soon. The GJF functionality has been moved to its own fix style
:doc:`fix gjf <fix_gjf>` and it is strongly recommended to use that
fix instead.
LAMMPS shell
------------

View File

@ -21,9 +21,14 @@ can be invoked via the *dpd/tstat* pair style:
* :doc:`fix nvt/sllod <fix_nvt_sllod>`
* :doc:`fix temp/berendsen <fix_temp_berendsen>`
* :doc:`fix temp/csvr <fix_temp_csvr>`
* :doc:`fix ffl <fix_ffl>`
* :doc:`fix gjf <fix_gjf>`
* :doc:`fix gld <fix_gld>`
* :doc:`fix gle <fix_gle>`
* :doc:`fix langevin <fix_langevin>`
* :doc:`fix temp/rescale <fix_temp_rescale>`
* :doc:`pair_style dpd/tstat <pair_dpd>`
* :doc:`pair_style dpd/ext/tstat <pair_dpd_ext>`
:doc:`Fix nvt <fix_nh>` only thermostats the translational velocity of
particles. :doc:`Fix nvt/sllod <fix_nvt_sllod>` also does this,
@ -82,10 +87,10 @@ that:
.. note::
Only the nvt fixes perform time integration, meaning they update
Not all thermostat fixes perform time integration, meaning they update
the velocities and positions of particles due to forces and velocities
respectively. The other thermostat fixes only adjust velocities; they
do NOT perform time integration updates. Thus they should be used in
do NOT perform time integration updates. Thus, they should be used in
conjunction with a constant NVE integration fix such as these:
* :doc:`fix nve <fix_nve>`

View File

@ -256,6 +256,7 @@ accelerated styles exist.
* :doc:`flow/gauss <fix_flow_gauss>` - Gaussian dynamics for constant mass flux
* :doc:`freeze <fix_freeze>` - freeze atoms in a granular simulation
* :doc:`gcmc <fix_gcmc>` - grand canonical insertions/deletions
* :doc:`gjf <fix_gjf>` - statistically correct Langevin temperature control using the GJ methods
* :doc:`gld <fix_gld>` - generalized Langevin dynamics integrator
* :doc:`gle <fix_gle>` - generalized Langevin equation thermostat
* :doc:`gravity <fix_gravity>` - add gravity to atoms in a granular simulation

207
doc/src/fix_gjf.rst Normal file
View File

@ -0,0 +1,207 @@
.. index:: fix gjf
fix gjf command
========================
Syntax
""""""
.. code-block:: LAMMPS
fix ID group-ID gjf Tstart Tstop damp seed keyword values ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* gjf = style name of this fix command
* Tstart,Tstop = desired temperature at start/end of run (temperature units)
* Tstart can be a variable (see below)
* damp = damping parameter (time units)
* seed = random number seed to use for white noise (positive integer)
* zero or more keyword/value pairs may be appended
* keyword = *vel* or *method*
.. parsed-literal::
*vel* value = *vfull* or *vhalf*
*vfull* = use on-site velocity
*vhalf* = use half-step velocity
*method* value = *1-8*
*1-8* = choose one of the many GJ formulations
*7* = requires input of additional scalar between 0 and 1
Examples
""""""""
.. code-block:: LAMMPS
fix 3 boundary gjf 10.0 10.0 1.0 699483
fix 1 all gjf 10.0 100.0 100.0 48279 vel vfull method 4
fix 2 all gjf 10.0 10.0 1.0 26488 method 7 0.95
Description
"""""""""""
Apply a Langevin thermostat as described in :ref:`(Gronbech-Jensen-2020) <Gronbech-Jensen-2020>`
to a group of atoms which models an interaction with a background
implicit solvent. As described in the papers cited below, the GJ methods
provide exact diffusion, drift, and Boltzmann sampling for linear systems for
any time step within the stability limit. The purpose of this set of methods
is therefore to significantly improve statistical accuracy at longer time steps
compared to other thermostats.
The current implementation provides the user with the option to output
the velocity in one of two forms: *vfull* or *vhalf*. The option *vhalf*
outputs the 2GJ half-step velocity given in :ref:`Gronbech Jensen/Gronbech-Jensen
<Gronbech-Jensen-2019>`; for linear systems, this velocity is shown to not
have any statistical errors for any stable time step. The option *vfull*
outputs the on-site velocity given in :ref:`Gronbech-Jensen/Farago
<Gronbech-Jensen-Farago>`; this velocity is shown to be systematically lower
than the target temperature by a small amount, which grows
quadratically with the timestep. An overview of statistically correct Boltzmann
and Maxwell-Boltzmann sampling of true on-site and true half-step velocities is
given in :ref:`Gronbech-Jensen-2020 <Gronbech-Jensen-2020>`.
This fix allows the use of several GJ methods as listed in :ref:`Gronbech-Jensen-2020 <Gronbech-Jensen-2020>`.
The GJ-VII method is described in :ref:`Finkelstein <Finkelstein>` and GJ-VIII
is described in :ref:`Gronbech-Jensen-2024 <Gronbech-Jensen-2024>`.
The implementation follows the splitting form provided in Eqs. (24) and (25)
in :ref:`Gronbech-Jensen-2024 <Gronbech-Jensen-2024>`, including the application
of Gaussian noise values, per the description in
:ref:`Gronbech-Jensen-2023 <Gronbech-Jensen-2023>`.
.. note::
Unlike the :doc:`fix langevin <fix_langevin>` command which performs force
modifications only, this fix performs thermostatting and time integration.
Thus you no longer need a separate time integration fix, like :doc:`fix nve <fix_nve>`.
See the :doc:`Howto thermostat <Howto_thermostat>` page for
a discussion of different ways to compute temperature and perform
thermostatting.
The desired temperature at each timestep is a ramped value during the
run from *Tstart* to *Tstop*\ .
*Tstart* can be specified as an equal-style or atom-style
:doc:`variable <variable>`. In this case, the *Tstop* setting is
ignored. If the value is a variable, it should be specified as
v_name, where name is the variable name. In this case, the variable
will be evaluated each timestep, and its value used to determine the
target temperature.
Equal-style variables can specify formulas with various mathematical
functions, and include :doc:`thermo_style <thermo_style>` command
keywords for the simulation box parameters and timestep and elapsed
time. Thus it is easy to specify a time-dependent temperature.
Atom-style variables can specify the same formulas as equal-style
variables but can also include per-atom values, such as atom
coordinates. Thus it is easy to specify a spatially-dependent
temperature with optional time-dependence as well.
Like other fixes that perform thermostatting, this fix can be used
with :doc:`compute commands <compute>` that remove a "bias" from the
atom velocities. E.g. to apply the thermostat only to atoms within a
spatial :doc:`region <region>`, or to remove the center-of-mass
velocity from a group of atoms, or to remove the x-component of
velocity from the calculation.
This is not done by default, but only if the :doc:`fix_modify
<fix_modify>` command is used to assign a temperature compute to this
fix that includes such a bias term. See the doc pages for individual
:doc:`compute temp commands <compute>` to determine which ones include
a bias.
The *damp* parameter is specified in time units and determines how
rapidly the temperature is relaxed. For example, a value of 100.0 means
to relax the temperature in a timespan of (roughly) 100 time units
(:math:`\tau` or fs or ps - see the :doc:`units <units>` command). The
damp factor can be thought of as inversely related to the viscosity of
the solvent. I.e. a small relaxation time implies a high-viscosity
solvent and vice versa. See the discussion about :math:`\gamma` and
viscosity in the documentation for the :doc:`fix viscous <fix_viscous>`
command for more details.
The random # *seed* must be a positive integer. A Marsaglia random
number generator is used. Each processor uses the input seed to
generate its own unique seed and its own stream of random numbers.
Thus the dynamics of the system will not be identical on two runs on
different numbers of processors.
----------
The keyword/value option pairs are used in the following ways.
The keyword *vel* determines which velocity is used to determine
quantities of interest in the simulation.
The keyword *method* selects one of the eight GJ-methods implemented in LAMMPS.
----------
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`.
Because the state of the random number generator is not saved in restart files,
this means you cannot do "exact" restarts with this fix, where the simulation
continues on the same as if no restart had taken place. However, in a
statistical sense, a restarted simulation should produce the same behavior.
Additionally, the GJ methods implement noise exclusively within each time step
(unlike the BBK thermostat of the fix-langevin). The restart is done with
either vfull or vhalf velocity output for as long as the choice of vfull/vhalf
is the same for the simulation as it is in the restart file.
The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
fix. You can use it to assign a temperature :doc:`compute <compute>`
you have defined to this fix which will be used in its thermostatting
procedure, as described above. For consistency, the group used by
this fix and by the compute should be the same.
This fix can ramp its target temperature over multiple runs, using the
*start* and *stop* keywords of the :doc:`run <run>` command. See the
:doc:`run <run>` command for details of how to do this.
This fix is not invoked during :doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix is not compatible with run_style respa. It is not compatible with
accelerated packages such as KOKKOS.
Related commands
""""""""""""""""
:doc:`fix langevin <fix_langevin>`, :doc:`fix nvt <fix_nh>`
Default
"""""""
The option defaults are vel = vhalf, method = 1.
----------
.. _Gronbech-Jensen-2020:
**(Gronbech-Jensen-2020)** Gronbech-Jensen, Mol Phys 118, e1662506 (2020).
.. _Gronbech-Jensen-2019:
**(Gronbech Jensen/Gronbech-Jensen)** Gronbech Jensen and Gronbech-Jensen, Mol Phys, 117, 2511 (2019)
.. _Gronbech-Jensen-Farago:
**(Gronbech-Jensen/Farago)** Gronbech-Jensen and Farago, Mol Phys, 111, 983 (2013).
.. _Finkelstein:
**(Finkelstein)** Finkelstein, Cheng, Florin, Seibold, Gronbech-Jensen, J. Chem. Phys., 155, 18 (2021)
.. _Gronbech-Jensen-2024:
**(Gronbech-Jensen-2024)** Gronbech-Jensen, J. Stat. Phys. 191, 137 (2024).
.. _Gronbech-Jensen-2023:
**(Gronbech-Jensen-2023)** Gronbech-Jensen, J. Stat. Phys. 190, 96 (2023).

View File

@ -56,7 +56,7 @@ Examples
Description
"""""""""""
Apply a Langevin thermostat as described in :ref:`(Schneider) <Schneider1>`
Apply a Langevin thermostat as described in :ref:`(Bruenger) <Bruenger1>`
to a group of atoms which models an interaction with a background
implicit solvent. Used with :doc:`fix nve <fix_nve>`, this command
performs Brownian dynamics (BD), since the total force on each atom
@ -241,6 +241,13 @@ to zero by subtracting off an equal part of it from each atom in the
group. As a result, the center-of-mass of a system with zero initial
momentum will not drift over time.
.. deprecated:: TDB
The *gjf* keyword in fix langevin is deprecated and will be removed
soon. The GJF functionality has been moved to its own fix style
:doc:`fix gjf <fix_gjf>` and it is strongly recommended to use that
fix instead.
The keyword *gjf* can be used to run the :ref:`Gronbech-Jensen/Farago
<Gronbech-Jensen>` time-discretization of the Langevin model. As
described in the papers cited below, the purpose of this method is to
@ -324,14 +331,16 @@ types, tally = no, zero = no, gjf = no.
----------
.. _Bruenger1:
**(Bruenger)** Bruenger, Brooks, and Karplus, Chem. Phys. Lett. 105, 495 (1982).
[Previously attributed to Schneider and Stoll, Phys. Rev. B 17, 1302 (1978).
Implementation remains unchanged.]
.. _Dunweg1:
**(Dunweg)** Dunweg and Paul, Int J of Modern Physics C, 2, 817-27 (1991).
.. _Schneider1:
**(Schneider)** Schneider and Stoll, Phys Rev B, 17, 1302 (1978).
.. _Gronbech-Jensen:
**(Gronbech-Jensen)** Gronbech-Jensen and Farago, Mol Phys, 111, 983

View File

@ -1200,6 +1200,7 @@ filesystem
filesystems
Fily
Fincham
Finkelstein
Fint
fingerprintconstants
fingerprintsperelement
@ -1354,6 +1355,7 @@ Gingold
Gissinger
github
Giusti
GJ
gjf
gjwagne
gl
@ -3473,6 +3475,7 @@ sectoring
sed
Seddon
segmental
Seibold
Seifert
Seleson
sellerio
@ -4092,9 +4095,11 @@ versa
Verstraelen
ves
vf
vfull
vflag
vflow
vfrac
vhalf
vhi
vibrational
Vij

47
examples/gjf/README Normal file
View File

@ -0,0 +1,47 @@
LAMMPS GJ THERMOSTAT EXAMPLE
Required LAMMPS packages: EXTRA-FIX, MOLECULE, EXTRA-PAIR
This directory contains the ingredients to run an NVT simulation using the
GJ thermostats.
Example:
NP=4 #number of processors
mpirun -np $NP lmp_mpi -in.gjf.vhalf
Compared to other thermostats, the GJ thermostat allows for larger timesteps
with the correct Boltzmann statistics. A comparison using averaged properties
from this example's input file is shown below. 'X' denotes a failed simulation.
The theoretical value for KE is 1.1168 eV.
POTENTIAL ENERGY (eV)
| Δt || 0.01 | 0.05 | 0.10 | 0.11 | 0.12 | 0.13 | 0.14 |
|===================||========|========|========|========|========|========|========|
| gjf half || -55.11 | -55.11 | -55.11 | -55.11 | -55.11 | -55.10 | -55.07 |
| gjf full || -55.11 | -55.11 | -55.11 | -55.11 | -55.11 | -55.10 | -55.07 |
| langevin || -55.11 | -55.07 | -54.87 | -54.79 | -54.65 | X | X |
| nvt (Nose-Hoover) || -55.14 | -55.07 | -54.90 | -54.84 | -54.76 | X | X |
|-------------------||--------|--------|--------|--------|--------|--------|--------|
KINETIC ENERGY (eV)
| Δt || 0.01 | 0.05 | 0.10 | 0.11 | 0.12 | 0.13 | 0.14 |
|===================||========|========|========|========|========|========|========|
| gjf half || 1.117 | 1.116 | 1.119 | 1.119 | 1.123 | 1.136 | 1.170 |
| gjf full || 1.116 | 1.071 | 0.938 | 0.898 | 0.858 | 0.817 | 0.780 |
| langevin || 1.110 | 1.113 | 1.121 | 1.129 | 1.157 | X | X |
| nvt (Nose-Hoover) || 1.083 | 1.109 | 1.112 | 1.113 | 1.114 | X | X |
|-------------------||--------|--------|--------|--------|--------|--------|--------|
Script Commands:
--
fix lang all gjf 10 10 1 26488
--
fix lang all gjf 10 10 1 26488 vel vfull
--
fix nve all nve
fix lang all langevin 10 10 1 26488
--
fix noho all nvt temp 10 10 1
--

View File

@ -1,13 +0,0 @@
# LAMMPS GJF-2GJ THERMOSTAT EXAMPLE
## GJF-2GJ THERMOSTAT
This directory contains the ingredients to run an NVT simulation using the GJF-2GJ thermostat.
Example:
```
NP=4 #number of processors
mpirun -np $NP lmp_mpi -in.gjf.vhalf
```
## Required LAMMPS packages: MOLECULE package

View File

@ -1,23 +1,25 @@
# GJF-2GJ thermostat
# GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
include ff-argon.lmp
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
fix lang all langevin 10 10 1 26488 gjf vfull
fix nve all nve
compute myKE all ke
compute myPE all pe
fix lang all gjf 10 10 1 26488 vel vfull method 1
thermo 200
run 5000
fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out
thermo 2000
run 35000

View File

@ -1,23 +1,25 @@
# GJF-2GJ thermostat
# GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
include ff-argon.lmp
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
fix lang all langevin 10 10 1 26488 gjf vhalf
fix nve all nve
compute myKE all ke
compute myPE all pe
fix lang all gjf 10 10 1 26488
thermo 200
run 5000
fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out
thermo 2000
run 35000

View File

@ -1,125 +0,0 @@
LAMMPS (19 Sep 2019)
using 1 OpenMP thread(s) per MPI task
# GJF-2GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
orthogonal box = (0 0 0) to (32.146 32.146 32.146)
1 by 1 by 1 MPI processor grid
reading atoms ...
864 atoms
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000150019 secs
read_data CPU = 0.001946 secs
include ff-argon.lmp
#############################
#Atoms types - mass - charge#
#############################
#@ 1 atom types #!THIS LINE IS NECESSARY DON'T SPEND HOURS FINDING THAT OUT!#
variable Ar equal 1
#############
#Atom Masses#
#############
mass ${Ar} 39.903
mass 1 39.903
###########################
#Pair Potentials - Tersoff#
###########################
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
fix lang all langevin 10 10 1 26488 gjf vfull
fix nve all nve
thermo 200
run 5000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.94072
ghost atom cutoff = 6.94072
binsize = 3.47036, bins = 10 10 10
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cubic, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.875 | 6.875 | 6.875 Mbytes
Step Temp E_pair E_mol TotEng Press
0 11.080223 -56.207655 0 -54.97164 37.215524
200 8.2588471 -55.073602 0 -54.152316 339.80416
400 8.1427292 -55.072244 0 -54.16391 338.91883
600 8.7595618 -55.066739 0 -54.089596 344.25426
800 8.550633 -55.148315 0 -54.194479 318.9385
1000 8.5394337 -55.125709 0 -54.173122 326.59471
1200 8.565973 -55.114892 0 -54.159345 328.5193
1400 8.2092914 -55.109233 0 -54.193475 329.56161
1600 8.209495 -55.138161 0 -54.22238 321.39971
1800 8.4039924 -55.13355 0 -54.196072 322.64214
2000 8.4548937 -55.062994 0 -54.119838 343.29888
2200 8.3775139 -55.13364 0 -54.199116 323.63744
2400 8.537332 -55.163702 0 -54.21135 315.62864
2600 8.672488 -55.112054 0 -54.144625 330.1106
2800 8.3000218 -55.147275 0 -54.221396 318.73112
3000 8.3552421 -55.135164 0 -54.203124 323.53075
3200 8.4126798 -55.135753 0 -54.197306 321.48817
3400 8.4986413 -55.135408 0 -54.187372 323.42951
3600 8.38431 -55.103932 0 -54.16865 330.68929
3800 8.8262454 -55.103648 0 -54.119067 332.97779
4000 7.9658136 -55.120402 0 -54.231803 324.9595
4200 8.2265544 -55.129011 0 -54.211327 323.87069
4400 8.1253738 -55.153089 0 -54.246691 316.304
4600 8.2010823 -55.124053 0 -54.20921 325.98402
4800 8.5512149 -55.075877 0 -54.121976 338.30137
5000 8.4737659 -55.158604 0 -54.213343 316.22418
Loop time of 2.73236 on 1 procs for 5000 steps with 864 atoms
Performance: 15810.507 ns/day, 0.002 hours/ns, 1829.920 timesteps/s
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.4262 | 1.4262 | 1.4262 | 0.0 | 52.20
Bond | 0.00042836 | 0.00042836 | 0.00042836 | 0.0 | 0.02
Neigh | 0.12819 | 0.12819 | 0.12819 | 0.0 | 4.69
Comm | 0.058611 | 0.058611 | 0.058611 | 0.0 | 2.15
Output | 0.00047283 | 0.00047283 | 0.00047283 | 0.0 | 0.02
Modify | 1.0924 | 1.0924 | 1.0924 | 0.0 | 39.98
Other | | 0.02605 | | | 0.95
Nlocal: 864 ave 864 max 864 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1593 ave 1593 max 1593 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 18143 ave 18143 max 18143 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 18143
Ave neighs/atom = 20.9988
Ave special neighs/atom = 0
Neighbor list builds = 158
Dangerous builds = 5
Total wall time: 0:00:02

View File

@ -1,125 +0,0 @@
LAMMPS (19 Sep 2019)
using 1 OpenMP thread(s) per MPI task
# GJF-2GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
orthogonal box = (0 0 0) to (32.146 32.146 32.146)
1 by 2 by 2 MPI processor grid
reading atoms ...
864 atoms
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000556268 secs
read_data CPU = 0.003817 secs
include ff-argon.lmp
#############################
#Atoms types - mass - charge#
#############################
#@ 1 atom types #!THIS LINE IS NECESSARY DON'T SPEND HOURS FINDING THAT OUT!#
variable Ar equal 1
#############
#Atom Masses#
#############
mass ${Ar} 39.903
mass 1 39.903
###########################
#Pair Potentials - Tersoff#
###########################
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
fix lang all langevin 10 10 1 26488 gjf vfull
fix nve all nve
thermo 200
run 5000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.94072
ghost atom cutoff = 6.94072
binsize = 3.47036, bins = 10 10 10
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cubic, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.808 | 6.808 | 6.808 Mbytes
Step Temp E_pair E_mol TotEng Press
0 11.080228 -56.207655 0 -54.971639 37.215541
200 8.4818184 -55.127334 0 -54.181174 324.96159
400 8.5960916 -55.09236 0 -54.133453 334.83136
600 8.1607556 -55.073136 0 -54.162791 339.035
800 8.8350489 -55.133382 0 -54.147819 324.48149
1000 8.5692704 -55.118463 0 -54.162548 327.26328
1200 8.4174147 -55.126297 0 -54.187322 324.4248
1400 8.6362603 -55.123075 0 -54.159688 326.7798
1600 8.222512 -55.153799 0 -54.236565 317.8147
1800 8.324523 -55.116698 0 -54.188085 327.35373
2000 7.9615959 -55.155825 0 -54.267697 315.37215
2200 8.495968 -55.083943 0 -54.136205 336.67775
2400 7.7926986 -55.044816 0 -54.175529 344.87758
2600 8.1551351 -55.069404 0 -54.159687 339.60901
2800 8.2593599 -55.084151 0 -54.162807 336.54935
3000 8.2860869 -55.110296 0 -54.185971 328.99074
3200 8.4074534 -55.123576 0 -54.185712 326.06823
3400 8.6694364 -55.128925 0 -54.161836 324.67512
3600 8.5718984 -55.129861 0 -54.173653 325.20586
3800 8.508102 -55.099093 0 -54.150001 333.91437
4000 8.2966658 -55.117782 0 -54.192276 327.13516
4200 8.7641728 -55.135792 0 -54.158136 324.00844
4400 8.8827909 -55.096369 0 -54.10548 335.08467
4600 8.7666577 -55.127213 0 -54.149279 326.15539
4800 8.6670762 -55.163395 0 -54.19657 316.48383
5000 8.1893094 -55.073756 0 -54.160226 337.95271
Loop time of 0.870594 on 4 procs for 5000 steps with 864 atoms
Performance: 49621.267 ns/day, 0.000 hours/ns, 5743.202 timesteps/s
96.5% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.33582 | 0.35125 | 0.3724 | 2.3 | 40.35
Bond | 0.00030267 | 0.00031316 | 0.00033538 | 0.0 | 0.04
Neigh | 0.034246 | 0.03479 | 0.035904 | 0.4 | 4.00
Comm | 0.15068 | 0.17419 | 0.19191 | 3.6 | 20.01
Output | 0.00044776 | 0.00054703 | 0.00083177 | 0.0 | 0.06
Modify | 0.27679 | 0.28079 | 0.28849 | 0.9 | 32.25
Other | | 0.02871 | | | 3.30
Nlocal: 216 ave 216 max 216 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 888.75 ave 899 max 876 min
Histogram: 1 0 1 0 0 0 0 0 0 2
Neighs: 4536 ave 4737 max 4335 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 18144
Ave neighs/atom = 21
Ave special neighs/atom = 0
Neighbor list builds = 178
Dangerous builds = 11
Total wall time: 0:00:00

View File

@ -1,125 +0,0 @@
LAMMPS (19 Sep 2019)
using 1 OpenMP thread(s) per MPI task
# GJF-2GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
orthogonal box = (0 0 0) to (32.146 32.146 32.146)
1 by 1 by 1 MPI processor grid
reading atoms ...
864 atoms
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000147804 secs
read_data CPU = 0.00194898 secs
include ff-argon.lmp
#############################
#Atoms types - mass - charge#
#############################
#@ 1 atom types #!THIS LINE IS NECESSARY DON'T SPEND HOURS FINDING THAT OUT!#
variable Ar equal 1
#############
#Atom Masses#
#############
mass ${Ar} 39.903
mass 1 39.903
###########################
#Pair Potentials - Tersoff#
###########################
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
fix lang all langevin 10 10 1 26488 gjf vhalf
fix nve all nve
thermo 200
run 5000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.94072
ghost atom cutoff = 6.94072
binsize = 3.47036, bins = 10 10 10
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cubic, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.5 | 6.5 | 6.5 Mbytes
Step Temp E_pair E_mol TotEng Press
0 11.080223 -56.207655 0 -54.97164 37.215524
200 9.8808568 -55.073602 0 -53.971378 345.62207
400 9.8712816 -55.072244 0 -53.971088 345.11889
600 10.528988 -55.066739 0 -53.892214 350.60093
800 10.167171 -55.148315 0 -54.014152 324.73679
1000 10.029026 -55.125709 0 -54.006956 331.93766
1200 10.194424 -55.114892 0 -53.977688 334.36032
1400 9.3473846 -55.109233 0 -54.066518 333.64378
1600 9.7774071 -55.138161 0 -54.047477 327.02358
1800 9.9814275 -55.13355 0 -54.020107 328.30017
2000 10.2515 -55.062994 0 -53.919424 349.74304
2200 9.8126922 -55.13364 0 -54.039019 328.78521
2400 10.044314 -55.163702 0 -54.043244 321.03397
2600 10.543316 -55.112054 0 -53.935932 336.82099
2800 9.7874375 -55.147275 0 -54.055472 324.06626
3000 9.7703821 -55.135164 0 -54.045263 328.60665
3200 10.141958 -55.135753 0 -54.004402 327.69084
3400 10.160576 -55.135408 0 -54.00198 329.39063
3600 10.044652 -55.103932 0 -53.983436 336.64469
3800 10.662403 -55.103648 0 -53.914241 339.56382
4000 9.2921047 -55.120402 0 -54.083853 329.71671
4200 9.8744553 -55.129011 0 -54.027501 329.78147
4400 9.4085964 -55.153089 0 -54.103546 320.90673
4600 9.5463801 -55.124053 0 -54.05914 330.80941
4800 10.223884 -55.075877 0 -53.935387 344.30099
5000 9.6243338 -55.158604 0 -54.084996 320.3511
Loop time of 2.29551 on 1 procs for 5000 steps with 864 atoms
Performance: 18819.358 ns/day, 0.001 hours/ns, 2178.166 timesteps/s
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.4393 | 1.4393 | 1.4393 | 0.0 | 62.70
Bond | 0.0004441 | 0.0004441 | 0.0004441 | 0.0 | 0.02
Neigh | 0.12136 | 0.12136 | 0.12136 | 0.0 | 5.29
Comm | 0.059342 | 0.059342 | 0.059342 | 0.0 | 2.59
Output | 0.00046968 | 0.00046968 | 0.00046968 | 0.0 | 0.02
Modify | 0.64937 | 0.64937 | 0.64937 | 0.0 | 28.29
Other | | 0.02522 | | | 1.10
Nlocal: 864 ave 864 max 864 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1593 ave 1593 max 1593 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 18143 ave 18143 max 18143 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 18143
Ave neighs/atom = 20.9988
Ave special neighs/atom = 0
Neighbor list builds = 158
Dangerous builds = 5
Total wall time: 0:00:02

View File

@ -1,125 +0,0 @@
LAMMPS (19 Sep 2019)
using 1 OpenMP thread(s) per MPI task
# GJF-2GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
orthogonal box = (0 0 0) to (32.146 32.146 32.146)
1 by 2 by 2 MPI processor grid
reading atoms ...
864 atoms
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000315903 secs
read_data CPU = 0.0653752 secs
include ff-argon.lmp
#############################
#Atoms types - mass - charge#
#############################
#@ 1 atom types #!THIS LINE IS NECESSARY DON'T SPEND HOURS FINDING THAT OUT!#
variable Ar equal 1
#############
#Atom Masses#
#############
mass ${Ar} 39.903
mass 1 39.903
###########################
#Pair Potentials - Tersoff#
###########################
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
fix lang all langevin 10 10 1 26488 gjf vhalf
fix nve all nve
thermo 200
run 5000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.94072
ghost atom cutoff = 6.94072
binsize = 3.47036, bins = 10 10 10
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cubic, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.433 | 6.433 | 6.433 Mbytes
Step Temp E_pair E_mol TotEng Press
0 11.080228 -56.207655 0 -54.971639 37.215541
200 9.8046716 -55.127334 0 -54.033608 329.70647
400 10.174622 -55.09236 0 -53.957366 340.49331
600 9.9812299 -55.073136 0 -53.959714 345.56477
800 10.512874 -55.133382 0 -53.960655 330.4996
1000 9.9587885 -55.118463 0 -54.007545 332.24728
1200 10.236607 -55.126297 0 -53.984388 330.94998
1400 10.134679 -55.123075 0 -53.992537 332.15441
1600 9.8934078 -55.153799 0 -54.050174 323.80795
1800 10.064966 -55.116698 0 -53.993936 333.59644
2000 9.6736107 -55.155825 0 -54.076719 321.5129
2200 10.264537 -55.083943 0 -53.938918 343.02135
2400 9.5640032 -55.044816 0 -53.977937 351.23099
2600 9.6581077 -55.069404 0 -53.992028 344.99996
2800 9.9622575 -55.084151 0 -53.972846 342.6574
3000 9.8724909 -55.110296 0 -54.009005 334.68094
3200 10.032027 -55.123576 0 -54.004488 331.89534
3400 10.221132 -55.128925 0 -53.988742 330.24082
3600 10.085802 -55.129861 0 -54.004774 330.63601
3800 10.098545 -55.099093 0 -53.972585 339.61905
4000 10.000257 -55.117782 0 -54.002238 333.24569
4200 10.20477 -55.135792 0 -53.997435 329.17565
4400 10.545132 -55.096369 0 -53.920044 341.04725
4600 10.376108 -55.127213 0 -53.969743 331.92825
4800 10.247392 -55.163395 0 -54.020283 322.15219
5000 9.7753102 -55.073756 0 -53.983305 343.64146
Loop time of 1.19785 on 4 procs for 5000 steps with 864 atoms
Performance: 36064.674 ns/day, 0.001 hours/ns, 4174.152 timesteps/s
88.6% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.36387 | 0.38652 | 0.44086 | 5.1 | 32.27
Bond | 0.00028847 | 0.00030833 | 0.000338 | 0.0 | 0.03
Neigh | 0.033934 | 0.034959 | 0.036917 | 0.6 | 2.92
Comm | 0.39292 | 0.47821 | 0.52198 | 7.3 | 39.92
Output | 0.00050343 | 0.0012343 | 0.0023338 | 1.9 | 0.10
Modify | 0.1605 | 0.17963 | 0.19457 | 2.9 | 15.00
Other | | 0.117 | | | 9.77
Nlocal: 216 ave 216 max 216 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 888.75 ave 899 max 876 min
Histogram: 1 0 1 0 0 0 0 0 0 2
Neighs: 4536 ave 4737 max 4335 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 18144
Ave neighs/atom = 21
Ave special neighs/atom = 0
Neighbor list builds = 178
Dangerous builds = 11
Total wall time: 0:00:01

View File

@ -0,0 +1,193 @@
LAMMPS (2 Apr 2025 - Development - d4867ab55e-modified)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:99)
using 1 OpenMP thread(s) per MPI task
# GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
Reading data file ...
orthogonal box = (0 0 0) to (32.146 32.146 32.146)
1 by 1 by 1 MPI processor grid
reading atoms ...
864 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.007 seconds
include ff-argon.lmp
#############################
#Atoms types - mass - charge#
#############################
#@ 1 atom types #!THIS LINE IS NECESSARY DON'T SPEND HOURS FINDING THAT OUT!#
variable Ar equal 1
#############
#Atom Masses#
#############
mass ${Ar} 39.903
mass 1 39.903
###########################
#Pair Potentials - Tersoff#
###########################
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
compute myKE all ke
compute myPE all pe
fix lang all gjf 10 10 1 26488 vel vfull method 1
run 5000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- GJ methods: doi:10.1080/00268976.2019.1662506
@Article{gronbech-jensen_complete_2020,
title = {Complete set of stochastic Verlet-type thermostats for correct Langevin simulations},
volume = {118},
number = {8},
url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1662506},
doi = {10.1080/00268976.2019.1662506},
journal = {Molecular Physics},
author = {Grønbech-Jensen, Niels},
year = {2020}
}
- GJ-I vfull method: doi:10.1080/00268976.2012.760055
@Article{gronbech-jensen_simple_2013,
title = {A simple and effective Verlet-type algorithm for simulating Langevin dynamics},
volume = {111},
url = {http://www.tandfonline.com/doi/abs/10.1080/00268976.2012.760055},
doi = {10.1080/00268976.2012.760055},
pages = {983-991},
number = {8},
journal = {Molecular Physics},
author = {Grønbech-Jensen, Niels and Farago, Oded},
year = {2013}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.9407173
ghost atom cutoff = 6.9407173
binsize = 3.4703587, bins = 10 10 10
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cubic, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.481 | 6.481 | 6.481 Mbytes
Step Temp E_pair E_mol TotEng Press
0 10 -56.207652 0 -55.092137 33.341103
5000 8.4535562 -55.150518 0 -54.207511 318.20862
Loop time of 2.26831 on 1 procs for 5000 steps with 864 atoms
Performance: 19044.977 ns/day, 0.001 hours/ns, 2204.280 timesteps/s, 1.904 Matom-step/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.2802 | 1.2802 | 1.2802 | 0.0 | 56.44
Bond | 0.00051213 | 0.00051213 | 0.00051213 | 0.0 | 0.02
Neigh | 0.27007 | 0.27007 | 0.27007 | 0.0 | 11.91
Comm | 0.057527 | 0.057527 | 0.057527 | 0.0 | 2.54
Output | 6.3876e-05 | 6.3876e-05 | 6.3876e-05 | 0.0 | 0.00
Modify | 0.63364 | 0.63364 | 0.63364 | 0.0 | 27.93
Other | | 0.02635 | | | 1.16
Nlocal: 864 ave 864 max 864 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1593 ave 1593 max 1593 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 18143 ave 18143 max 18143 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 18143
Ave neighs/atom = 20.998843
Ave special neighs/atom = 0
Neighbor list builds = 258
Dangerous builds = 0
fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out
thermo 2000
run 35000
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 6.481 | 6.481 | 6.481 Mbytes
Step Temp E_pair E_mol TotEng Press
5000 8.4535562 -55.150518 0 -54.207511 318.20862
6000 8.4899401 -55.108242 0 -54.161176 331.10703
8000 8.3618893 -55.092171 0 -54.15939 334.11831
10000 8.8684311 -55.100316 0 -54.111029 334.09931
12000 8.4339192 -55.07343 0 -54.132614 340.00487
14000 8.072393 -55.115121 0 -54.214633 327.98965
16000 8.3420289 -55.077813 0 -54.147247 337.74926
18000 8.3803911 -55.12201 0 -54.187164 326.10485
20000 8.4676985 -55.176339 0 -54.231754 311.57092
22000 8.8560138 -55.110505 0 -54.122603 330.66179
24000 8.3187826 -55.120592 0 -54.192619 327.01148
26000 8.0327666 -55.116664 0 -54.220596 326.25179
28000 8.3672169 -55.130413 0 -54.197037 324.2368
30000 8.1669275 -55.057678 0 -54.146645 344.9168
32000 8.3819314 -55.08989 0 -54.154873 335.45317
34000 8.109088 -55.17222 0 -54.267639 310.83717
36000 8.3048574 -55.079475 0 -54.153056 338.04291
38000 8.8708544 -55.108991 0 -54.119434 330.70097
40000 8.4012779 -55.080817 0 -54.143642 338.54326
Loop time of 18.9699 on 1 procs for 35000 steps with 864 atoms
Performance: 15941.040 ns/day, 0.002 hours/ns, 1845.028 timesteps/s, 1.594 Matom-step/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 11.593 | 11.593 | 11.593 | 0.0 | 61.11
Bond | 0.0041801 | 0.0041801 | 0.0041801 | 0.0 | 0.02
Neigh | 2.2671 | 2.2671 | 2.2671 | 0.0 | 11.95
Comm | 0.42339 | 0.42339 | 0.42339 | 0.0 | 2.23
Output | 0.00062204 | 0.00062204 | 0.00062204 | 0.0 | 0.00
Modify | 4.4976 | 4.4976 | 4.4976 | 0.0 | 23.71
Other | | 0.1839 | | | 0.97
Nlocal: 864 ave 864 max 864 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1592 ave 1592 max 1592 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 18144 ave 18144 max 18144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 18144
Ave neighs/atom = 21
Ave special neighs/atom = 0
Neighbor list builds = 2122
Dangerous builds = 0
Total wall time: 0:00:21

View File

@ -0,0 +1,193 @@
LAMMPS (2 Apr 2025 - Development - d4867ab55e-modified)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:99)
using 1 OpenMP thread(s) per MPI task
# GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
Reading data file ...
orthogonal box = (0 0 0) to (32.146 32.146 32.146)
1 by 2 by 2 MPI processor grid
reading atoms ...
864 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.002 seconds
read_data CPU = 0.015 seconds
include ff-argon.lmp
#############################
#Atoms types - mass - charge#
#############################
#@ 1 atom types #!THIS LINE IS NECESSARY DON'T SPEND HOURS FINDING THAT OUT!#
variable Ar equal 1
#############
#Atom Masses#
#############
mass ${Ar} 39.903
mass 1 39.903
###########################
#Pair Potentials - Tersoff#
###########################
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
compute myKE all ke
compute myPE all pe
fix lang all gjf 10 10 1 26488 vel vfull method 1
run 5000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- GJ methods: doi:10.1080/00268976.2019.1662506
@Article{gronbech-jensen_complete_2020,
title = {Complete set of stochastic Verlet-type thermostats for correct Langevin simulations},
volume = {118},
number = {8},
url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1662506},
doi = {10.1080/00268976.2019.1662506},
journal = {Molecular Physics},
author = {Grønbech-Jensen, Niels},
year = {2020}
}
- GJ-I vfull method: doi:10.1080/00268976.2012.760055
@Article{gronbech-jensen_simple_2013,
title = {A simple and effective Verlet-type algorithm for simulating Langevin dynamics},
volume = {111},
url = {http://www.tandfonline.com/doi/abs/10.1080/00268976.2012.760055},
doi = {10.1080/00268976.2012.760055},
pages = {983-991},
number = {8},
journal = {Molecular Physics},
author = {Grønbech-Jensen, Niels and Farago, Oded},
year = {2013}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.9407173
ghost atom cutoff = 6.9407173
binsize = 3.4703587, bins = 10 10 10
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cubic, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.427 | 6.427 | 6.427 Mbytes
Step Temp E_pair E_mol TotEng Press
0 10 -56.207652 0 -55.092137 33.341103
5000 7.946377 -55.076514 0 -54.190084 337.31999
Loop time of 2.0998 on 4 procs for 5000 steps with 864 atoms
Performance: 20573.405 ns/day, 0.001 hours/ns, 2381.181 timesteps/s, 2.057 Matom-step/s
65.2% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.53641 | 0.54389 | 0.54721 | 0.6 | 25.90
Bond | 0.00056487 | 0.0006195 | 0.00068462 | 0.0 | 0.03
Neigh | 0.10567 | 0.1086 | 0.11128 | 0.7 | 5.17
Comm | 0.96913 | 0.97758 | 0.98191 | 0.5 | 46.56
Output | 0.00025213 | 0.00025642 | 0.00026405 | 0.0 | 0.01
Modify | 0.25061 | 0.25105 | 0.25172 | 0.1 | 11.96
Other | | 0.2178 | | | 10.37
Nlocal: 216 ave 216 max 216 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 884.75 ave 885 max 884 min
Histogram: 1 0 0 0 0 0 0 0 0 3
Neighs: 4536 ave 4737 max 4335 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 18144
Ave neighs/atom = 21
Ave special neighs/atom = 0
Neighbor list builds = 273
Dangerous builds = 0
fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out
thermo 2000
run 35000
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 6.428 | 6.428 | 6.428 Mbytes
Step Temp E_pair E_mol TotEng Press
5000 7.946377 -55.076514 0 -54.190084 337.31999
6000 8.2565866 -55.129244 0 -54.208209 324.57967
8000 7.9942397 -55.101417 0 -54.209648 331.24127
10000 8.5413968 -55.083292 0 -54.130486 337.82599
12000 8.3682078 -55.090905 0 -54.157419 335.08066
14000 8.5082065 -55.085051 0 -54.135948 336.2765
16000 8.1944037 -55.090733 0 -54.176635 334.03786
18000 8.2607106 -55.030131 0 -54.108637 352.49892
20000 8.1154691 -55.104072 0 -54.198779 330.14203
22000 8.5592601 -55.152019 0 -54.197221 318.03507
24000 8.3182914 -55.115242 0 -54.187324 328.46084
26000 8.3691375 -55.125275 0 -54.191685 325.43673
28000 8.531632 -55.107097 0 -54.155381 331.42771
30000 8.1102222 -55.099011 0 -54.194304 332.04678
32000 8.5558571 -55.077016 0 -54.122598 339.87746
34000 8.4213946 -55.097068 0 -54.157649 333.34935
36000 8.0936615 -55.152202 0 -54.249342 316.20169
38000 7.999652 -55.048407 0 -54.156034 345.07945
40000 8.6699753 -55.087634 0 -54.120485 337.23709
Loop time of 17.6726 on 4 procs for 35000 steps with 864 atoms
Performance: 17111.263 ns/day, 0.001 hours/ns, 1980.470 timesteps/s, 1.711 Matom-step/s
65.4% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.0739 | 5.1178 | 5.1689 | 1.5 | 28.96
Bond | 0.0043764 | 0.004688 | 0.0051706 | 0.4 | 0.03
Neigh | 0.83797 | 0.85506 | 0.87554 | 1.8 | 4.84
Comm | 6.816 | 6.8932 | 6.9215 | 1.7 | 39.00
Output | 0.0043624 | 0.0045336 | 0.004998 | 0.4 | 0.03
Modify | 3.3008 | 3.3033 | 3.3066 | 0.1 | 18.69
Other | | 1.494 | | | 8.45
Nlocal: 216 ave 222 max 210 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 905.5 ave 911 max 899 min
Histogram: 1 1 0 0 0 0 0 0 0 2
Neighs: 4535.75 ave 4837 max 4218 min
Histogram: 1 0 0 1 0 0 1 0 0 1
Total # of neighbors = 18143
Ave neighs/atom = 20.998843
Ave special neighs/atom = 0
Neighbor list builds = 2140
Dangerous builds = 0
Total wall time: 0:00:19

View File

@ -0,0 +1,192 @@
LAMMPS (2 Apr 2025 - Development - d4867ab55e-modified)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:99)
using 1 OpenMP thread(s) per MPI task
# GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
Reading data file ...
orthogonal box = (0 0 0) to (32.146 32.146 32.146)
1 by 1 by 1 MPI processor grid
reading atoms ...
864 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.010 seconds
include ff-argon.lmp
#############################
#Atoms types - mass - charge#
#############################
#@ 1 atom types #!THIS LINE IS NECESSARY DON'T SPEND HOURS FINDING THAT OUT!#
variable Ar equal 1
#############
#Atom Masses#
#############
mass ${Ar} 39.903
mass 1 39.903
###########################
#Pair Potentials - Tersoff#
###########################
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
compute myKE all ke
compute myPE all pe
fix lang all gjf 10 10 1 26488
run 5000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- GJ methods: doi:10.1080/00268976.2019.1662506
@Article{gronbech-jensen_complete_2020,
title = {Complete set of stochastic Verlet-type thermostats for correct Langevin simulations},
volume = {118},
number = {8},
url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1662506},
doi = {10.1080/00268976.2019.1662506},
journal = {Molecular Physics},
author = {Grønbech-Jensen, Niels},
year = {2020}
}
- GJ-I vhalf method: doi:10.1080/00268976.2019.1570369
@Article{jensen_accurate_2019,
title = {Accurate configurational and kinetic statistics in discrete-time Langevin systems},
volume = {117},
url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1570369},
doi = {10.1080/00268976.2019.1570369},
number = {18},
journal = {Molecular Physics},
author = {Jensen, Lucas Frese Grønbech and Grønbech-Jensen, Niels},
year = {2019}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.9407173
ghost atom cutoff = 6.9407173
binsize = 3.4703587, bins = 10 10 10
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cubic, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.481 | 6.481 | 6.481 Mbytes
Step Temp E_pair E_mol TotEng Press
0 10 -56.207652 0 -55.092137 33.341103
5000 9.7731898 -55.150518 0 -54.060304 322.94195
Loop time of 2.28421 on 1 procs for 5000 steps with 864 atoms
Performance: 18912.438 ns/day, 0.001 hours/ns, 2188.940 timesteps/s, 1.891 Matom-step/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.2715 | 1.2715 | 1.2715 | 0.0 | 55.66
Bond | 0.00057126 | 0.00057126 | 0.00057126 | 0.0 | 0.03
Neigh | 0.27008 | 0.27008 | 0.27008 | 0.0 | 11.82
Comm | 0.057938 | 0.057938 | 0.057938 | 0.0 | 2.54
Output | 6.1954e-05 | 6.1954e-05 | 6.1954e-05 | 0.0 | 0.00
Modify | 0.658 | 0.658 | 0.658 | 0.0 | 28.81
Other | | 0.0261 | | | 1.14
Nlocal: 864 ave 864 max 864 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1593 ave 1593 max 1593 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 18143 ave 18143 max 18143 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 18143
Ave neighs/atom = 20.998843
Ave special neighs/atom = 0
Neighbor list builds = 258
Dangerous builds = 0
fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out
thermo 2000
run 35000
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 6.481 | 6.481 | 6.481 Mbytes
Step Temp E_pair E_mol TotEng Press
5000 9.7731898 -55.150518 0 -54.060304 322.94195
6000 10.024842 -55.108242 0 -53.989956 336.6125
8000 10.118994 -55.092171 0 -53.963382 340.42078
10000 10.541359 -55.100316 0 -53.924412 340.09986
12000 10.023234 -55.07343 0 -53.955323 345.70551
14000 9.5912018 -55.115121 0 -54.045208 333.43739
16000 9.9450498 -55.077813 0 -53.968428 343.49906
18000 10.113744 -55.12201 0 -53.993806 332.32214
20000 9.9345204 -55.176339 0 -54.068128 316.83219
22000 10.585719 -55.110505 0 -53.929652 336.86599
24000 10.024757 -55.120592 0 -54.002315 333.13056
26000 9.7787474 -55.116664 0 -54.02583 332.51437
28000 9.6092087 -55.130413 0 -54.058491 328.69165
30000 9.8245787 -55.057678 0 -53.961731 350.86255
32000 10.066994 -55.08989 0 -53.966902 341.49724
34000 9.5677059 -55.17222 0 -54.104928 316.06902
36000 9.7252627 -55.079475 0 -53.994608 343.13769
38000 10.438984 -55.108991 0 -53.944506 336.32562
40000 10.238268 -55.080817 0 -53.938723 345.13228
Loop time of 19.138 on 1 procs for 35000 steps with 864 atoms
Performance: 15801.041 ns/day, 0.002 hours/ns, 1828.824 timesteps/s, 1.580 Matom-step/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 11.568 | 11.568 | 11.568 | 0.0 | 60.44
Bond | 0.0042372 | 0.0042372 | 0.0042372 | 0.0 | 0.02
Neigh | 2.2577 | 2.2577 | 2.2577 | 0.0 | 11.80
Comm | 0.42841 | 0.42841 | 0.42841 | 0.0 | 2.24
Output | 0.00060128 | 0.00060128 | 0.00060128 | 0.0 | 0.00
Modify | 4.694 | 4.694 | 4.694 | 0.0 | 24.53
Other | | 0.1852 | | | 0.97
Nlocal: 864 ave 864 max 864 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1592 ave 1592 max 1592 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 18144 ave 18144 max 18144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 18144
Ave neighs/atom = 21
Ave special neighs/atom = 0
Neighbor list builds = 2122
Dangerous builds = 0
Total wall time: 0:00:21

View File

@ -0,0 +1,192 @@
LAMMPS (2 Apr 2025 - Development - d4867ab55e-modified)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:99)
using 1 OpenMP thread(s) per MPI task
# GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
Reading data file ...
orthogonal box = (0 0 0) to (32.146 32.146 32.146)
1 by 2 by 2 MPI processor grid
reading atoms ...
864 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.002 seconds
read_data CPU = 0.015 seconds
include ff-argon.lmp
#############################
#Atoms types - mass - charge#
#############################
#@ 1 atom types #!THIS LINE IS NECESSARY DON'T SPEND HOURS FINDING THAT OUT!#
variable Ar equal 1
#############
#Atom Masses#
#############
mass ${Ar} 39.903
mass 1 39.903
###########################
#Pair Potentials - Tersoff#
###########################
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
compute myKE all ke
compute myPE all pe
fix lang all gjf 10 10 1 26488
run 5000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- GJ methods: doi:10.1080/00268976.2019.1662506
@Article{gronbech-jensen_complete_2020,
title = {Complete set of stochastic Verlet-type thermostats for correct Langevin simulations},
volume = {118},
number = {8},
url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1662506},
doi = {10.1080/00268976.2019.1662506},
journal = {Molecular Physics},
author = {Grønbech-Jensen, Niels},
year = {2020}
}
- GJ-I vhalf method: doi:10.1080/00268976.2019.1570369
@Article{jensen_accurate_2019,
title = {Accurate configurational and kinetic statistics in discrete-time Langevin systems},
volume = {117},
url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1570369},
doi = {10.1080/00268976.2019.1570369},
number = {18},
journal = {Molecular Physics},
author = {Jensen, Lucas Frese Grønbech and Grønbech-Jensen, Niels},
year = {2019}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.9407173
ghost atom cutoff = 6.9407173
binsize = 3.4703587, bins = 10 10 10
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cubic, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.427 | 6.427 | 6.427 Mbytes
Step Temp E_pair E_mol TotEng Press
0 10 -56.207652 0 -55.092137 33.341103
5000 9.3726166 -55.076514 0 -54.030985 342.43571
Loop time of 2.11818 on 4 procs for 5000 steps with 864 atoms
Performance: 20394.822 ns/day, 0.001 hours/ns, 2360.512 timesteps/s, 2.039 Matom-step/s
63.1% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.53987 | 0.54922 | 0.56044 | 1.2 | 25.93
Bond | 0.00058281 | 0.00063674 | 0.00075153 | 0.0 | 0.03
Neigh | 0.10821 | 0.10912 | 0.11017 | 0.2 | 5.15
Comm | 0.96075 | 0.97484 | 0.98645 | 1.1 | 46.02
Output | 0.00026318 | 0.00026575 | 0.00027192 | 0.0 | 0.01
Modify | 0.26142 | 0.2634 | 0.26465 | 0.2 | 12.44
Other | | 0.2207 | | | 10.42
Nlocal: 216 ave 216 max 216 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 884.75 ave 885 max 884 min
Histogram: 1 0 0 0 0 0 0 0 0 3
Neighs: 4536 ave 4737 max 4335 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 18144
Ave neighs/atom = 21
Ave special neighs/atom = 0
Neighbor list builds = 273
Dangerous builds = 0
fix energies all ave/time 1 20000 20000 c_myKE c_myPE #file ave.out
thermo 2000
run 35000
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Per MPI rank memory allocation (min/avg/max) = 6.428 | 6.428 | 6.428 Mbytes
Step Temp E_pair E_mol TotEng Press
5000 9.3726166 -55.076514 0 -54.030985 342.43571
6000 9.6911866 -55.129244 0 -54.048177 329.72537
8000 9.7296551 -55.101417 0 -54.016059 337.46595
10000 10.098808 -55.083292 0 -53.956755 343.4122
12000 10.114344 -55.090905 0 -53.962635 341.3438
14000 10.230012 -55.085051 0 -53.943878 342.45237
16000 9.5989709 -55.090733 0 -54.019954 339.07584
18000 10.016071 -55.030131 0 -53.912824 358.79514
20000 9.7197057 -55.104072 0 -54.019824 335.89619
22000 9.959647 -55.152019 0 -54.041005 323.05805
24000 10.075138 -55.115242 0 -53.991345 334.76239
26000 10.227192 -55.125275 0 -53.984416 332.10131
28000 10.177109 -55.107097 0 -53.971825 337.32979
30000 9.521036 -55.099011 0 -54.036925 337.10716
32000 10.265633 -55.077016 0 -53.93187 346.01018
34000 10.173978 -55.097068 0 -53.962146 339.63562
36000 9.6032778 -55.152202 0 -54.080942 321.61646
38000 9.8802995 -55.048407 0 -53.946245 351.82506
40000 10.372288 -55.087634 0 -53.93059 343.34304
Loop time of 17.867 on 4 procs for 35000 steps with 864 atoms
Performance: 16925.013 ns/day, 0.001 hours/ns, 1958.914 timesteps/s, 1.693 Matom-step/s
65.3% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.0932 | 5.1683 | 5.2256 | 2.5 | 28.93
Bond | 0.0044473 | 0.0048347 | 0.0058137 | 0.8 | 0.03
Neigh | 0.85262 | 0.8601 | 0.87438 | 0.9 | 4.81
Comm | 6.8164 | 6.8981 | 6.9859 | 2.6 | 38.61
Output | 0.0046884 | 0.0047093 | 0.0047322 | 0.0 | 0.03
Modify | 3.4107 | 3.4186 | 3.4248 | 0.3 | 19.13
Other | | 1.512 | | | 8.47
Nlocal: 216 ave 222 max 210 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 905.5 ave 911 max 899 min
Histogram: 1 1 0 0 0 0 0 0 0 2
Neighs: 4535.75 ave 4837 max 4218 min
Histogram: 1 0 0 1 0 0 1 0 0 1
Total # of neighbors = 18143
Ave neighs/atom = 20.998843
Ave special neighs/atom = 0
Neighbor list builds = 2140
Dangerous builds = 0
Total wall time: 0:00:21

2
src/.gitignore vendored
View File

@ -875,6 +875,8 @@
/fix_freeze.h
/fix_gcmc.cpp
/fix_gcmc.h
/fix_gjf.cpp
/fix_gjf.h
/fix_gld.cpp
/fix_gld.h
/fix_gle.cpp

687
src/EXTRA-FIX/fix_gjf.cpp Normal file
View File

@ -0,0 +1,687 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Tim Linke & Niels Gronbech-Jensen (UC Davis)
------------------------------------------------------------------------- */
#include "fix_gjf.h"
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "compute.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "input.h"
#include "memory.h"
#include "modify.h"
#include "random_mars.h"
#include "respa.h"
#include "update.h"
#include "variable.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
enum { NOBIAS, BIAS };
enum { CONSTANT, EQUAL, ATOM };
static const char cite_gjf[] =
"GJ methods: doi:10.1080/00268976.2019.1662506\n\n"
"@Article{gronbech-jensen_complete_2020,\n"
"title = {Complete set of stochastic Verlet-type thermostats for correct Langevin simulations},\n"
"volume = {118},\n"
"number = {8},\n"
"url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1662506},\n"
"doi = {10.1080/00268976.2019.1662506},\n"
"journal = {Molecular Physics},\n"
"author = {Grønbech-Jensen, Niels},\n"
"year = {2020}\n"
"}\n\n";
static const char cite_gjf_7[] =
"GJ-VII method: doi:10.1063/5.0066008\n\n"
"@Article{finkelstein_2021,\n"
"title = {Bringing discrete-time Langevin splitting methods into agreement with thermodynamics},\n"
"volume = {155},\n"
"number = {18},\n"
"url = {https://doi.org/10.1063/5.0066008},\n"
"doi = {10.1063/5.0066008},\n"
"journal = {J. Chem. Phys.},\n"
"author = {Finkelstein, Joshua and Cheng, Chungho and Fiorin, Giacomo and Seibold, Benjamin and Grønbech-Jensen, Niels},\n"
"year = {2021},\n"
"pages = {184104}\n"
"}\n\n";
static const char cite_gjf_8[] =
"GJ-VIII method: doi:10.1007/s10955-024-03345-1\n\n"
"@Article{gronbech_jensen_2024,\n"
"title = {On the Definition of Velocity in Discrete-Time, Stochastic Langevin Simulations},\n"
"volume = {191},\n"
"number = {10},\n"
"url = {https://doi.org/10.1007/s10955-024-03345-1},\n"
"doi = {10.1007/s10955-024-03345-1},\n"
"journal = {J. Stat. Phys.},\n"
"author = {Gronbech-Jensen, Niels},\n"
"year = {2024},\n"
"pages = {137}\n"
"}\n\n";
static const char cite_gjf_vhalf[] =
"GJ-I vhalf method: doi:10.1080/00268976.2019.1570369\n\n"
"@Article{jensen_accurate_2019,\n"
"title = {Accurate configurational and kinetic statistics in discrete-time Langevin systems},\n"
"volume = {117},\n"
"url = {https://www.tandfonline.com/doi/full/10.1080/00268976.2019.1570369},\n"
"doi = {10.1080/00268976.2019.1570369},\n"
"number = {18},\n"
"journal = {Molecular Physics},\n"
"author = {Jensen, Lucas Frese Grønbech and Grønbech-Jensen, Niels},\n"
"year = {2019}\n"
"}\n\n";
static const char cite_gjf_vfull[] =
"GJ-I vfull method: doi:10.1080/00268976.2012.760055\n\n"
"@Article{gronbech-jensen_simple_2013,\n"
"title = {A simple and effective Verlet-type algorithm for simulating Langevin dynamics},\n"
"volume = {111},\n"
"url = {http://www.tandfonline.com/doi/abs/10.1080/00268976.2012.760055},\n"
"doi = {10.1080/00268976.2012.760055},\n"
"pages = {983-991},\n"
"number = {8},\n"
"journal = {Molecular Physics},\n"
"author = {Grønbech-Jensen, Niels and Farago, Oded},\n"
"year = {2013}\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
FixGJF::FixGJF(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), tstr(nullptr), tforce(nullptr), lv(nullptr), id_temp(nullptr), random(nullptr)
{
if (lmp->citeme) lmp->citeme->add(cite_gjf);
if (narg < 7) error->all(FLERR, "Illegal fix gjf command");
time_integrate = 1;
global_freq = 1;
nevery = 1;
if (utils::strmatch(arg[3], "^v_")) {
tstr = utils::strdup(arg[3] + 2);
} else {
t_start = utils::numeric(FLERR, arg[3], false, lmp);
t_target = t_start;
tstyle = CONSTANT;
}
t_stop = utils::numeric(FLERR, arg[4], false, lmp);
t_period = utils::numeric(FLERR, arg[5], false, lmp);
seed = utils::inumeric(FLERR, arg[6], false, lmp);
if (t_period <= 0.0) error->all(FLERR, "Fix gjf period must be > 0.0");
if (seed <= 0) error->all(FLERR, "Illegal fix gjf command");
// initialize Marsaglia RNG with processor-unique seed
random = new RanMars(lmp, seed + comm->me);
int GJmethods = 8; // number of currently implemented GJ methods
maxatom = 0;
// optional args
// per default, use half step and GJ-I
osflag = 0;
GJmethod = 1;
lv_allocated = 0;
int iarg = 7;
while (iarg < narg) {
if (strcmp(arg[iarg], "vel") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix gjf command");
if (strcmp(arg[iarg + 1], "vfull") == 0) {
osflag = 1;
} else if (strcmp(arg[iarg + 1], "vhalf") == 0) {
osflag = 0;
} else
error->all(FLERR, "Illegal fix gjf command");
iarg += 2;
} else if (strcmp(arg[iarg], "method") == 0) {
GJmethod = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
if (GJmethod == 7) {
if (iarg + 3 > narg) error->all(FLERR, "Illegal fix gjf command for GJ-VII");
gjfc2 = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
if (gjfc2 < 0 || gjfc2 > 1) error->all(FLERR, "Choice of c2 in GJ-VII must be 0≤c2≤1");
iarg += 3;
if (lmp->citeme) lmp->citeme->add(cite_gjf_7);
}
else {
if (iarg + 2 > narg) error->all(FLERR, "Illegal fix gjf command");
if (GJmethod < 0 || GJmethod > GJmethods) error->all(FLERR, "Invalid GJ method choice in gjf command");
if (GJmethod == 8) if (lmp->citeme) lmp->citeme->add(cite_gjf_8);
iarg += 2;
}
} else
error->all(FLERR, "Illegal fix gjf command");
}
if (GJmethod == 1 && osflag == 0) if (lmp->citeme) lmp->citeme->add(cite_gjf_vhalf);
if (GJmethod == 1 && osflag == 1) if (lmp->citeme) lmp->citeme->add(cite_gjf_vfull);
// set temperature = nullptr, user can override via fix_modify if wants bias
id_temp = nullptr;
temperature = nullptr;
lv = nullptr;
tforce = nullptr;
// setup atom-based array for lv
// register with Atom class
// no need to set peratom_flag, b/c data is for internal use only
FixGJF::grow_arrays(atom->nmax);
atom->add_callback(Atom::GROW);
// initialize lv to onsite velocity
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
lv[i][0] = 0.0;
lv[i][1] = 0.0;
lv[i][2] = 0.0;
}
}
/* ---------------------------------------------------------------------- */
FixGJF::~FixGJF()
{
if (copymode) return;
delete random;
delete[] tstr;
delete[] id_temp;
memory->destroy(tforce);
memory->destroy(lv);
if (modify->get_fix_by_id(id)) atom->delete_callback(id, Atom::GROW);
}
/* ---------------------------------------------------------------------- */
int FixGJF::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
if (!osflag) mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixGJF::init()
{
if (id_temp) {
temperature = modify->get_compute_by_id(id_temp);
if (!temperature) {
error->all(FLERR, "Temperature compute ID {} for fix {} does not exist", id_temp, style);
} else {
if (temperature->tempflag == 0)
error->all(FLERR, "Compute ID {} for fix {} does not compute temperature", id_temp, style);
}
}
// check variable
if (tstr) {
tvar = input->variable->find(tstr);
if (tvar < 0) error->all(FLERR, "Variable name {} for fix gjf does not exist", tstr);
if (input->variable->equalstyle(tvar))
tstyle = EQUAL;
else if (input->variable->atomstyle(tvar))
tstyle = ATOM;
else
error->all(FLERR, "Variable {} for fix gjf is invalid style", tstr);
}
if (utils::strmatch(update->integrate_style, "^respa")) {
error->all(FLERR, "Fix gjf and run style respa are not compatible");
}
if (temperature && temperature->tempbias)
tbiasflag = BIAS;
else
tbiasflag = NOBIAS;
// Complete set of thermostats is given in Gronbech-Jensen, Molecular Physics, 118 (2020)
switch (GJmethod) {
case 1:
gjfc2 = (1.0 - update->dt / 2.0 / t_period) / (1.0 + update->dt / 2.0 / t_period);
break;
case 2:
gjfc2 = exp(-update->dt / t_period);
break;
case 3:
gjfc2 = 1.0 - update->dt / t_period;
break;
case 4:
gjfc2 = ( sqrt(1.0 + 4.0 * (update->dt / t_period) ) - 1.0 ) / ( 2.0 * update->dt / t_period );
break;
case 5:
gjfc2 = 1.0 / (1.0 + update->dt / t_period);
break;
case 6:
gjfc2 = (1.0 / (1.0 + update->dt / 2.0 / t_period)) * (1.0 / (1.0 + update->dt / 2.0 / t_period));
break;
case 7: // provided in Finkelstein (2021)
update->dt = (1.0 + gjfc2) / (1.0 - gjfc2) * log(gjfc2) * log(gjfc2) * 0.5 * t_period;
break;
case 8: // provided in Gronbech-Jensen (2024)
gjfc2 = sqrt( (update->dt / t_period) * (update->dt / t_period) + 1.0 ) - update->dt / t_period;
break;
case 0:
gjfc2 = 0.0;
break;
default:
error->all(FLERR, "Fix gjf method not found");
break;
}
gjfc1 = (1.0 + gjfc2) / 2.0;
gjfc3 = (1.0 - gjfc2) * t_period / update->dt;
}
/* ----------------------------------------------------------------------
integrate position and velocity according to the GJ methods
in Grønbech-Jensen, J Stat Phys 191, 137 (2024). The general workflow is
1. GJ Initial Integration
2. Force Update
3. GJ Final Integration
4. Velocity Choice in end_of_step()
------------------------------------------------------------------------- */
void FixGJF::initial_integrate(int /* vflag */)
{
// This function provides the integration of the GJ formulation 24 a-e
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double fran[3];
double boltz = force->boltz;
double dt = update->dt;
double mvv2e = force->mvv2e;
double ftm2v = force->ftm2v;
double dtf = 0.5 * dt * ftm2v;
double dtfm;
double c1sqrt = sqrt(gjfc1);
double c3sqrt = sqrt(gjfc3);
double csq = sqrt(gjfc3 / gjfc1);
double m, beta;
// If user elected vhalf, v needs to be reassigned to onsite velocity for integration
if (!osflag && lv_allocated) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
// lv is Eq. 24f from previous time step
v[i][0] = lv[i][0];
v[i][1] = lv[i][1];
v[i][2] = lv[i][2];
}
}
compute_target();
if (tbiasflag) temperature->compute_scalar();
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
m = rmass[i];
beta = tsqrt * sqrt(2.0*dt*m*boltz/t_period/mvv2e) / ftm2v;
fran[0] = beta*random->gaussian();
fran[1] = beta*random->gaussian();
fran[2] = beta*random->gaussian();
// First integration delivers Eq. 24a and 24b:
dtfm = dtf / m;
v[i][0] += csq * dtfm * f[i][0];
v[i][1] += csq * dtfm * f[i][1];
v[i][2] += csq * dtfm * f[i][2];
x[i][0] += 0.5 * csq * dt * v[i][0];
x[i][1] += 0.5 * csq * dt * v[i][1];
x[i][2] += 0.5 * csq * dt * v[i][2];
if (tbiasflag) temperature->remove_bias(i, v[i]);
// Calculate Eq. 24c:
lv[i][0] = c1sqrt*v[i][0] + ftm2v * (c3sqrt / (2.0 * m)) * fran[0];
lv[i][1] = c1sqrt*v[i][1] + ftm2v * (c3sqrt / (2.0 * m)) * fran[1];
lv[i][2] = c1sqrt*v[i][2] + ftm2v * (c3sqrt / (2.0 * m)) * fran[2];
// Calculate Eq. 24d
v[i][0] = (gjfc2 / c1sqrt) * lv[i][0] + ftm2v * csq * (0.5 / m) * fran[0];
v[i][1] = (gjfc2 / c1sqrt) * lv[i][1] + ftm2v * csq * (0.5 / m) * fran[1];
v[i][2] = (gjfc2 / c1sqrt) * lv[i][2] + ftm2v * csq * (0.5 / m) * fran[2];
if (tbiasflag) temperature->restore_bias(i, v[i]);
if (tbiasflag) temperature->restore_bias(i, lv[i]);
// Calculate Eq. 24e. Final integrator then calculates Eq. 24f after force update.
x[i][0] += 0.5 * csq * dt * v[i][0];
x[i][1] += 0.5 * csq * dt * v[i][1];
x[i][2] += 0.5 * csq * dt * v[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
m = mass[type[i]];
beta = tsqrt * sqrt(2.0*dt*m*boltz/t_period/mvv2e) / ftm2v;
fran[0] = beta*random->gaussian();
fran[1] = beta*random->gaussian();
fran[2] = beta*random->gaussian();
// First integration delivers Eq. 24a and 24b:
dtfm = dtf / m;
v[i][0] += csq * dtfm * f[i][0];
v[i][1] += csq * dtfm * f[i][1];
v[i][2] += csq * dtfm * f[i][2];
x[i][0] += 0.5 * csq * dt * v[i][0];
x[i][1] += 0.5 * csq * dt * v[i][1];
x[i][2] += 0.5 * csq * dt * v[i][2];
if (tbiasflag) temperature->remove_bias(i, v[i]);
// Calculate Eq. 24c:
lv[i][0] = c1sqrt*v[i][0] + ftm2v * (c3sqrt / (2.0 * m)) * fran[0];
lv[i][1] = c1sqrt*v[i][1] + ftm2v * (c3sqrt / (2.0 * m)) * fran[1];
lv[i][2] = c1sqrt*v[i][2] + ftm2v * (c3sqrt / (2.0 * m)) * fran[2];
// Calculate Eq. 24d
v[i][0] = (gjfc2 / c1sqrt) * lv[i][0] + ftm2v * csq * (0.5 / m) * fran[0];
v[i][1] = (gjfc2 / c1sqrt) * lv[i][1] + ftm2v * csq * (0.5 / m) * fran[1];
v[i][2] = (gjfc2 / c1sqrt) * lv[i][2] + ftm2v * csq * (0.5 / m) * fran[2];
if (tbiasflag) temperature->restore_bias(i, v[i]);
if (tbiasflag) temperature->restore_bias(i, lv[i]);
// Calculate Eq. 24e. Final integrator then calculates Eq. 24f after force update.
x[i][0] += 0.5 * csq * dt * v[i][0];
x[i][1] += 0.5 * csq * dt * v[i][1];
x[i][2] += 0.5 * csq * dt * v[i][2];
}
}
}
}
void FixGJF::final_integrate()
{
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
double dtfm;
double dtf = 0.5 * update->dt * force->ftm2v;
double csq = sqrt(gjfc3 / gjfc1);
// Calculate Eq. 24f.
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += csq * dtfm * f[i][0];
v[i][1] += csq * dtfm * f[i][1];
v[i][2] += csq * dtfm * f[i][2];
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += csq * dtfm * f[i][0];
v[i][1] += csq * dtfm * f[i][1];
v[i][2] += csq * dtfm * f[i][2];
}
}
lv_allocated = 1;
}
/* ----------------------------------------------------------------------
set current t_target and t_sqrt
------------------------------------------------------------------------- */
void FixGJF::compute_target()
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
// if variable temp, evaluate variable, wrap with clear/add
// reallocate tforce array if necessary
if (tstyle == CONSTANT) {
t_target = t_start + delta * (t_stop-t_start);
tsqrt = sqrt(t_target);
} else {
modify->clearstep_compute();
if (tstyle == EQUAL) {
t_target = input->variable->compute_equal(tvar);
if (t_target < 0.0)
error->one(FLERR, "Fix gjf variable returned negative temperature");
tsqrt = sqrt(t_target);
} else {
if (atom->nmax > maxatom) {
maxatom = atom->nmax;
memory->destroy(tforce);
memory->create(tforce,maxatom,"gjf:tforce");
}
input->variable->compute_atom(tvar,igroup,tforce,1,0);
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (tforce[i] < 0.0)
error->one(FLERR, "Fix gjf variable returned negative temperature");
}
modify->addstep_compute(update->ntimestep + 1);
}
}
/* ----------------------------------------------------------------------
select velocity for GJ
------------------------------------------------------------------------- */
void FixGJF::end_of_step()
{
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// After the final integrator delivers 24f, either the on-site or half-step
// velocity is used in remaining simulation tasks, depending on user input
double tmp[3];
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
// v is Eq. 24f
tmp[0] = v[i][0];
tmp[1] = v[i][1];
tmp[2] = v[i][2];
// Move on with half-step velocity
v[i][0] = lv[i][0];
v[i][1] = lv[i][1];
v[i][2] = lv[i][2];
// store Eq. 24f in lv for next timestep
lv[i][0] = tmp[0];
lv[i][1] = tmp[1];
lv[i][2] = tmp[2];
}
}
// clang-format on
/* ---------------------------------------------------------------------- */
void FixGJF::reset_target(double t_new)
{
t_target = t_start = t_stop = t_new;
}
/* ---------------------------------------------------------------------- */
void FixGJF::reset_dt()
{
// Complete set of thermostats is given in Gronbech-Jensen, Molecular Physics, 118 (2020)
switch (GJmethod) {
case 1:
gjfc2 = (1.0 - update->dt / 2.0 / t_period) / (1.0 + update->dt / 2.0 / t_period);
break;
case 2:
gjfc2 = exp(-update->dt / t_period);
break;
case 3:
gjfc2 = 1.0 - update->dt / t_period;
break;
case 4:
gjfc2 = ( sqrt(1.0 + 4.0 * (update->dt / t_period) ) - 1.0 ) / ( 2.0 * update->dt / t_period );
break;
case 5:
gjfc2 = 1.0 / (1.0 + update->dt / t_period);
break;
case 6:
gjfc2 = (1.0 / (1.0 + update->dt / 2.0 / t_period)) * (1.0 / (1.0 + update->dt / 2.0 / t_period));
break;
case 7: // provided in Finkelstein (2021)
update->dt = (1.0 + gjfc2) / (1.0 - gjfc2) * log(gjfc2) * log(gjfc2) * 0.5 * t_period;
break;
case 8: // provided in Gronbech-Jensen (2024)
gjfc2 = sqrt( (update->dt / t_period)*(update->dt / t_period) + 1.0 ) - update->dt / t_period;
break;
case 0:
gjfc2 = 0.0;
break;
default:
error->all(FLERR, "Fix gjf method not found");
break;
}
gjfc1 = (1.0 + gjfc2) / 2.0;
gjfc3 = (1.0 - gjfc2) * t_period / update->dt;
}
/* ---------------------------------------------------------------------- */
int FixGJF::modify_param(int narg, char **arg)
{
if (strcmp(arg[0], "temp") == 0) {
if (narg < 2) utils::missing_cmd_args(FLERR, "fix_modify", error);
delete[] id_temp;
id_temp = utils::strdup(arg[1]);
temperature = modify->get_compute_by_id(id_temp);
if (!temperature)
error->all(FLERR, "Could not find fix_modify temperature compute ID: {}", id_temp);
if (temperature->tempflag == 0)
error->all(FLERR, "Fix_modify temperature compute {} does not compute temperature", id_temp);
if (temperature->igroup != igroup && comm->me == 0)
error->warning(FLERR, "Group for fix_modify temp != fix group: {} vs {}",
group->names[igroup], group->names[temperature->igroup]);
return 2;
}
return 0;
}
/* ----------------------------------------------------------------------
extract thermostat properties
------------------------------------------------------------------------- */
void *FixGJF::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str, "t_target") == 0) { return &t_target; }
return nullptr;
}
/* ----------------------------------------------------------------------
memory usage of tally array
------------------------------------------------------------------------- */
double FixGJF::memory_usage()
{
double bytes = 0.0;
bytes += (double) atom->nmax * 3 * sizeof(double);
if (tforce) bytes += (double) atom->nmax * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
allocate atom-based array for lv
------------------------------------------------------------------------- */
void FixGJF::grow_arrays(int nmax)
{
memory->grow(lv, nmax, 3, "fix_gjf:lv");
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
------------------------------------------------------------------------- */
void FixGJF::copy_arrays(int i, int j, int /*delflag*/)
{
lv[j][0] = lv[i][0];
lv[j][1] = lv[i][1];
lv[j][2] = lv[i][2];
}
/* ----------------------------------------------------------------------
pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */
int FixGJF::pack_exchange(int i, double *buf)
{
int n = 0;
buf[n++] = lv[i][0];
buf[n++] = lv[i][1];
buf[n++] = lv[i][2];
return n;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */
int FixGJF::unpack_exchange(int nlocal, double *buf)
{
int n = 0;
lv[nlocal][0] = buf[n++];
lv[nlocal][1] = buf[n++];
lv[nlocal][2] = buf[n++];
return n;
}

68
src/EXTRA-FIX/fix_gjf.h Normal file
View File

@ -0,0 +1,68 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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(gjf,FixGJF);
// clang-format on
#else
#ifndef LMP_FIX_GJF_H
#define LMP_FIX_GJF_H
#include "fix.h"
namespace LAMMPS_NS {
class FixGJF : public Fix {
public:
FixGJF(class LAMMPS *, int, char **);
~FixGJF() override;
int setmask() override;
void init() override;
void initial_integrate(int) override;
void final_integrate() override;
void end_of_step() override;
void reset_target(double) override;
void reset_dt() override;
int modify_param(int, char **) override;
double memory_usage() override;
void *extract(const char *, int &) override;
void grow_arrays(int) override;
void copy_arrays(int, int, int) override;
int pack_exchange(int, double *) override;
int unpack_exchange(int, double *) override;
protected:
int osflag, tbiasflag, GJmethod, maxatom, lv_allocated;
double t_start, t_stop, t_period, t_target, tsqrt;
double gjfc1, gjfc2, gjfc3;
int tstyle, tvar;
char *tstr;
double *tforce;
double **lv; //half step velocity
char *id_temp;
class Compute *temperature;
class RanMars *random;
int seed;
void compute_target();
};
} // namespace LAMMPS_NS
#endif
#endif

View File

@ -170,6 +170,9 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
// no need to set peratom_flag, b/c data is for internal use only
if (gjfflag) {
if (comm->me == 0)
error->warning(FLERR, "The GJF formulation in fix {} is deprecated and will be removed soon. "
"\nPlease use fix gjf instead: https://docs.lammps.org/fix_gjf.html", style);
FixLangevin::grow_arrays(atom->nmax);
atom->add_callback(Atom::GROW);