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.