## 2 Structural Features and Methodology
### Computational Details  
In this work, we propose a new procedure that makes calculations 'computationally tractable' in the complex microporous zeolitic systems. Computationally tractable in this context means having reliable accuracy with excellent computational efficiency. The number of atoms in the unit cell of a clinoptilolite makes the choice of CP2K [29] natural and rational because this open-source software package employs the hybrid Gaussian and plane wave (GPW) approach [30]. That is, a plane wave basis is used to represent the density, and the atom-centered Gaussian orbitals calculate the wave functions. [PERSON] et al. [31] reported that this approach would show even better accuracy if the Goedecker-Teter-Hutter (GTH) pseudopotentials were employed. Those pseudopotentials [32] are complete, solved analytically, transferable, fully relativistic, norm-conserving, based on all-electron calculations, and, like the GPW approach, can be integrated on a real space grid.  
In the present study, first-principles DFT [33; 34] calculations in the framework of dispersion correction [35] were carried out using the GPW approach and GTH pseudopotentials, as implemented in the CP2K/Quickstep program [31]. Atomic positions and cell parameters were optimized using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) [36] minimization scheme. Ionic relaxations were stopped when forces on all atoms were less than \(10^{-5}\) Ha bohr\({}^{-1}\), SCF accuracy was \(10^{-10}\), and displacements were below \(10^{-4}\) bohr. The wave function optimization was carried out using diagonalization with the Broyden mixing approach for the density-mixing. For other systems, the orbital transformation (OT) method [37] was reported to speed up convergence; however, diagonalization with Broyden mixing showed more robust behavior in the following system. Diagonalization solves the non-linear eigenvalue problem and allows for the molecular orbitals to be occupied from the lowest-energy orbital up; hence, the current input parameters can be transferred to metallic systems just by including smearing.  
CP2K divides the space of a periodic unit cell by a set of equidistant grid points; the quality of plane wave expansion within this equidistant grid (computational box) is determined by the plane wave cutoff. The utilization of multigrids stands out as one of the primary factors enhancing the effective mapping of the Kohn-Sham (KS) density matrix and basis functions to the expansion coefficients. Essentially, multigrids are grids with varying mesh sizes; one can think of them as finite small elements used to partition the space of interest, as well as being utilized to calculate densities or matrix elements. After setting up multigrids, CP2K needs to map the Gaussians onto those grids, with the relative cutoff being the foremost determinant of which Gaussian will be assigned to which multigrid  
Figure 5: The possible cation sites M1, M2, and M3 for Na\({}^{+}\), Ca\({}^{2+}\), and K\({}^{+}\), distributed over channels A, B, and C, respectively, in one of our models.  
level. In essence, PW cutoff and relative cutoff play a very important role in the accuracy of the calculation due to their effect on the integration grid. For this, one should and must conduct a convergence test to minimize the total energy. This procedure is thus system-dependent and, for the current system, 350 Ry provided optimal energy minimization and convergence performance, as confirmed by several energy-volume curves. The relative cutoff was set to 70 Ry. The scripts and codes required for such calculation are available upon readers' request.