# MOF CALF-15 and fix gcmc

As a special student (I already have my PhD from 13 years ago), I’m exempted from assignments and exams so you’re not helping me for a grade. However, I’m doing an assignment on MOF CALF-15 in LAMMPS instead of DLPOLY for self-learning. Also to demonstrate to the prof he should be teaching using LAMMPS instead of obsolete and barely supported DLPOLY in future years.

Here’s the assignment question edited for brevity:

3. In this question we will calculate the CO2 gas adsorption isotherm for a Zn based MOF, named CALF-15.

We will be using a GCMC code developed [… for DL_POLY …]. In the directory /scratch/CHM4390/Assignment3/Grad/Q3 you will find the input files needed to start a GCMC calculation for CO2 guests in the MOF framework. In these GCMC simulations an appropriate force field has already been set up using conventional Lennard-Jones parameters for the MOF framework atoms. However, we will modify it using electrostatic potential fitted charges for CO2.

Using the CO2 geometry from the FIELD file, using WebMO and Gaussian to evaluate the CHELPG ESP charges at the B3LYP/cc-pVTZ level of theory. To do this, build a CO2 molecule in the WebMO builder with the bond distances from the FIELD file and select ‘cc-pVTZ’ as the ‘Basis Set’ type. (Do not optimize the geometry)

[…] Run GCMC simulations at 0.05, 0.1, 0.2, 0.3, 0.4, 0.5 0.6, 0.7, 0.8, 0.9 and 1.0 bar pressure. […]

a. [5] Plot the number of guest molecules as a function of the MC step for the 1 bar simulation. After how many time steps would you say that the number of guest molecules has equilibrated?

b. [5] In this question you will visualize the probability distribution resulting from the simulation at 1 bar. This will provide an indication of the location of the CO2 binding sites in the material. […]

c. [10] Construct a gas adsorption isotherm from a series of GCMC simulations at the various pressures by plotting the amount of CO2 adsorbed in (mmol/g of MOF) as a function of the pressure in bar. Include the statistical errors in the uptake […].

HINTS:
- Once the GCMC calculation is completed […] the OUTPUT file will contain the average number of guest molecules and statistical error estimates.
- The experimental isotherm data at 273K can be found in the file “experimental.dat”.
- The chemical formula of the simulation cell is: Zn16N64C48O32H48

(i) are alphanumeric type labels supported only in input file with pair_coeff but not in PairIJ Coeffs section of read_data file ? if so, what section(s)/line(s) of /src/read_data.cpp would need to be modified to implement this feature, ReadData::parse_coeffs() and/or ReadData::pairIJcoeffs() ? i prefer keeping as much as possible in .data to have a clean readable .in.

ERROR: Expected integer parameter instead of ‘Nc’ in input script or data file (src/read_data.cpp:2461)

(ii) is there a way from lammps to generate a lammpstrj that can be visualized in vmd given a changing number of atoms due to fix gcmc ?

Info) /Applications/VMD 1.9.4a57-arm64-Rev12.app/Contents/vmd/plugins/MACOSXARM64/molfile
[…]
vmd >
[…]
lammpsplugin) Too many atoms in timestep. 211 vs. 208

(iii) temperature is not thermostated at 273.15 and goes to zero quickly even with fix nvt and tfac_insert ?

``````Step          Temp          Press          PotEng         KinEng        Density        Atoms     v_numguests
0   273.15         1.9530805e+08  5572497.7      168.54111      0.6523943            208   0
10   4.3463547e-15  1.8997265e+08  5542949.6      2.6818212e-15  0.66171809           211   3
20   5.1290856e-29  1.8997319e+08  5542949.6      3.1647878e-29  0.66171809           211   3
30   6.2723694e-43  1.8997369e+08  5542949.5      3.8702256e-43  0.66171809           211   3
40   2.0555804e-56  1.9510588e+08  5592474.7      1.26835e-56    0.66171809           211   3
50   5.7490945e-76  1.942311e+08   5592474.7      3.5473504e-76  0.66171809           211   3
60   5.4997385e-97  1.8997665e+08  5542949.3      3.3934909e-97  0.66171809           211   3
70   1.8184378e-109  1.908419e+08   5592474.2      1.1220265e-109  0.67104189           214   6
80   9.4816341e-128  1.8997566e+08  5542950.5      5.8504307e-128  0.66171809           211   3
90   4.7229926e-142  1.899571e+08   5542977.3      2.9142172e-142  0.66171809           211   3
100   1.8440604e-158  1.899722e+08   5542949.4      1.1378363e-158  0.66171809           211   3
110   1.0407035e-170  1.9518321e+08  5592945.5      6.4214286e-171  0.66171809           211   3
120   0              1.9366016e+08  5543044.2      0              0.66171809           211   3
130   0              1.9009203e+08  5543394.9      0              0.66171809           211   3
140   0              1.9000212e+08  5542952.2      0              0.66171809           211   3
150   0              1.9001354e+08  5542950.1      0              0.66171809           211   3
160   0              1.9001261e+08  5542948.5      0              0.66171809           211   3
170   0              1.9088027e+08  5592471.3      0              0.67104189           214   6

ERROR: Lost atoms: original 214 current 211 (src/thermo.cpp:490)

Last command: run 20000
``````

(iv) cont. from (iii), then i get

ERROR: Lost atoms: original 214 current 211 (src/thermo.cpp:490)

DLPOLY files:
CONFIG (15.9 KB)
CONTROL (731 Bytes)
experimental.dat (2.0 KB)
FIELD (3.6 KB)

LAMMPS files:
calf15.in (7.6 KB)
calf15.data (7.9 KB)
CO2.txt (248 Bytes)

calf15.data was created from DLPOLY files using mathematica code:

``````SetDirectory[NotebookDirectory[]];
atomTypes =
"!egrep '^..\\s+[0-9.]+\\s+' q3/_dlpoly/FIELD", {Word, Number,
Number, String}];
types = Union[atomTypes[[All, {1, 2}]]];
triazole = atomTypes[[4 ;; 12]];
oxalate = atomTypes[[13 ;; 18]];
zinc = atomTypes[[{19}]];
"!egrep '^..\\s+..\\s+lj' q3/_dlpoly/FIELD", {Word, Word, Word,
Number, Number}];
"!tail -n +6 q3/_dlpoly/CONFIG", {Word, Number, Number, Number,
Number}];

id = 0;
stream = OpenWrite["q3/calf15.data"];

WriteString[stream,
StringJoin[
"# CHM8173 - A3 - Q3 - LAMMPS data file for CALF-15\n\n\t",
ToString@Length@atoms, " atoms\n\t10 atom types\n\n",
ToString@
StringForm["\t`` `` xlo xhi\n", Min@atoms[[All, 3]] - 1,
Max@atoms[[All, 3]] + 1],
ToString@
StringForm["\t`` `` ylo yhi\n", Min@atoms[[All, 4]] - 1,
Max@atoms[[All, 4]] + 1],
ToString@
StringForm["\t`` `` zlo zhi\n", Min@atoms[[All, 5]] - 1,
Max@atoms[[All, 5]] + 1],
"\t0.0 0.0 0.0 xy xz yz\n\nAtom Type Labels\n\n",
Table[
ToString@StringForm["\t`` ``\n", i, types[[i, 1]]], {i,
Length@types}],
"\nMasses\n\n",
Table[
ToString@StringForm["\t`` ``\n", i, types[[i, 2]]], {i,
Length@types}],
(*"\nPairIJ Coeffs\n\n",
StringJoin@Table[ToString@StringForm["\t`` `` `` ``\n",Sequence@@
vdw[[i,{1,2,4,5}]]],{i,Length@vdw}],*)"\nAtoms\n\n",
Flatten@
Join[
Table[
ToString@
StringForm["\t`` `` `` `` `` `` ``\n", (id++) + 1, i,
atoms[[id, 1]], triazole[[j, 3]],
Sequence @@ atoms[[id, 3 ;;]]], {i, 16}, {j, 9}],
Table[
ToString@
StringForm["\t`` `` `` `` `` `` ``\n", (id++) + 1, i + 16,
atoms[[id, 1]], oxalate[[j, 3]],
Sequence @@ atoms[[id, 3 ;;]]], {i, 8}, {j, 6}],
Table[
ToString@
StringForm["\t`` `` `` `` `` `` ``\n", (id++) + 1, i + 24,
atoms[[id, 1]], zinc[[j, 3]],
Sequence @@ atoms[[id, 3 ;;]]], {i, 16}, {j, 1}]]
]];

Close[stream];
``````

Shouldn’t your force field use partial charges? And if it does, how could you possibly expect correct results using the `lj/cut` pair style??

Parsing “Coeffs” sections in the data file is handled by the same (style class) functions as when issuing the corresponding input file coeff commands.

I consider this less optimal and it specifically kills most of the benefits of using type labels. You can keep your input file just as clean by putting your coeff commands into separate file and read it with the include command. This way you have more flexibility since the data file is no longer fully tied to the force field and adjustment are more easily made in the include file(s) or by selecting between multiple of them.