Merge pull request #2230 from evoyiatzis/widom

Implementation of Widom insertions through a new fix widom command
This commit is contained in:
Axel Kohlmeyer 2020-07-18 11:40:04 -04:00 committed by GitHub
commit 211beaee48
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
18 changed files with 21449 additions and 35 deletions

View File

@ -244,3 +244,4 @@ OPT.
* :doc:`wall/region <fix_wall_region>`
* :doc:`wall/region/ees <fix_wall_ees>`
* :doc:`wall/srd <fix_wall_srd>`
* :doc:`widom <fix_widom>`

View File

@ -387,6 +387,7 @@ accelerated styles exist.
* :doc:`wall/region <fix_wall_region>` - use region surface as wall
* :doc:`wall/region/ees <fix_wall_ees>` - use region surface as wall for ellipsoidal particles
* :doc:`wall/srd <fix_wall_srd>` - slip/no-slip wall for SRD particles
* :doc:`widom <fix_widom>` - Widom insertions of atoms or molecules
Restrictions
""""""""""""

View File

@ -68,7 +68,7 @@ Description
This fix performs grand canonical Monte Carlo (GCMC) exchanges of
atoms or molecules with an imaginary ideal gas
reservoir at the specified T and chemical potential (mu) as discussed
in :ref:`(Frenkel) <Frenkel>`. It also
in :ref:`(Frenkel) <Frenkel2>`. It also
attempts Monte Carlo (MC) moves (translations and molecule
rotations) within the simulation cell or
region. If used with the :doc:`fix nvt <fix_nh>`
@ -481,7 +481,7 @@ listed above.
----------
.. _Frenkel:
.. _Frenkel2:
**(Frenkel)** Frenkel and Smit, Understanding Molecular Simulation,
Academic Press, London, 2002.

227
doc/src/fix_widom.rst Normal file
View File

@ -0,0 +1,227 @@
.. index:: fix widom
fix widom command
=================
Syntax
""""""
.. parsed-literal::
fix ID group-ID widom N M type seed T keyword values ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* widom = style name of this fix command
* N = invoke this fix every N steps
* M = number of Widom insertions to attempt every N steps
* type = atom type for inserted atoms (must be 0 if mol keyword used)
* seed = random # seed (positive integer)
* T = temperature of the system (temperature units)
* zero or more keyword/value pairs may be appended to args
.. parsed-literal::
keyword = *mol*\ , *region*\ , *full_energy*, *charge*\ , *intra_energy*
*mol* value = template-ID
template-ID = ID of molecule template specified in a separate :doc:`molecule <molecule>` command
*region* value = region-ID
region-ID = ID of region where Widom insertions are allowed
*full_energy* = compute the entire system energy when performing Widom insertions
*charge* value = charge of inserted atoms (charge units)
*intra_energy* value = intramolecular energy (energy units)
Examples
""""""""
.. code-block:: LAMMPS
fix 2 gas widom 1 50000 1 19494 2.0
fix 3 water widom 1000 100 0 29494 300.0 mol h2omol full_energy
Description
"""""""""""
This fix performs Widom insertions of atoms or molecules at the given
temperature as discussed in :ref:`(Frenkel) <Frenkel1>`. Specific uses
include computation of Henry constants of small molecules in microporous
materials or amorphous systems.
Every N timesteps the fix attempts M number of Widom insertions of atoms
or molecules.
If the *mol* keyword is used, only molecule insertions are performed.
Conversely, if the *mol* keyword is not used, only atom insertions are
performed.
This command may optionally use the *region* keyword to define an
insertion volume. The specified region must have been previously
defined with a :doc:`region <region>` command. It must be defined with
side = *in*\ . Insertion attempts occur only within the specified
region. For non-rectangular regions, random trial points are generated
within the rectangular bounding box until a point is found that lies
inside the region. If no valid point is generated after 1000 trials, no
insertion is performed. If an attempted insertion places the atom or
molecule center-of-mass outside the specified region, a new attempted
insertion is generated. This process is repeated until the atom or
molecule center-of-mass is inside the specified region.
Note that neighbor lists are re-built every timestep that this fix is
invoked, so you should not set N to be too small. See the :doc:`neighbor
<neighbor>` command for details.
When an atom or molecule is to be inserted, its coordinates are chosen
at a random position within the current simulation cell or region.
Relative coordinates for atoms in a molecule are taken from the
template molecule provided by the user. The center of mass of the
molecule is placed at the insertion point. The orientation of the
molecule is chosen at random by rotating about this point.
Individual atoms are inserted, unless the *mol* keyword is used. It
specifies a *template-ID* previously defined using the :doc:`molecule
<molecule>` command, which reads a file that defines the molecule. The
coordinates, atom types, charges, etc., as well as any bonding and
special neighbor information for the molecule can be specified in the
molecule file. See the :doc:`molecule <molecule>` command for details.
The only settings required to be in this file are the coordinates and
types of atoms in the molecule.
If you wish to insert molecules via the *mol* keyword, that will have
their bonds or angles constrained via SHAKE, use the *shake* keyword,
specifying as its value the ID of a separate :doc:`fix shake
<fix_shake>` command which also appears in your input script.
Note that fix widom does not use configurational bias MC or any other
kind of sampling of intramolecular degrees of freedom. Inserted
molecules can have different orientations, but they will all have the
same intramolecular configuration, which was specified in the molecule
command input.
For atoms, inserted particles have the specified atom type. For
molecules, they use the same atom types as in the template molecule
supplied by the user.
The excess chemical potential mu_ex is defined as:
.. math::
\mu_{ex} = -kT \ln(<\exp(-(U_{N+1}-U_{N})/{kT})>)
where *k* is Boltzman's constant, *T* is the user-specified temperature,
U_N and U_{N+1} is the potential energy of the system with N and N+1
particles.
The *full_energy* option means that the fix calculates the total
potential energy of the entire simulated system, instead of just the
energy of the part that is changed. By default, this option is off, in
which case only partial energies are computed to determine the energy
difference due to the proposed change.
The *full_energy* option is needed for systems with complicated
potential energy calculations, including the following:
* long-range electrostatics (kspace)
* many-body pair styles
* hybrid pair styles
* eam pair styles
* tail corrections
* need to include potential energy contributions from other fixes
In these cases, LAMMPS will automatically apply the *full_energy*
keyword and issue a warning message.
When the *mol* keyword is used, the *full_energy* option also includes
the intramolecular energy of inserted and deleted molecules, whereas
this energy is not included when *full_energy* is not used. If this is
not desired, the *intra_energy* keyword can be used to define an amount
of energy that is subtracted from the final energy when a molecule is
inserted, and subtracted from the initial energy when a molecule is
deleted. For molecules that have a non-zero intramolecular energy, this
will ensure roughly the same behavior whether or not the *full_energy*
option is used.
Some fixes have an associated potential energy. Examples of such fixes
include: :doc:`efield <fix_efield>`, :doc:`gravity <fix_gravity>`,
:doc:`addforce <fix_addforce>`, :doc:`restrain <fix_restrain>`, and
:doc:`wall fixes <fix_wall>`. For that energy to be included in the
total potential energy of the system (the quantity used when performing
Widom insertions), you MUST enable the :doc:`fix_modify <fix_modify>`
*energy* option for that fix. The doc pages for individual :doc:`fix
<fix>` commands specify if this should be done.
Use the *charge* option to insert atoms with a user-specified point
charge. Note that doing so will cause the system to become non-neutral.
LAMMPS issues a warning when using long-range electrostatics (kspace)
with non-neutral systems. See the :doc:`compute group/group
<compute_group_group>` documentation for more details about simulating
non-neutral systems with kspace on.
**Restart, fix_modify, output, run start/stop, minimize info:**
This fix writes the state of the fix to :doc:`binary restart files
<restart>`. This includes information about the random number generator
seed, the next timestep for Widom insertions etc. See the
:doc:`read_restart <read_restart>` command for info on how to re-specify
a fix in an input script that reads a restart file, so that the
operation of the fix continues in an uninterrupted fashion.
.. note::
For this to work correctly, the timestep must **not** be changed
after reading the restart with :doc:`reset_timestep
<reset_timestep>`. The fix will try to detect it and stop with an
error.
None of the :doc:`fix_modify <fix_modify>` options are relevant to this
fix.
This fix computes a global vector of length 3, which can be accessed by
various :doc:`output commands <Howto_output>`. The vector values are
the following global cumulative quantities:
* 1 = average excess chemical potential on each timestep
* 2 = average difference in potential energy on each timestep
* 3 = volume of the insertion region
The vector values calculated by this fix are "extensive".
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix is part of the MC package. It is only enabled if LAMMPS was
built with that package. See the :doc:`Build package <Build_package>`
doc page for more info.
Do not set "neigh_modify once yes" or else this fix will never be
called. Reneighboring is **required**.
Can be run in parallel, but aspects of the GCMC part will not scale well
in parallel. Only usable for 3D simulations.
Related commands
""""""""""""""""
:doc:`fix gcmc <fix_gcmc>`
:doc:`fix atom/swap <fix_atom_swap>`,
:doc:`neighbor <neighbor>`,
:doc:`fix deposit <fix_deposit>`, :doc:`fix evaporate <fix_evaporate>`,
Default
"""""""
The option defaults are mol = no, intra_energy = 0.0 and full_energy =
no, except for the situations where full_energy is required, as listed
above.
----------
.. _Frenkel1:
**(Frenkel)** Frenkel and Smit, Understanding Molecular Simulation,
Academic Press, London, 2002.

View File

@ -3250,6 +3250,8 @@ whitesmoke
whitespace
Wi
Wicaksono
Widom
widom
Wijk
Wikipedia
wildcard

View File

@ -75,7 +75,7 @@ eim: NaCl using the EIM potential
ellipse: ellipsoidal particles in spherical solvent, 2d system
flow: Couette and Poiseuille flow in a 2d channel
friction: frictional contact of spherical asperities between 2d surfaces
gcmc: Grand Canonical Monte Carlo (GCMC) via the fix gcmc command
gcmc: Grand Canonical MC with fix gcmc, Widom insertion with fix widom
gjf: use of fix langevin Gronbech-Jensen/Farago option
granregion: use of fix wall/region/gran as boundary on granular particles
hugoniostat: Hugoniostat shock dynamics

17322
examples/gcmc/data.spce Normal file

File diff suppressed because it is too large Load Diff

2019
examples/gcmc/data.widom.lj Normal file

File diff suppressed because it is too large Load Diff

28
examples/gcmc/in.widom.lj Normal file
View File

@ -0,0 +1,28 @@
# Kob and Andersen model Phys. Rev. E 51, 4626 (1995)
units lj
atom_style atomic
pair_style lj/cut 2.5
pair_modify shift yes
read_data data.widom.lj
pair_coeff 1 1 1.0 1.0 2.5
pair_coeff 1 2 1.5 0.8 2.0
pair_coeff 2 2 0.5 0.88 2.2
neighbor 0.3 bin
neigh_modify delay 0 every 5 check yes
fix mywidom all widom 10 100000 2 29494 0.75
fix 1 all langevin 0.75 0.75 0.1 48279 zero yes
fix 2 all nve
timestep 0.002
thermo_style custom step temp pe etotal press vol density f_mywidom[1] f_mywidom[2] f_mywidom[3]
thermo 10
run 100

View File

@ -0,0 +1,38 @@
units real
dimension 3
boundary p p p
atom_style full
pair_style lj/cut/coul/long 10.0
bond_style harmonic
angle_style harmonic
read_data data.spce
molecule h2omol H2O.txt
### Flexible SPC/E Potential Parameters ###
### Zhang et al., Fluid Phase Equilibria, 262 (2007) 210-216 ###
pair_coeff 1 1 0.1502629 3.1169
pair_coeff 1 2 0.0341116368 2.04845
pair_coeff 2 2 0.00774378 0.98
bond_coeff 1 176.864 0.9611
angle_coeff 1 42.1845 109.4712
kspace_style pppm 1.0e-4
fix mywidom all widom 10 20 0 29494 298 mol h2omol
fix 2 all nvt temp 298.0 298.0 100.0
neighbor 2.0 bin
neigh_modify delay 10 every 2 check yes
#run variables
timestep 0.5
thermo 10
thermo_style custom step etotal pe temp press vol density f_mywidom[1] f_mywidom[2] f_mywidom[3]
run 100

View File

@ -0,0 +1,89 @@
LAMMPS (30 Jun 2020)
using 1 OpenMP thread(s) per MPI task
# Kob and Andersen model Phys. Rev. E 51, 4626 (1995)
units lj
atom_style atomic
pair_style lj/cut 2.5
pair_modify shift yes
read_data data.widom.lj
orthogonal box = (0.0000000 0.0000000 0.0000000) to (9.4000000 9.4000000 9.4000000)
1 by 1 by 1 MPI processor grid
reading atoms ...
1000 atoms
reading velocities ...
1000 velocities
read_data CPU = 0.003 seconds
pair_coeff 1 1 1.0 1.0 2.5
pair_coeff 1 2 1.5 0.8 2.0
pair_coeff 2 2 0.5 0.88 2.2
neighbor 0.3 bin
neigh_modify delay 0 every 5 check yes
fix mywidom all widom 10 100000 2 29494 0.75
fix 1 all langevin 0.75 0.75 0.1 48279 zero yes
fix 2 all nve
timestep 0.002
thermo_style custom step temp pe etotal press vol density f_mywidom[1] f_mywidom[2] f_mywidom[3]
thermo 10
run 100
Neighbor list info ...
update every 5 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 3.181 | 3.181 | 3.181 Mbytes
Step Temp PotEng TotEng Press Volume Density f_mywidom[1] f_mywidom[2] f_mywidom[3]
0 0.75021245 -6.4204308 -5.2962374 7.2962696 830.584 1.2039721 0 0 830.584
10 0.7358936 -6.4405082 -5.3377717 7.1699962 830.584 1.2039721 -3.8577501 171.3429 830.584
20 0.75426414 -6.4267946 -5.2965298 7.2833985 830.584 1.2039721 -4.0708206 227.63895 830.584
30 0.72947489 -6.4064078 -5.3132896 7.3872583 830.584 1.2039721 -4.4304803 367.7146 830.584
40 0.73504751 -6.4248725 -5.3234038 7.2927069 830.584 1.2039721 -4.1904189 266.99373 830.584
50 0.76497439 -6.4352472 -5.2889331 7.3046861 830.584 1.2039721 -3.8628472 172.51133 830.584
60 0.75752861 -6.4147051 -5.2795485 7.4416 830.584 1.2039721 -3.5355467 111.5042 830.584
70 0.77775078 -6.4210798 -5.2556202 7.4473703 830.584 1.2039721 -3.4754802 102.92223 830.584
80 0.80937104 -6.4320008 -5.2191583 7.4121087 830.584 1.2039721 -3.9287513 188.35625 830.584
90 0.76321255 -6.4203633 -5.2766893 7.4307727 830.584 1.2039721 -4.2257529 279.87337 830.584
100 0.74561743 -6.4010576 -5.2837499 7.52907 830.584 1.2039721 -3.6817835 135.5099 830.584
Loop time of 25.8264 on 1 procs for 100 steps with 1000 atoms
Performance: 669.082 tau/day, 3.872 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.08186 | 0.08186 | 0.08186 | 0.0 | 0.32
Neigh | 0.023613 | 0.023613 | 0.023613 | 0.0 | 0.09
Comm | 0.0053532 | 0.0053532 | 0.0053532 | 0.0 | 0.02
Output | 0.00037837 | 0.00037837 | 0.00037837 | 0.0 | 0.00
Modify | 25.715 | 25.715 | 25.715 | 0.0 | 99.57
Other | | 0.0005643 | | | 0.00
Nlocal: 1000.00 ave 1000 max 1000 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 3049.00 ave 3049 max 3049 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 46176.0 ave 46176 max 46176 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 46176
Ave neighs/atom = 46.176000
Neighbor list builds = 10
Dangerous builds = 0
Total wall time: 0:00:25

View File

@ -0,0 +1,89 @@
LAMMPS (30 Jun 2020)
using 1 OpenMP thread(s) per MPI task
# Kob and Andersen model Phys. Rev. E 51, 4626 (1995)
units lj
atom_style atomic
pair_style lj/cut 2.5
pair_modify shift yes
read_data data.widom.lj
orthogonal box = (0.0000000 0.0000000 0.0000000) to (9.4000000 9.4000000 9.4000000)
1 by 2 by 2 MPI processor grid
reading atoms ...
1000 atoms
reading velocities ...
1000 velocities
read_data CPU = 0.002 seconds
pair_coeff 1 1 1.0 1.0 2.5
pair_coeff 1 2 1.5 0.8 2.0
pair_coeff 2 2 0.5 0.88 2.2
neighbor 0.3 bin
neigh_modify delay 0 every 5 check yes
fix mywidom all widom 10 100000 2 29494 0.75
fix 1 all langevin 0.75 0.75 0.1 48279 zero yes
fix 2 all nve
timestep 0.002
thermo_style custom step temp pe etotal press vol density f_mywidom[1] f_mywidom[2] f_mywidom[3]
thermo 10
run 100
Neighbor list info ...
update every 5 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 7 7 7
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 3.13 | 3.13 | 3.131 Mbytes
Step Temp PotEng TotEng Press Volume Density f_mywidom[1] f_mywidom[2] f_mywidom[3]
0 0.75021245 -6.4204308 -5.2962374 7.2962696 830.584 1.2039721 0 0 830.584
10 0.74279259 -6.4332442 -5.3201696 7.2071344 830.584 1.2039721 -3.4462258 98.984938 830.584
20 0.73016272 -6.4222911 -5.3281423 7.2714238 830.584 1.2039721 -4.4887329 397.41346 830.584
30 0.74416342 -6.429139 -5.3140101 7.2713375 830.584 1.2039721 -4.3313574 322.19056 830.584
40 0.73295457 -6.4370942 -5.3387618 7.227091 830.584 1.2039721 -4.7419281 557.00309 830.584
50 0.76914235 -6.4473959 -5.2948361 7.2179148 830.584 1.2039721 -3.1794069 69.352982 830.584
60 0.74099083 -6.4433012 -5.3329265 7.204973 830.584 1.2039721 -3.5231065 109.66994 830.584
70 0.74514671 -6.4288463 -5.312244 7.2771269 830.584 1.2039721 -1.0154832 3.8727995 830.584
80 0.72787097 -6.4457574 -5.3550427 7.1130709 830.584 1.2039721 -3.6691709 133.2501 830.584
90 0.78452846 -6.5034376 -5.3278217 6.8238659 830.584 1.2039721 -2.0798339 16.008373 830.584
100 0.77188499 -6.487313 -5.3306433 6.9133194 830.584 1.2039721 -2.066579 15.727938 830.584
Loop time of 3.37748 on 4 procs for 100 steps with 1000 atoms
Performance: 5116.239 tau/day, 29.608 timesteps/s
98.7% 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.02025 | 0.020952 | 0.021322 | 0.3 | 0.62
Neigh | 0.0059071 | 0.0060568 | 0.0061543 | 0.1 | 0.18
Comm | 0.016558 | 0.059554 | 0.094282 | 12.2 | 1.76
Output | 0.00022173 | 0.00042927 | 0.0010471 | 0.0 | 0.01
Modify | 3.2552 | 3.2902 | 3.3335 | 1.7 | 97.42
Other | | 0.0003003 | | | 0.01
Nlocal: 250.000 ave 256 max 242 min
Histogram: 1 0 0 0 1 0 0 0 1 1
Nghost: 1666.00 ave 1670 max 1659 min
Histogram: 1 0 0 0 0 0 1 0 0 2
Neighs: 11538.0 ave 11832 max 11091 min
Histogram: 1 0 0 0 0 0 1 1 0 1
Total # of neighbors = 46152
Ave neighs/atom = 46.152000
Neighbor list builds = 10
Dangerous builds = 0
Total wall time: 0:00:03

View File

@ -0,0 +1,133 @@
LAMMPS (30 Jun 2020)
using 1 OpenMP thread(s) per MPI task
units real
dimension 3
boundary p p p
atom_style full
pair_style lj/cut/coul/long 10.0
bond_style harmonic
angle_style harmonic
read_data data.spce
orthogonal box = (-0.031613 -0.023523 -0.085255) to (43.234352 44.939753 42.306533)
1 by 1 by 1 MPI processor grid
reading atoms ...
8640 atoms
scanning bonds ...
2 = max bonds/atom
scanning angles ...
1 = max angles/atom
reading bonds ...
5760 bonds
reading angles ...
2880 angles
2 = max # of 1-2 neighbors
1 = max # of 1-3 neighbors
1 = max # of 1-4 neighbors
2 = max # of special neighbors
special bonds CPU = 0.008 seconds
read_data CPU = 0.028 seconds
molecule h2omol H2O.txt
Read molecule template h2omol:
1 molecules
3 atoms with max type 2
2 bonds with max type 1
1 angles with max type 1
0 dihedrals with max type 0
0 impropers with max type 0
### Flexible SPC/E Potential Parameters ###
### Zhang et al., Fluid Phase Equilibria, 262 (2007) 210-216 ###
pair_coeff 1 1 0.1502629 3.1169
pair_coeff 1 2 0.0341116368 2.04845
pair_coeff 2 2 0.00774378 0.98
bond_coeff 1 176.864 0.9611
angle_coeff 1 42.1845 109.4712
kspace_style pppm 1.0e-4
fix mywidom all widom 10 20 0 29494 298 mol h2omol
fix 2 all nvt temp 298.0 298.0 100.0
neighbor 2.0 bin
neigh_modify delay 10 every 2 check yes
#run variables
timestep 0.5
thermo 10
thermo_style custom step etotal pe temp press vol density f_mywidom[1] f_mywidom[2] f_mywidom[3]
run 100
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:330)
G vector (1/distance) = 0.2690183
grid = 24 24 24
stencil order = 5
estimated absolute RMS force accuracy = 0.024843102
estimated relative force accuracy = 7.4814263e-05
using double precision FFTW3
3d grid and FFT values/proc = 29791 13824
WARNING: Fix Widom using full_energy option (src/MC/fix_widom.cpp:297)
0 atoms in group FixWidom:widom_exclusion_group:mywidom
0 atoms in group FixWidom:rotation_gas_atoms:mywidom
WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (src/neighbor.cpp:487)
Neighbor list info ...
update every 2 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 8 8 8
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 29.66 | 29.66 | 29.66 Mbytes
Step TotEng PotEng Temp Press Volume Density f_mywidom[1] f_mywidom[2] f_mywidom[3]
0 -29703.973 -29703.973 0 -4764.5901 82468.116 1.9140713 0 0 82468.116
10 -29712.131 -31110.041 54.285179 -3154.4423 82468.116 1.9140713 241.93348 3.7366217e-178 82468.116
20 -29711.939 -32614.429 112.71273 -4216.1592 82468.116 1.9140713 16.095006 1.5716469e-12 82468.116
30 -29688.142 -32368.506 104.08688 -4093.6515 82468.116 1.9140713 5.7862327 5.7086352e-05 82468.116
40 -29662.343 -32252.144 100.57005 -1458.5339 82468.116 1.9140713 126.68071 1.2467216e-93 82468.116
50 -29646.78 -32837.635 123.91081 -4607.1155 82468.116 1.9140713 74.622397 1.8790479e-55 82468.116
60 -29628.968 -33001.229 130.9554 -4589.5296 82468.116 1.9140713 3.6575433 0.0020780497 82468.116
70 -29602.78 -32816.28 124.79023 -3082.1133 82468.116 1.9140713 13.983097 5.561247e-11 82468.116
80 -29577.552 -33141.454 138.39742 -6332.8138 82468.116 1.9140713 41.98931 1.6075608e-31 82468.116
90 -29550.865 -33792.115 164.70094 -4607.6419 82468.116 1.9140713 68.690681 4.2082269e-51 82468.116
100 -29515.107 -34052.782 176.21207 -3609.5709 82468.116 1.9140713 41.090597 7.3326206e-31 82468.116
Loop time of 163.407 on 1 procs for 100 steps with 8640 atoms
Performance: 0.026 ns/day, 907.819 hours/ns, 0.612 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 8.5495 | 8.5495 | 8.5495 | 0.0 | 5.23
Bond | 0.031981 | 0.031981 | 0.031981 | 0.0 | 0.02
Kspace | 2.3995 | 2.3995 | 2.3995 | 0.0 | 1.47
Neigh | 5.0542 | 5.0542 | 5.0542 | 0.0 | 3.09
Comm | 0.051965 | 0.051965 | 0.051965 | 0.0 | 0.03
Output | 0.0018802 | 0.0018802 | 0.0018802 | 0.0 | 0.00
Modify | 147.31 | 147.31 | 147.31 | 0.0 | 90.15
Other | | 0.003614 | | | 0.00
Nlocal: 8640.00 ave 8640 max 8640 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 23499.0 ave 23499 max 23499 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 3.27380e+06 ave 3.2738e+06 max 3.2738e+06 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 3273800
Ave neighs/atom = 378.91204
Ave special neighs/atom = 2.0000000
Neighbor list builds = 220
Dangerous builds = 0
Total wall time: 0:02:44

View File

@ -493,11 +493,7 @@ void FixGCMC::init()
}
}
if (full_flag) {
char *id_pe = (char *) "thermo_pe";
int ipe = modify->find_compute(id_pe);
c_pe = modify->compute[ipe];
}
if (full_flag) c_pe = modify->compute[modify->find_compute("thermo_pe")];
int *type = atom->type;
@ -580,18 +576,12 @@ void FixGCMC::init()
// skip if already exists from previous init()
if (full_flag && !exclusion_group_bit) {
char **group_arg = new char*[4];
// create unique group name for atoms to be excluded
int len = strlen(id) + 30;
group_arg[0] = new char[len];
sprintf(group_arg[0],"FixGCMC:gcmc_exclusion_group:%s",id);
group_arg[1] = (char *) "subtract";
group_arg[2] = (char *) "all";
group_arg[3] = (char *) "all";
group->assign(4,group_arg);
exclusion_group = group->find(group_arg[0]);
auto group_id = std::string("FixGCMC:gcmc_exclusion_group:") + id;
group->assign(group_id + " subtract all all");
exclusion_group = group->find(group_id);
if (exclusion_group == -1)
error->all(FLERR,"Could not find fix gcmc exclusion group ID");
exclusion_group_bit = group->bitmask[exclusion_group];
@ -603,34 +593,25 @@ void FixGCMC::init()
char **arg = new char*[narg];;
arg[0] = (char *) "exclude";
arg[1] = (char *) "group";
arg[2] = group_arg[0];
arg[2] = (char *) group_id.c_str();
arg[3] = (char *) "all";
neighbor->modify_params(narg,arg);
delete [] group_arg[0];
delete [] group_arg;
delete [] arg;
}
// create a new group for temporary use with selected molecules
if (exchmode == EXCHMOL || movemode == MOVEMOL) {
char **group_arg = new char*[3];
// create unique group name for atoms to be rotated
int len = strlen(id) + 30;
group_arg[0] = new char[len];
sprintf(group_arg[0],"FixGCMC:rotation_gas_atoms:%s",id);
group_arg[1] = (char *) "molecule";
char digits[12];
sprintf(digits,"%d",-1);
group_arg[2] = digits;
group->assign(3,group_arg);
molecule_group = group->find(group_arg[0]);
auto group_id = std::string("FixGCMC:rotation_gas_atoms:") + id;
group->assign(group_id + " molecule -1");
molecule_group = group->find(group_id);
if (molecule_group == -1)
error->all(FLERR,"Could not find fix gcmc rotation group ID");
molecule_group_bit = group->bitmask[molecule_group];
molecule_group_inversebit = molecule_group_bit ^ ~0;
delete [] group_arg[0];
delete [] group_arg;
}
// get all of the needed molecule data if exchanging

1227
src/MC/fix_widom.cpp Normal file

File diff suppressed because it is too large Load Diff

257
src/MC/fix_widom.h Normal file
View File

@ -0,0 +1,257 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(widom,FixWidom)
#else
#ifndef LMP_FIX_WIDOM_H
#define LMP_FIX_WIDOM_H
#include <cstdio>
#include "fix.h"
namespace LAMMPS_NS {
class FixWidom : public Fix {
public:
FixWidom(class LAMMPS *, int, char **);
~FixWidom();
int setmask();
void init();
void pre_exchange();
void attempt_atomic_insertion();
void attempt_molecule_insertion();
void attempt_atomic_insertion_full();
void attempt_molecule_insertion_full();
double energy(int, int, tagint, double *);
double molecule_energy(tagint);
double energy_full();
void update_gas_atoms_list();
double compute_vector(int);
double memory_usage();
void write_restart(FILE *);
void restart(char *);
void grow_molecule_arrays(int);
private:
int molecule_group,molecule_group_bit;
int molecule_group_inversebit;
int exclusion_group,exclusion_group_bit;
int nwidom_type,nevery,seed;
int ninsertions;
int ngas; // # of gas atoms on all procs
int ngas_local; // # of gas atoms on this proc
int exchmode; // exchange ATOM or MOLECULE
int movemode; // move ATOM or MOLECULE
int regionflag; // 0 = anywhere in box, 1 = specific region
int iregion; // widom region
char *idregion; // widom region id
bool charge_flag; // true if user specified atomic charge
bool full_flag; // true if doing full system energy calculations
int natoms_per_molecule; // number of atoms in each inserted molecule
int nmaxmolatoms; // number of atoms allocated for molecule arrays
double ave_widom_chemical_potential;
int widom_nmax;
int max_region_attempts;
double gas_mass;
double insertion_temperature;
double beta,sigma,volume;
double charge;
double xlo,xhi,ylo,yhi,zlo,zhi;
double region_xlo,region_xhi,region_ylo,region_yhi,region_zlo,region_zhi;
double region_volume;
double energy_stored; // full energy of old/current configuration
double *sublo,*subhi;
int *local_gas_list;
double **cutsq;
double **molcoords;
double *molq;
imageint *molimage;
imageint imagezero;
double energy_intra;
class Pair *pair;
class RanPark *random_equal;
class Atom *model_atom;
class Molecule **onemols;
int imol,nmol;
int triclinic; // 0 = orthog box, 1 = triclinic
class Compute *c_pe;
void options(int, char **);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix Widom does not (yet) work with atom_style template
Self-explanatory.
E: Fix Widom region does not support a bounding box
Not all regions represent bounded volumes. You cannot use
such a region with the fix Widom command.
E: Fix Widom region cannot be dynamic
Only static regions can be used with fix Widom.
E: Fix Widom region extends outside simulation box
Self-explanatory.
E: Fix Widom molecule must have coordinates
The defined molecule does not specify coordinates.
E: Fix Widom molecule must have atom types
The defined molecule does not specify atom types.
E: Atom type must be zero in fix Widom mol command
Self-explanatory.
E: Fix Widom molecule has charges, but atom style does not
Self-explanatory.
E: Fix Widom molecule template ID must be same as atom_style template ID
When using atom_style template, you cannot insert molecules that are
not in that template.
E: Fix Widom atom has charge, but atom style does not
Self-explanatory.
UNDOCUMENTED
E: Molecule template ID for fix Widom does not exist
Self-explanatory.
W: Molecule template for fix Widom has multiple molecules
The fix Widom command will only create molecules of a single type,
i.e. the first molecule in the template.
E: Region ID for fix Widom does not exist
Self-explanatory.
W: Fix Widom using full_energy option
Fix Widom has automatically turned on the full_energy option since it
is required for systems like the one specified by the user. User input
included one or more of the following: kspace, a hybrid
pair style, an eam pair style, tail correction,
or no "single" function for the pair style.
E: Invalid atom type in fix Widom command
The atom type specified in the Widom command does not exist.
E: Fix Widom cannot exchange individual atoms belonging to a molecule
This is an error since you should not delete only one atom of a
molecule. The user has specified atomic (non-molecular) gas
exchanges, but an atom belonging to a molecule could be deleted.
E: All mol IDs should be set for fix Widom group atoms
The molecule flag is on, yet not all molecule ids in the fix group
have been set to non-zero positive values by the user. This is an
error since all atoms in the fix Widom group are eligible for deletion,
rotation, and translation and therefore must have valid molecule ids.
E: Fix Widom molecule command requires that atoms have molecule attributes
Should not choose the Widom molecule feature if no molecules are being
simulated. The general molecule flag is off, but Widom's molecule flag
is on.
E: Cannot use fix Widom in a 2d simulation
Fix Widom is set up to run in 3d only. No 2d simulations with fix Widom
are allowed.
E: Could not find fix Widom exclusion group ID
Self-explanatory.
E: Illegal fix Widom gas mass <= 0
The computed mass of the designated gas molecule or atom type was less
than or equal to zero.
E: Cannot do Widom on atoms in atom_modify first group
This is a restriction due to the way atoms are organized in a list to
enable the atom_modify first command.
W: Fix Widom is being applied to the default group all
This is allowed, but it will result in Monte Carlo moves
being performed on all the atoms in the system, which is
often not what is intended.
E: Could not find specified fix Widom group ID
Self-explanatory.
E: fix Widom does currently not support full_energy option with molecules on more than 1 MPI process.
UNDOCUMENTED
E: Fix Widom put atom outside box
This should not normally happen. Contact the developers.
E: Fix Widom ran out of available molecule IDs
See the setting for tagint in the src/lmptype.h file.
E: Fix Widom ran out of available atom IDs
See the setting for tagint in the src/lmptype.h file.
E: Too many total atoms
See the setting for bigint in the src/lmptype.h file.
*/

View File

@ -584,10 +584,10 @@ void Group::create(char *name, int *flag)
return group index if name matches existing group, -1 if no such group
------------------------------------------------------------------------- */
int Group::find(const char *name)
int Group::find(const std::string &name)
{
for (int igroup = 0; igroup < MAX_GROUP; igroup++)
if (names[igroup] && strcmp(name,names[igroup]) == 0) return igroup;
if (names[igroup] && (name == names[igroup])) return igroup;
return -1;
}

View File

@ -33,7 +33,7 @@ class Group : protected Pointers {
void assign(int, char **); // assign atoms to a group
void assign(const std::string &); // convenience function
void create(char *, int *); // add flagged atoms to a group
int find(const char *); // lookup name in list of groups
int find(const std::string &); // lookup name in list of groups
int find_or_create(const char *); // lookup name or create new group
void write_restart(FILE *);
void read_restart(FILE *);