A worked out GCMC exercise is available from a 2-day ReaxFF workshop.


The GCMC code in ADF-ReaxFF needs the following input files to run:

  • control_MC : The GCMC control file, which holds MC settings and the atoms/molecules to insert/move/delete
  • control : A reaxff control file, in which only a small number of parameters is of interest
  • ffield : A standard reaxff forcefield file
  • geo : A geometry file, preferably in biograph format (code not yet tested with xyz)
  • iopt : Text file that should only contain a “5” (without the quotes)
  • insertData_MC : (optional file) Table used when restarting GCMC simulations

Also, the current version of the GCMC code can only run in serial, so please set the NSCM environment variable to 1 (type “export NSCM=1” (without quotes) in the shell before starting a GCMC reaxff run).

the control_MC file

Lines in the control_MC file that start with ! or # will be ignored, so those can be used for comments. Empty lines are also ignored, so feel free to leave some in the file. Lines with keywords should have their value in the first 8 columns, followed by a couple of spaces (at least 1), followed by the 6-character keyword. The order is free, except for the nmols keyword, which should be the last one. The nmols keyword signals the parser that the next section of control_MC will define X new MC Molecule Types.

This is an example for the control_MC file:

! GCMC control file example
      0   iensmb !select MC ensemble (0=Mu-NVT with fixed volume, 1=Mu-NPT with variable volume)
   5000   niter  !number of MC iterations to do this simulation
      0   nstart !start the iteration counter with an offset, usefull for restarts to avoid double files
  300.0   mctemp !Temperature of the simulation, affects acceptance rate for steps that increase the energy
    0.0   mcpres !NPT pressure in GPa (set to zero for incompressible solid systems unless at very high pressures)
    3.0   rmaxpl !Max radius for atom placement on insert/displace move
    1.2   rminpl !Min radius for atom placement on insert/displace move
   2000   nmctry !Maximum number of trials allowed when inserting or moving a molecule. If the
!                !  rmaxpl and rminpl variables are very strict, this number needs to be large
      1   igcfac !include GC prefactor in probabilities? 0 = no 1 = yes
      0   ivol   !select MC volume calculation technique:
!                !  0: vvacu needed! volume = total volume - occupied volume - specified vacuum volume (vvacu)
!                !  1:               volume = total cell volume
!                !  2: vacc needed!  volume = specified accessible volume (vacc)
!                !  3:               volume = total cell volume - occupied volume
!                !  4: vacc needed!  volume = specified accessible volume (vacc) - occupied volume
  435.0   vacc   !if ivol=2 or ivol=4, specify Vacc in angrstoms^3
    0.0   vvacu  !if ivol=0 specify non-accessible (vacuum) volume Vvacu in angrstoms^3
   0.25   ivlim  !volume change limit (value between between 0 and 1, Vnew = ((1+ivlim)*V1)
      1   resopt !write restart files: 0=no, 1=yes
      1   resfrq !frequency of writing restart files (MC code only writes files if the MC move is accepted)
      0   debug  !print debug output if set to 1, print even more debug output when set to 2
      5   nmols  !Number of MC molecule types, must match the number of molecule blocks that follow!
! Molecule Specific Data: C2H2 example
!    This part is fixed format!
!    We need cmpot on line 1,
!    possibly followed by the noinsr on line two,
!    and forced to be ended with nmatom on line 2 or 3, followed by nmatom lines of coordinates.
! the coordinates are FIXED FORMAT! (24d.15,1x,A2) x,y,z (24 wide, 15 after decimal), 1 space, 2chars symbol)
 -75.00   cmpot   !chemical potential of molecule
      4   nmatom  !number of atoms in molecule
      12.180480000000000       0.421696000000000       1.316689000000000 C
      13.124731000000000       0.376902000000000       0.568360000000000 C
      11.349475000000000       0.459560000000000       1.988208000000000 H
      13.957314000000000       0.335843000000000      -0.101000000000000 H

!Molecule Specific Data: C2H4 example
 -75.00   cmpot   !chemical potential of molecule
      6   nmatom  !number of atoms in molecule
      13.989222000000000       0.405391000000000       1.000150000000000 C
      13.316784000000000       0.399646000000000       0.885795000000000 C
      11.494513000000000       0.461837000000000       1.970612000000000 H
      11.335219000000000       0.353577000000000       0.129581000000000 H
      13.811701000000000       0.340224000000000      -0.084000000000000 H
      13.970561000000000       0.453325000000000       1.756236000000000 H

!Molecule Specific Data: H2O Example
 -75.00   cmpot   !chemical potential of molecule
      1   noinsr  !setting this to 1 disables insert/deletion moves. If it is set to 1 for all types, the ensamble becomes NVT/NPT
      3   nmatom  !number of atoms in molecule
      39.996720000000000      40.747660000000000      40.512210000000000 H
      40.000210000000000      39.999520000000000      39.934730000000000 O
      40.000030000000000      39.259880000000000      40.523700000000000 H

!Molecule Specific Data: H2 Example
 -75.00   cmpot   !chemical potential of molecule
      2   nmatom  !number of atoms in molecule
       5.025812000000000      0.0000000000000000       0.000000000000000 H
       5.774188000000000      0.0000000000000000       0.000000000000000 H

!Molecule Specific Data: Single atom example
 -75.00   cmpot   !chemical potential of molecule
      1   nmatom  !number of atoms in molecule
       0.000000000000000      0.0000000000000000       0.000000000000000 H

The Molecule Specific Data blocks define the molecules (or atoms) that can be inserted/moved/deleted with the MC code. The atoms named here should of course be in the forcefield files, and the coordinates should form a reasonable structure. The MC code uses these coordinates during the insertion step by giving them a random rotation, followed by a random translation to generate a random position of the molecule inside the box. Currently, there is no check to make sure the molecule stays inside the boundries of the box, the code only checks that the rmaxpl/rminpl values are satisfied. If you plan on inserting large molecules, make sure there is enough room in the rmaxpl value, otherwise the code will stop with an error message.

The chemical potential (cmpot) keyword

The cmpot keyword sets the chemical potential of the molecule (or atom) reservoir, and is employed when calculating the Boltzmann accept/reject criteria after a MC move is executed. This value can be derived from first principles using statistical mechanics, or equivalently, it can be determined from thermochemical tables available in literature sources.

For example, the proper chemical potential for a GCMC simulation in which single oxygen atoms are exchanged with a reservoir of O2 gas, should equal 1/2 the chemical potential of O2 at the temperature and pressure of the reservoir:

cmpot = Mu_O(T,P) = 1/2*Mu_O2(T,P) = 1/2 * [Mu_ref(T,P_ref) + kT*Log(P/Pref) - E_diss]

where the reference chemical potential [Mu_ref(T,P_ref)] is the experimentally determined chemical potential of O2 at T and Pref, kT*Log(P/Pref) is the pressure correction to the free energy, and E_diss is the dissociation energy of the O2 molecule.

The no insert (noinsr) keyword

The noinsr keyword tells the GCMC code to keep the number of molecules/atoms of this type fixed. It will thus disable Insert/Delete moves on this type, meaning it can only do a displacement move, or volume change move (if the iensmb keyword is set to 1).

the control file

The control file is a regular reaxff control file and it influences the minimization step after an MC trial move. Because of this, only a small number of the reaxff keywords are used during the GCMC simulation.

An example of the control file:

# some of the parameters that influence the minimization step in the GCMC code
      1 icentr     Put the center of mass at the center of the cube
      1 igeofo     0:xyz-input geometry 1: Biograf input geometry 2: xmol-input geometry
2.50000 endmm      End point criterium for MM energy minimisation
    500 imaxit     Maximum number of iterations
      0 icelop     Optimize cell parameters 0=no 1=yes
1.00050 celopt     Cell parameter change
      0 imaxmo     In this case: 0: POLAK_RIBIERE Conj.Grad method, 1: Limited-memory BFGS method

The code has been tested with various imaxit and endmm values, the other options have not been fully tested. Other reaxff keywords might also influence the minimization procedure, but those are best left to their default settings

the ffield file

The ffield file should be an normal reaxff forcefield file, as described in the reaxff documentation by A. van Duin (visit the documentation section on the SCM website to obtain this document).

the geo file

The GCMC code has been tested with biograph input files, but other input formats might work. The details of this file are also described in the original reaxff documentation by van Duin.

the iopt file

The iopt file is a text file with a single digit inside that selects the execution mode of the reaxff code. To run the GCMC code, this file should contain a “5” (without the quotes).

the insertData_MC file

The GCMC code can insert multiple atom/molecule types in a single simulation, so it needs to keep track of what atom belongs to which insert. This information is automatically stored and updated when insertion/deletion/moving of atoms or molecules during the simulation, but is by default unknown of the atoms of the starting geometry. The GCMC code will therefore by default not modify the atoms in the original input in the MC trial moves (keep in mind that they can move around during the minimization step). the insertData_MC file can be used to tell the GCMC code what atoms in the geo file belong to which molecule.

An example of the insertData_MC file:

#   atomNumber  MCInsMolType  MCInsertNmbr
      30       1       1
      40       2       1
      46       2       1
      47       1       2
      48       1       3

This example specifies 4 molecules/atoms that are modifiable by the GCMC code, belonging to 2 different GCMC molecules/atoms that are defined in the control_MC file. The first “molecule” in the control_MC file should thus consist of a single atom (if this doesn’t match, the code will most likely crash!). It was inserted three times (atom 30, 47 and 48) The second molecule has two atoms, and was inserted once.

The atoms do not have a fixed order, and not all atoms have to be defined. If an atom is not appointed to a certain MCInsMolType and MCInsertNmbr, if will simply not be modified during the MC moves. The insertData_MCXXXXXX files generated by the restart option of the code can be directly used as valid insertData_MC files, just remove the digits from the filename and replace the geo file with the corresponding geo_MCXXXXXX file.