diff --git a/src/constants/element.rs b/src/constants/element.rs
new file mode 100644
index 0000000000000000000000000000000000000000..227874225d58102fb9f96f9d35b03c309fbd29b8
--- /dev/null
+++ b/src/constants/element.rs
@@ -0,0 +1,381 @@
+use lazy_static::lazy_static;
+use std::collections::HashMap;
+
+pub const SPECIES_NAME: [&str; 118] = ["H", "He",
+        "Li","Be","B", "C","N", "O", "F", "Ne",
+        "Na","Mg","Al","Si","P", "S", "Cl","Ar",
+        "K", "Ca","Sc","Ti","V", "Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr",
+        "Rb","Sr","Y", "Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I", "Xe",
+        "Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W", "Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn",
+        "Fr","Ra","Ac","Th","Pa","U", "Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs","Mt","Ds","Rg","Cn","Nh","Fl","Mc","Lv","Ts","Og"
+        ];
+
+// IUPAC 2021 standard atomic weights (Prohaska et al., Pure Appl. Chem., 2022)
+// Table 1 column 7: abridged to 5 significant figures; interval elements use conventional single values
+// 14 interval elements: H, Li, B, C, N, O, Mg, Si, S, Cl, Ar, Br, Tl, Pb
+// Tc, Pm, Po-At, Rn-Ac, Np and beyond: most stable isotope mass
+//   AME2020 / NUBASE2020 via IUPAC 2021 Table 2
+
+pub const MASS_CHARGE: [(f64,f64);118] = [
+    (1.0080,1.0),       (4.0026,2.0),
+
+    (6.94,3.0),         (9.0122,4.0),
+    (10.81,5.0),        (12.011,6.0),       (14.007,7.0),       (15.999,8.0),
+    (18.998,9.0),       (20.180,10.0),
+
+    (22.990,11.0),      (24.305,12.0),
+    (26.982,13.0),      (28.085,14.0),      (30.974,15.0),      (32.06,16.0),
+    (35.45,17.0),       (39.95,18.0),
+
+    (39.098,19.0),      (40.078,20.0),
+    (44.956,21.0),      (47.867,22.0),      (50.942,23.0),      (51.996,24.0),
+    (54.938,25.0),      (55.845,26.0),      (58.933,27.0),      (58.693,28.0),
+    (63.546,29.0),      (65.38,30.0),
+
+    (69.723,31.0),      (72.630,32.0),
+    (74.922,33.0),      (78.971,34.0),      (79.904,35.0),      (83.798,36.0),
+
+    (85.468,37.0),      (87.62,38.0),       (88.906,39.0),      (91.224,40.0),
+    (92.906,41.0),      (95.95,42.0),       (97.90721,43.0),    (101.07,44.0),
+    (102.91,45.0),      (106.42,46.0),      (107.87,47.0),      (112.41,48.0),
+
+    (114.82,49.0),      (118.71,50.0),      (121.76,51.0),      (127.60,52.0),
+    (126.90,53.0),      (131.29,54.0),
+
+    (132.91,55.0),      (137.33,56.0),
+
+    (138.91,57.0),      (140.12,58.0),      (140.91,59.0),      (144.24,60.0),
+    (144.91276,61.0),   (150.36,62.0),      (151.96,63.0),      (157.25,64.0),
+    (158.93,65.0),      (162.50,66.0),      (164.93,67.0),      (167.26,68.0),
+    (168.93,69.0),      (173.05,70.0),
+
+    (174.97,71.0),      (178.49,72.0),
+    (180.95,73.0),      (183.84,74.0),      (186.21,75.0),      (190.23,76.0),
+    (192.22,77.0),      (195.08,78.0),      (196.97,79.0),      (200.59,80.0),
+
+    (204.38,81.0),      (207.2,82.0),       (208.98,83.0),      (208.98243,84.0),
+    (209.98715,85.0),   (222.01758,86.0),
+
+    (223.01973,87.0),   (226.02541,88.0),
+    (227.02775,89.0),   (232.04,90.0),      (231.04,91.0),      (238.03,92.0),
+
+    (237.04817,93.0),   (244.06420,94.0),   (243.06138,95.0),   (247.07035,96.0),
+    (247.07031,97.0),   (251.07959,98.0),   (252.08298,99.0),   (257.09511,100.0),
+    (258.09843,101.0),  (259.10100,102.0),
+
+    (262.10962,103.0),  (267.12179,104.0),
+    (268.12567,105.0),  (271.13378,106.0),  (270.13337,107.0),  (269.13365,108.0),
+    (278.15649,109.0),  (281.16455,110.0),  (282.16934,111.0),  (285.17723,112.0),
+
+    (286.18246,113.0),  (289.19052,114.0),  (289.19397,115.0),  (293.20458,116.0),
+    (294.21084,117.0),  (294.21398,118.0),
+];
+
+pub const ELEM1ST: [&str;2]  = ["H", "He"];
+pub const ELEM2ND: [&str;8]  = ["Li","Be","B", "C", "N", "O", "F", "Ne"];
+pub const ELEM3RD: [&str;8]  = ["Na","Mg","Al","Si","P", "S", "Cl","Ar"];
+pub const ELEM4TH: [&str;18] = ["K", "Ca","Sc","Ti","V", "Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr"];
+pub const ELEM5TH: [&str;18] = ["Rb","Sr","Y", "Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I", "Xe"];
+pub const ELEM6TH: [&str;32] = ["Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W", "Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn"];
+pub const ELEM7TH: [&str;32] = ["Fr","Ra","Ac","Th","Pa","U", "Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs","Mt","Ds","Rg","Cn","Nh","Fl","Mc","Lv","Ts","Og"];
+
+pub const ELEMTMS: [&str;40] = [
+    "Sc","Ti","V", "Cr","Mn","Fe","Co","Ni","Cu","Zn",
+    "Y", "Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd",
+    "La","Hf","Ta","W", "Re","Os","Ir","Pt","Au","Hg",
+    "Ac","Rf","Db","Sg","Bh","Hs","Mt","Ds","Rg","Cn"
+];
+
+
+lazy_static!{
+    pub static ref SPECIES_INFO: HashMap<&'static str, &'static (f64,f64)> = {
+        let mut m = HashMap::new();
+        SPECIES_NAME.iter().zip(MASS_CHARGE.iter()).for_each(|(name,info)| {
+            m.insert(*name,info);
+        });
+        m
+    };
+}
+
+pub const ATOM_CONFIGURATION: [[usize; 4]; 119] = [
+    [ 0, 0, 0, 0],     //  0  GHOST
+    [ 1, 0, 0, 0],     //  1  H
+    [ 2, 0, 0, 0],     //  2  He
+    [ 3, 0, 0, 0],     //  3  Li
+    [ 4, 0, 0, 0],     //  4  Be
+    [ 4, 1, 0, 0],     //  5  B
+    [ 4, 2, 0, 0],     //  6  C
+    [ 4, 3, 0, 0],     //  7  N
+    [ 4, 4, 0, 0],     //  8  O
+    [ 4, 5, 0, 0],     //  9  F
+    [ 4, 6, 0, 0],     // 10  Ne
+    [ 5, 6, 0, 0],     // 11  Na
+    [ 6, 6, 0, 0],     // 12  Mg
+    [ 6, 7, 0, 0],     // 13  Al
+    [ 6, 8, 0, 0],     // 14  Si
+    [ 6, 9, 0, 0],     // 15  P
+    [ 6,10, 0, 0],     // 16  S
+    [ 6,11, 0, 0],     // 17  Cl
+    [ 6,12, 0, 0],     // 18  Ar
+    [ 7,12, 0, 0],     // 19  K
+    [ 8,12, 0, 0],     // 20  Ca
+    [ 8,12, 1, 0],     // 21  Sc
+    [ 8,12, 2, 0],     // 22  Ti
+    [ 8,12, 3, 0],     // 23  V
+    [ 7,12, 5, 0],     // 24  Cr
+    [ 8,12, 5, 0],     // 25  Mn
+    [ 8,12, 6, 0],     // 26  Fe
+    [ 8,12, 7, 0],     // 27  Co
+    [ 8,12, 8, 0],     // 28  Ni
+    [ 7,12,10, 0],     // 29  Cu
+    [ 8,12,10, 0],     // 30  Zn
+    [ 8,13,10, 0],     // 31  Ga
+    [ 8,14,10, 0],     // 32  Ge
+    [ 8,15,10, 0],     // 33  As
+    [ 8,16,10, 0],     // 34  Se
+    [ 8,17,10, 0],     // 35  Br
+    [ 8,18,10, 0],     // 36  Kr
+    [ 9,18,10, 0],     // 37  Rb
+    [10,18,10, 0],     // 38  Sr
+    [10,18,11, 0],     // 39  Y
+    [10,18,12, 0],     // 40  Zr
+    [ 9,18,14, 0],     // 41  Nb
+    [ 9,18,15, 0],     // 42  Mo
+    [10,18,15, 0],     // 43  Tc
+    [ 9,18,17, 0],     // 44  Ru
+    [ 9,18,18, 0],     // 45  Rh
+    [ 8,18,20, 0],     // 46  Pd
+    [ 9,18,20, 0],     // 47  Ag
+    [10,18,20, 0],     // 48  Cd
+    [10,19,20, 0],     // 49  In
+    [10,20,20, 0],     // 50  Sn
+    [10,21,20, 0],     // 51  Sb
+    [10,22,20, 0],     // 52  Te
+    [10,23,20, 0],     // 53  I
+    [10,24,20, 0],     // 54  Xe
+    [11,24,20, 0],     // 55  Cs
+    [12,24,20, 0],     // 56  Ba
+    [12,24,21, 0],     // 57  La
+    [12,24,21, 1],     // 58  Ce
+    [12,24,20, 3],     // 59  Pr
+    [12,24,20, 4],     // 60  Nd
+    [12,24,20, 5],     // 61  Pm
+    [12,24,20, 6],     // 62  Sm
+    [12,24,20, 7],     // 63  Eu
+    [12,24,21, 7],     // 64  Gd
+    [12,24,21, 8],     // 65  Tb
+    [12,24,20,10],     // 66  Dy
+    [12,24,20,11],     // 67  Ho
+    [12,24,20,12],     // 68  Er
+    [12,24,20,13],     // 69  Tm
+    [12,24,20,14],     // 70  Yb
+    [12,24,21,14],     // 71  Lu
+    [12,24,22,14],     // 72  Hf
+    [12,24,23,14],     // 73  Ta
+    [12,24,24,14],     // 74  W
+    [12,24,25,14],     // 75  Re
+    [12,24,26,14],     // 76  Os
+    [12,24,27,14],     // 77  Ir
+    [11,24,29,14],     // 78  Pt
+    [11,24,30,14],     // 79  Au
+    [12,24,30,14],     // 80  Hg
+    [12,25,30,14],     // 81  Tl
+    [12,26,30,14],     // 82  Pb
+    [12,27,30,14],     // 83  Bi
+    [12,28,30,14],     // 84  Po
+    [12,29,30,14],     // 85  At
+    [12,30,30,14],     // 86  Rn
+    [13,30,30,14],     // 87  Fr
+    [14,30,30,14],     // 88  Ra
+    [14,30,31,14],     // 89  Ac
+    [14,30,32,14],     // 90  Th
+    [14,30,31,16],     // 91  Pa
+    [14,30,31,17],     // 92  U
+    [14,30,31,18],     // 93  Np
+    [14,30,30,20],     // 94  Pu
+    [14,30,30,21],     // 95  Am
+    [14,30,31,21],     // 96  Cm
+    [14,30,31,22],     // 97  Bk
+    [14,30,30,24],     // 98  Cf
+    [14,30,30,25],     // 99  Es
+    [14,30,30,26],     //100  Fm
+    [14,30,30,27],     //101  Md
+    [14,30,30,28],     //102  No
+    [14,30,31,28],     //103  Lr
+    [14,30,32,28],     //104  Rf
+    [14,30,33,28],     //105  Db
+    [14,30,34,28],     //106  Sg
+    [14,30,35,28],     //107  Bh
+    [14,30,36,28],     //108  Hs
+    [14,30,37,28],     //109  Mt
+    [14,30,38,28],     //110  Ds
+    [14,30,39,28],     //111  Rg
+    [14,30,40,28],     //112  Cn
+    [14,31,40,28],     //113  Nh
+    [14,32,40,28],     //114  Fl
+    [14,33,40,28],     //115  Mc
+    [14,34,40,28],     //116  Lv
+    [14,35,40,28],     //117  Ts
+    [14,36,40,28],     //118  Og
+];
+
+// SAD and ECP configuration
+pub const S_SHELL: [f64; 1] = [2.0];
+pub const P_SHELL: [f64; 3] = [2.0, 2.0, 2.0];
+pub const D_SHELL: [f64; 5] = [2.0, 2.0, 2.0, 2.0, 2.0];
+pub const F_SHELL: [f64; 7] = [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0];
+pub const XE_SHELL: [f64; 27] =
+//   1s   2s   2p   2p   2p   3s   3p   3p   3p   4s   3d   3d   3d   3d   3d   4p   4p   4p   5s   4d   4d   4d   4d   4d   5p   5p   5p  
+    [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0];
+pub const KR_SHELL: [f64; 18] =
+//   1s   2s   2p   2p   2p   3s   3p   3p   3p   4s   3d   3d   3d   3d   3d   4p   4p   4p  
+    [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0];
+
+//                                    1s   2s   2p   3s   3p   4s   3d    4p   5s   4d    5p   4f    5d    6s   6p 
+pub const NELE_IN_SHELLS: [f64; 15] = [2.0, 2.0, 6.0, 2.0, 6.0, 2.0, 10.0, 6.0, 2.0, 10.0, 6.0, 14.0, 10.0, 2.0, 6.0];
+
+lazy_static!{
+    pub static ref ATOMIC_RADII: HashMap<&'static str, &'static f64> = {
+    /// 获取元素周期表中所有元素的原子半径（单位：Å）
+    /// 数据来源：Clementi-Raimondi半径、共价半径和范德华半径的综合
+        let mut radii = HashMap::new();
+        
+        // 第1周期
+        radii.insert( "H", &0.53);  // 氢
+        radii.insert("He", &0.31);  // 氦
+
+        // 第2周期
+        radii.insert("Li", &1.67);  // 锂
+        radii.insert("Be", &1.12);  // 铍
+        radii.insert( "B", &0.87);  // 硼
+        radii.insert( "C", &0.67);  // 碳
+        radii.insert( "N", &0.56);  // 氮
+        radii.insert( "O", &0.48);  // 氧
+        radii.insert( "F", &0.42);  // 氟
+        radii.insert("Ne", &0.38);  // 氖
+
+        // 第3周期
+        radii.insert("Na", &1.90);  // 钠
+        radii.insert("Mg", &1.45);  // 镁
+        radii.insert("Al", &1.18);  // 铝
+        radii.insert("Si", &1.11);  // 硅
+        radii.insert( "P", &0.98);  // 磷
+        radii.insert( "S", &0.88);  // 硫
+        radii.insert("Cl", &0.79);  // 氯
+        radii.insert("Ar", &0.71);  // 氩
+
+        // 第4周期
+        radii.insert( "K", &2.43);  // 钾
+        radii.insert("Ca", &1.94);  // 钙
+        radii.insert("Sc", &1.84);  // 钪
+        radii.insert("Ti", &1.76);  // 钛
+        radii.insert( "V", &1.71);  // 钒
+        radii.insert("Cr", &1.66);  // 铬
+        radii.insert("Mn", &1.61);  // 锰
+        radii.insert("Fe", &1.56);  // 铁
+        radii.insert("Co", &1.52);  // 钴
+        radii.insert("Ni", &1.49);  // 镍
+        radii.insert("Cu", &1.45);  // 铜
+        radii.insert("Zn", &1.42);  // 锌
+        radii.insert("Ga", &1.36);  // 镓
+        radii.insert("Ge", &1.25);  // 锗
+        radii.insert("As", &1.14);  // 砷
+        radii.insert("Se", &1.03);  // 硒
+        radii.insert("Br", &0.94);  // 溴
+        radii.insert("Kr", &0.88);  // 氪
+
+        // 第5周期
+        radii.insert("Rb", &2.65);  // 铷
+        radii.insert("Sr", &2.19);  // 锶
+        radii.insert( "Y", &2.12);  // 钇
+        radii.insert("Zr", &2.06);  // 锆
+        radii.insert("Nb", &1.98);  // 铌
+        radii.insert("Mo", &1.90);  // 钼
+        radii.insert("Tc", &1.83);  // 锝
+        radii.insert("Ru", &1.78);  // 钌
+        radii.insert("Rh", &1.73);  // 铑
+        radii.insert("Pd", &1.69);  // 钯
+        radii.insert("Ag", &1.65);  // 银
+        radii.insert("Cd", &1.61);  // 镉
+        radii.insert("In", &1.56);  // 铟
+        radii.insert("Sn", &1.45);  // 锡
+        radii.insert("Sb", &1.33);  // 锑
+        radii.insert("Te", &1.23);  // 碲
+        radii.insert( "I", &1.15);  // 碘
+        radii.insert("Xe", &1.08);  // 氙
+
+        // 第6周期
+        radii.insert("Cs", &2.98);  // 铯
+        radii.insert("Ba", &2.53);  // 钡
+        radii.insert("La", &2.50);  // 镧
+        radii.insert("Ce", &2.48);  // 铈
+        radii.insert("Pr", &2.47);  // 镨
+        radii.insert("Nd", &2.45);  // 钕
+        radii.insert("Pm", &2.43);  // 钷
+        radii.insert("Sm", &2.42);  // 钐
+        radii.insert("Eu", &2.40);  // 铕
+        radii.insert("Gd", &2.38);  // 钆
+        radii.insert("Tb", &2.37);  // 铽
+        radii.insert("Dy", &2.35);  // 镝
+        radii.insert("Ho", &2.33);  // 钬
+        radii.insert("Er", &2.32);  // 铒
+        radii.insert("Tm", &2.30);  // 铥
+        radii.insert("Yb", &2.28);  // 镱
+        radii.insert("Lu", &2.27);  // 镥
+        radii.insert("Hf", &2.25);  // 铪
+        radii.insert("Ta", &2.20);  // 钽
+        radii.insert( "W", &2.10);  // 钨
+        radii.insert("Re", &2.05);  // 铼
+        radii.insert("Os", &2.00);  // 锇
+        radii.insert("Ir", &1.97);  // 铱
+        radii.insert("Pt", &1.92);  // 铂
+        radii.insert("Au", &1.87);  // 金
+        radii.insert("Hg", &1.75);  // 汞
+        radii.insert("Tl", &1.70);  // 铊
+        radii.insert("Pb", &1.54);  // 铅
+        radii.insert("Bi", &1.43);  // 铋
+        radii.insert("Po", &1.35);  // 钋
+        radii.insert("At", &1.27);  // 砹
+        radii.insert("Rn", &1.20);  // 氡
+
+        // 第7周期
+        radii.insert("Fr", &3.00);  // 钫
+        radii.insert("Ra", &2.70);  // 镭
+        radii.insert("Ac", &2.60);  // 锕
+        radii.insert("Th", &2.50);  // 钍
+        radii.insert("Pa", &2.40);  // 镤
+        radii.insert( "U",  &2.30); // 铀
+        radii.insert("Np", &2.30);  // 镎
+        radii.insert("Pu", &2.30);  // 钚
+        radii.insert("Am", &2.30);  // 镅
+        radii.insert("Cm", &2.30);  // 锔
+        radii.insert("Bk", &2.30);  // 锫
+        radii.insert("Cf", &2.30);  // 锎
+        radii.insert("Es", &2.30);  // 锿
+        radii.insert("Fm", &2.30);  // 镄
+        radii.insert("Md", &2.30);  // 钔
+        radii.insert("No", &2.30);  // 锘
+        radii.insert("Lr", &2.30);  // 铹
+        radii.insert("Rf", &2.30);  // 卢瑟福
+        radii.insert("Db", &2.30);  // 𨧀
+        radii.insert("Sg", &2.30);  // 𨭎
+        radii.insert("Bh", &2.30);  // 𨨏
+        radii.insert("Hs", &2.30);  // 𨭆
+        radii.insert("Mt", &2.30);  // 䥑
+        radii.insert("Ds", &2.30);  // 鐽
+        radii.insert("Rg", &2.30);  // 錀
+        radii.insert("Cn", &2.30);  // 鎶
+        radii.insert("Nh", &2.30);  // 鉨
+        radii.insert("Fl", &2.30);  // 鈇
+        radii.insert("Mc", &2.30);  // 鏌
+        radii.insert("Lv", &2.30);  // 鉝
+        radii.insert("Ts", &2.30);  // 鿬
+        radii.insert("Og", &2.30);  // 鿫
+
+        // 添加一些常见的同位素和特殊表示
+        radii.insert("D", &0.53);   // 氘
+        radii.insert("T", &0.53);   // 氚
+        
+        radii
+    };
+}
\ No newline at end of file
diff --git a/src/constants/mod.rs b/src/constants/mod.rs
index d92d9922716884d4b8df79150c4c5b7e4c72d044..08be236d13d00149ac0f262e4c34c60eaf9da90e 100644
--- a/src/constants/mod.rs
+++ b/src/constants/mod.rs
@@ -2,127 +2,24 @@ pub mod vsap;
 pub mod c2s;
 pub mod solvent;
 mod cartesian_gto;
+pub mod element;
 
-use std::collections::HashMap;
+// use std::collections::HashMap;
 
-use lazy_static::lazy_static;
+// use lazy_static::lazy_static;
 
 pub use crate::constants::vsap::*;
 pub use crate::constants::c2s::*;
 pub use crate::constants::cartesian_gto::*;
+pub use crate::constants::element::*;
 
-//struct Matrix3x3 {
-//    size
-//
-//}
-
-pub const SPECIES_NAME: [&str; 118] = ["H", "He",
-        "Li","Be","B", "C","N", "O", "F", "Ne",
-        "Na","Mg","Al","Si","P", "S", "Cl","Ar",
-        "K", "Ca","Sc","Ti","V", "Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr",
-        "Rb","Sr","Y", "Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I", "Xe",
-        "Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W", "Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn",
-        "Fr","Ra","Ac","Th","Pa","U", "Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs","Mt","Ds","Rg","Cn","Nh","Fl","Mc","Lv","Ts","Og"
-        ];
-
-// IUPAC 2021 standard atomic weights (Prohaska et al., Pure Appl. Chem., 2022)
-// Table 1 column 7: abridged to 5 significant figures; interval elements use conventional single values
-// 14 interval elements: H, Li, B, C, N, O, Mg, Si, S, Cl, Ar, Br, Tl, Pb
-// Tc, Pm, Po-At, Rn-Ac, Np and beyond: most stable isotope mass
-//   AME2020 / NUBASE2020 via IUPAC 2021 Table 2
-
-pub const MASS_CHARGE: [(f64,f64);118] = [
-    (1.0080,1.0),       (4.0026,2.0),
-
-    (6.94,3.0),         (9.0122,4.0),
-    (10.81,5.0),        (12.011,6.0),       (14.007,7.0),       (15.999,8.0),
-    (18.998,9.0),       (20.180,10.0),
-
-    (22.990,11.0),      (24.305,12.0),
-    (26.982,13.0),      (28.085,14.0),      (30.974,15.0),      (32.06,16.0),
-    (35.45,17.0),       (39.95,18.0),
-
-    (39.098,19.0),      (40.078,20.0),
-    (44.956,21.0),      (47.867,22.0),      (50.942,23.0),      (51.996,24.0),
-    (54.938,25.0),      (55.845,26.0),      (58.933,27.0),      (58.693,28.0),
-    (63.546,29.0),      (65.38,30.0),
-
-    (69.723,31.0),      (72.630,32.0),
-    (74.922,33.0),      (78.971,34.0),      (79.904,35.0),      (83.798,36.0),
-
-    (85.468,37.0),      (87.62,38.0),       (88.906,39.0),      (91.224,40.0),
-    (92.906,41.0),      (95.95,42.0),       (97.90721,43.0),    (101.07,44.0),
-    (102.91,45.0),      (106.42,46.0),      (107.87,47.0),      (112.41,48.0),
-
-    (114.82,49.0),      (118.71,50.0),      (121.76,51.0),      (127.60,52.0),
-    (126.90,53.0),      (131.29,54.0),
-
-    (132.91,55.0),      (137.33,56.0),
-
-    (138.91,57.0),      (140.12,58.0),      (140.91,59.0),      (144.24,60.0),
-    (144.91276,61.0),   (150.36,62.0),      (151.96,63.0),      (157.25,64.0),
-    (158.93,65.0),      (162.50,66.0),      (164.93,67.0),      (167.26,68.0),
-    (168.93,69.0),      (173.05,70.0),
-
-    (174.97,71.0),      (178.49,72.0),
-    (180.95,73.0),      (183.84,74.0),      (186.21,75.0),      (190.23,76.0),
-    (192.22,77.0),      (195.08,78.0),      (196.97,79.0),      (200.59,80.0),
-
-    (204.38,81.0),      (207.2,82.0),       (208.98,83.0),      (208.98243,84.0),
-    (209.98715,85.0),   (222.01758,86.0),
-
-    (223.01973,87.0),   (226.02541,88.0),
-    (227.02775,89.0),   (232.04,90.0),      (231.04,91.0),      (238.03,92.0),
-
-    (237.04817,93.0),   (244.06420,94.0),   (243.06138,95.0),   (247.07035,96.0),
-    (247.07031,97.0),   (251.07959,98.0),   (252.08298,99.0),   (257.09511,100.0),
-    (258.09843,101.0),  (259.10100,102.0),
-
-    (262.10962,103.0),  (267.12179,104.0),
-    (268.12567,105.0),  (271.13378,106.0),  (270.13337,107.0),  (269.13365,108.0),
-    (278.15649,109.0),  (281.16455,110.0),  (282.16934,111.0),  (285.17723,112.0),
 
-    (286.18246,113.0),  (289.19052,114.0),  (289.19397,115.0),  (293.20458,116.0),
-    (294.21084,117.0),  (294.21398,118.0),
-];
-
-pub const ELEM1ST: [&str;2]  = ["H", "He"];
-pub const ELEM2ND: [&str;8]  = ["Li","Be","B", "C", "N", "O", "F", "Ne"];
-pub const ELEM3RD: [&str;8]  = ["Na","Mg","Al","Si","P", "S", "Cl","Ar"];
-pub const ELEM4TH: [&str;18] = ["K", "Ca","Sc","Ti","V", "Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr"];
-pub const ELEM5TH: [&str;18] = ["Rb","Sr","Y", "Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I", "Xe"];
-pub const ELEM6TH: [&str;32] = ["Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W", "Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn"];
-pub const ELEM7TH: [&str;32] = ["Fr","Ra","Ac","Th","Pa","U", "Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs","Mt","Ds","Rg","Cn","Nh","Fl","Mc","Lv","Ts","Og"];
-
-pub const ELEMTMS: [&str;40] = [
-    "Sc","Ti","V", "Cr","Mn","Fe","Co","Ni","Cu","Zn",
-    "Y", "Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd",
-    "La","Hf","Ta","W", "Re","Os","Ir","Pt","Au","Hg",
-    "Ac","Rf","Db","Sg","Bh","Hs","Mt","Ds","Rg","Cn"
-];
-
-
-lazy_static!{
-    pub static ref SPECIES_INFO: HashMap<&'static str, &'static (f64,f64)> = {
-        let mut m = HashMap::new();
-        SPECIES_NAME.iter().zip(MASS_CHARGE.iter()).for_each(|(name,info)| {
-            m.insert(*name,info);
-        });
-        m
-    };
-}
 
 // use for the inverse (sqrt inverse) of the auxiliary coulomb matrix
 pub const AUXBAS_THRESHOLD: f64 = 1.0e-10;
 pub const INVERSE_THRESHOLD: f64 = 1.0e-10;
 pub const SQRT_THRESHOLD: f64 = 1.0e-10;
 
-pub const CM: f64 = 8065.541;
-pub const ANG:f64 = 0.5291772083;
-pub const EV: f64 = 27.2113845;
-pub const FQ: f64 = 1822.888;
-pub const E:  f64 = std::f64::consts::E;
-pub const PI: f64 = std::f64::consts::PI;
 
 pub const E5: f64 = 1.0e5;
 pub const E6: f64 = 1.0e6;
@@ -136,127 +33,7 @@ pub struct DMatrix<const L: usize> {
     indicing:[usize;2],
     data: [f64; L]
 }
-pub const ATOM_CONFIGURATION: [[usize; 4]; 119] = [
-    [ 0, 0, 0, 0],     //  0  GHOST
-    [ 1, 0, 0, 0],     //  1  H
-    [ 2, 0, 0, 0],     //  2  He
-    [ 3, 0, 0, 0],     //  3  Li
-    [ 4, 0, 0, 0],     //  4  Be
-    [ 4, 1, 0, 0],     //  5  B
-    [ 4, 2, 0, 0],     //  6  C
-    [ 4, 3, 0, 0],     //  7  N
-    [ 4, 4, 0, 0],     //  8  O
-    [ 4, 5, 0, 0],     //  9  F
-    [ 4, 6, 0, 0],     // 10  Ne
-    [ 5, 6, 0, 0],     // 11  Na
-    [ 6, 6, 0, 0],     // 12  Mg
-    [ 6, 7, 0, 0],     // 13  Al
-    [ 6, 8, 0, 0],     // 14  Si
-    [ 6, 9, 0, 0],     // 15  P
-    [ 6,10, 0, 0],     // 16  S
-    [ 6,11, 0, 0],     // 17  Cl
-    [ 6,12, 0, 0],     // 18  Ar
-    [ 7,12, 0, 0],     // 19  K
-    [ 8,12, 0, 0],     // 20  Ca
-    [ 8,12, 1, 0],     // 21  Sc
-    [ 8,12, 2, 0],     // 22  Ti
-    [ 8,12, 3, 0],     // 23  V
-    [ 7,12, 5, 0],     // 24  Cr
-    [ 8,12, 5, 0],     // 25  Mn
-    [ 8,12, 6, 0],     // 26  Fe
-    [ 8,12, 7, 0],     // 27  Co
-    [ 8,12, 8, 0],     // 28  Ni
-    [ 7,12,10, 0],     // 29  Cu
-    [ 8,12,10, 0],     // 30  Zn
-    [ 8,13,10, 0],     // 31  Ga
-    [ 8,14,10, 0],     // 32  Ge
-    [ 8,15,10, 0],     // 33  As
-    [ 8,16,10, 0],     // 34  Se
-    [ 8,17,10, 0],     // 35  Br
-    [ 8,18,10, 0],     // 36  Kr
-    [ 9,18,10, 0],     // 37  Rb
-    [10,18,10, 0],     // 38  Sr
-    [10,18,11, 0],     // 39  Y
-    [10,18,12, 0],     // 40  Zr
-    [ 9,18,14, 0],     // 41  Nb
-    [ 9,18,15, 0],     // 42  Mo
-    [10,18,15, 0],     // 43  Tc
-    [ 9,18,17, 0],     // 44  Ru
-    [ 9,18,18, 0],     // 45  Rh
-    [ 8,18,20, 0],     // 46  Pd
-    [ 9,18,20, 0],     // 47  Ag
-    [10,18,20, 0],     // 48  Cd
-    [10,19,20, 0],     // 49  In
-    [10,20,20, 0],     // 50  Sn
-    [10,21,20, 0],     // 51  Sb
-    [10,22,20, 0],     // 52  Te
-    [10,23,20, 0],     // 53  I
-    [10,24,20, 0],     // 54  Xe
-    [11,24,20, 0],     // 55  Cs
-    [12,24,20, 0],     // 56  Ba
-    [12,24,21, 0],     // 57  La
-    [12,24,21, 1],     // 58  Ce
-    [12,24,20, 3],     // 59  Pr
-    [12,24,20, 4],     // 60  Nd
-    [12,24,20, 5],     // 61  Pm
-    [12,24,20, 6],     // 62  Sm
-    [12,24,20, 7],     // 63  Eu
-    [12,24,21, 7],     // 64  Gd
-    [12,24,21, 8],     // 65  Tb
-    [12,24,20,10],     // 66  Dy
-    [12,24,20,11],     // 67  Ho
-    [12,24,20,12],     // 68  Er
-    [12,24,20,13],     // 69  Tm
-    [12,24,20,14],     // 70  Yb
-    [12,24,21,14],     // 71  Lu
-    [12,24,22,14],     // 72  Hf
-    [12,24,23,14],     // 73  Ta
-    [12,24,24,14],     // 74  W
-    [12,24,25,14],     // 75  Re
-    [12,24,26,14],     // 76  Os
-    [12,24,27,14],     // 77  Ir
-    [11,24,29,14],     // 78  Pt
-    [11,24,30,14],     // 79  Au
-    [12,24,30,14],     // 80  Hg
-    [12,25,30,14],     // 81  Tl
-    [12,26,30,14],     // 82  Pb
-    [12,27,30,14],     // 83  Bi
-    [12,28,30,14],     // 84  Po
-    [12,29,30,14],     // 85  At
-    [12,30,30,14],     // 86  Rn
-    [13,30,30,14],     // 87  Fr
-    [14,30,30,14],     // 88  Ra
-    [14,30,31,14],     // 89  Ac
-    [14,30,32,14],     // 90  Th
-    [14,30,31,16],     // 91  Pa
-    [14,30,31,17],     // 92  U
-    [14,30,31,18],     // 93  Np
-    [14,30,30,20],     // 94  Pu
-    [14,30,30,21],     // 95  Am
-    [14,30,31,21],     // 96  Cm
-    [14,30,31,22],     // 97  Bk
-    [14,30,30,24],     // 98  Cf
-    [14,30,30,25],     // 99  Es
-    [14,30,30,26],     //100  Fm
-    [14,30,30,27],     //101  Md
-    [14,30,30,28],     //102  No
-    [14,30,31,28],     //103  Lr
-    [14,30,32,28],     //104  Rf
-    [14,30,33,28],     //105  Db
-    [14,30,34,28],     //106  Sg
-    [14,30,35,28],     //107  Bh
-    [14,30,36,28],     //108  Hs
-    [14,30,37,28],     //109  Mt
-    [14,30,38,28],     //110  Ds
-    [14,30,39,28],     //111  Rg
-    [14,30,40,28],     //112  Cn
-    [14,31,40,28],     //113  Nh
-    [14,32,40,28],     //114  Fl
-    [14,33,40,28],     //115  Mc
-    [14,34,40,28],     //116  Lv
-    [14,35,40,28],     //117  Ts
-    [14,36,40,28],     //118  Og
-];
+
 
 // =========== libcint ===================================
 // for the bas index - libcint
@@ -300,184 +77,32 @@ pub const NUC_FRAC_CHARGE: i32 = 3;
 
 pub const ENV_PRT_START: usize = 20;
 
-
-// SAD and ECP configuration
-pub const S_SHELL: [f64; 1] = [2.0];
-pub const P_SHELL: [f64; 3] = [2.0, 2.0, 2.0];
-pub const D_SHELL: [f64; 5] = [2.0, 2.0, 2.0, 2.0, 2.0];
-pub const F_SHELL: [f64; 7] = [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0];
-pub const XE_SHELL: [f64; 27] =
-//   1s   2s   2p   2p   2p   3s   3p   3p   3p   4s   3d   3d   3d   3d   3d   4p   4p   4p   5s   4d   4d   4d   4d   4d   5p   5p   5p  
-    [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0];
-pub const KR_SHELL: [f64; 18] =
-//   1s   2s   2p   2p   2p   3s   3p   3p   3p   4s   3d   3d   3d   3d   3d   4p   4p   4p  
-    [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0];
-
-//                                    1s   2s   2p   3s   3p   4s   3d    4p   5s   4d    5p   4f    5d    6s   6p 
-pub const NELE_IN_SHELLS: [f64; 15] = [2.0, 2.0, 6.0, 2.0, 6.0, 2.0, 10.0, 6.0, 2.0, 10.0, 6.0, 14.0, 10.0, 2.0, 6.0];
-
+// math, physics
+// NOTE: these constants come from several different CODATA releases (see per-line sources).
+pub const EV: f64 = 27.2113845;                 // Hartree energy in eV, CODATA 2002
+pub const HARTREE2KCAL: f64 = 627.509451;       // Hartree -> kcal/mol; CODATA 2002 (Hartree energy in J / (kcal * Avogadro))
+pub const HARTREE2WAVENUMBER: f64 = 219474.63;  // Hartree -> cm^-1 (hartree-inverse meter relationship / 100); CODATA (~2.194746313705e7 m^-1)
+pub const FQ: f64 = 1822.8884861920776;         // u/m_e ratio (= 1 / electron mass in u 5.48579909070e-4), CODATA 2014
+pub const E:  f64 = std::f64::consts::E;         // math constant (std); year-insensitive
+pub const PI: f64 = std::f64::consts::PI;        // math constant (std); year-insensitive
 //
-pub const LIGHT_SPEED: f64 = 137.03599967994;   // http://physics.nist.gov/cgi-bin/cuu/Value?alph
+pub const LIGHT_SPEED: f64 = 137.03599967994;   // inverse fine-structure constant, CODATA 2006; http://physics.nist.gov/cgi-bin/cuu/Value?alph
 // BOHR = .529 177 210 92(17) e-10m  // http://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0
+// source: CODATA 2010 (https://physics.nist.gov/cuu/Constants/ArchiveASCII/allascii_2010.txt)
 pub const BOHR: f64 = 0.52917721092;  // Angstroms
 pub const BOHR_SI: f64 = BOHR * 1e-10;
 
-pub const G_ELECTRON: f64 = 2.00231930436182;  // http://physics.nist.gov/cgi-bin/cuu/Value?gem
-pub const E_MASS: f64 = 9.10938356e-31;         // kg https://physics.nist.gov/cgi-bin/cuu/Value?me
-pub const AVOGADRO: f64 = 6.022140857e23;       // https://physics.nist.gov/cgi-bin/cuu/Value?na
-pub const PLANCK: f64 = 6.626070040e-34;        // J*s http://physics.nist.gov/cgi-bin/cuu/Value?h
-pub const BOLTZMANN: f64 = 1.380649e-23;        // J/K https://physics.nist.gov/cgi-bin/cuu/Value?k
-pub const CLIGHT_CMS: f64 = 2.99792458e10;      // speed of light, cm/s
-pub const R_GAS: f64 = BOLTZMANN * AVOGADRO;    // J/(mol*K) ideal gas constant
-pub const E_CHARGE: f64 = 1.6021766208e-19;
-pub const DEBYE:f64 = 3.335641e-30;            // C*m = 1e-18/LIGHT_SPEED_SI https://cccbdb.nist.gov/debye.asp
-pub const AU2DEBYE:f64 = E_CHARGE * BOHR*1e-10 / DEBYE; // 2.541746
+pub const G_ELECTRON: f64 = 2.00231930436182;  // CODATA 2014; http://physics.nist.gov/cgi-bin/cuu/Value?gem
+pub const E_MASS: f64 = 9.10938356e-31;         // kg, CODATA 2014; https://physics.nist.gov/cgi-bin/cuu/Value?me
+pub const AVOGADRO: f64 = 6.022140857e23;       // CODATA 2014; https://physics.nist.gov/cgi-bin/cuu/Value?na
+pub const PLANCK: f64 = 6.626070040e-34;        // J*s, CODATA 2014; http://physics.nist.gov/cgi-bin/cuu/Value?h
+pub const BOLTZMANN: f64 = 1.380649e-23;        // J/K, CODATA 2018 (exact, SI 2019); https://physics.nist.gov/cgi-bin/cuu/Value?k
+pub const CLIGHT_CMS: f64 = 2.99792458e10;      // speed of light, cm/s; exact by SI definition (year-insensitive)
+pub const R_GAS: f64 = BOLTZMANN * AVOGADRO;    // J/(mol*K) ideal gas constant; derived
+pub const E_CHARGE: f64 = 1.6021766208e-19;    // C, CODATA 2014
+pub const DEBYE:f64 = 3.335641e-30;            // C*m = 1e-18/LIGHT_SPEED_SI; defined via c, year-insensitive; https://cccbdb.nist.gov/debye.asp
+pub const AU2DEBYE:f64 = E_CHARGE * BOHR*1e-10 / DEBYE; // 2.541746; derived
 
 
 pub const MPI_CHUNK:usize = 134217728; // around 1 GB
 
-lazy_static!{
-    pub static ref ATOMIC_RADII: HashMap<&'static str, &'static f64> = {
-    /// 获取元素周期表中所有元素的原子半径（单位：Å）
-    /// 数据来源：Clementi-Raimondi半径、共价半径和范德华半径的综合
-        let mut radii = HashMap::new();
-        
-        // 第1周期
-        radii.insert( "H", &0.53);  // 氢
-        radii.insert("He", &0.31);  // 氦
-
-        // 第2周期
-        radii.insert("Li", &1.67);  // 锂
-        radii.insert("Be", &1.12);  // 铍
-        radii.insert( "B", &0.87);  // 硼
-        radii.insert( "C", &0.67);  // 碳
-        radii.insert( "N", &0.56);  // 氮
-        radii.insert( "O", &0.48);  // 氧
-        radii.insert( "F", &0.42);  // 氟
-        radii.insert("Ne", &0.38);  // 氖
-
-        // 第3周期
-        radii.insert("Na", &1.90);  // 钠
-        radii.insert("Mg", &1.45);  // 镁
-        radii.insert("Al", &1.18);  // 铝
-        radii.insert("Si", &1.11);  // 硅
-        radii.insert( "P", &0.98);  // 磷
-        radii.insert( "S", &0.88);  // 硫
-        radii.insert("Cl", &0.79);  // 氯
-        radii.insert("Ar", &0.71);  // 氩
-
-        // 第4周期
-        radii.insert( "K", &2.43);  // 钾
-        radii.insert("Ca", &1.94);  // 钙
-        radii.insert("Sc", &1.84);  // 钪
-        radii.insert("Ti", &1.76);  // 钛
-        radii.insert( "V", &1.71);  // 钒
-        radii.insert("Cr", &1.66);  // 铬
-        radii.insert("Mn", &1.61);  // 锰
-        radii.insert("Fe", &1.56);  // 铁
-        radii.insert("Co", &1.52);  // 钴
-        radii.insert("Ni", &1.49);  // 镍
-        radii.insert("Cu", &1.45);  // 铜
-        radii.insert("Zn", &1.42);  // 锌
-        radii.insert("Ga", &1.36);  // 镓
-        radii.insert("Ge", &1.25);  // 锗
-        radii.insert("As", &1.14);  // 砷
-        radii.insert("Se", &1.03);  // 硒
-        radii.insert("Br", &0.94);  // 溴
-        radii.insert("Kr", &0.88);  // 氪
-
-        // 第5周期
-        radii.insert("Rb", &2.65);  // 铷
-        radii.insert("Sr", &2.19);  // 锶
-        radii.insert( "Y", &2.12);  // 钇
-        radii.insert("Zr", &2.06);  // 锆
-        radii.insert("Nb", &1.98);  // 铌
-        radii.insert("Mo", &1.90);  // 钼
-        radii.insert("Tc", &1.83);  // 锝
-        radii.insert("Ru", &1.78);  // 钌
-        radii.insert("Rh", &1.73);  // 铑
-        radii.insert("Pd", &1.69);  // 钯
-        radii.insert("Ag", &1.65);  // 银
-        radii.insert("Cd", &1.61);  // 镉
-        radii.insert("In", &1.56);  // 铟
-        radii.insert("Sn", &1.45);  // 锡
-        radii.insert("Sb", &1.33);  // 锑
-        radii.insert("Te", &1.23);  // 碲
-        radii.insert( "I", &1.15);  // 碘
-        radii.insert("Xe", &1.08);  // 氙
-
-        // 第6周期
-        radii.insert("Cs", &2.98);  // 铯
-        radii.insert("Ba", &2.53);  // 钡
-        radii.insert("La", &2.50);  // 镧
-        radii.insert("Ce", &2.48);  // 铈
-        radii.insert("Pr", &2.47);  // 镨
-        radii.insert("Nd", &2.45);  // 钕
-        radii.insert("Pm", &2.43);  // 钷
-        radii.insert("Sm", &2.42);  // 钐
-        radii.insert("Eu", &2.40);  // 铕
-        radii.insert("Gd", &2.38);  // 钆
-        radii.insert("Tb", &2.37);  // 铽
-        radii.insert("Dy", &2.35);  // 镝
-        radii.insert("Ho", &2.33);  // 钬
-        radii.insert("Er", &2.32);  // 铒
-        radii.insert("Tm", &2.30);  // 铥
-        radii.insert("Yb", &2.28);  // 镱
-        radii.insert("Lu", &2.27);  // 镥
-        radii.insert("Hf", &2.25);  // 铪
-        radii.insert("Ta", &2.20);  // 钽
-        radii.insert( "W", &2.10);  // 钨
-        radii.insert("Re", &2.05);  // 铼
-        radii.insert("Os", &2.00);  // 锇
-        radii.insert("Ir", &1.97);  // 铱
-        radii.insert("Pt", &1.92);  // 铂
-        radii.insert("Au", &1.87);  // 金
-        radii.insert("Hg", &1.75);  // 汞
-        radii.insert("Tl", &1.70);  // 铊
-        radii.insert("Pb", &1.54);  // 铅
-        radii.insert("Bi", &1.43);  // 铋
-        radii.insert("Po", &1.35);  // 钋
-        radii.insert("At", &1.27);  // 砹
-        radii.insert("Rn", &1.20);  // 氡
-
-        // 第7周期
-        radii.insert("Fr", &3.00);  // 钫
-        radii.insert("Ra", &2.70);  // 镭
-        radii.insert("Ac", &2.60);  // 锕
-        radii.insert("Th", &2.50);  // 钍
-        radii.insert("Pa", &2.40);  // 镤
-        radii.insert( "U",  &2.30); // 铀
-        radii.insert("Np", &2.30);  // 镎
-        radii.insert("Pu", &2.30);  // 钚
-        radii.insert("Am", &2.30);  // 镅
-        radii.insert("Cm", &2.30);  // 锔
-        radii.insert("Bk", &2.30);  // 锫
-        radii.insert("Cf", &2.30);  // 锎
-        radii.insert("Es", &2.30);  // 锿
-        radii.insert("Fm", &2.30);  // 镄
-        radii.insert("Md", &2.30);  // 钔
-        radii.insert("No", &2.30);  // 锘
-        radii.insert("Lr", &2.30);  // 铹
-        radii.insert("Rf", &2.30);  // 卢瑟福
-        radii.insert("Db", &2.30);  // 𨧀
-        radii.insert("Sg", &2.30);  // 𨭎
-        radii.insert("Bh", &2.30);  // 𨨏
-        radii.insert("Hs", &2.30);  // 𨭆
-        radii.insert("Mt", &2.30);  // 䥑
-        radii.insert("Ds", &2.30);  // 鐽
-        radii.insert("Rg", &2.30);  // 錀
-        radii.insert("Cn", &2.30);  // 鎶
-        radii.insert("Nh", &2.30);  // 鉨
-        radii.insert("Fl", &2.30);  // 鈇
-        radii.insert("Mc", &2.30);  // 鏌
-        radii.insert("Lv", &2.30);  // 鉝
-        radii.insert("Ts", &2.30);  // 鿬
-        radii.insert("Og", &2.30);  // 鿫
-
-        // 添加一些常见的同位素和特殊表示
-        radii.insert("D", &0.53);   // 氘
-        radii.insert("T", &0.53);   // 氚
-        
-        radii
-    };
-}
\ No newline at end of file
diff --git a/src/dft/gen_grids/atom.rs b/src/dft/gen_grids/atom.rs
index eb7b74435d881825ca4ddc7e323c33e07ab4724e..ac294f9f04395758c6973b1373e42814062ab9c1 100644
--- a/src/dft/gen_grids/atom.rs
+++ b/src/dft/gen_grids/atom.rs
@@ -11,6 +11,7 @@ use super::lebedev;
 use super::prune::*;
 use super::radial;
 use super::parameters::LEBEDEV_NGRID;
+use crate::constants::BOHR;
 
 /// Generate atomic grids based on given arguments. <br>
 /// Download basis set information from www.basissetexchange.org to determine alpha_max and alpha_min.
@@ -121,7 +122,7 @@ pub fn atom_grid(
     //println!("radial num = {}", default_radial_num(proton_charges[center_index] as usize));
     // factors match DIRAC code
     //println!("rs = {:?}, w = {:?}", rs, weights_radial);
-    let rb = bragg::get_bragg_angstrom(proton_charges[center_index]) / (5.0 * 0.529177249);
+    let rb = bragg::get_bragg_angstrom(proton_charges[center_index]) / (5.0 * BOHR);
 
     let mut coordinates = Vec::new();
     let mut weights = Vec::new();
diff --git a/src/dft/gen_grids/atom_old.rs b/src/dft/gen_grids/atom_old.rs
index 87a6ef95aa3000a338c36dcb35dd134cd0b83d46..24b68f05ecf60fea3aa486f2b8c18bd66a55d9f8 100644
--- a/src/dft/gen_grids/atom_old.rs
+++ b/src/dft/gen_grids/atom_old.rs
@@ -91,7 +91,7 @@ pub fn atom_grid(
     //println!("radial num = {}", default_radial_num(proton_charges[center_index] as usize));
     // factors match DIRAC code
     //println!("rs = {:?}, w = {:?}", rs, weights_radial);
-    let rb = bragg::get_bragg_angstrom(proton_charges[center_index]) / (5.0 * 0.529177249);
+    let rb = bragg::get_bragg_angstrom(proton_charges[center_index]) / (5.0 * crate::constants::BOHR);
 
     let mut coordinates = Vec::new();
     let mut weights = Vec::new();
diff --git a/src/dft/gen_grids/parameters.rs b/src/dft/gen_grids/parameters.rs
index 265a7dfdc94e980fc3e25fb6a127ff80e5687aa1..897af95671e6e9ee8441f4fcaa41b5a706bce382 100644
--- a/src/dft/gen_grids/parameters.rs
+++ b/src/dft/gen_grids/parameters.rs
@@ -7,9 +7,6 @@ pub static SG1RADII: [f64; 19] = [0.0, 1.0000, 0.5882,   //H, He
 3.0769, 2.0513, 1.5385, 1.2308, 1.0256, 0.8791, 0.7692, 0.6838,  //2nd Period
 4.0909, 3.1579, 2.5714, 2.1687, 1.8750, 1.6514, 1.4754, 1.3333];  //3rd Period
 
-/// Bohr radius.
-pub const BOHR: f64 = 0.52917721092;
-
 /// Bragg radii for elements.
 pub static BRAGG0: [f64; 131] = [0.0,  // Ghost atom
 0.35,                                     1.40,             // 1s
diff --git a/src/dft/gen_grids/prune.rs b/src/dft/gen_grids/prune.rs
index 89004b4358cd37b0cfbd60188196bce0d5b2b48e..34baac0ba1880cb5c0e7caa166b40a3038aa2f7c 100644
--- a/src/dft/gen_grids/prune.rs
+++ b/src/dft/gen_grids/prune.rs
@@ -7,9 +7,10 @@
 
 use num_traits::{ToPrimitive};
 use tensors::MatrixFull;
-//use super::parameters::{SG1RADII, BOHR, BRAGG0, LEBEDEV_NGRID};
+//use super::parameters::{SG1RADII, BRAGG0, LEBEDEV_NGRID};
 use crate::{dft::Grids, utilities::balancing};
-use super::{parameters::{SG1RADII, BOHR, BRAGG0, LEBEDEV_NGRID}, atom::default_angular_num};
+use crate::constants::BOHR;
+use super::{parameters::{SG1RADII, BRAGG0, LEBEDEV_NGRID}, atom::default_angular_num};
 
 /// Standard Grid 1 according to _P. M. W. Gill, B. G. Johnson, J. A. Pople. Chemical Physics Letters 209, 506-512 (1993)_.<br>
 /// Reference can be found [here](https://doi.org/10.1016/0009-2614(93)80125-9).
diff --git a/src/dft/gen_grids/radial.rs b/src/dft/gen_grids/radial.rs
index 4a6e57b7bb735d82a27ff2d9347e85a2517c7a32..25f1205cd98131461e3c9212fc2980d71e6899d8 100644
--- a/src/dft/gen_grids/radial.rs
+++ b/src/dft/gen_grids/radial.rs
@@ -28,7 +28,8 @@ use std::f64::consts::PI;
 use super::bragg;
 use super::bse;
 use super::parameters;
-use super::parameters::{BOHR, BRAGG0};
+use super::parameters::BRAGG0;
+use crate::constants::BOHR;
 use statrs::function::gamma;
 
 /// Krack-Koster radial grid according to _M. Krack, A. M. Köster. The Journal of Chemical Physics 108, 3226-3234 (1998)_, eqs. 9-13.<br>
diff --git a/src/external_field/num_dipole.rs b/src/external_field/num_dipole.rs
index 244f4a83e3a7ba8e2a2611acc02098abebb936e2..9ebb6e7698d7dbace1d5935e96b56700c95d2ed3 100644
--- a/src/external_field/num_dipole.rs
+++ b/src/external_field/num_dipole.rs
@@ -71,8 +71,8 @@ mod debug {
      print_level =          2
      job_type =             "numerical dipole"
      xc =                   "hf"
-     basis_path =           "basis-set-pool/def2-TZVP"
-     auxbas_path =          "basis-set-pool/def2-SVP-JKFIT"
+     basis_path =           "def2-TZVP"
+     auxbas_path =          "def2-universal-JKFIT"
      guessfile =            "none"
      charge =               0.0
      spin =                 1.0
@@ -102,8 +102,8 @@ mod debug {
 [ctrl]
      print_level =          2
      xc =                   "hf"
-     basis_path =           "basis-set-pool/def2-TZVP"
-     auxbas_path =          "basis-set-pool/def2-SV(P)-JKFIT"
+     basis_path =           "def2-TZVP"
+     auxbas_path =          "def2-SV(P)-JKFIT"
      guessfile =            "none"
      charge =               0.0
      spin =                 1.0
diff --git a/src/fileop/chkfile.rs b/src/fileop/chkfile.rs
index aea7d1929d07220cc77c4c08a369a0487dd3bd1e..d06f274d8be786459ce22cf0ea26895f5c1fcb99 100644
--- a/src/fileop/chkfile.rs
+++ b/src/fileop/chkfile.rs
@@ -7,6 +7,7 @@ use crate::scf_io::{SCF, SCFType};
 use crate::geom_io::get_mass_charge;
 use crate::basis_io::Basis4Elem;
 use rest_libcint::CintType;
+use crate::constants::BOHR;
 
 pub fn write_scf_attribute<T>(group: &hdf5::Group, dataset_name: &str, value: &[T]) 
 where 
@@ -183,7 +184,7 @@ pub fn save_overlap(scf_data: &SCF) {
 }
 
 pub fn save_geometry(scf_data: &SCF) {
-    let ang = crate::constants::ANG;
+    let ang = BOHR;
     let chkfile= &scf_data.mol.ctrl.chkfile;
     let path = Path::new(chkfile);
     let file = if path.exists() {
diff --git a/src/geom_io/mod.rs b/src/geom_io/mod.rs
index 7143f60fa337829cfd88b098f21983d82c6fcfee..6c280d5da4787f0dbd9c6066829f44bad81a8eac 100644
--- a/src/geom_io/mod.rs
+++ b/src/geom_io/mod.rs
@@ -14,7 +14,7 @@ use serde_json::Value;
 //use tensors::Tensors;
 
 use crate::basis_io::Basis4Elem;
-use crate::constants::{ANG, ATOMIC_RADII, MASS_CHARGE, SPECIES_INFO, SPECIES_NAME};
+use crate::constants::{BOHR, ATOMIC_RADII, MASS_CHARGE, SPECIES_INFO, SPECIES_NAME};
 use crate::external_field::ExtField;
 mod pyrest_geom_io;
 
@@ -343,7 +343,7 @@ impl GeomCell {
     }
     pub fn geom_update(&mut self, new_position:&[f64], unit: GeomUnit) {
         let factor = match unit {
-            GeomUnit::Angstrom => ANG,
+            GeomUnit::Angstrom => BOHR,
             GeomUnit::Bohr => 1.0,
         };
         if self.position.data.len() != new_position.len() {
@@ -423,7 +423,7 @@ impl GeomCell {
         let mut tmp_pos_tensor = MatrixFull::from_vec(tmp_size, tmp_pos).unwrap();
         if let GeomUnit::Angstrom = unit {
             // To store the geometry position in "Bohr" according to the convention of quantum chemistry. 
-            tmp_pos_tensor.self_multiple(ANG.powf(-1.0));
+            tmp_pos_tensor.self_multiple(BOHR.powf(-1.0));
         };
         Ok((tmp_ele, tmp_fix, tmp_pos_tensor, tmp_nfree))
     }
@@ -457,7 +457,7 @@ impl GeomCell {
         let mut tmp_lat = unsafe{MatrixFull::from_vec_unchecked([3,3],tmp_vec)};
         if let GeomUnit::Angstrom = unit {
             // To store the lattice vector in "Bohr" according to the convention of quantum chemistry. 
-            tmp_lat.self_multiple(ANG.powf(-1.0));
+            tmp_lat.self_multiple(BOHR.powf(-1.0));
         };
         Ok(tmp_lat)
         //if frac_bool {
@@ -534,7 +534,7 @@ impl GeomCell {
         let mut tmp_pos_tensor = MatrixFull::from_vec(tmp_size, tmp_pos).unwrap();
         if let GeomUnit::Angstrom = unit {
             // To store the geometry position in "Bohr" according to the convention of quantum chemistry. 
-            tmp_pos_tensor.self_multiple(ANG.powf(-1.0));
+            tmp_pos_tensor.self_multiple(BOHR.powf(-1.0));
         };
         Ok((tmp_ele, tmp_fix, tmp_pos_tensor, tmp_nfree))
     }
@@ -606,7 +606,7 @@ impl GeomCell {
         let mut tmp_pos_tensor = MatrixFull::from_vec(tmp_size, tmp_pos).unwrap();
         if let GeomUnit::Angstrom = unit {
             // To store the geometry position in "Bohr" according to the convention of quantum chemistry. 
-            tmp_pos_tensor.self_multiple(ANG.powf(-1.0));
+            tmp_pos_tensor.self_multiple(BOHR.powf(-1.0));
         };
         Ok((tmp_ele, tmp_fix, tmp_pos_tensor, tmp_nfree))
 
@@ -700,7 +700,7 @@ impl GeomCell {
             let mut tmp_pos_tensor = MatrixFull::from_vec(tmp_size, tmp_bs_pos).unwrap();
             if let GeomUnit::Angstrom = unit {
                 // To store the geometry position in "Bohr" according to the convention of quantum chemistry. 
-                tmp_pos_tensor.self_multiple(ANG.powf(-1.0));
+                tmp_pos_tensor.self_multiple(BOHR.powf(-1.0));
             };
             Some((tmp_bs_ele, tmp_pos_tensor))
         };
@@ -712,7 +712,7 @@ impl GeomCell {
             let mut tmp_pos_tensor = MatrixFull::from_vec(tmp_size, tmp_pc_pos).unwrap();
             if let GeomUnit::Angstrom = unit {
                 // To store the geometry position in "Bohr" according to the convention of quantum chemistry. 
-                tmp_pos_tensor.self_multiple(ANG.powf(-1.0));
+                tmp_pos_tensor.self_multiple(BOHR.powf(-1.0));
             };
             Some((tmp_pc_chg, tmp_pos_tensor))
         };
@@ -724,7 +724,7 @@ impl GeomCell {
             let mut tmp_pos_tensor = MatrixFull::from_vec(tmp_size, tmp_ep_pos).unwrap();
             if let GeomUnit::Angstrom = unit {
                 // To store the geometry position in "Bohr" according to the convention of quantum chemistry. 
-                tmp_pos_tensor.self_multiple(ANG.powf(-1.0));
+                tmp_pos_tensor.self_multiple(BOHR.powf(-1.0));
             };
             Some((tmp_ep_pth, tmp_pos_tensor))
         };
@@ -738,7 +738,7 @@ impl GeomCell {
     }
 
     pub fn to_xyz(&self, filename: String) {
-        let ang = crate::constants::ANG;
+        let ang = BOHR;
         let mut input = fs::File::create(&filename).unwrap();
         write!(input, "{}\n\n", self.elem.len());
         self.position.iter_columns_full().zip(self.elem.iter()).for_each(|(pos, elem)| {
@@ -747,7 +747,7 @@ impl GeomCell {
     }
 
     pub fn formated_geometry(&self) -> String {
-        let ang = crate::constants::ANG;
+        let ang = BOHR;
         let mut input = String::new();
         //write!(input, "{}\n\n", self.elem.len());
         self.position.iter_columns_full().zip(self.elem.iter()).for_each(|(pos, elem)| {
@@ -808,7 +808,7 @@ impl GeomCell {
         let rj = self.position.iter_column(j);
         let mut dd = ri.zip(rj)
             .fold(0.0,|acc,(ri,rj)| acc + (ri-rj).powf(2.0)).sqrt();
-        dd *= crate::constants::ANG.powf(-1.0);
+        dd *= BOHR.powf(-1.0);
         dd
     }
 
diff --git a/src/grad/mod.rs b/src/grad/mod.rs
index d341dc17b462ac9b70eff9fdce8f65f44ae2e18a..787c25ffd731807d2af3c11ee63a191397f2cd1b 100644
--- a/src/grad/mod.rs
+++ b/src/grad/mod.rs
@@ -8,7 +8,7 @@ pub mod traits;
 use std::io::{self, Write};
 
 use crate::main_driver::{collect_total_energy, performance_essential_calculations};
-use crate::{constants::{ANG, EV}, scf_io::{initialize_scf, scf_without_build, SCF}, utilities};
+use crate::{constants::{BOHR, EV}, scf_io::{initialize_scf, scf_without_build, SCF}, utilities};
 use crate::mpi_io::{MPIData, MPIOperator};
 use tensors::MatrixFull;
 
@@ -132,7 +132,7 @@ pub fn formated_force(force: &MatrixFull<f64>, elem: &Vec<String>) -> String {
 pub fn formated_force_ev(force: &MatrixFull<f64>, elem: &Vec<String>) -> String {
     let mut output = String::new();
     force.iter_columns_full().zip(elem.iter()).for_each(|(force, elem)| {
-        output  = format!("{}{:3}{:16.8}{:16.8}{:16.8}\n", output, elem, force[0]*EV/ANG,force[1]*EV/ANG,force[2]*EV/ANG);
+        output  = format!("{}{:3}{:16.8}{:16.8}{:16.8}\n", output, elem, force[0]*EV/BOHR,force[1]*EV/BOHR,force[2]*EV/BOHR);
     });
 
     output
diff --git a/src/hessian/rhf.rs b/src/hessian/rhf.rs
index b2cb7f983248807aacc8b951b1f66d9ce93a0d9b..965b369dc6f95356d2042a010fe66f6fcc85f0b8 100644
--- a/src/hessian/rhf.rs
+++ b/src/hessian/rhf.rs
@@ -10,6 +10,7 @@ use tensors::matrix_blas_lapack::_power_rayon_for_symmetric_matrix;
 use tensors::MatrixFull;
 
 use crate::utilities::rstsr_util::*;
+use crate::constants::{FQ, HARTREE2WAVENUMBER};
 
 /// Routing for each optimizable G-term in calc_ej_ek().
 #[derive(Clone, Copy, PartialEq, Debug)]
@@ -4247,7 +4248,7 @@ pub fn compute_frequencies_from_hessian(
     // Convert eigenvalues to frequencies in cm⁻¹
     // ω² = λ (in a.u.), ω = √λ (a.u.), ν = ω/(2π) (a.u.)
     // 1 Hartree = 219474.63 cm⁻¹
-    let conv = 219474.63 / 1822.8885f64.sqrt();  // = 5140.49
+    let conv = HARTREE2WAVENUMBER / FQ.sqrt();  // = 5140.49
     let mut freqs = vec![0.0; n3];
     for i in 0..n3 {
         let lambda = eigvals[i];
diff --git a/src/lib_rint/mod.rs b/src/lib_rint/mod.rs
index a3394cc8f838b7a42fa7286bff6dbc0dfefd20a5..2ea4c7f41d5cc510485574789caf616c1c4b8993 100644
--- a/src/lib_rint/mod.rs
+++ b/src/lib_rint/mod.rs
@@ -1,6 +1,6 @@
 use crate::basis_io::etb::{etb_gen_for_atom_list, get_etb_elem};
 use crate::basis_io::Basis4Elem;
-use crate::constants::{AUXBAS_THRESHOLD, SQRT_THRESHOLD};
+use crate::constants::{AUXBAS_THRESHOLD, BOHR, SQRT_THRESHOLD};
 use crate::geom_io::GeomCell;
 use crate::molecule_io::Molecule;
 use crate::scf_io::{
@@ -11728,7 +11728,7 @@ position = [
             scf_without_build(&mut scf_data, &None);
             let cov_a = fragment_dipole_fluctuation_tensor(&scf_data, &[0]);
             let cov_b = fragment_dipole_fluctuation_tensor(&scf_data, &[1]);
-            let distance_bohr = distance_ang / crate::constants::ANG;
+            let distance_bohr = distance_ang / BOHR;
             let x = fragment_connected_dipole_x(cov_a, cov_b, distance_bohr);
             println!(
                 "occ-closure dipole X_disp(He2): R={distance_ang:.3} Ang ({distance_bohr:.6} bohr) X={x:.16e} logR={:.8} log|X|={:.8}",
diff --git a/src/main_driver.rs b/src/main_driver.rs
index 3f4b4edb29637a0a5ce034056821139e817b9702..b64a868a4835628cede2458bff77841b93641ad9 100644
--- a/src/main_driver.rs
+++ b/src/main_driver.rs
@@ -13,7 +13,7 @@ use num_traits::Pow;
 use pyo3::prelude::*;
 //use autocxx::prelude::*;
 use crate::ctrl_io::JobType;
-use crate::constants::{ANG, AU2DEBYE};
+use crate::constants::{BOHR, AU2DEBYE};
 use crate::scf_io::{scf_without_build, SCFType, SCF};
 use tensors::{MathMatrix, MatrixFull};
 use tensors::matrix_blas_lapack::_dsyevd;
@@ -180,7 +180,7 @@ pub fn main_driver() -> anyhow::Result<()> {
                 if scf_data.mol.ctrl.print_level>0 {
                     println!("Geometry optimization invoked");
                 }
-                let displace = scf_data.mol.ctrl.nforce_displacement/ANG;
+                let displace = scf_data.mol.ctrl.nforce_displacement/BOHR;
 
                 let mut position = scf_data.mol.geom.position.iter().map(|x| *x).collect::<Vec<f64>>();
                 lbfgs().minimize(
@@ -637,7 +637,7 @@ fn eval_force(scf_data: &mut SCF, time_mark: &mut utilities::TimeRecords, mpi_op
         if scf_data.mol.ctrl.print_level > 1 {
             println!("Gradient evaluation using numerical differentiation");
         }
-        let displace = scf_data.mol.ctrl.nforce_displacement / ANG;
+        let displace = scf_data.mol.ctrl.nforce_displacement / BOHR;
         let (energy, nforce) = numerical_force(&scf_data, displace, &mpi_operator);
         println!("------ Output gradient [a.u.] ------");
         println!("{}", formated_force(&nforce, &scf_data.mol.geom.elem));
@@ -755,7 +755,7 @@ fn eval_normal_modes(
     let num_atoms = scf_data.mol.geom.nfree;
     let dim = num_atoms * 3;
     let displace_ang = scf_data.mol.ctrl.nhessian_displacement;
-    let displace = displace_ang / ANG; // convert Angstrom to Bohr
+    let displace = displace_ang / BOHR; // convert Angstrom to Bohr
 
     if scf_data.mol.ctrl.print_level > 0 {
         println!("");
@@ -954,7 +954,6 @@ mod geometric_pyo3_impl {
         pyo3::prepare_freethreaded_python();
         
         let elem = scf_data.mol.geom.elem.iter().map(|x| x.as_str()).collect::<Vec<&str>>();
-        const BOHR: f64 = crate::constants::BOHR;
         let xyz = scf_data.mol.geom.position.data.iter().map(|x| x * BOHR).collect::<Vec<f64>>();
         //let xyz = scf_data.mol.geom.position.iter().map(|x| *x).collect::<Vec<f64>>();
         let xyzs = vec![xyz];
diff --git a/src/post_scf_analysis/mod.rs b/src/post_scf_analysis/mod.rs
index 360109fbbe9049a5f7e23fc6eaf177280d21d136..7f9753a08a4b0094346894d6321e7d81db32f344 100644
--- a/src/post_scf_analysis/mod.rs
+++ b/src/post_scf_analysis/mod.rs
@@ -10,7 +10,7 @@ pub mod spin_correction;
 use rest_libcint::prelude::rest_libcint_wrapper::int1e_r;
 use tensors::{MathMatrix, RIFull};
 
-use crate::constants::{ANG, AU2DEBYE};
+use crate::constants::{BOHR, AU2DEBYE};
 use crate::grad::{formated_force, formated_force_ev, numerical_force};
 use crate::mpi_io::MPIOperator;
 use crate::ri_pt2::sbge2::{close_shell_sbge2_rayon, open_shell_sbge2_rayon};
@@ -130,7 +130,7 @@ pub fn post_scf_output(scf_data: &SCF, mpi_operator: &Option<MPIOperator>) {
             }
         } else if output_type.eq("num_force") {
             let displace = match scf_data.mol.geom.unit {
-                crate::geom_io::GeomUnit::Angstrom => scf_data.mol.ctrl.nforce_displacement/ANG,
+                crate::geom_io::GeomUnit::Angstrom => scf_data.mol.ctrl.nforce_displacement/BOHR,
                 crate::geom_io::GeomUnit::Bohr => scf_data.mol.ctrl.nforce_displacement,
             };
             let (energy, num_force) = numerical_force(scf_data, displace, mpi_operator);
diff --git a/src/post_scf_analysis/rrs_pbc.rs b/src/post_scf_analysis/rrs_pbc.rs
index 83c865503847dfda9f971101304dabb9ba3a1d17..1b79252076ecd18c218d776e3bffa87ef039173c 100644
--- a/src/post_scf_analysis/rrs_pbc.rs
+++ b/src/post_scf_analysis/rrs_pbc.rs
@@ -7,7 +7,7 @@ use std::f64::consts::PI;
 use rayon::prelude::*;
 use approx::abs_diff_eq;
 use std::collections::HashMap;
-use crate::constants::ANG;
+use crate::constants::{BOHR, EV};
 use crate::geom_io::{GeomUnit,get_charge};
 use crate::molecule_io::Molecule;
 
@@ -52,12 +52,12 @@ pub fn rrs_pbc_output(scf_data: &SCF, unit_cell_elem: Vec<String>, index_map: Ha
     }
     let mut homo_lumo_gap = vec![0.0;tot_k_points];
     homo_lumo_gap.iter_mut().zip(pbc_result.iter()).for_each(|(h,p)| {
-        *h = (p[(cell_charge/2.0) as usize] - p[(cell_charge/2.0) as usize - 1])*27.2113863;
+        *h = (p[(cell_charge/2.0) as usize] - p[(cell_charge/2.0) as usize - 1])*EV;
     });
     let min_hlg = homo_lumo_gap.iter().filter(|&&x| !x.is_nan()).min_by(|a, b| a.partial_cmp(b).unwrap());
     println!("HOCO-LUCO gap: {:?} eV",min_hlg.unwrap());
     let elec = scf_data.mol.num_elec[0];
-    let hlg = (scf_data.eigenvalues[0][(elec/2.0) as usize] - scf_data.eigenvalues[0][(elec/2.0) as usize - 1])*27.2113863;
+    let hlg = (scf_data.eigenvalues[0][(elec/2.0) as usize] - scf_data.eigenvalues[0][(elec/2.0) as usize - 1])*EV;
     println!("HOMO-LUMO gap: {:?} eV",hlg);
 }
 
@@ -126,7 +126,7 @@ pub fn rrs_pbc_match(mol: &Molecule) -> (HashMap<(i32,i32,i32),Vec<usize>>, Vec<
     let vec_rrs_pbc_vec = rrs_pbc_vec.data;
     let mut tot_pbc_vec = vec![vec![0.0;3];3]; //[3,3]
     let factor = match unit {
-        GeomUnit::Angstrom => ANG,
+        GeomUnit::Angstrom => BOHR,
         GeomUnit::Bohr => 1.0,
     };
     for i in 0..rrs_pbc_dim {
@@ -194,7 +194,7 @@ pub fn rrs_pbc_new(scf: &SCF,index_map: HashMap<(i32,i32,i32),Vec<usize>>) -> (V
     let rrs_pbc_dim = scf.mol.geom.pbc_dim;
     let unit = scf.mol.geom.unit.clone();
     let factor = match unit {
-        GeomUnit::Angstrom => ANG,
+        GeomUnit::Angstrom => BOHR,
         GeomUnit::Bohr => 1.0,
     };
     let mut vec_rrs_pbc_vec = rrs_pbc_vec.data.clone();
diff --git a/src/ri_bse/dipoles.rs b/src/ri_bse/dipoles.rs
index 408eadbfb926c2072f1ae6f4ac17ec3563c39085..f18e51312183fb01a26f50939dd1d5f473b3de00 100644
--- a/src/ri_bse/dipoles.rs
+++ b/src/ri_bse/dipoles.rs
@@ -3,7 +3,7 @@ use crate::mpi_io::MPIOperator;
 use std::path::Path;
 use rest_libcint::prelude::rest_libcint_wrapper::int1e_r;
 use tensors::{MathMatrix, MatrixFull, RIFull,MatrixFullSlice};
-use crate::constants::{ANG, AU2DEBYE, SPECIES_INFO};
+use crate::constants::{AU2DEBYE, SPECIES_INFO};
 use itertools::Itertools;
 use crate::ri_gw::get_occupation_parameters;
 use rest_tensors::matrix::matrix_blas_lapack::{_dgemv};
diff --git a/src/ri_bse/dynamicbse.rs b/src/ri_bse/dynamicbse.rs
index 6b32cc4173463b8a8dcc68d61ee3798b1c60e2af..4ee720608fe9c959c56932028db7a6b8ac3c9b98 100644
--- a/src/ri_bse/dynamicbse.rs
+++ b/src/ri_bse/dynamicbse.rs
@@ -17,6 +17,7 @@ use num::Complex;
 use rand::Rng;
 use std::time::Instant;
 use tensors::{MathMatrix, MatrixFull};
+use crate::constants::EV;
 use rest_tensors::matrix::matrix_blas_lapack::{
     _dgemm_full, _dgemm_scaled, _dgeev,
 };
@@ -833,7 +834,7 @@ pub fn dynamic_bse_main(scf_data: &SCF, qp_ctrl: &QuasiParticle) {
     println!("  #       Excitation energy (eV)    ‖T(λ)x‖");
     println!("  ───     ─────────────────────    ──────────");
     for k in 0..result.n_found {
-        let lam_ev = result.eigenvalues[k] * 27.2114;
+        let lam_ev = result.eigenvalues[k] * EV;
         println!("  {:>3}     {:>12.6} eV             {:>9.2e}",
             k, lam_ev, result.residuals[k]);
     }
@@ -842,7 +843,7 @@ pub fn dynamic_bse_main(scf_data: &SCF, qp_ctrl: &QuasiParticle) {
         for k in 0..result.n_found {
             let xv: Vec<f64> = (0..n).map(|i| result.eigenvectors[[i, k]]).collect();
             println!("\n  Excitation #{}: λ = {:.6} Ha = {:.6} eV",
-                k, result.eigenvalues[k], result.eigenvalues[k] * 27.2114);
+                k, result.eigenvalues[k], result.eigenvalues[k] * EV);
             super::leading_components(&xv, occ_size, vir_size);
         }
     }
diff --git a/src/ri_bse/nonlinbse.rs b/src/ri_bse/nonlinbse.rs
index 89af13a0ab1ef30712f23cd6f81de1128e04b3cb..e51222349204fafcecb6d1b892cdc78def649918 100644
--- a/src/ri_bse/nonlinbse.rs
+++ b/src/ri_bse/nonlinbse.rs
@@ -23,6 +23,7 @@ use rand::Rng;
 use rayon::prelude::*;
 use std::time::Instant;
 use tensors::{MathMatrix, MatrixFull};
+use crate::constants::EV;
 use rest_tensors::matrix::matrix_blas_lapack::{
     _dgemm_full, _dgemm_scaled, _dgeev,
 };
@@ -1218,7 +1219,7 @@ pub fn nlfeast_bse_main(scf_data: &SCF, qp_ctrl: &QuasiParticle) {
     println!("  #       Excitation energy (eV)    ‖T(λ)x‖");
     println!("  ───     ─────────────────────    ──────────");
     for k in 0..result.n_found {
-        let lam_ev = result.eigenvalues[k] * 27.2114;
+        let lam_ev = result.eigenvalues[k] * EV;
         println!("  {:>3}     {:>12.6} eV             {:>9.2e}",
             k, lam_ev, result.residuals[k]);
     }
@@ -1227,7 +1228,7 @@ pub fn nlfeast_bse_main(scf_data: &SCF, qp_ctrl: &QuasiParticle) {
         for k in 0..result.n_found {
             let xv: Vec<f64> = (0..n).map(|i| result.eigenvectors[[i, k]]).collect();
             println!("\n  Excitation #{}: λ = {:.6} Ha = {:.6} eV",
-                k, result.eigenvalues[k], result.eigenvalues[k] * 27.2114);
+                k, result.eigenvalues[k], result.eigenvalues[k] * EV);
             super::leading_components(&xv, occ_size, vir_size);
         }
     }
@@ -1659,7 +1660,7 @@ pub fn nlfeast_dynamical_bse_main(scf_data: &SCF, qp_ctrl: &QuasiParticle) {
     println!("  #       Excitation energy (eV)    ‖T(λ)x‖");
     println!("  ───     ─────────────────────    ──────────");
     for k in 0..result.n_found {
-        let lam_ev = result.eigenvalues[k] * 27.2114;
+        let lam_ev = result.eigenvalues[k] * EV;
         println!("  {:>3}     {:>12.6} eV             {:>9.2e}", k, lam_ev, result.residuals[k]);
     }
 }
diff --git a/src/ri_bse/response.rs b/src/ri_bse/response.rs
index 758b150743170804adf76912b418152948dafd56..7b7afbb3278d46b20b1134054575558f2cf43dad 100644
--- a/src/ri_bse/response.rs
+++ b/src/ri_bse/response.rs
@@ -3,7 +3,7 @@ use crate::mpi_io::MPIOperator;
 use std::path::Path;
 use rest_libcint::rest_libcint_wrapper::int1e_r;
 use tensors::{MathMatrix, MatrixFull, RIFull,MatrixFullSlice};
-use crate::constants::{ANG, AU2DEBYE, SPECIES_INFO};
+use crate::constants::{AU2DEBYE, SPECIES_INFO};
 use itertools::Itertools;
 use crate::ri_gw::get_occupation_parameters;
 use std::f64::consts::SQRT_2;
diff --git a/src/ri_gw/mod.rs b/src/ri_gw/mod.rs
index f5e010ed8f254bb1b7a5b66b474560f3efe8ca40..1c1b947423bf90980c5a3ef6a837875b3671e4d7 100644
--- a/src/ri_gw/mod.rs
+++ b/src/ri_gw/mod.rs
@@ -1,5 +1,5 @@
 //use std::simd::num;
-use crate::constants::PI;
+use crate::constants::{EV, PI};
 use itertools::Itertools;
 use std::ops::Range;
 use crate::utilities;
@@ -141,7 +141,6 @@ pub fn gw_calculations(scf_data:&mut SCF,num_freq:usize,vxc_nn:&Vec<f64>,cancel_
     let v_matrix=v_matrix(&scf_data,&ri_mat);
     let (start_mo,num_state,occ_size,vir_size,homo,lumo)=get_occupation_parameters(scf_data,'Y');
     let (start_mo,num_state_cutoff,occ_size,vir_size_cutoff,homo,lumo)=get_occupation_parameters(scf_data,'N');
-    let hatree=27.2113863;
     let spin_channel=scf_data.mol.ctrl.spin_channel;
     let mut quasiparticle_energies_g:Vec<f64>=Vec::new();
     let mut quasiparticle_energies_w:Vec<f64>=Vec::new();
@@ -1483,7 +1482,7 @@ pub fn generate_real_axis_vchiv(
 
     let grid_label = if grid_type == 1 { "quadratic (power-law)" } else { "linear" };
     println!("Low-rank contour: Maximum real frequency needed for v*chi*v = {:.6} Ha = {:.6} eV",
-             de_max, de_max * 27.2113863);
+             de_max, de_max * EV);
     println!("Low-rank contour: Computing sqrt(v)*chi*sqrt(v) at {} real-axis grid points ({})",
              nomega_chi_real, grid_label);
     if print_level > 2 {
@@ -1508,7 +1507,7 @@ pub fn generate_real_axis_vchiv(
         };
         if print_level > 2 {
             println!("[DEBUG] Real-axis freq {} / {} : omega = {:.10} Ha = {:.6} eV  (t={:.6})",
-                     i_omega + 1, nomega_chi_real, omega_real, omega_real * 27.2113863,
+                     i_omega + 1, nomega_chi_real, omega_real, omega_real * EV,
                      (i_omega as f64) / ((nomega_chi_real - 1) as f64));
         }
         let lr = low_rank_vchi_vsqrt(
diff --git a/src/ri_tddft/response.rs b/src/ri_tddft/response.rs
index 01ab240484ab5aa21b222bd9aa3206aa2ec7d092..c573e40422fffa4af23bd95fdb1a37bd7d637df1 100644
--- a/src/ri_tddft/response.rs
+++ b/src/ri_tddft/response.rs
@@ -20,6 +20,7 @@ use std::fs::OpenOptions;
 use std::io::Write;
 
 use rest_tensors::MatrixFull;
+use crate::constants::EV;
 use rest_tensors::matrix::matrix_blas_lapack::{_dsolve, _dgemm_full};
 
 use crate::ri_bse::response::{
@@ -222,7 +223,7 @@ fn print_response_tddft_results(
     let pol_im: f64 = density_imag.iter().zip(mu_z_vec.iter()).map(|(d, m)| d * m).sum();
 
     println!("  --- Response TDDFT Results ---");
-    println!("    Frequency ω = {:.8} Ha ({:.4} eV)", omega, omega * 27.2114);
+    println!("    Frequency ω = {:.8} Ha ({:.4} eV)", omega, omega * EV);
     println!("    Lifetime  γ = {:.8} Ha", gamma);
     println!("    Re[α_zz(ω)] = {:.12e} a.u.", pol_re);
     println!("    Im[α_zz(ω)] = {:.12e} a.u.", pol_im);
@@ -685,7 +686,7 @@ pub fn response_tddft(scf: &mut SCF) -> Result<(), String> {
     println!("  Method: Response (frequency-domain)");
     println!("  Spin: {}", if xlet == 'S' { "Singlet" } else if xlet == 'T' { "Triplet" } else { "Generic" });
     println!("  occ_size={}, vir_size={}, dim={}", occ_size, vir_size, dim);
-    println!("  ω = {:.8} Ha ({:.4} eV)", omega, omega * 27.2114);
+    println!("  ω = {:.8} Ha ({:.4} eV)", omega, omega * EV);
     println!("  γ = {:.8} Ha", gamma);
 
     // Prepare fxc data
diff --git a/src/solvent/smd_cds.rs b/src/solvent/smd_cds.rs
index c2d60443ce937e8478e1665a8425c72b2d2cc87f..aa54da59cda80c99ebdaaa3ccbe9deb6d8ea1fa7 100644
--- a/src/solvent/smd_cds.rs
+++ b/src/solvent/smd_cds.rs
@@ -62,8 +62,8 @@
 //!
 //! ## Unit Conversions
 //!
-//! - `TO_ANGS = 0.529177` Bohr → Å
-//! - `TO_KCAL = 627.509` Hartree → kcal/mol
+//! - `TO_ANGS = BOHR` Bohr → Å
+//! - `TO_KCAL = HARTREE2KCAL` Hartree → kcal/mol
 //! - Internal computation in Å and kcal/mol; public API input/output in Bohr and Hartree.
 //!
 //! ## References
@@ -74,6 +74,7 @@
 
 use std::f64::consts::PI;
 use tensors::MatrixFull;
+use crate::constants::{BOHR, HARTREE2KCAL};
 
 // ============================================================================
 //  Debug printing helper (controlled by env var REST_CDS_DEBUG=1)
@@ -1880,8 +1881,8 @@ fn cds_eg(
     coords: &[[f64; 3]], atomic_numbers: &[usize],
     sigma: &[f64; 151], hsigma: &[f64; 151], rad: &[f64],
 ) -> (f64, f64, Vec<[f64; 3]>) {
-    const TO_ANGS: f64 = 0.52917724924;    // Bohr → Å
-    const TO_KCAL: f64 = 627.509451;       // Hartree → kcal/mol
+    const TO_ANGS: f64 = BOHR;    // Bohr → Å
+    const TO_KCAL: f64 = HARTREE2KCAL;       // Hartree → kcal/mol
 
     // ---- Build rlio (pairwise distances, Å) and urlio (unit vectors) ----
     let ncot = nat * (nat + 1) / 2;
@@ -2135,12 +2136,12 @@ pub fn compute_cds(
         cds_print_scalar("gcds_kcal (CDS energy, kcal/mol)", gcds_kcal);
         cds_print_scalar("tarea (total SASA, A^2)", tarea);
         cds_print_grad("dcds (CDS gradient, Hartree/Bohr)", &dcds);
-        const TO_KCAL: f64 = 627.509451;
+        const TO_KCAL: f64 = HARTREE2KCAL;
         cds_print_scalar("gcds_hartree (CDS energy, Hartree)", gcds_kcal / TO_KCAL);
     }
 
     // kcal/mol → Hartree (gradient already in Hartree/Bohr from cds_eg)
-    const TO_KCAL: f64 = 627.509451;
+    const TO_KCAL: f64 = HARTREE2KCAL;
     let gcds = gcds_kcal / TO_KCAL;
 
     (gcds, tarea, dcds)
diff --git a/src/x2c/mod.rs b/src/x2c/mod.rs
index 7f54f3ddb12fe96fd5281b162326f70763f4abec..7cd75f571344c3ac3c7108614f79f69840f1c109 100644
--- a/src/x2c/mod.rs
+++ b/src/x2c/mod.rs
@@ -6,7 +6,7 @@ use rest_libcint::prelude::*;
 use rest_libcint::CINTR2CDATA;
 use rest_libcint_wrapper::*;
 use crate::basis_io::BasInfo;
-use crate::constants::ENV_PRT_START;
+use crate::constants::{ENV_PRT_START, LIGHT_SPEED};
 use crate::molecule_io::{Molecule, get_basis_name};
 use serde::{Deserialize, Serialize};
 
@@ -168,7 +168,7 @@ impl Molecule {
 
     
     pub fn generate_sfx2c_hamiltonian(&self) -> MatrixUpper<f64> {
-        let light_speed: f64 = 137.03599967994;
+        let light_speed: f64 = LIGHT_SPEED;
 
         let (xmol, contr_coeff) = self.decontract_basis();
         let n_contracted = self.num_basis;