LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands

fix gcmc command

Syntax: 

fix ID group-ID gcmc N X M type seed T mu displace keyword values ...

Examples: 

fix 2 gas gcmc 10 1000 1000 2 29494 298.0 -0.5 0.01
fix 3 water gcmc 10 100 100 0 3456543 3.0 -2.5 0.1 mol my_one_water maxangle 180 full_energy
fix 4 my_gas gcmc 1 10 10 1 123456543 300.0 -12.5 1.0 region disk

Description:

This fix performs grand canonical Monte Carlo (GCMC) exchanges of atoms or molecules of the given type with an imaginary ideal gas reservoir at the specified T and chemical potential (mu) as discussed in (Frenkel). If used with the fix nvt command, simulations in the grand canonical ensemble (muVT, constant chemical potential, constant volume, and constant temperature) can be performed. Specific uses include computing isotherms in microporous materials, or computing vapor-liquid coexistence curves.

Perform up to X exchanges (on average) of gas atoms or molecules of the given type between the simulation domain and the imaginary reservoir every N timesteps. Also perform M (on average) Monte Carlo displacements or rotations (for molecules) of gas of the given type within the simulation domain. M should typically be chosen to be approximately equal to the expected number of gas atoms or molecules of the given type within the domain, which will result in roughly one MC translation per atom or molecule per MC cycle.

For MC moves of molecular gasses, rotations and translations are each attempted with 50% probability. For MC moves of atomic gasses, translations are attempted 100% of the time. For MC exchanges of either molecular or atomic gasses, deletions and insertions are each attempted with 50% probability.

All inserted particles are assigned to two groups: the default group "all" and the group specified in the fix gcmc command (which can also be "all"). If inserted particles are individual atoms, they are assigned the specified atom type. If they are molecules, the type of each atom in the inserted molecule is specified in the file read by the molecule command, and those values are added to the specified atom type. E.g. if type = 2, and the file specifies atom types 1,2,3, then the inserted molecule will have atom types 3,4,5.

This fix cannot be used to perform MC insertions of gas atoms or molecules other than the exchanged type, but MC deletions, translations, and rotations can be performed on any atom/molecule in the fix group. All atoms in the simulation domain can be moved using regular time integration displacements, e.g. via fix_nvt, resulting in a hybrid GCMC+MD simulation. A smaller-than-usual timestep size may be needed when running such a hybrid simulation, especially if the inserted molecules are not well equilibrated.

This command may optionally use the region keyword to define an exchange and move volume. The specified region must have been previously defined with a region command. It must be defined with side = in. Insertion attempts occur only within the specified region. Move and deletion attempt candidates are selected from gas atoms or molecules within the region. If no candidate can be found within the specified region after randomly selecting candidates 1000 times, the move or deletion attempt is considered a failure. Moves must start and end with the atom or molecule center-of-mass within the specified region. Moves are not otherwise attempted.

If used with fix_nvt, the temperature of the imaginary reservoir, T, should be set to be equivalent to the target temperature used in fix_nvt. Otherwise, the imaginary reservoir will not be in thermal equilibrium with the simulation domain.

Note that neighbor lists are re-built every timestep that this fix is invoked, so you should not set N to be too small. However, periodic rebuilds are necessary in order to avoid dangerous rebuilds and missed interactions. Specifically, avoid performing so many MC displacements per timestep that atoms can move beyond the neighbor list skin distance. See the neighbor command for details.

When an atom or molecule is to be inserted, its center-of-mass coordinates are chosen as a random position within the current simulation domain, and new atom velocities are randomly chosen from the specified temperature distribution given by T. Relative coordinates for atoms in a molecule are taken from the template molecule provided by the user. A random initial rotation is used in the case of molecule insertions.

Individual atoms are inserted, unless the mol keyword is used. It specifies a template-ID previously defined using the molecule command, which reads a file that defines the molecule. The coordinates, atom types, charges, etc, as well as any bond/angle/etc and special neighbor information for the molecule can be specified in the molecule file. See the molecule command for details. The only settings required to be in this file are the coordinates and types of atoms in the molecule.

When not using the mol keyword, you should ensure you do not delete only a portion of a molecule (only some of its atoms), or LAMMPS will soon generate an error when it tries to find those atoms. LAMMPS will warn you if any of the atoms eligible for deletion have a non-zero molecule ID, but does not check for this at the time of deletion.

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 fix shake command which also appears in your input script.

Optionally, users may specify the maximum rotation angle for molecular rotations using the maxangle keyword and specifying the angle in degrees. The specified angle will apply to all three Euler angles used internally to define the rotation matrix for molecular rotations. The max angle can be set to zero, but rotations will be pointless. Note that the default is ten degrees for each Euler angle. Also note that fix GCMC does not do configurational bias MC. In other words, while inserted molecules have different rotations, they all have the same molecular configuration, which was specified in the molecule command input.

For atomic gasses, inserted atoms have the specified atom type, but deleted atoms are any atoms that have been inserted or that belong to the user-specified fix group. For molecular gasses, exchanged molecules use the same atom types as in the template molecule supplied by the user. In both cases, exchanged atoms/molecules are assigned to two groups: the default group "all" and the group specified in the fix gcmc command (which can also be "all").

The gas reservoir pressure can be specified using the pressure keyword, in which case the user-specified chemical potential is ignored. For non-ideal gas reservoirs, the user may also specify the fugacity coefficient using the fugacity_coeff keyword.

The full_energy option means that fix GCMC will compute the total potential energy of the entire simulated system. The total system energy before and after the proposed GCMC move is then used in the Metropolis criterion to determine whether or not to accept the proposed GCMC move. By default, this option is off, in which case only partial energies are computed to determine the difference in energy that would be caused by the proposed GCMC move.

The full_energy option is needed for systems with complicated potential energy calculations, including the following:

Some fixes have an associated potential energy. Examples of such fixes include: efield, gravity, addforce, langevin, restrain, temp/berendsen, temp/rescale, and wall fixes. For that energy to be included in the total potential energy of the system (the quantity used when performing GCMC moves), you MUST enable the fix_modify energy option for that fix. The doc pages for individual 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 compute_group_group documentation for more details about simulating non-neutral systems with kspace on.

Use of this fix typically will cause the number of atoms to fluctuate, therefore, you will want to use the compute_modify command to insure that the current number of atoms is used as a normalizing factor each time temperature is computed. Here is the necessary command: 

compute_modify thermo_temp dynamic yes

If LJ units are used, note that a value of 0.18292026 is used by this fix as the reduced value for Planck's constant. This value was derived from LJ paramters for argon, where h* = h/sqrt(sigma^2 * epsilon * mass), sigma = 3.429 angstroms, epsilon/k = 121.85 K, and mass = 39.948 amu.

Restart, fix_modify, output, run start/stop, minimize info:

This fix writes the state of the fix to binary restart files. This includes information about the random number generator seed, the next timestep for MC exchanges, etc. See the 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.

None of the fix_modify options are relevant to this fix.

This fix computes a global vector of length 8, which can be accessed by various output commands. The vector values are the following global cumulative quantities:

The vector values calculated by this fix are "extensive".

No parameter of this fix can be used with the start/stop keywords of the run command. This fix is not invoked during energy minimization.

Restrictions:

This fix is part of the MC package. It is only enabled if LAMMPS was built with that package. See the Making LAMMPS section 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.

Note that very lengthy simulations involving insertions/deletions of billions of gas molecules may run out of atom or molecule IDs and trigger an error, so it is better to run multiple shorter-duration simulations. Likewise, very large molecules have not been tested and may turn out to be problematic.

Use of multiple fix gcmc commands in the same input script can be problematic if using a template molecule. The issue is that the user-referenced template molecule in the second fix gcmc command may no longer exist since it might have been deleted by the first fix gcmc command. An existing template molecule will need to be referenced by the user for each subsequent fix gcmc command.

Related commands:

fix_atom_swap, fix_nvt, neighbor, fix_deposit, fix_evaporate, delete_atoms

Default:

The option defaults are mol = no, maxangle = 10, full_energy = no.

(Frenkel) Frenkel and Smit, Understanding Molecular Simulation, Academic Press, London, 2002.