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.
reading atom labelmap …
ERROR: Expected integer parameter instead of ‘Nc’ in input script or data file (src/read_data.cpp:2461)
Last command: read_data calf15.data
(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 =
ReadList[
"!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}]];
vdw = ReadList[
"!egrep '^..\\s+..\\s+lj' q3/_dlpoly/FIELD", {Word, Word, Word,
Number, Number}];
atoms = ReadList[
"!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];