This is an old revision of the document!
What is GEMC
In a Gibbs Ensemble Monte Carlo (GEMC) simulation, 2 boxes are utilized containing a vapor and liquid phase. In order to equilibrate the system, a molecule is allowed to perform 3 types of moves:
1)Translations, rotations, and conformational changes
2)Volume exchanges
3)Particle Swaps.
These particles are swapped between boxes to equilibrate the chemical potential, volume moves equilibrate pressure, and the rest of the moves within a box are performed to maintain thermal equilibrium. The main advantage of the GEMC simulation is that coexisting phases can be observed without a physical interface using a unified partition function. Thus, we used GEMC simulations to determined vapor-liquid coexistence curves for a system.
Files required to run GEMC
In order to run GEMC certain input files are needed; the two main input files for each box gemc_nvt_box1.inp.zipgemc_nvt_box2.inp.zip, a topology filetopology_atoms_wat.psf.zip for the particular component, and a bias file that contain information of the approximate potential for that component.bias_template.inp.zip. For example, we have provided a sample files for 64 water molecules.
Sample of input files
Starting with the input gem_nvt_box1.inp file first, we note that the location of the basis set file, as well as the potential file, is declared( in this case these two files are present in the current working directory). FORCE_EVAL initializes the parameters need to calculate the energy and forces to describe your system. The Quickstep method is used in order to declare the use of electronic structure methods.
&FORCE_EVAL METHOD Quickstep &DFT BASIS_SET_FILE_NAME GTH_BASIS_SETS POTENTIAL_FILE_NAME POTENTIAL
The SCF section in the code below generates atomic density.
&MGRID CUTOFF 280 &END MGRID &QS &END QS &SCF SCF_GUESS ATOMIC &END SCF
One can increase the SCF iteration by including the
MAX_SCF
command. Adding
EPS_SCF
declares the expected SCF convergence. Continuing down the code of the input file, the section shown below details which functional we intend on using, in our case BLYP. We also use the DFTD2 dispersion correction. Additionally, the XC_DERIV function applies derivatives.
&XC &XC_FUNCTIONAL BLYP &END XC_FUNCTIONAL &VDW_POTENTIAL POTENTIAL_TYPE PAIR_POTENTIAL &PAIR_POTENTIAL R_CUTOFF 40.0 TYPE DFTD2 REFERENCE_FUNCTIONAL BLYP &END PAIR_POTENTIAL &END VDW_POTENTIAL &XC_GRID XC_DERIV SPLINE2 XC_SMOOTH_RHO NONE &END XC_GRID &END XC
The code snippet below declares the cell and cell ref size.
&CELL ABC 13.7151207699 13.7151207699 13.7151207699 &CELL_REF ABC 13.7151207699 13.7151207699 13.7151207699 &END CELL_REF
After the above section of code, there is a listing of coordinates for each atom. After this, the section
&KIND H BASIS_SET TZV2P-GTH POTENTIAL GTH-BLYP-q1 &END KIND
declares the basis set intended to be used for the simulation.
&MOTION
&MC ENSEMBLE GEMC_NVT TEMPERATURE 398.0 IPRINT 1 LBIAS yes LSTOP yes NMOVES 8 NSWAPMOVES 640 NSTEP 5 PRESSURE 1.013 RESTART no BOX2_FILE_NAME GEMC_NVT_box1.inp RESTART_FILE_NAME mc_restart_2
It should be noted that the condition 'LBIAS yes' must be true, as we pre sample moves in a different potential. The LSTOP function determines whether the simulation increment will be in cycles (yes), or in steps (no). The NSTEP feature gives the number of MC cycles in a particular simulation run, and should be adjusted according to the length of the simulation. The line 'RESTART no' should only be set to 'no' for the initial run, then switched to 'yes' after the first simulation run is complete. The line BOX2_FILE_NAME GEMC_NVT_box2.inp gives the file name of the input for Box2 and uses it as a reference such that the two input files are read together. Box 2 has a similar line that references Box 1, shown below
BOX2_FILE_NAME GEMC_NVT_box1.inp
Sample of output files
Sample input files for 64 H2O molecules are provided and explained below. As the input file for Box 1 and Box2 has already been explained above, we will not go into further detailgemc_nvt_box1.inp.zipgemc_nvt_box2.inp.zip. We include a topology file,topology_atoms_wat.psf.zip. This file is molecule specific and will consequently change accordingly. We additionally include a sample output file,output4.out.zip. Information used to calculate the density at the end of the run is towards the end of the file and looks like the following:
| BOX 1 | ------------------------------------------------ ******************************************************************************** Average Energy [Hartrees]: -1085.04223205 Average number of molecules: 63.00000000 Average Volume [angstroms**3]: 2579.876452 -------------------------------------------------------------------------------- Quickstep Moves Attempted Accepted Percent 4 1 25.000 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- Move Data for Molecule Type 1 -------------------------------------------------------------------------------- Conformational Moves Attempted Accepted Percent 9 2 22.222 Bond Changes Attempted Accepted Percent 5 2 40.000 -------------------------------------------------------------------------------- Angle Changes Attempted Accepted Percent 4 0 0.000 -------------------------------------------------------------------------------- Conformational Moves Rejected BecauseBox Was Empty: 0 ------------------------------------------------------------------------------- Translation Moves Attempted Accepted Percent 12 1 8.333 -------------------------------------------------------------------------------- Rotation Moves Attempted Accepted Percent 11 1 9.091 -------------------------------------------------------------------------------- Biased Move Data -------------------------------------------------------------------------------- Bond Changes Attempted Accepted Percent 5 5 100.000 -------------------------------------------------------------------------------- Angle Changes Attempted Accepted Percent 4 4 100.000 -------------------------------------------------------------------------------- Translation Moves Attempted Accepted Percent 12 8 66.667 -------------------------------------------------------------------------------- Rotation Moves Attempted Accepted Percent 11 5 45.455 -------------------------------------------------------------------------------- ********************************************************************************
Additionally, other relevant information is included in this file, such as the percentage of accepted moves.