diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 35c3804969..e80f01b87a 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -77,6 +77,7 @@ OPT. * :doc:`flow/gauss ` * :doc:`freeze (k) ` * :doc:`gcmc ` + * :doc:`gjf ` * :doc:`gld ` * :doc:`gle ` * :doc:`gravity (ko) ` diff --git a/doc/src/Commands_removed.rst b/doc/src/Commands_removed.rst index 847fe76456..c92d3d9064 100644 --- a/doc/src/Commands_removed.rst +++ b/doc/src/Commands_removed.rst @@ -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 ` and it is strongly recommended to use that +fix instead. + + LAMMPS shell ------------ diff --git a/doc/src/Howto_thermostat.rst b/doc/src/Howto_thermostat.rst index fe53fff540..bda3bc6cb4 100644 --- a/doc/src/Howto_thermostat.rst +++ b/doc/src/Howto_thermostat.rst @@ -21,9 +21,14 @@ can be invoked via the *dpd/tstat* pair style: * :doc:`fix nvt/sllod ` * :doc:`fix temp/berendsen ` * :doc:`fix temp/csvr ` +* :doc:`fix ffl ` +* :doc:`fix gjf ` +* :doc:`fix gld ` +* :doc:`fix gle ` * :doc:`fix langevin ` * :doc:`fix temp/rescale ` * :doc:`pair_style dpd/tstat ` +* :doc:`pair_style dpd/ext/tstat ` :doc:`Fix nvt ` only thermostats the translational velocity of particles. :doc:`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 ` diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 3b1bc4a75c..cdc0e786ce 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -256,6 +256,7 @@ accelerated styles exist. * :doc:`flow/gauss ` - Gaussian dynamics for constant mass flux * :doc:`freeze ` - freeze atoms in a granular simulation * :doc:`gcmc ` - grand canonical insertions/deletions +* :doc:`gjf ` - statistically correct Langevin temperature control using the GJ methods * :doc:`gld ` - generalized Langevin dynamics integrator * :doc:`gle ` - generalized Langevin equation thermostat * :doc:`gravity ` - add gravity to atoms in a granular simulation diff --git a/doc/src/fix_gjf.rst b/doc/src/fix_gjf.rst new file mode 100644 index 0000000000..f55a811ab7 --- /dev/null +++ b/doc/src/fix_gjf.rst @@ -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 ` 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) ` +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 +`; 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 +`; 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 `. + +This fix allows the use of several GJ methods as listed in :ref:`Gronbech-Jensen-2020 `. +The GJ-VII method is described in :ref:`Finkelstein ` and GJ-VIII +is described in :ref:`Gronbech-Jensen-2024 `. +The implementation follows the splitting form provided in Eqs. (24) and (25) +in :ref:`Gronbech-Jensen-2024 `, including the application +of Gaussian noise values, per the description in +:ref:`Gronbech-Jensen-2023 `. + + +.. note:: + + Unlike the :doc:`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 `. + +See the :doc:`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 `. 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 ` 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 ` that remove a "bias" from the +atom velocities. E.g. to apply the thermostat only to atoms within a +spatial :doc:`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 +` 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 ` 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 ` 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 ` +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 `. +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 ` *temp* option is supported by this +fix. You can use it to assign a temperature :doc:`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 ` command. See the +:doc:`run ` command for details of how to do this. + +This fix is not invoked during :doc:`energy minimization `. + +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 `, :doc:`fix nvt ` + +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). diff --git a/doc/src/fix_langevin.rst b/doc/src/fix_langevin.rst index 30e4c48270..b8be16e16d 100644 --- a/doc/src/fix_langevin.rst +++ b/doc/src/fix_langevin.rst @@ -56,7 +56,7 @@ Examples Description """"""""""" -Apply a Langevin thermostat as described in :ref:`(Schneider) ` +Apply a Langevin thermostat as described in :ref:`(Bruenger) ` to a group of atoms which models an interaction with a background implicit solvent. Used with :doc:`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 ` and it is strongly recommended to use that +fix instead. + The keyword *gjf* can be used to run the :ref:`Gronbech-Jensen/Farago ` 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 diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 4bd99e4c80..c81be91cbe 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -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 diff --git a/examples/gjf/README b/examples/gjf/README new file mode 100644 index 0000000000..f7c955080f --- /dev/null +++ b/examples/gjf/README @@ -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 +-- \ No newline at end of file diff --git a/examples/gjf/README.md b/examples/gjf/README.md deleted file mode 100644 index e6886cb2dd..0000000000 --- a/examples/gjf/README.md +++ /dev/null @@ -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 diff --git a/examples/gjf/in.gjf.vfull b/examples/gjf/in.gjf.vfull index 40512ac37a..fad6df4a9b 100644 --- a/examples/gjf/in.gjf.vfull +++ b/examples/gjf/in.gjf.vfull @@ -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 \ No newline at end of file diff --git a/examples/gjf/in.gjf.vhalf b/examples/gjf/in.gjf.vhalf index 63fb8bd467..2f5c482928 100644 --- a/examples/gjf/in.gjf.vhalf +++ b/examples/gjf/in.gjf.vhalf @@ -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 \ No newline at end of file diff --git a/examples/gjf/log.15Oct19.gjf.vfull.g++.1 b/examples/gjf/log.15Oct19.gjf.vfull.g++.1 deleted file mode 100644 index e3e9cce124..0000000000 --- a/examples/gjf/log.15Oct19.gjf.vfull.g++.1 +++ /dev/null @@ -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 diff --git a/examples/gjf/log.15Oct19.gjf.vfull.g++.4 b/examples/gjf/log.15Oct19.gjf.vfull.g++.4 deleted file mode 100644 index 95caed5dc9..0000000000 --- a/examples/gjf/log.15Oct19.gjf.vfull.g++.4 +++ /dev/null @@ -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 diff --git a/examples/gjf/log.15Oct19.gjf.vhalf.g++.1 b/examples/gjf/log.15Oct19.gjf.vhalf.g++.1 deleted file mode 100644 index a87b20a887..0000000000 --- a/examples/gjf/log.15Oct19.gjf.vhalf.g++.1 +++ /dev/null @@ -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 diff --git a/examples/gjf/log.15Oct19.gjf.vhalf.g++.4 b/examples/gjf/log.15Oct19.gjf.vhalf.g++.4 deleted file mode 100644 index a70a67a89c..0000000000 --- a/examples/gjf/log.15Oct19.gjf.vhalf.g++.4 +++ /dev/null @@ -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 diff --git a/examples/gjf/log.2Apr25.gjf.vfull.g++.1 b/examples/gjf/log.2Apr25.gjf.vfull.g++.1 new file mode 100644 index 0000000000..ca739b3a22 --- /dev/null +++ b/examples/gjf/log.2Apr25.gjf.vfull.g++.1 @@ -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 diff --git a/examples/gjf/log.2Apr25.gjf.vfull.g++.4 b/examples/gjf/log.2Apr25.gjf.vfull.g++.4 new file mode 100644 index 0000000000..93466e8dfe --- /dev/null +++ b/examples/gjf/log.2Apr25.gjf.vfull.g++.4 @@ -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 diff --git a/examples/gjf/log.2Apr25.gjf.vhalf.g++.1 b/examples/gjf/log.2Apr25.gjf.vhalf.g++.1 new file mode 100644 index 0000000000..569b13bad6 --- /dev/null +++ b/examples/gjf/log.2Apr25.gjf.vhalf.g++.1 @@ -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 diff --git a/examples/gjf/log.2Apr25.gjf.vhalf.g++.4 b/examples/gjf/log.2Apr25.gjf.vhalf.g++.4 new file mode 100644 index 0000000000..2b5e19e634 --- /dev/null +++ b/examples/gjf/log.2Apr25.gjf.vhalf.g++.4 @@ -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 diff --git a/src/.gitignore b/src/.gitignore index a25a884a8a..20d1bccc4a 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -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 diff --git a/src/EXTRA-FIX/fix_gjf.cpp b/src/EXTRA-FIX/fix_gjf.cpp new file mode 100644 index 0000000000..f602e6181d --- /dev/null +++ b/src/EXTRA-FIX/fix_gjf.cpp @@ -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 +#include + +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; +} diff --git a/src/EXTRA-FIX/fix_gjf.h b/src/EXTRA-FIX/fix_gjf.h new file mode 100644 index 0000000000..154f6543da --- /dev/null +++ b/src/EXTRA-FIX/fix_gjf.h @@ -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 diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp index 60a55bbbb4..5f754b6d2b 100644 --- a/src/fix_langevin.cpp +++ b/src/fix_langevin.cpp @@ -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);