To calculate the elastic constants, the crystal was strained in various directions, and the corresponding pressures (stresses) were collected. However, prior to apply any deformation, the structure has to be relaxed to its equilibrium shape and volume, and the atomic positions have to be fully optimized. As such, geometry optimization was conducted by energy minimization of the atomic system using the conjugate gradient (CG) algorithm, and was coupled with anisotropic box relaxation at zero pressure. Then, another energy minimization stage was performed but this time using the fire algorithm at a fixed volume and shape. Because the fire algorithm is more efficient than CG but cannot be coupled with box relaxation. For both minimization stages, the stopping tolerance for energy and force was set to 10 -15 kcal mol -1 Å -1 , while the maximum iterations and evaluations of force/energy were limited to 10 5 . To ensure that the structure is fully optimized, the two minimization stages were looped for 100 cycles while slightly disturbing the atoms between the cycles. At the end of geometry optimization stages, the optimized atomic positions and other simulation settings were saved to a restart-file. It worth noting that, for some of the MOFs under study, the angle and dihedral styles were hybrid, so such information cannot be stored together with the other force field settings in the restart file. Therefore, angle and dihedral coefficients must be provided separately whenever the restart-file is recalled. Next, positive and negative deformations (±1 % strain) were applied to the relaxed crystal, in various directions (x, y, z, yz, xz, and xy). After applying each deformation, the atomic positions were optimized through energy minimization using the fire algorithm (at fixed volume and shape), with the same above-mentioned setting. Finally, the yielded pressures (stresses) were collected and correspondingly the elastic constants were calculated.