<p>Dose calculation is the foundation of boron neutron capture therapy (BNCT). MagicDose, a dose calculation program for the BNCT treatment planning system, is developed based on the Monte Carlo method. First, the voxel phantom of the modified Snyder head with 16 and 8 mm is constructed, and the results from MagicDose and MCNP are presented as two-dimensional coordinate points (X_n, Y_n), comparing their relationship relative to the <i>y</i> = <i>x</i> linear function, while analyzing their respective calculation time. A modified Snyder head phantom with a tumor at three different spatial resolutions of 16, 8, and 1 mm was constructed, and the depth-dose rate curves and spatial distribution maps are analyzed. Finally, the patients’ head CT data were used for the application. The results indicate that the calculations from MagicDose and MCNP exhibit high consistency and demonstrate that MagicDose offers superior computational efficiency compared to MCNP, with improvements of approximately 31.24% and 28.65% at spatial resolutions of 16 and 8 mm, respectively. As the spatial resolution increased, the variability in the dose rate results decreased. The voxel size and number of threads are both inversely proportional to the calculation time. For the CT model, a voxel phantom with a spatial resolution of 1 mm × 1 mm × 1 mm is successfully constructed. The calculation results showed that the boron dose rate contribution significantly exceeds that of the other dose components, with the spatial distribution of the total relative biological effect dose rate clearly delineating the boundaries between the high- and low-dose rate regions. The above results verify the correctness of MagicDose, which also provides a reference for optimizing the design of voxel phantoms for clinical treatment.</p>