mopac_core 0.1.1

//! Extended d-orbital rotation, overlap, and two-center two-electron multipoles.
//!
//! Licensed under the Apache License, Version 2.0 (the "License").
//! Authentic mathematical translation of OpenMOPAC `mndod.F90`, `rotmat`, `charg`, `rijkl`, `reppd2`, `tx`, `rotatd`.

#![allow(clippy::too_many_arguments, clippy::needless_range_loop)]

use crate::constants::codata2018::{BOHR_RADIUS_ANGSTROMS as A0_BOHR, HARTREE_TO_EV as EV_HARTREE};

/// OpenMOPAC `indexd` 9x9 address mapping.
pub const INDEXD: [[usize; 9]; 9] = [
    [1, 2, 3, 4, 5, 6, 7, 8, 9],
    [2, 10, 11, 12, 13, 14, 15, 16, 17],
    [3, 11, 18, 19, 20, 21, 22, 23, 24],
    [4, 12, 19, 25, 26, 27, 28, 29, 30],
    [5, 13, 20, 26, 31, 32, 33, 34, 35],
    [6, 14, 21, 27, 32, 36, 37, 38, 39],
    [7, 15, 22, 28, 33, 37, 40, 41, 42],
    [8, 16, 23, 29, 34, 38, 41, 43, 44],
    [9, 17, 24, 30, 35, 39, 42, 44, 45],
];

/// OpenMOPAC standard pair index `indx` 9x9.
pub const INDX: [[usize; 9]; 9] = [
    [1, 2, 4, 7, 11, 16, 22, 29, 37],
    [2, 3, 5, 8, 12, 17, 23, 30, 38],
    [4, 5, 6, 9, 13, 18, 24, 31, 39],
    [7, 8, 9, 10, 14, 19, 25, 32, 40],
    [11, 12, 13, 14, 15, 20, 26, 33, 41],
    [16, 17, 18, 19, 20, 21, 27, 34, 42],
    [22, 23, 24, 25, 26, 27, 28, 35, 43],
    [29, 30, 31, 32, 33, 34, 35, 36, 44],
    [37, 38, 39, 40, 41, 42, 43, 44, 45],
];

pub const INDPP: [[usize; 3]; 3] = [[1, 4, 5], [4, 2, 6], [5, 6, 3]];

pub const INDDP: [[usize; 3]; 5] = [[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]];

pub const INDDD: [[usize; 5]; 5] = [
    [1, 6, 7, 9, 12],
    [6, 2, 8, 10, 13],
    [7, 8, 3, 11, 14],
    [9, 10, 11, 4, 15],
    [12, 13, 14, 15, 5],
];

/// OpenMOPAC `ind2(45, 45)` mapping pair-indices (ij, kl) to 1..491.
pub const IND2: [[usize; 45]; 45] = [
    [
        1, 2, 0, 0, 35, 0, 0, 0, 0, 3, 0, 0, 36, 0, 0, 0, 0, 4, 0, 0, 38, 0, 0, 0, 5, 0, 0, 40, 0,
        0, 37, 0, 0, 0, 0, 39, 0, 0, 0, 41, 0, 0, 42, 0, 43,
    ],
    [
        6, 7, 0, 0, 44, 0, 0, 0, 0, 8, 0, 0, 45, 0, 0, 0, 0, 9, 0, 0, 47, 0, 0, 0, 10, 0, 0, 49, 0,
        0, 46, 0, 0, 0, 0, 48, 0, 0, 0, 50, 0, 0, 51, 0, 52,
    ],
    [
        0, 0, 16, 0, 0, 63, 0, 0, 0, 0, 17, 0, 0, 64, 0, 0, 0, 0, 0, 62, 0, 0, 66, 0, 0, 0, 0, 0,
        0, 68, 0, 65, 0, 0, 0, 0, 0, 67, 0, 0, 0, 69, 0, 0, 0,
    ],
    [
        0, 0, 0, 25, 0, 0, 91, 0, 0, 0, 0, 26, 0, 0, 92, 0, 0, 0, 0, 0, 0, 0, 0, 96, 0, 90, 0, 0,
        94, 0, 0, 0, 93, 0, 0, 0, 0, 0, 97, 0, 95, 0, 0, 0, 0,
    ],
    [
        124, 125, 0, 0, 129, 0, 0, 0, 0, 126, 0, 0, 130, 0, 0, 0, 0, 127, 0, 0, 132, 0, 0, 0, 128,
        0, 0, 134, 0, 0, 131, 0, 0, 0, 0, 133, 0, 0, 0, 135, 0, 0, 136, 0, 137,
    ],
    [
        0, 0, 186, 0, 0, 189, 0, 0, 0, 0, 187, 0, 0, 190, 0, 0, 0, 0, 0, 188, 0, 0, 192, 0, 0, 0,
        0, 0, 0, 194, 0, 191, 0, 0, 0, 0, 0, 193, 0, 0, 0, 195, 0, 0, 0,
    ],
    [
        0, 0, 0, 257, 0, 0, 260, 0, 0, 0, 0, 258, 0, 0, 261, 0, 0, 0, 0, 0, 0, 0, 0, 265, 0, 259,
        0, 0, 263, 0, 0, 0, 262, 0, 0, 0, 0, 0, 266, 0, 264, 0, 0, 0, 0,
    ],
    [
        0, 0, 0, 0, 0, 0, 0, 341, 0, 0, 0, 0, 0, 0, 0, 342, 0, 335, 0, 0, 337, 0, 0, 0, 336, 0, 0,
        339, 0, 0, 0, 0, 0, 343, 0, 338, 0, 0, 0, 340, 0, 0, 0, 0, 0,
    ],
    [
        0, 0, 0, 0, 0, 0, 0, 0, 420, 0, 0, 0, 0, 0, 0, 0, 421, 0, 416, 0, 0, 418, 0, 0, 0, 0, 417,
        0, 0, 0, 0, 0, 0, 0, 422, 0, 419, 0, 0, 0, 0, 0, 0, 0, 0,
    ],
    [
        11, 12, 0, 0, 53, 0, 0, 0, 0, 13, 0, 0, 54, 0, 0, 0, 0, 14, 0, 0, 56, 0, 0, 0, 15, 0, 0,
        58, 0, 0, 55, 0, 0, 0, 0, 57, 0, 0, 0, 59, 0, 0, 60, 0, 61,
    ],
    [
        0, 0, 18, 0, 0, 71, 0, 0, 0, 0, 19, 0, 0, 72, 0, 0, 0, 0, 0, 70, 0, 0, 74, 0, 0, 0, 0, 0,
        0, 76, 0, 73, 0, 0, 0, 0, 0, 75, 0, 0, 0, 77, 0, 0, 0,
    ],
    [
        0, 0, 0, 27, 0, 0, 99, 0, 0, 0, 0, 28, 0, 0, 100, 0, 0, 0, 0, 0, 0, 0, 0, 104, 0, 98, 0, 0,
        102, 0, 0, 0, 101, 0, 0, 0, 0, 0, 105, 0, 103, 0, 0, 0, 0,
    ],
    [
        138, 139, 0, 0, 143, 0, 0, 0, 0, 140, 0, 0, 144, 0, 0, 0, 0, 141, 0, 0, 146, 0, 0, 0, 142,
        0, 0, 148, 0, 0, 145, 0, 0, 0, 0, 147, 0, 0, 0, 149, 0, 0, 150, 0, 151,
    ],
    [
        0, 0, 196, 0, 0, 199, 0, 0, 0, 0, 197, 0, 0, 200, 0, 0, 0, 0, 0, 198, 0, 0, 202, 0, 0, 0,
        0, 0, 0, 204, 0, 201, 0, 0, 0, 0, 0, 203, 0, 0, 0, 205, 0, 0, 0,
    ],
    [
        0, 0, 0, 267, 0, 0, 270, 0, 0, 0, 0, 268, 0, 0, 271, 0, 0, 0, 0, 0, 0, 0, 0, 275, 0, 269,
        0, 0, 273, 0, 0, 0, 272, 0, 0, 0, 0, 0, 276, 0, 274, 0, 0, 0, 0,
    ],
    [
        0, 0, 0, 0, 0, 0, 0, 350, 0, 0, 0, 0, 0, 0, 0, 351, 0, 344, 0, 0, 346, 0, 0, 0, 345, 0, 0,
        348, 0, 0, 0, 0, 0, 352, 0, 347, 0, 0, 0, 349, 0, 0, 0, 0, 0,
    ],
    [
        0, 0, 0, 0, 0, 0, 0, 0, 427, 0, 0, 0, 0, 0, 0, 0, 428, 0, 423, 0, 0, 425, 0, 0, 0, 0, 424,
        0, 0, 0, 0, 0, 0, 0, 429, 0, 426, 0, 0, 0, 0, 0, 0, 0, 0,
    ],
    [
        20, 21, 0, 0, 78, 0, 0, 85, 0, 22, 0, 0, 79, 0, 0, 86, 0, 23, 0, 0, 81, 0, 0, 0, 24, 0, 0,
        83, 0, 0, 80, 0, 0, 87, 0, 82, 0, 0, 0, 84, 0, 0, 88, 0, 89,
    ],
    [
        0, 0, 0, 0, 0, 0, 0, 0, 109, 0, 0, 0, 0, 0, 0, 0, 110, 0, 29, 0, 0, 107, 0, 0, 0, 0, 106,
        0, 0, 0, 0, 0, 0, 0, 111, 0, 108, 0, 0, 0, 0, 0, 0, 0, 0,
    ],
    [
        0, 0, 152, 0, 0, 155, 0, 0, 0, 0, 153, 0, 0, 156, 0, 0, 0, 0, 0, 154, 0, 0, 158, 0, 0, 0,
        0, 0, 0, 160, 0, 157, 0, 0, 0, 0, 0, 159, 0, 0, 0, 161, 0, 0, 0,
    ],
    [
        206, 207, 0, 0, 211, 0, 0, 218, 0, 208, 0, 0, 212, 0, 0, 219, 0, 209, 0, 0, 214, 0, 0, 0,
        210, 0, 0, 216, 0, 0, 213, 0, 0, 220, 0, 215, 0, 0, 0, 217, 0, 0, 221, 0, 222,
    ],
    [
        0, 0, 0, 0, 0, 0, 0, 0, 281, 0, 0, 0, 0, 0, 0, 0, 282, 0, 277, 0, 0, 279, 0, 0, 0, 0, 278,
        0, 0, 0, 0, 0, 0, 0, 283, 0, 280, 0, 0, 0, 0, 0, 0, 0, 0,
    ],
    [
        0, 0, 353, 0, 0, 356, 0, 0, 0, 0, 354, 0, 0, 357, 0, 0, 0, 0, 0, 355, 0, 0, 359, 0, 0, 0,
        0, 0, 0, 361, 0, 358, 0, 0, 0, 0, 0, 360, 0, 0, 0, 362, 0, 0, 0,
    ],
    [
        0, 0, 0, 430, 0, 0, 433, 0, 0, 0, 0, 431, 0, 0, 434, 0, 0, 0, 0, 0, 0, 0, 0, 438, 0, 432,
        0, 0, 436, 0, 0, 0, 435, 0, 0, 0, 0, 0, 439, 0, 437, 0, 0, 0, 0,
    ],
    [
        30, 31, 0, 0, 112, 0, 0, 119, 0, 32, 0, 0, 113, 0, 0, 120, 0, 33, 0, 0, 115, 0, 0, 0, 34,
        0, 0, 117, 0, 0, 114, 0, 0, 121, 0, 116, 0, 0, 0, 118, 0, 0, 122, 0, 123,
    ],
    [
        0, 0, 0, 162, 0, 0, 165, 0, 0, 0, 0, 163, 0, 0, 166, 0, 0, 0, 0, 0, 0, 0, 0, 170, 0, 164,
        0, 0, 168, 0, 0, 0, 167, 0, 0, 0, 0, 0, 171, 0, 169, 0, 0, 0, 0,
    ],
    [
        0, 0, 0, 0, 0, 0, 0, 0, 227, 0, 0, 0, 0, 0, 0, 0, 228, 0, 223, 0, 0, 225, 0, 0, 0, 0, 224,
        0, 0, 0, 0, 0, 0, 0, 229, 0, 226, 0, 0, 0, 0, 0, 0, 0, 0,
    ],
    [
        284, 285, 0, 0, 289, 0, 0, 296, 0, 286, 0, 0, 290, 0, 0, 297, 0, 287, 0, 0, 292, 0, 0, 0,
        288, 0, 0, 294, 0, 0, 291, 0, 0, 298, 0, 293, 0, 0, 0, 295, 0, 0, 299, 0, 300,
    ],
    [
        0, 0, 0, 363, 0, 0, 366, 0, 0, 0, 0, 364, 0, 0, 367, 0, 0, 0, 0, 0, 0, 0, 0, 371, 0, 365,
        0, 0, 369, 0, 0, 0, 368, 0, 0, 0, 0, 0, 372, 0, 370, 0, 0, 0, 0,
    ],
    [
        0, 0, 440, 0, 0, 443, 0, 0, 0, 0, 441, 0, 0, 444, 0, 0, 0, 0, 0, 442, 0, 0, 446, 0, 0, 0,
        0, 0, 0, 448, 0, 445, 0, 0, 0, 0, 0, 447, 0, 0, 0, 449, 0, 0, 0,
    ],
    [
        172, 173, 0, 0, 177, 0, 0, 0, 0, 174, 0, 0, 178, 0, 0, 0, 0, 175, 0, 0, 180, 0, 0, 0, 176,
        0, 0, 182, 0, 0, 179, 0, 0, 0, 0, 181, 0, 0, 0, 183, 0, 0, 184, 0, 185,
    ],
    [
        0, 0, 230, 0, 0, 233, 0, 0, 0, 0, 231, 0, 0, 234, 0, 0, 0, 0, 0, 232, 0, 0, 236, 0, 0, 0,
        0, 0, 0, 238, 0, 235, 0, 0, 0, 0, 0, 237, 0, 0, 0, 239, 0, 0, 0,
    ],
    [
        0, 0, 0, 301, 0, 0, 304, 0, 0, 0, 0, 302, 0, 0, 305, 0, 0, 0, 0, 0, 0, 0, 0, 309, 0, 303,
        0, 0, 307, 0, 0, 0, 306, 0, 0, 0, 0, 0, 310, 0, 308, 0, 0, 0, 0,
    ],
    [
        0, 0, 0, 0, 0, 0, 0, 379, 0, 0, 0, 0, 0, 0, 0, 380, 0, 373, 0, 0, 375, 0, 0, 0, 374, 0, 0,
        377, 0, 0, 0, 0, 0, 381, 0, 376, 0, 0, 0, 378, 0, 0, 0, 0, 0,
    ],
    [
        0, 0, 0, 0, 0, 0, 0, 0, 454, 0, 0, 0, 0, 0, 0, 0, 455, 0, 450, 0, 0, 452, 0, 0, 0, 0, 451,
        0, 0, 0, 0, 0, 0, 0, 456, 0, 453, 0, 0, 0, 0, 0, 0, 0, 0,
    ],
    [
        240, 241, 0, 0, 245, 0, 0, 252, 0, 242, 0, 0, 246, 0, 0, 253, 0, 243, 0, 0, 248, 0, 0, 0,
        244, 0, 0, 250, 0, 0, 247, 0, 0, 254, 0, 249, 0, 0, 0, 251, 0, 0, 255, 0, 256,
    ],
    [
        0, 0, 0, 0, 0, 0, 0, 0, 315, 0, 0, 0, 0, 0, 0, 0, 316, 0, 311, 0, 0, 313, 0, 0, 0, 0, 312,
        0, 0, 0, 0, 0, 0, 0, 317, 0, 314, 0, 0, 0, 0, 0, 0, 0, 0,
    ],
    [
        0, 0, 382, 0, 0, 385, 0, 0, 0, 0, 383, 0, 0, 386, 0, 0, 0, 0, 0, 384, 0, 0, 388, 0, 0, 0,
        0, 0, 0, 390, 0, 387, 0, 0, 0, 0, 0, 389, 0, 0, 0, 391, 0, 0, 0,
    ],
    [
        0, 0, 0, 457, 0, 0, 460, 0, 0, 0, 0, 458, 0, 0, 461, 0, 0, 0, 0, 0, 0, 0, 0, 465, 0, 459,
        0, 0, 463, 0, 0, 0, 462, 0, 0, 0, 0, 0, 466, 0, 464, 0, 0, 0, 0,
    ],
    [
        318, 319, 0, 0, 323, 0, 0, 330, 0, 320, 0, 0, 324, 0, 0, 331, 0, 321, 0, 0, 326, 0, 0, 0,
        322, 0, 0, 328, 0, 0, 325, 0, 0, 332, 0, 327, 0, 0, 0, 329, 0, 0, 333, 0, 334,
    ],
    [
        0, 0, 0, 392, 0, 0, 395, 0, 0, 0, 0, 393, 0, 0, 396, 0, 0, 0, 0, 0, 0, 0, 0, 400, 0, 394,
        0, 0, 398, 0, 0, 0, 397, 0, 0, 0, 0, 0, 401, 0, 399, 0, 0, 0, 0,
    ],
    [
        0, 0, 467, 0, 0, 470, 0, 0, 0, 0, 468, 0, 0, 471, 0, 0, 0, 0, 0, 469, 0, 0, 473, 0, 0, 0,
        0, 0, 0, 475, 0, 472, 0, 0, 0, 0, 0, 474, 0, 0, 0, 476, 0, 0, 0,
    ],
    [
        402, 403, 0, 0, 407, 0, 0, 0, 0, 404, 0, 0, 408, 0, 0, 0, 0, 405, 0, 0, 410, 0, 0, 0, 406,
        0, 0, 412, 0, 0, 409, 0, 0, 0, 0, 411, 0, 0, 0, 413, 0, 0, 414, 0, 415,
    ],
    [
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 477, 0,
    ],
    [
        478, 479, 0, 0, 483, 0, 0, 0, 0, 480, 0, 0, 484, 0, 0, 0, 0, 481, 0, 0, 486, 0, 0, 0, 482,
        0, 0, 488, 0, 0, 485, 0, 0, 0, 0, 487, 0, 0, 0, 489, 0, 0, 490, 0, 491,
    ],
];

/// OpenMOPAC `isym(492)` symmetry relationships between local frame integrals.
pub const ISYM: [i32; 492] = [
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 38, 39, 0, 42, 0, 0, 0, 0, 0, 47, 48, 0, 51, 0, 0, 0, 0, 0, 56, 57, 0,
    60, 0, 0, 0, 0, 0, 0, 66, 67, 0, 0, 0, 0, 0, 0, 74, 75, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 88,
    62, 63, 64, 65, -66, -67, 66, 67, 70, 71, 72, 73, -74, -75, 74, 75, 86, 86, 0, 85, 86, 87, 78,
    79, 80, 83, 84, 81, 82, -85, -86, -87, 88, 88, 0, 0, 0, 0, 127, 0, 0, 0, 0, 0, 132, 133, 0,
    136, 0, 0, 0, 0, 141, 0, 0, 0, 0, 0, 146, 147, 0, 150, 0, 0, 0, 0, 0, 0, 0, 0, 158, 159, 152,
    153, 154, 155, 156, 157, -158, -159, 158, 159, 0, 0, 0, 0, 175, 0, 0, 0, 0, 0, 180, 181, 0,
    184, 0, 0, 0, 0, 0, 0, 0, 0, 192, 193, 0, 0, 0, 0, 0, 0, 0, 0, 202, 203, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 221, 0, 219, 219, 0, 218, 219, 220, 0, 0, 0, 0, 0, 0, 0, 0, 236,
    237, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 255, 186, 187, 188, 189, 190, 191, -192,
    -193, 192, 193, 196, 197, 198, 199, 200, 201, -202, -203, 202, 203, 223, 219, 219, 226, 218,
    219, 220, 206, 207, 208, 210, 209, 211, 212, 213, 216, 217, 214, 215, -218, -219, -220, 221,
    221, 230, 231, 232, 233, 234, 235, -236, -237, 236, 237, 0, 253, 253, 0, 252, 253, 254, 240,
    241, 242, 244, 243, 245, 246, 247, 250, 251, 248, 249, -252, -253, -254, 255, 255, 0, -335, 0,
    0, -337, -338, 0, 337, 0, 223, -223, 219, 226, -219, -226, 218, 219, 220, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, -353, -354, -355, -356, -357, -358, 359, 360, -361, -362, 0, -373, 0, 0, -375, -376,
    0, 375, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -382, -383, -384, -385, -386, -387, 388, 389, -390,
    -391, 0, 0, 0, 0, 405, 0, 0, 0, 0, 0, 410, 411, 0, 0, 335, 337, 337, 338, 341, 337, 343, 223,
    219, 219, 226, 218, 219, 220, 353, 354, 355, 356, 357, 358, -361, -362, 359, 360, 353, 354,
    355, 356, 357, 358, 361, 362, 359, 360, 373, 375, 375, 376, 379, 375, 381, 382, 383, 384, 385,
    386, 387, -390, -391, 388, 389, 382, 383, 384, 385, 386, 387, 390, 391, 388, 389, 0, 402, 403,
    404, 405, 405, 407, 408, 409, 410, 411, 410, 411, 415, 414,
];

/// OpenMOPAC `ch(45, 3, 5)` multipole conversion coefficients (m mapped to m+2).
pub const CH: [[[f64; 5]; 3]; 45] = [
    [
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 1.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 1.15470054, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 1.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 1.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [1.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 1.33333333, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 1.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 1.15470054, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 1.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, -0.66666667, 0.0, 1.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [1.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, -0.57735027, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 1.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, -0.66666667, 0.0, -1.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, -0.57735027, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, -1.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 1.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 1.33333333, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.57735027, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.57735027, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, -1.15470054],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [-1.15470054, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.66666667, 0.0, 1.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [1.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 1.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 1.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.66666667, 0.0, -1.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, -1.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 1.0, 0.0],
    ],
    [
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, -1.33333333, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
    ],
    [
        [0.0, 0.0, 1.0, 0.0, 0.0],
        [0.0, 0.0, 0.0, 0.0, 0.0],
        [0.0, 0.0, -1.33333333, 0.0, 0.0],
    ],
];

/// 3D Spatial Frame Rotation Engine for p- and d-orbitals.
///
/// Implements exact SO(3) representation transformations matching OpenMOPAC `rotmat.F90`.
#[derive(Debug, Clone, Copy)]
pub struct DOrbitalRotation3D {
    pub r_bohr: f64,
    /// 3x3 rotation matrix for p orbitals (px, py, pz)
    pub p: [[f64; 3]; 3],
    /// 5x5 rotation matrix for d orbitals (dx2-y2, dxz, dz2, dyz, dxy)
    pub d: [[f64; 5]; 5],
    /// Tensor product S-P: sp(3, 3)
    pub sp: [[f64; 3]; 3],
    /// Tensor product P-P: pp(6, 3, 3)
    pub pp: [[[f64; 3]; 3]; 6],
    /// Tensor product S-D: sd(5, 5)
    pub sd: [[f64; 5]; 5],
    /// Tensor product D-P: dp(15, 5, 3)
    pub dp: [[[f64; 3]; 5]; 15],
    /// Tensor product D-D: d_d(15, 5, 5)
    pub d_d: [[[f64; 5]; 5]; 15],
}

impl DOrbitalRotation3D {
    /// Construct rotation from displacement vector pointing from atom A to atom B.
    pub fn new(dx_ang: f64, dy_ang: f64, dz_ang: f64, r_ang: f64) -> Self {
        let r_bohr = r_ang / A0_BOHR;
        let small = 1.0e-7;
        let pt5sq3 = 3.0f64.sqrt() * 0.5;

        let b = dx_ang * dx_ang + dy_ang * dy_ang;
        let sqb = b.sqrt();
        let sb = sqb / r_ang;

        let (ca, sa, cb) = if sb > small {
            (dx_ang / sqb, dy_ang / sqb, dz_ang / r_ang)
        } else {
            let (ca_val, cb_val) = if dz_ang < 0.0 {
                (-1.0, -1.0)
            } else if dz_ang > 0.0 {
                (1.0, 1.0)
            } else {
                (0.0, 0.0)
            };
            (ca_val, 0.0, cb_val)
        };

        // p(3, 3) matrix
        let mut p = [[0.0f64; 3]; 3];
        p[0][0] = ca * sb;
        p[1][0] = ca * cb;
        p[2][0] = -sa;
        p[0][1] = sa * sb;
        p[1][1] = sa * cb;
        p[2][1] = ca;
        p[0][2] = cb;
        p[1][2] = -sb;
        p[2][2] = 0.0;

        let c2a = 2.0 * ca * ca - 1.0;
        let c2b = 2.0 * cb * cb - 1.0;
        let s2a = 2.0 * sa * ca;
        let s2b = 2.0 * sb * cb;

        // d(5, 5) matrix
        let mut d = [[0.0f64; 5]; 5];
        d[0][0] = pt5sq3 * c2a * sb * sb;
        d[1][0] = 0.5 * c2a * s2b;
        d[2][0] = -s2a * sb;
        d[3][0] = c2a * (cb * cb + 0.5 * sb * sb);
        d[4][0] = -s2a * cb;

        d[0][1] = pt5sq3 * ca * s2b;
        d[1][1] = ca * c2b;
        d[2][1] = -sa * cb;
        d[3][1] = -0.5 * ca * s2b;
        d[4][1] = sa * sb;

        d[0][2] = cb * cb - 0.5 * sb * sb;
        d[1][2] = -pt5sq3 * s2b;
        d[2][2] = 0.0;
        d[3][2] = pt5sq3 * sb * sb;
        d[4][2] = 0.0;

        d[0][3] = pt5sq3 * sa * s2b;
        d[1][3] = sa * c2b;
        d[2][3] = ca * cb;
        d[3][3] = -0.5 * sa * s2b;
        d[4][3] = -ca * sb;

        d[0][4] = pt5sq3 * s2a * sb * sb;
        d[1][4] = 0.5 * s2a * s2b;
        d[2][4] = c2a * sb;
        d[3][4] = s2a * (cb * cb + 0.5 * sb * sb);
        d[4][4] = c2a * cb;

        // sp = p
        let sp = p;

        // pp(6, 3, 3)
        let mut pp = [[[0.0f64; 3]; 3]; 6];
        for k in 0..3 {
            pp[0][k][k] = p[k][0] * p[k][0];
            pp[1][k][k] = p[k][1] * p[k][1];
            pp[2][k][k] = p[k][2] * p[k][2];
            pp[3][k][k] = p[k][0] * p[k][1];
            pp[4][k][k] = p[k][0] * p[k][2];
            pp[5][k][k] = p[k][1] * p[k][2];
            for prev in 0..k {
                pp[0][k][prev] = 2.0 * p[k][0] * p[prev][0];
                pp[1][k][prev] = 2.0 * p[k][1] * p[prev][1];
                pp[2][k][prev] = 2.0 * p[k][2] * p[prev][2];
                pp[3][k][prev] = p[k][0] * p[prev][1] + p[k][1] * p[prev][0];
                pp[4][k][prev] = p[k][0] * p[prev][2] + p[k][2] * p[prev][0];
                pp[5][k][prev] = p[k][1] * p[prev][2] + p[k][2] * p[prev][1];
            }
        }

        // sd = d
        let sd = d;

        // dp(15, 5, 3)
        let mut dp = [[[0.0f64; 3]; 5]; 15];
        for k in 0..5 {
            for col in 0..3 {
                dp[0][k][col] = d[k][0] * p[col][0];
                dp[1][k][col] = d[k][0] * p[col][1];
                dp[2][k][col] = d[k][0] * p[col][2];
                dp[3][k][col] = d[k][1] * p[col][0];
                dp[4][k][col] = d[k][1] * p[col][1];
                dp[5][k][col] = d[k][1] * p[col][2];
                dp[6][k][col] = d[k][2] * p[col][0];
                dp[7][k][col] = d[k][2] * p[col][1];
                dp[8][k][col] = d[k][2] * p[col][2];
                dp[9][k][col] = d[k][3] * p[col][0];
                dp[10][k][col] = d[k][3] * p[col][1];
                dp[11][k][col] = d[k][3] * p[col][2];
                dp[12][k][col] = d[k][4] * p[col][0];
                dp[13][k][col] = d[k][4] * p[col][1];
                dp[14][k][col] = d[k][4] * p[col][2];
            }
        }

        // d_d(15, 5, 5)
        let mut d_d = [[[0.0f64; 5]; 5]; 15];
        for k in 0..5 {
            d_d[0][k][k] = d[k][0] * d[k][0];
            d_d[1][k][k] = d[k][1] * d[k][1];
            d_d[2][k][k] = d[k][2] * d[k][2];
            d_d[3][k][k] = d[k][3] * d[k][3];
            d_d[4][k][k] = d[k][4] * d[k][4];
            d_d[5][k][k] = d[k][0] * d[k][1];
            d_d[6][k][k] = d[k][0] * d[k][2];
            d_d[7][k][k] = d[k][1] * d[k][2];
            d_d[8][k][k] = d[k][0] * d[k][3];
            d_d[9][k][k] = d[k][1] * d[k][3];
            d_d[10][k][k] = d[k][2] * d[k][3];
            d_d[11][k][k] = d[k][0] * d[k][4];
            d_d[12][k][k] = d[k][1] * d[k][4];
            d_d[13][k][k] = d[k][2] * d[k][4];
            d_d[14][k][k] = d[k][3] * d[k][4];
            for prev in 0..k {
                d_d[0][k][prev] = 2.0 * d[k][0] * d[prev][0];
                d_d[1][k][prev] = 2.0 * d[k][1] * d[prev][1];
                d_d[2][k][prev] = 2.0 * d[k][2] * d[prev][2];
                d_d[3][k][prev] = 2.0 * d[k][3] * d[prev][3];
                d_d[4][k][prev] = 2.0 * d[k][4] * d[prev][4];
                d_d[5][k][prev] = d[k][0] * d[prev][1] + d[k][1] * d[prev][0];
                d_d[6][k][prev] = d[k][0] * d[prev][2] + d[k][2] * d[prev][0];
                d_d[7][k][prev] = d[k][1] * d[prev][2] + d[k][2] * d[prev][1];
                d_d[8][k][prev] = d[k][0] * d[prev][3] + d[k][3] * d[prev][0];
                d_d[9][k][prev] = d[k][1] * d[prev][3] + d[k][3] * d[prev][1];
                d_d[10][k][prev] = d[k][2] * d[prev][3] + d[k][3] * d[prev][2];
                d_d[11][k][prev] = d[k][0] * d[prev][4] + d[k][4] * d[prev][0];
                d_d[12][k][prev] = d[k][1] * d[prev][4] + d[k][4] * d[prev][1];
                d_d[13][k][prev] = d[k][2] * d[prev][4] + d[k][4] * d[prev][2];
                d_d[14][k][prev] = d[k][3] * d[prev][4] + d[k][4] * d[prev][3];
            }
        }

        Self {
            r_bohr,
            p,
            d,
            sp,
            pp,
            sd,
            dp,
            d_d,
        }
    }

    /// Check strict SO(3) orthonormality of p (3x3) and d (5x5) rotation matrices.
    /// Returns (||P P^T - I_3||_inf, ||D D^T - I_5||_inf).
    pub fn check_orthonormality(&self) -> (f64, f64) {
        let mut err_p = 0.0f64;
        for i in 0..3 {
            for j in 0..3 {
                let mut dot = 0.0;
                for k in 0..3 {
                    dot += self.p[i][k] * self.p[j][k];
                }
                let target = if i == j { 1.0 } else { 0.0 };
                err_p = err_p.max((dot - target).abs());
            }
        }

        let mut err_d = 0.0f64;
        for i in 0..5 {
            for j in 0..5 {
                let mut dot = 0.0;
                for k in 0..5 {
                    dot += self.d[i][k] * self.d[j][k];
                }
                let target = if i == j { 1.0 } else { 0.0 };
                err_d = err_d.max((dot - target).abs());
            }
        }

        (err_p, err_d)
    }

    /// Compute Casimir invariant preservation: ||D v||^2 - ||v||^2 for any 5D spherical vector.
    pub fn casimir_invariance(&self, v: &[f64; 5]) -> f64 {
        let mut norm_in_sq = 0.0;
        for &x in v {
            norm_in_sq += x * x;
        }

        let mut norm_out_sq = 0.0;
        for i in 0..5 {
            let mut dot = 0.0;
            for j in 0..5 {
                dot += self.d[i][j] * v[j];
            }
            norm_out_sq += dot * dot;
        }

        (norm_out_sq - norm_in_sq).abs()
    }
}

/// Evaluation of the two point-charge interaction configurations (MNDO/d).
/// Direct port of OpenMOPAC `charg` (`mndod.F90` lines 2160-2284).
pub fn charg(r: f64, l1: usize, l2: usize, m: usize, da: f64, db: f64, add: f64) -> f64 {
    let mut val = 0.0f64;
    if l1 == 0 && l2 == 0 {
        val = 1.0 / (r * r + add).sqrt();
    } else if l1 == 1 && l2 == 0 {
        val = (-1.0 / ((r + da) * (r + da) + add).sqrt()
            + 1.0 / ((r - da) * (r - da) + add).sqrt())
            * 0.5;
    } else if l1 == 0 && l2 == 1 {
        val = (1.0 / ((r + db) * (r + db) + add).sqrt() - 1.0 / ((r - db) * (r - db) + add).sqrt())
            * 0.5;
    } else if l1 == 1 && l2 == 1 && m == 0 {
        let dzdz = 1.0 / ((r + da - db) * (r + da - db) + add).sqrt()
            + 1.0 / ((r - da + db) * (r - da + db) + add).sqrt()
            - 1.0 / ((r - da - db) * (r - da - db) + add).sqrt()
            - 1.0 / ((r + da + db) * (r + da + db) + add).sqrt();
        val = dzdz * 0.25;
    } else if l1 == 1 && l2 == 1 && m == 1 {
        let dxdx = 2.0 / (r * r + (da - db) * (da - db) + add).sqrt()
            - 2.0 / (r * r + (da + db) * (da + db) + add).sqrt();
        val = dxdx * 0.25;
    } else if l1 == 0 && l2 == 2 {
        let qqzz = 1.0 / ((r - db) * (r - db) + add).sqrt() - 2.0 / (r * r + db * db + add).sqrt()
            + 1.0 / ((r + db) * (r + db) + add).sqrt();
        val = qqzz * 0.25;
    } else if l1 == 2 && l2 == 0 {
        let qzzq = 1.0 / ((r - da) * (r - da) + add).sqrt() - 2.0 / (r * r + da * da + add).sqrt()
            + 1.0 / ((r + da) * (r + da) + add).sqrt();
        val = qzzq * 0.25;
    } else if l1 == 1 && l2 == 2 && m == 0 {
        let dzqzz = 1.0 / ((r - da - db) * (r - da - db) + add).sqrt()
            - 2.0 / ((r - da) * (r - da) + db * db + add).sqrt()
            + 1.0 / ((r + db - da) * (r + db - da) + add).sqrt()
            - 1.0 / ((r - db + da) * (r - db + da) + add).sqrt()
            + 2.0 / ((r + da) * (r + da) + db * db + add).sqrt()
            - 1.0 / ((r + da + db) * (r + da + db) + add).sqrt();
        val = dzqzz * 0.125;
    } else if l1 == 2 && l2 == 1 && m == 0 {
        let qzzdz = -1.0 / ((r - da - db) * (r - da - db) + add).sqrt()
            + 2.0 / ((r - db) * (r - db) + da * da + add).sqrt()
            - 1.0 / ((r + da - db) * (r + da - db) + add).sqrt()
            + 1.0 / ((r - da + db) * (r - da + db) + add).sqrt()
            - 2.0 / ((r + db) * (r + db) + da * da + add).sqrt()
            + 1.0 / ((r + da + db) * (r + da + db) + add).sqrt();
        val = qzzdz * 0.125;
    } else if l1 == 2 && l2 == 2 && m == 0 {
        let zzzz = 1.0 / ((r - da - db) * (r - da - db) + add).sqrt()
            + 1.0 / ((r + da + db) * (r + da + db) + add).sqrt()
            + 1.0 / ((r - da + db) * (r - da + db) + add).sqrt()
            + 1.0 / ((r + da - db) * (r + da - db) + add).sqrt()
            - 2.0 / ((r - da) * (r - da) + db * db + add).sqrt()
            - 2.0 / ((r - db) * (r - db) + da * da + add).sqrt()
            - 2.0 / ((r + da) * (r + da) + db * db + add).sqrt()
            - 2.0 / ((r + db) * (r + db) + da * da + add).sqrt()
            + 2.0 / (r * r + (da - db) * (da - db) + add).sqrt()
            + 2.0 / (r * r + (da + db) * (da + db) + add).sqrt();
        let xyxy = 4.0 / (r * r + (da - db) * (da - db) + add).sqrt()
            + 4.0 / (r * r + (da + db) * (da + db) + add).sqrt()
            - 8.0 / (r * r + da * da + db * db + add).sqrt();
        val = zzzz * 0.0625 - xyxy * 0.015625;
    } else if l1 == 1 && l2 == 2 && m == 1 {
        let ab = db / 2.0f64.sqrt();
        let dxqxz = -2.0 / ((r - ab) * (r - ab) + (da - ab) * (da - ab) + add).sqrt()
            + 2.0 / ((r + ab) * (r + ab) + (da - ab) * (da - ab) + add).sqrt()
            + 2.0 / ((r - ab) * (r - ab) + (da + ab) * (da + ab) + add).sqrt()
            - 2.0 / ((r + ab) * (r + ab) + (da + ab) * (da + ab) + add).sqrt();
        val = dxqxz * 0.125;
    } else if l1 == 2 && l2 == 1 && m == 1 {
        let aa = da / 2.0f64.sqrt();
        let qxzdx = -2.0 / ((r + aa) * (r + aa) + (aa - db) * (aa - db) + add).sqrt()
            + 2.0 / ((r - aa) * (r - aa) + (aa - db) * (aa - db) + add).sqrt()
            + 2.0 / ((r + aa) * (r + aa) + (aa + db) * (aa + db) + add).sqrt()
            - 2.0 / ((r - aa) * (r - aa) + (aa + db) * (aa + db) + add).sqrt();
        val = qxzdx * 0.125;
    } else if l1 == 2 && l2 == 2 && m == 1 {
        let aa = da / 2.0f64.sqrt();
        let ab = db / 2.0f64.sqrt();
        let qxzqxz = 2.0 / ((r + aa - ab) * (r + aa - ab) + (aa - ab) * (aa - ab) + add).sqrt()
            - 2.0 / ((r + aa + ab) * (r + aa + ab) + (aa - ab) * (aa - ab) + add).sqrt()
            - 2.0 / ((r - aa - ab) * (r - aa - ab) + (aa - ab) * (aa - ab) + add).sqrt()
            + 2.0 / ((r - aa + ab) * (r - aa + ab) + (aa - ab) * (aa - ab) + add).sqrt()
            - 2.0 / ((r + aa - ab) * (r + aa - ab) + (aa + ab) * (aa + ab) + add).sqrt()
            + 2.0 / ((r + aa + ab) * (r + aa + ab) + (aa + ab) * (aa + ab) + add).sqrt()
            + 2.0 / ((r - aa - ab) * (r - aa - ab) + (aa + ab) * (aa + ab) + add).sqrt()
            - 2.0 / ((r - aa + ab) * (r - aa + ab) + (aa + ab) * (aa + ab) + add).sqrt();
        val = qxzqxz * 0.0625;
    } else if l1 == 2 && l2 == 2 && m == 2 {
        let xyxy = 4.0 / (r * r + (da - db) * (da - db) + add).sqrt()
            + 4.0 / (r * r + (da + db) * (da + db) + add).sqrt()
            - 8.0 / (r * r + da * da + db * db + add).sqrt();
        val = xyxy * 0.0625;
    }
    val
}

/// Compute diatomic multipole interaction element in local frame.
/// Direct port of OpenMOPAC `rijkl` (`mndod.F90` lines 1000-1093).
pub fn rijkl(
    po_a: &[f64; 10],
    ddp_a: &[f64; 7],
    po_b: &[f64; 10],
    ddp_b: &[f64; 7],
    ij_pair: usize, // 1-based index 1..45
    kl_pair: usize, // 1-based index 1..45
    li: usize,
    lj: usize,
    lk: usize,
    ll: usize,
    ic: usize,
    r_bohr: f64,
) -> f64 {
    let l1min = (li as isize - lj as isize).unsigned_abs().min(2);
    let l1max = (li + lj).min(2);
    let lij = INDX[li][lj]; // 1-based pair index: 1..6

    let l2min = (lk as isize - ll as isize).unsigned_abs().min(2);
    let l2max = (lk + ll).min(2);
    let lkl = INDX[lk][ll]; // 1-based pair index: 1..6

    let mut sum = 0.0f64;

    for l1 in l1min..=l1max {
        let (dij, pij) = if l1 == 0 {
            let p = match lij {
                1 => {
                    if ic == 1 {
                        po_a[9]
                    } else {
                        po_a[1]
                    }
                }
                3 => po_a[7],
                6 => po_a[8],
                _ => po_a[1],
            };
            (0.0, p)
        } else {
            (ddp_a[lij], po_a[lij])
        };

        for l2 in l2min..=l2max {
            let (dkl, pkl) = if l2 == 0 {
                let p = match lkl {
                    1 => {
                        if ic == 2 {
                            po_b[9]
                        } else {
                            po_b[1]
                        }
                    }
                    3 => po_b[7],
                    6 => po_b[8],
                    _ => po_b[1],
                };
                (0.0, p)
            } else {
                (ddp_b[lkl], po_b[lkl])
            };

            let add = (pij + pkl) * (pij + pkl);
            let lmin = l1.min(l2);
            let mut s1 = 0.0f64;

            for m in -(lmin as isize)..=(lmin as isize) {
                let m_idx = (m + 2) as usize;
                let ccc = CH[ij_pair - 1][l1][m_idx] * CH[kl_pair - 1][l2][m_idx];
                if ccc.abs() < 1e-15 {
                    continue;
                }
                let mm = m.unsigned_abs();
                s1 += charg(r_bohr, l1, l2, mm, dij, dkl, add) * ccc;
            }
            sum += s1;
        }
    }
    sum
}

pub const MET: [usize; 45] = [
    1, 2, 3, 2, 3, 3, 2, 3, 3, 3, 4, 5, 5, 5, 6, 4, 5, 5, 5, 6, 6, 4, 5, 5, 5, 6, 6, 6, 4, 5, 5, 5,
    6, 6, 6, 6, 4, 5, 5, 5, 6, 6, 6, 6, 6,
];

pub const IPOS: [usize; 34] = [
    1, 5, 11, 12, 12, 2, 6, 13, 14, 14, 3, 8, 16, 18, 18, 7, 15, 10, 20, 4, 9, 17, 19, 21, 7, 15,
    10, 20, 22, 4, 9, 17, 21, 19,
];

pub const LORB: [usize; 9] = [0, 1, 1, 1, 2, 2, 2, 2, 2];

/// Compute 491 local-frame two-center two-electron repulsion integrals.
/// Direct port of OpenMOPAC `reppd2` (`mndod.F90` lines 846-997).
pub fn compute_reppd2(
    norb_a: usize,
    norb_b: usize,
    r_bohr: f64,
    ri: &[f64; 22],
    po_a: &[f64; 10],
    ddp_a: &[f64; 7],
    po_b: &[f64; 10],
    ddp_b: &[f64; 7],
    tore_a: f64,
    tore_b: f64,
    rep: &mut [f64; 492],
    core: &mut [[f64; 2]; 10],
    feather: bool,
) {
    let ev = EV_HARTREE;
    for (m, &pos) in IPOS.iter().enumerate() {
        rep[m + 1] = ri[pos - 1];
    }

    let has_d_a = norb_a >= 9;
    let has_d_b = norb_b >= 9;

    if has_d_a || has_d_b {
        let lasti = norb_a;
        let lastk = norb_b;

        for i in 1..=lasti {
            let li = LORB[i - 1];
            for j in 1..=i {
                let lj = LORB[j - 1];
                let ij = INDEXD[i - 1][j - 1];
                let coul = i == j;

                for k in 1..=lastk {
                    let lk = LORB[k - 1];
                    for l in 1..=k {
                        let ll = LORB[l - 1];
                        let kl = INDEXD[k - 1][l - 1];

                        let numb = IND2[ij - 1][kl - 1];
                        if numb <= 34 || numb > 491 {
                            continue;
                        }

                        let nold = ISYM[numb];
                        if nold >= 35 {
                            rep[numb] = rep[nold as usize];
                        } else if nold <= -35 {
                            rep[numb] = -rep[(-nold) as usize];
                        } else if nold == 0 {
                            let mut val =
                                rijkl(po_a, ddp_a, po_b, ddp_b, ij, kl, li, lj, lk, ll, 0, r_bohr)
                                    * ev;
                            if feather {
                                let r_angstrom = r_bohr * A0_BOHR;
                                let const_val = if r_angstrom < 7.0 {
                                    let dr = r_angstrom - 7.0;
                                    1.0 - (-dr * dr * 0.22).exp()
                                } else {
                                    0.0
                                };
                                let point =
                                    crate::constants::codata2018::EV_ANGSTROM_FACTOR / r_angstrom;
                                let coulomb = coul && k == l;
                                if coulomb {
                                    val = val * const_val + (1.0 - const_val) * point;
                                } else {
                                    val *= const_val;
                                }
                            }
                            rep[numb] = val;
                        }
                    }
                }
            }
        }

        for c in 4..10 {
            core[c][0] = 0.0;
            core[c][1] = 0.0;
        }

        if has_d_b {
            let ij = INDEXD[0][0]; // (s, s) on center A
                                   // <S S | D S>
            let mut kl = INDEXD[4][0];
            core[4][1] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, ij, kl, 0, 0, 2, 0, 1, r_bohr) * ev * tore_a;
            // <S S | D P>
            kl = INDEXD[4][1];
            core[5][1] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, ij, kl, 0, 0, 2, 1, 1, r_bohr) * ev * tore_a;
            // <S S | D D>
            kl = INDEXD[4][4];
            core[6][1] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, ij, kl, 0, 0, 2, 2, 1, r_bohr) * ev * tore_a;
            // <S S | D+ P+>
            kl = INDEXD[5][2];
            core[7][1] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, ij, kl, 0, 0, 2, 1, 1, r_bohr) * ev * tore_a;
            // <S S | D+ D+>
            kl = INDEXD[5][5];
            core[8][1] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, ij, kl, 0, 0, 2, 2, 1, r_bohr) * ev * tore_a;
            // <S S | D# D#>
            kl = INDEXD[7][7];
            core[9][1] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, ij, kl, 0, 0, 2, 2, 1, r_bohr) * ev * tore_a;
        }

        if has_d_a {
            let ij = INDEXD[0][0]; // (s, s) on center B
                                   // <D S | S S>
            let mut kl = INDEXD[4][0];
            core[4][0] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, kl, ij, 2, 0, 0, 0, 2, r_bohr) * ev * tore_b;
            // <D P | S S>
            kl = INDEXD[4][1];
            core[5][0] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, kl, ij, 2, 1, 0, 0, 2, r_bohr) * ev * tore_b;
            // <D D | S S>
            kl = INDEXD[4][4];
            core[6][0] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, kl, ij, 2, 2, 0, 0, 2, r_bohr) * ev * tore_b;
            // <D+ P+ | S S>
            kl = INDEXD[5][2];
            core[7][0] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, kl, ij, 2, 1, 0, 0, 2, r_bohr) * ev * tore_b;
            // <D+ D+ | S S>
            kl = INDEXD[5][5];
            core[8][0] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, kl, ij, 2, 2, 0, 0, 2, r_bohr) * ev * tore_b;
            // <D# D# | S S>
            kl = INDEXD[7][7];
            core[9][0] =
                -rijkl(po_a, ddp_a, po_b, ddp_b, kl, ij, 2, 2, 0, 0, 2, r_bohr) * ev * tore_b;
        }
    }
}

/// First transformation step of two-center two-electron integrals.
/// Direct port of OpenMOPAC `tx` (`mndod.F90` lines 1870-1989).
pub fn tx(
    norb_a: usize,
    norb_b: usize,
    rep: &[f64; 492],
    logv: &mut [[bool; 45]; 45],
    v: &mut [[f64; 45]; 45],
    rot: &DOrbitalRotation3D,
) {
    let limkl = INDX[norb_b - 1][norb_b - 1];
    for row in logv.iter_mut() {
        row[..limkl].fill(false);
    }
    for row in v.iter_mut() {
        row[..limkl].fill(0.0);
    }

    for i1 in 1..=norb_a {
        for j1 in 1..=i1 {
            let ij = INDEXD[i1 - 1][j1 - 1];

            for k1 in 1..=norb_b {
                for l1 in 1..=k1 {
                    let kl = INDEXD[k1 - 1][l1 - 1];
                    let nd = IND2[ij - 1][kl - 1];
                    if nd == 0 {
                        continue;
                    }
                    let wrepp = rep[nd];
                    let ll = INDX[k1 - 1][l1 - 1];
                    let mm = MET[ll - 1];

                    match mm {
                        1 => {
                            v[ij - 1][0] = wrepp;
                        }
                        2 => {
                            let k = k1 - 2;
                            v[ij - 1][1] += rot.sp[k][0] * wrepp;
                            v[ij - 1][3] += rot.sp[k][1] * wrepp;
                            v[ij - 1][6] += rot.sp[k][2] * wrepp;
                        }
                        3 => {
                            let k = k1 - 2;
                            let l = l1 - 2;
                            v[ij - 1][2] += rot.pp[0][k][l] * wrepp;
                            v[ij - 1][5] += rot.pp[1][k][l] * wrepp;
                            v[ij - 1][9] += rot.pp[2][k][l] * wrepp;
                            v[ij - 1][4] += rot.pp[3][k][l] * wrepp;
                            v[ij - 1][7] += rot.pp[4][k][l] * wrepp;
                            v[ij - 1][8] += rot.pp[5][k][l] * wrepp;
                        }
                        4 => {
                            let k = k1 - 5;
                            v[ij - 1][10] += rot.sd[k][0] * wrepp;
                            v[ij - 1][15] += rot.sd[k][1] * wrepp;
                            v[ij - 1][21] += rot.sd[k][2] * wrepp;
                            v[ij - 1][28] += rot.sd[k][3] * wrepp;
                            v[ij - 1][36] += rot.sd[k][4] * wrepp;
                        }
                        5 => {
                            let k = k1 - 5;
                            let l = l1 - 2;
                            v[ij - 1][11] += rot.dp[0][k][l] * wrepp;
                            v[ij - 1][12] += rot.dp[1][k][l] * wrepp;
                            v[ij - 1][13] += rot.dp[2][k][l] * wrepp;
                            v[ij - 1][16] += rot.dp[3][k][l] * wrepp;
                            v[ij - 1][17] += rot.dp[4][k][l] * wrepp;
                            v[ij - 1][18] += rot.dp[5][k][l] * wrepp;
                            v[ij - 1][22] += rot.dp[6][k][l] * wrepp;
                            v[ij - 1][23] += rot.dp[7][k][l] * wrepp;
                            v[ij - 1][24] += rot.dp[8][k][l] * wrepp;
                            v[ij - 1][29] += rot.dp[9][k][l] * wrepp;
                            v[ij - 1][30] += rot.dp[10][k][l] * wrepp;
                            v[ij - 1][31] += rot.dp[11][k][l] * wrepp;
                            v[ij - 1][37] += rot.dp[12][k][l] * wrepp;
                            v[ij - 1][38] += rot.dp[13][k][l] * wrepp;
                            v[ij - 1][39] += rot.dp[14][k][l] * wrepp;
                        }
                        6 => {
                            let k = k1 - 5;
                            let l = l1 - 5;
                            v[ij - 1][14] += rot.d_d[0][k][l] * wrepp;
                            v[ij - 1][20] += rot.d_d[1][k][l] * wrepp;
                            v[ij - 1][27] += rot.d_d[2][k][l] * wrepp;
                            v[ij - 1][35] += rot.d_d[3][k][l] * wrepp;
                            v[ij - 1][44] += rot.d_d[4][k][l] * wrepp;
                            v[ij - 1][19] += rot.d_d[5][k][l] * wrepp;
                            v[ij - 1][25] += rot.d_d[6][k][l] * wrepp;
                            v[ij - 1][26] += rot.d_d[7][k][l] * wrepp;
                            v[ij - 1][32] += rot.d_d[8][k][l] * wrepp;
                            v[ij - 1][33] += rot.d_d[9][k][l] * wrepp;
                            v[ij - 1][34] += rot.d_d[10][k][l] * wrepp;
                            v[ij - 1][40] += rot.d_d[11][k][l] * wrepp;
                            v[ij - 1][41] += rot.d_d[12][k][l] * wrepp;
                            v[ij - 1][42] += rot.d_d[13][k][l] * wrepp;
                            v[ij - 1][43] += rot.d_d[14][k][l] * wrepp;
                        }
                        _ => {}
                    }
                }
            }
            for kl in 0..limkl {
                if v[ij - 1][kl].abs() > 1.0e-15 {
                    logv[ij - 1][kl] = true;
                }
            }
        }
    }
}

/// Second transformation step of two-center two-electron integrals.
/// Rotates first center from `v` into `ww: [f64; 2025]`.
/// Direct port of OpenMOPAC `rotatd` (`mndod.F90` lines 1183-1350).
pub fn rotatd_step2(
    norb_a: usize,
    norb_b: usize,
    v: &[[f64; 45]; 45],
    logv: &[[bool; 45]; 45],
    rot: &DOrbitalRotation3D,
    ww: &mut [f64; 2025],
) {
    let limkl = INDX[norb_b - 1][norb_b - 1];
    let limij = INDX[norb_a - 1][norb_a - 1];
    ww[..limij * limkl].fill(0.0);

    for i1 in 1..=norb_a {
        for j1 in 1..=i1 {
            let ij = INDEXD[i1 - 1][j1 - 1];
            let jj = INDX[i1 - 1][j1 - 1];
            let mm = MET[jj - 1];

            for k in 1..=norb_b {
                for l in 1..=k {
                    let kl = INDX[k - 1][l - 1];
                    if !logv[ij - 1][kl - 1] {
                        continue;
                    }
                    let wrepp = v[ij - 1][kl - 1];

                    let indw = |i_orb: usize, j_orb: usize| -> usize {
                        (INDX[i_orb - 1][j_orb - 1] - 1) * limkl + (kl - 1)
                    };

                    match mm {
                        1 => {
                            let iw = indw(1, 1);
                            ww[iw] = wrepp;
                        }
                        2 => {
                            for i in 1..=3 {
                                let iw = indw(i + 1, 1);
                                ww[iw] += rot.sp[i1 - 2][i - 1] * wrepp;
                            }
                        }
                        3 => {
                            for i in 1..=3 {
                                let cc = rot.pp[0][i1 - 2][j1 - 2];
                                let iw = indw(i + 1, i + 1);
                                ww[iw] += cc * wrepp;
                                let iminus = i - 1;
                                if iminus == 0 {
                                    continue;
                                }
                                for j in 1..=iminus {
                                    let cc = rot.pp[i + j][i1 - 2][j1 - 2];
                                    let iw = indw(i + 1, j + 1);
                                    ww[iw] += cc * wrepp;
                                }
                            }
                        }
                        4 => {
                            for i in 1..=5 {
                                let iw = indw(i + 4, 1);
                                ww[iw] += rot.sd[i1 - 5][i - 1] * wrepp;
                            }
                        }
                        5 => {
                            for i in 1..=5 {
                                for j in 1..=3 {
                                    let iw = indw(i + 4, j + 1);
                                    let ij1 = 3 * (i - 1) + j;
                                    ww[iw] += rot.dp[ij1 - 1][i1 - 5][j1 - 2] * wrepp;
                                }
                            }
                        }
                        6 => {
                            for i in 1..=5 {
                                let cc = rot.d_d[i - 1][i1 - 5][j1 - 5];
                                let iw = indw(i + 4, i + 4);
                                ww[iw] += cc * wrepp;
                                let iminus = i - 1;
                                if iminus == 0 {
                                    continue;
                                }
                                for j in 1..=iminus {
                                    let ij1 = INDDD[i - 1][j - 1];
                                    let cc = rot.d_d[ij1 - 1][i1 - 5][j1 - 5];
                                    let iw = indw(i + 4, j + 4);
                                    ww[iw] += cc * wrepp;
                                }
                            }
                        }
                        _ => {}
                    }
                }
            }
        }
    }
}

/// Apply PM7 specific d-orbital coulomb corrections to ww and cored.
/// Direct port of OpenMOPAC `rotatd` lines 1264-1350.
pub fn apply_pm7_d_correction(
    norb_a: usize,
    norb_b: usize,
    has_d_a: bool,
    has_d_b: bool,
    ww: &mut [f64; 2025],
    cored: &mut [[f64; 2]; 10],
) {
    let k_stride = INDX[norb_b - 1][norb_b - 1];

    if has_d_a {
        if norb_b > 1 {
            let mut sum = 0.0;
            for i in 5..=9 {
                let j = k_stride * ((i * (i + 1)) / 2 - 1);
                sum += ww[j + 2] + ww[j + 5] + ww[j + 9];
            }
            sum = ww[0] - sum / 15.0;
            for i in 5..=9 {
                let j = k_stride * ((i * (i + 1)) / 2 - 1);
                for l in 2..=4 {
                    let idx = (l * (l + 1)) / 2 + j - 1;
                    ww[idx] += sum;
                }
            }
            let sum_core = cored[0][0] - (cored[2][0] + 2.0 * cored[3][0]) / 3.0;
            cored[2][0] += sum_core;
            cored[3][0] += sum_core;
        }

        let mut sum = 0.0;
        for i in 5..=9 {
            sum += ww[k_stride * ((i * (i + 1)) / 2 - 1)];
        }
        sum = ww[0] - sum / 5.0;
        for i in 5..=9 {
            let idx = k_stride * ((i * (i + 1)) / 2 - 1);
            ww[idx] += sum;
        }
        let sum_core = cored[0][0] - (cored[6][0] + 2.0 * cored[8][0] + 2.0 * cored[9][0]) / 5.0;
        cored[6][0] += sum_core;
        cored[8][0] += sum_core;
        cored[9][0] += sum_core;
    }

    if has_d_b {
        if has_d_a {
            let mut sum = 0.0;
            for i in 5..=9 {
                let j = 45 * ((i * (i + 1)) / 2 - 1);
                sum += ww[j + 14] + ww[j + 20] + ww[j + 27] + ww[j + 35] + ww[j + 44];
            }
            sum = ww[0] - sum / 25.0;
            for i in 5..=9 {
                let j = 45 * ((i * (i + 1)) / 2 - 1);
                for k in 5..=9 {
                    let idx = (k * (k + 1)) / 2 + j - 1;
                    ww[idx] += sum;
                }
            }
        }
        if norb_a > 1 {
            let mut sum = 0.0;
            for i in 2..=4 {
                let j = 45 * ((i * (i + 1)) / 2 - 1);
                sum += ww[j + 14] + ww[j + 20] + ww[j + 27] + ww[j + 35] + ww[j + 44];
            }
            sum = ww[0] - sum / 15.0;
            for i in 2..=4 {
                let j = 45 * ((i * (i + 1)) / 2 - 1);
                for k in 5..=9 {
                    let idx = (k * (k + 1)) / 2 + j - 1;
                    ww[idx] += sum;
                }
            }
            let sum_core = cored[0][1] - (cored[2][1] + 2.0 * cored[3][1]) / 3.0;
            cored[2][1] += sum_core;
            cored[3][1] += sum_core;
        }
        let sum = ww[0] - (ww[14] + ww[20] + ww[27] + ww[35] + ww[44]) / 5.0;
        for k in 5..=9 {
            let idx = (k * (k + 1)) / 2 - 1;
            ww[idx] += sum;
        }
        let sum_core = cored[0][1] - (cored[6][1] + 2.0 * cored[8][1] + 2.0 * cored[9][1]) / 5.0;
        cored[6][1] += sum_core;
        cored[8][1] += sum_core;
        cored[9][1] += sum_core;
    }
}

/// Accumulate nuclear stabilization terms `elenuc` into one-electron matrix `h` for atom pair A and B.
/// Direct port of OpenMOPAC `elenuc` (`mndod.F90` lines 2367-2488).
pub fn accumulate_elenuc_block(
    norb_a: usize,
    norb_b: usize,
    orb_start_a: usize,
    orb_start_b: usize,
    cored: &[[f64; 2]; 10],
    rot: &DOrbitalRotation3D,
    h: &mut crate::types::AlignedMatrix<f64>,
) {
    // Atom A: center n=0
    for i in 0..norb_a {
        for j in 0..=i {
            let val = if i == 0 && j == 0 {
                cored[0][0]
            } else if i < 4 && j == 0 {
                rot.sp[0][i - 1] * cored[1][0]
            } else if i < 4 && j < 4 {
                let ipp = INDPP[i - 1][j - 1] - 1;
                cored[2][0] * rot.pp[ipp][0][0]
                    + cored[3][0] * (rot.pp[ipp][1][1] + rot.pp[ipp][2][2])
            } else if i >= 4 && j == 0 {
                rot.sd[0][i - 4] * cored[4][0]
            } else if i >= 4 && j < 4 {
                let idp = INDDP[i - 4][j - 1] - 1;
                cored[5][0] * rot.dp[idp][0][0]
                    + cored[7][0] * (rot.dp[idp][1][1] + rot.dp[idp][2][2])
            } else {
                let idd = INDDD[i - 4][j - 4] - 1;
                cored[6][0] * rot.d_d[idd][0][0]
                    + cored[8][0] * (rot.d_d[idd][1][1] + rot.d_d[idd][2][2])
                    + cored[9][0] * (rot.d_d[idd][3][3] + rot.d_d[idd][4][4])
            };
            let idx_i = orb_start_a + i;
            let idx_j = orb_start_a + j;
            let cur = h.get(idx_i, idx_j);
            h.set(idx_i, idx_j, cur + val);
            if i != j {
                h.set(idx_j, idx_i, cur + val);
            }
        }
    }

    // Atom B: center n=1
    for i in 0..norb_b {
        for j in 0..=i {
            let val = if i == 0 && j == 0 {
                cored[0][1]
            } else if i < 4 && j == 0 {
                rot.sp[0][i - 1] * cored[1][1]
            } else if i < 4 && j < 4 {
                let ipp = INDPP[i - 1][j - 1] - 1;
                cored[2][1] * rot.pp[ipp][0][0]
                    + cored[3][1] * (rot.pp[ipp][1][1] + rot.pp[ipp][2][2])
            } else if i >= 4 && j == 0 {
                rot.sd[0][i - 4] * cored[4][1]
            } else if i >= 4 && j < 4 {
                let idp = INDDP[i - 4][j - 1] - 1;
                cored[5][1] * rot.dp[idp][0][0]
                    + cored[7][1] * (rot.dp[idp][1][1] + rot.dp[idp][2][2])
            } else {
                let idd = INDDD[i - 4][j - 4] - 1;
                cored[6][1] * rot.d_d[idd][0][0]
                    + cored[8][1] * (rot.d_d[idd][1][1] + rot.d_d[idd][2][2])
                    + cored[9][1] * (rot.d_d[idd][3][3] + rot.d_d[idd][4][4])
            };
            let idx_i = orb_start_b + i;
            let idx_j = orb_start_b + j;
            let cur = h.get(idx_i, idx_j);
            h.set(idx_i, idx_j, cur + val);
            if i != j {
                h.set(idx_j, idx_i, cur + val);
            }
        }
    }
}

/// Principal quantum numbers n for s, p shells (elements 1..=107).
pub const III: [usize; 108] = [
    0, // 0
    1, 1, // H, He (1..2)
    2, 2, 2, 2, 2, 2, 2, 2, // Li..Ne (3..10)
    3, 3, 3, 3, 3, 3, 3, 3, // Na..Ar (11..18)
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, // K..Kr (19..36)
    5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, // Rb..Xe (37..54)
    6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
    6, // Cs..Rn (55..86)
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
];

/// Principal quantum numbers n for d shells (elements 1..=107).
pub const IIID: [usize; 108] = [
    0, // 0
    3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
    3, // 1..30 (3d)
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, // 31..48 (4d)
    5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
    5, // 49..80 (5d)
    6, 6, 6, 6, 6, 6, // 81..86 (6d)
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
];

/// INTIJ table (243 entries) mapping pair (ij) for wstore.
pub const INTIJ: [usize; 243] = [
    1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6,
    7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12,
    12, 12, 13, 13, 13, 13, 13, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16,
    16, 17, 17, 17, 17, 17, 18, 18, 18, 19, 19, 19, 19, 19, 20, 20, 20, 20, 20, 21, 21, 21, 21, 21,
    21, 21, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 24, 24, 24,
    24, 24, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 28, 28,
    28, 28, 28, 29, 29, 29, 29, 29, 30, 30, 30, 31, 31, 31, 31, 31, 32, 32, 32, 32, 32, 33, 33, 33,
    33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 37,
    37, 37, 37, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 40, 40, 40, 41, 42, 42, 42, 42, 42, 43, 43,
    43, 43, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 45, 45, 45, 45, 45,
];

/// INTKL table (243 entries) mapping pair (kl) for wstore.
pub const INTKL: [usize; 243] = [
    15, 21, 28, 36, 45, 12, 19, 23, 39, 11, 15, 21, 22, 26, 28, 36, 45, 13, 24, 32, 38, 34, 37, 43,
    11, 15, 21, 22, 26, 28, 36, 45, 17, 25, 31, 16, 20, 27, 44, 29, 33, 35, 42, 15, 21, 22, 28, 36,
    45, 3, 6, 11, 21, 26, 36, 2, 12, 19, 23, 39, 4, 13, 24, 32, 38, 14, 17, 31, 1, 3, 6, 10, 15,
    21, 22, 28, 36, 45, 8, 16, 20, 27, 44, 7, 14, 17, 25, 31, 18, 30, 40, 2, 12, 19, 23, 39, 8, 16,
    20, 27, 44, 1, 3, 6, 10, 11, 15, 21, 22, 26, 28, 36, 45, 3, 6, 10, 15, 21, 22, 28, 36, 45, 2,
    12, 19, 23, 39, 4, 13, 24, 32, 38, 7, 17, 25, 31, 3, 6, 11, 21, 26, 36, 8, 16, 20, 27, 44, 1,
    3, 6, 10, 15, 21, 22, 28, 36, 45, 9, 29, 33, 35, 42, 18, 30, 40, 7, 14, 17, 25, 31, 4, 13, 24,
    32, 38, 9, 29, 33, 35, 42, 5, 34, 37, 43, 9, 29, 33, 35, 42, 1, 3, 6, 10, 11, 15, 21, 22, 26,
    28, 36, 45, 5, 34, 37, 43, 4, 13, 24, 32, 38, 2, 12, 19, 23, 39, 18, 30, 40, 41, 9, 29, 33, 35,
    42, 5, 34, 37, 43, 8, 16, 20, 27, 44, 1, 3, 6, 10, 15, 21, 22, 28, 36, 45,
];

/// INTREP table (243 entries) mapping to repd(1..52) for wstore.
pub const INTREP: [usize; 243] = [
    1, 1, 1, 1, 1, 3, 3, 8, 3, 9, 6, 6, 12, 14, 13, 7, 6, 15, 8, 3, 3, 11, 9, 14, 17, 6, 7, 12, 18,
    13, 6, 6, 3, 2, 3, 9, 11, 10, 11, 9, 16, 10, 11, 7, 6, 4, 5, 6, 7, 9, 17, 19, 32, 22, 40, 3,
    33, 34, 27, 46, 15, 33, 28, 41, 47, 35, 35, 42, 1, 6, 6, 7, 29, 38, 22, 31, 38, 51, 9, 19, 32,
    21, 32, 3, 35, 33, 24, 34, 35, 35, 35, 3, 34, 33, 26, 34, 11, 32, 44, 37, 49, 1, 6, 7, 6, 32,
    38, 29, 21, 39, 30, 38, 38, 12, 12, 4, 22, 21, 19, 20, 21, 22, 8, 27, 26, 25, 27, 8, 28, 25,
    26, 27, 2, 24, 23, 24, 14, 18, 22, 39, 48, 45, 10, 21, 37, 36, 37, 1, 13, 13, 5, 31, 30, 20,
    29, 30, 31, 9, 19, 40, 21, 32, 35, 35, 35, 3, 42, 34, 24, 33, 3, 41, 26, 33, 34, 16, 40, 44,
    43, 50, 11, 44, 32, 39, 10, 21, 43, 36, 37, 1, 7, 6, 6, 40, 38, 38, 21, 45, 30, 29, 38, 9, 32,
    19, 22, 3, 47, 27, 34, 33, 3, 46, 34, 27, 33, 35, 35, 35, 52, 11, 32, 50, 37, 44, 14, 39, 22,
    48, 11, 32, 49, 37, 44, 1, 6, 6, 7, 51, 38, 22, 31, 38, 29,
];

/// Factorials fx(1..30) where fx[i] = (i-1)!, and Pascal triangle b(1..30, 1..30) where b[i][j] = binom(i-1, j-1).
pub fn init_fbx() -> ([f64; 31], [[f64; 31]; 31]) {
    let mut fx = [0.0f64; 31];
    let mut b = [[0.0f64; 31]; 31];
    fx[1] = 1.0;
    for i in 2..=30 {
        fx[i] = fx[i - 1] * (i - 1) as f64;
    }
    for n in 1..=30 {
        b[n][1] = 1.0;
    }
    for i in 2..=30 {
        for j in 2..=i {
            b[i][j] = b[i - 1][j - 1] + b[i - 1][j];
        }
    }
    (fx, b)
}

/// Calculate radial part of one-center two-electron integrals (Slater-Condon parameter).
/// Direct port of OpenMOPAC `rsc` (`mndod.F90` lines 1570-1638).
pub fn rsc(
    k: usize,
    na: usize,
    ea: f64,
    nb: usize,
    eb: f64,
    nc: usize,
    ec: f64,
    nd: usize,
    ed: f64,
    fx: &[f64; 31],
    b: &[[f64; 31]; 31],
) -> f64 {
    let aea = ea.ln();
    let aeb = eb.ln();
    let aec = ec.ln();
    let aed = ed.ln();
    let nab = na + nb;
    let ncd = nc + nd;
    let ecd = ec + ed;
    let eab = ea + eb;
    let e = ecd + eab;
    let n = nab + ncd;
    let ae = e.ln();
    let a2 = 2.0f64.ln();
    let acd = ecd.ln();
    let aab = eab.ln();
    let ff = fx[n] / (fx[2 * na + 1] * fx[2 * nb + 1] * fx[2 * nc + 1] * fx[2 * nd + 1]).sqrt();
    let c = EV_HARTREE
        * ff
        * (na as f64 * aea
            + nb as f64 * aeb
            + nc as f64 * aec
            + nd as f64 * aed
            + 0.5 * (aea + aeb + aec + aed)
            + a2 * (n as f64 + 2.0)
            - ae * (n as f64))
            .exp();
    let mut s0 = 1.0 / e;
    let mut s1 = 0.0;
    let mut s2 = 0.0;
    let m = ncd - k;
    for i in 1..=m {
        s0 = s0 * e / ecd;
        s1 += s0 * (b[ncd - k][i] - b[ncd + k + 1][i]) / b[n][i];
    }
    let m1 = m + 1;
    let m2 = ncd + k + 1;
    for i in m1..=m2 {
        s0 = s0 * e / ecd;
        s2 += s0 * b[m2][i] / b[n][i];
    }
    let s3 = (ae * (n as f64) - acd * (m2 as f64) - aab * ((nab - k) as f64)).exp() / b[n][m2];
    c * (s1 - s2 + s3)
}

/// Compute the 12 Slater-Condon parameters for an element with d-orbitals.
/// Direct port of OpenMOPAC `scprm` (`mndod.F90` lines 1641-1680).
pub fn scprm(
    ni: u8,
    zsn: f64,
    zpn: f64,
    zdn: f64,
    fx: &[f64; 31],
    b: &[[f64; 31]; 31],
) -> (f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64) {
    let ns = III[ni as usize];
    let nd = IIID[ni as usize];
    let es = zsn;
    let ep = zpn;
    let ed = zdn;
    let r016 = rsc(0, ns, es, ns, es, nd, ed, nd, ed, fx, b);
    let r036 = rsc(0, ns, ep, ns, ep, nd, ed, nd, ed, fx, b);
    let r066 = rsc(0, nd, ed, nd, ed, nd, ed, nd, ed, fx, b);
    let r155 = rsc(1, ns, ep, nd, ed, ns, ep, nd, ed, fx, b);
    let r125 = rsc(1, ns, es, ns, ep, ns, ep, nd, ed, fx, b);
    let r244 = rsc(2, ns, es, nd, ed, ns, es, nd, ed, fx, b);
    let r236 = rsc(2, ns, ep, ns, ep, nd, ed, nd, ed, fx, b);
    let r266 = rsc(2, nd, ed, nd, ed, nd, ed, nd, ed, fx, b);
    let r234 = rsc(2, ns, ep, ns, ep, ns, es, nd, ed, fx, b);
    let r246 = rsc(2, ns, es, nd, ed, nd, ed, nd, ed, fx, b);
    let r355 = rsc(3, ns, ep, nd, ed, ns, ep, nd, ed, fx, b);
    let r466 = rsc(4, nd, ed, nd, ed, nd, ed, nd, ed, fx, b);
    (
        r066, r266, r466, r016, r244, r036, r236, r155, r355, r125, r234, r246,
    )
}

/// Compute the 52 repd terms from the Slater-Condon parameters.
/// Direct port of OpenMOPAC `inighd` (`mndod.F90` lines 62-163).
pub fn inighd(
    mut r016: f64,
    r066: f64,
    mut r244: f64,
    r266: f64,
    r466: f64,
    r036: f64,
    r236: f64,
    r155: f64,
    r355: f64,
    r125: f64,
    r234: f64,
    r246: f64,
    f0sd_override: f64,
    g2sd_override: f64,
    repd: &mut [f64; 53],
) {
    let s3 = 3.0f64.sqrt();
    let s5 = 5.0f64.sqrt();
    let s15 = 15.0f64.sqrt();

    if f0sd_override > 0.001 {
        r016 = f0sd_override;
    }
    if g2sd_override > 0.001 {
        r244 = g2sd_override;
    }

    repd[1] = r016;
    repd[2] = 2.0 / (3.0 * s5) * r125;
    repd[3] = 1.0 / s15 * r125;
    repd[4] = 2.0 / (5.0 * s5) * r234;
    repd[5] = r036 + 4.0 / 35.0 * r236;
    repd[6] = r036 + 2.0 / 35.0 * r236;
    repd[7] = r036 - 4.0 / 35.0 * r236;
    repd[8] = -1.0 / (3.0 * s5) * r125;
    repd[9] = (3.0 / 125.0f64).sqrt() * r234;
    repd[10] = s3 / 35.0 * r236;
    repd[11] = 3.0 / 35.0 * r236;
    repd[12] = -1.0 / (5.0 * s5) * r234;
    repd[13] = r036 - 2.0 / 35.0 * r236;
    repd[14] = -2.0 * s3 / 35.0 * r236;
    repd[15] = -repd[3];
    repd[16] = -repd[11];
    repd[17] = -repd[9];
    repd[18] = -repd[14];
    repd[19] = 1.0 / 5.0 * r244;
    repd[20] = 2.0 / (7.0 * s5) * r246;
    repd[21] = repd[20] / 2.0;
    repd[22] = -repd[20];
    repd[23] = 4.0 / 15.0 * r155 + 27.0 / 245.0 * r355;
    repd[24] = 2.0 * s3 / 15.0 * r155 - 9.0 * s3 / 245.0 * r355;
    repd[25] = 1.0 / 15.0 * r155 + 18.0 / 245.0 * r355;
    repd[26] = -s3 / 15.0 * r155 + 12.0 * s3 / 245.0 * r355;
    repd[27] = -s3 / 15.0 * r155 - 3.0 * s3 / 245.0 * r355;
    repd[28] = -repd[27];
    repd[29] = r066 + 4.0 / 49.0 * r266 + 4.0 / 49.0 * r466;
    repd[30] = r066 + 2.0 / 49.0 * r266 - 24.0 / 441.0 * r466;
    repd[31] = r066 - 4.0 / 49.0 * r266 + 6.0 / 441.0 * r466;
    repd[32] = (3.0 / 245.0f64).sqrt() * r246;
    repd[33] = 1.0 / 5.0 * r155 + 24.0 / 245.0 * r355;
    repd[34] = 1.0 / 5.0 * r155 - 6.0 / 245.0 * r355;
    repd[35] = 3.0 / 49.0 * r355;
    repd[36] = 1.0 / 49.0 * r266 + 30.0 / 441.0 * r466;
    repd[37] = s3 / 49.0 * r266 - 5.0 * s3 / 441.0 * r466;
    repd[38] = r066 - 2.0 / 49.0 * r266 - 4.0 / 441.0 * r466;
    repd[39] = -2.0 * s3 / 49.0 * r266 + 10.0 * s3 / 441.0 * r466;
    repd[40] = -repd[32];
    repd[41] = -repd[34];
    repd[42] = -repd[35];
    repd[43] = -repd[37];
    repd[44] = 3.0 / 49.0 * r266 + 20.0 / 441.0 * r466;
    repd[45] = -repd[39];
    repd[46] = 1.0 / 5.0 * r155 - 3.0 / 35.0 * r355;
    repd[47] = -repd[46];
    repd[48] = 4.0 / 49.0 * r266 + 15.0 / 441.0 * r466;
    repd[49] = 3.0 / 49.0 * r266 - 5.0 / 147.0 * r466;
    repd[50] = -repd[49];
    repd[51] = r066 + 4.0 / 49.0 * r266 - 34.0 / 441.0 * r466;
    repd[52] = 35.0 / 441.0 * r466;
}

/// Auxiliary overlap integral moment for multipole charge separations.
/// Equation (7) of Thiel & Voityuk (1996) / OpenMOPAC `aijl` (`mndod.F90` lines 2146-2157).
pub fn aijl(z1: f64, z2: f64, n1: usize, n2: usize, l: usize, fx: &[f64; 31]) -> f64 {
    let zz = z1 + z2 + 1.0e-20;
    fx[n1 + n2 + l + 1] / (fx[2 * n1 + 1] * fx[2 * n2 + 1]).sqrt()
        * (2.0 * z1 / zz).powi(n1 as i32)
        * (2.0 * z1 / zz).sqrt()
        * (2.0 * z2 / zz).powi(n2 as i32)
        * (2.0 * z2 / zz).sqrt()
        * (2.0f64.powi(l as i32) / zz.powi(l as i32))
}

/// Calculate multipole moments aij for an element.
/// Direct port of OpenMOPAC `aijm` (`mndod.F90` lines 2097-2142).
pub fn aijm(ni: u8, zs: f64, zp: f64, zd: f64, has_d: bool, fx: &[f64; 31]) -> [f64; 7] {
    let mut aij = [0.0f64; 7];
    if ni < 3 {
        return aij;
    }
    let zz = zs * zp;
    if zz < 0.01 {
        return aij;
    }
    let nsp = III[ni as usize];
    aij[2] = aijl(zs, zp, nsp, nsp, 1, fx);
    aij[3] = aijl(zp, zp, nsp, nsp, 2, fx);
    if has_d {
        let nd = IIID[ni as usize];
        aij[4] = aijl(zs, zd, nsp, nd, 2, fx);
        aij[5] = aijl(zp, zd, nsp, nd, 1, fx);
        aij[6] = aijl(zd, zd, nd, nd, 2, fx);
    }
    aij
}

/// Determine additive terms rho = poij for two-center two-electron integrals.
/// Golden-section root finder reproducing one-center terms.
/// Direct port of OpenMOPAC `poij` (`mndod.F90` lines 165-242).
pub fn poij(l: usize, d: f64, fg: f64) -> f64 {
    let ev = EV_HARTREE;
    if l == 0 {
        return 0.5 * ev / fg;
    }
    let dsq = d * d;
    let ev4 = ev * 0.25;
    let ev8 = ev / 8.0;
    let mut a1 = 0.1;
    let mut a2 = 5.0;
    const EPSIL: f64 = 1.0e-8;
    const G1: f64 = 0.382;
    const G2: f64 = 0.618;

    if l == 1 {
        for _ in 0..100 {
            let delta = a2 - a1;
            if delta < EPSIL {
                break;
            }
            let y1 = a1 + delta * G1;
            let y2 = a1 + delta * G2;
            let f1_val = ev4 * (1.0 / y1 - 1.0 / (y1 * y1 + dsq).sqrt()) - fg;
            let f2_val = ev4 * (1.0 / y2 - 1.0 / (y2 * y2 + dsq).sqrt()) - fg;
            if f1_val * f1_val < f2_val * f2_val {
                a2 = y2;
            } else {
                a1 = y1;
            }
        }
    } else if l == 2 {
        for _ in 0..100 {
            let delta = a2 - a1;
            if delta < EPSIL {
                break;
            }
            let y1 = a1 + delta * G1;
            let y2 = a1 + delta * G2;
            let f1_val = ev8
                * (1.0 / y1 - 2.0 / (y1 * y1 + dsq * 0.5).sqrt() + 1.0 / (y1 * y1 + dsq).sqrt())
                - fg;
            let f2_val = ev8
                * (1.0 / y2 - 2.0 / (y2 * y2 + dsq * 0.5).sqrt() + 1.0 / (y2 * y2 + dsq).sqrt())
                - fg;
            if f1_val * f1_val < f2_val * f2_val {
                a2 = y2;
            } else {
                a1 = y1;
            }
        }
    }
    (a1 + a2) * 0.5
}

/// Calculate charge separations ddp and additive terms po.
/// Direct port of OpenMOPAC `ddpo` (`mndod.F90` lines 2287-2363).
pub fn ddpo(
    ni: u8,
    gss: f64,
    hsp: f64,
    gpp: f64,
    gp2: f64,
    has_d: bool,
    aij: &[f64; 7],
    repd: &[f64; 53],
    ddp: &mut [f64; 7],
    po: &mut [f64; 10],
) {
    if gss > 0.1 {
        po[1] = poij(0, 1.0, gss);
    }
    if ni >= 3 {
        // SP
        let d_sp = aij[2] / 12.0f64.sqrt();
        ddp[2] = d_sp;
        po[2] = poij(1, d_sp, hsp);

        // PP
        po[7] = po[1];
        let d_pp = (aij[3] * 0.1).sqrt();
        let fg_pp = 0.5 * (gpp - gp2);
        ddp[3] = d_pp;
        po[3] = poij(2, d_pp, fg_pp);

        if has_d {
            // SD
            let da = (1.0 / 60.0f64).sqrt();
            let d_sd = (aij[4] * da).sqrt();
            let fg_sd = repd[19];
            ddp[4] = d_sd;
            po[4] = poij(2, d_sd, fg_sd);

            // PD
            let d_pd = aij[5] / 20.0f64.sqrt();
            let fg_pd = repd[23] - 1.8 * repd[35];
            ddp[5] = d_pd;
            po[5] = poij(1, d_pd, fg_pd);

            // DD
            let fg_dd = 0.2 * (repd[29] + 2.0 * repd[30] + 2.0 * repd[31]);
            if fg_dd > 1.0e-5 {
                po[8] = poij(0, 1.0, fg_dd);
            } else {
                po[8] = 1.0e5;
            }
            let d_dd = (aij[6] / 14.0).sqrt();
            let fg_dd2 = repd[44] - (20.0 / 35.0) * repd[52];
            ddp[6] = d_dd;
            po[6] = poij(2, d_dd, fg_dd2);
        }
    }
}

/// Complete set of precomputed parameters for an element with d-orbitals.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct DElementParams {
    pub z: u8,
    pub repd: [f64; 53],
    pub ddp: [f64; 7],
    pub po: [f64; 10],
}

/// Compute derived parameters for an element with d-orbitals.
pub fn compute_d_element_params(
    z: u8,
    zs: f64,
    zp: f64,
    zd: f64,
    gss: f64,
    _gsp: f64,
    gpp: f64,
    gp2: f64,
    hsp: f64,
    zsn: f64,
    zpn: f64,
    zdn: f64,
    f0sd_override: f64,
    g2sd_override: f64,
    pocord: f64,
) -> DElementParams {
    let (fx, b) = init_fbx();
    let (r066, r266, r466, r016, r244, r036, r236, r155, r355, r125, r234, r246) =
        scprm(z, zsn, zpn, zdn, &fx, &b);

    let mut repd = [0.0f64; 53];
    inighd(
        r016,
        r066,
        r244,
        r266,
        r466,
        r036,
        r236,
        r155,
        r355,
        r125,
        r234,
        r246,
        f0sd_override,
        g2sd_override,
        &mut repd,
    );

    let aij = aijm(z, zs, zp, zd, true, &fx);
    let mut ddp = [0.0f64; 7];
    let mut po = [0.0f64; 10];
    ddpo(z, gss, hsp, gpp, gp2, true, &aij, &repd, &mut ddp, &mut po);

    po[9] = po[1];
    if pocord > 1.0e-5 {
        po[9] = pocord;
    }

    DElementParams { z, repd, ddp, po }
}

/// Fill square matrix of one-center two-electron integrals (wstore).
/// Direct port of OpenMOPAC `wstore` (`mndod.F90` lines 2027-2095).
pub fn wstore(
    _ni: u8,
    norb: usize,
    gss: f64,
    gsp: f64,
    gpp: f64,
    gp2: f64,
    hsp: f64,
    repd: Option<&[f64; 53]>,
    w: &mut [f64],
) {
    let ilim = (norb * (norb + 1)) / 2;
    w[..ilim * ilim].fill(0.0);

    let idx = |i: usize, j: usize| -> usize { (i - 1) * ilim + (j - 1) };

    w[idx(1, 1)] = gss;
    if norb >= 4 {
        let ipx = 3;
        let ipy = 6;
        let ipz = 10;
        w[idx(ipx, 1)] = gsp;
        w[idx(ipy, 1)] = gsp;
        w[idx(ipz, 1)] = gsp;
        w[idx(1, ipx)] = gsp;
        w[idx(1, ipy)] = gsp;
        w[idx(1, ipz)] = gsp;
        w[idx(ipx, ipx)] = gpp;
        w[idx(ipy, ipy)] = gpp;
        w[idx(ipz, ipz)] = gpp;
        w[idx(ipy, ipx)] = gp2;
        w[idx(ipz, ipx)] = gp2;
        w[idx(ipz, ipy)] = gp2;
        w[idx(ipx, ipy)] = gp2;
        w[idx(ipx, ipz)] = gp2;
        w[idx(ipy, ipz)] = gp2;
        w[idx(2, 2)] = hsp;
        w[idx(4, 4)] = hsp;
        w[idx(7, 7)] = hsp;
        w[idx(5, 5)] = 0.5 * (gpp - gp2);
        w[idx(8, 8)] = 0.5 * (gpp - gp2);
        w[idx(9, 9)] = 0.5 * (gpp - gp2);

        if ilim > 10 {
            if let Some(rep) = repd {
                for i in 0..243 {
                    let ij = INTIJ[i];
                    let kl = INTKL[i];
                    let int_rep = INTREP[i];
                    w[idx(ij, kl)] = rep[int_rep];
                }
            }
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_sulfur_d_parameters_parity() {
        // PM7 parameters for Sulfur (Z = 16)
        let z = 16u8;
        let zs = 2.046153;
        let zp = 1.807678;
        let zd = 3.510309;
        let gss = 8.728478;
        let gsp = 6.483871;
        let gpp = 7.357401;
        let gp2 = 6.875448;
        let hsp = 3.012199;
        let zsn = 1.131343;
        let zpn = 0.823803;
        let zdn = 2.296065;

        let d_params = compute_d_element_params(
            z, zs, zp, zd, gss, gsp, gpp, gp2, hsp, zsn, zpn, zdn, 0.0, 0.0, 0.0,
        );

        // Empirical check against OpenMOPAC v23.2.5 printed values:
        // DD2 = 1.03469697
        // DD3 = 1.30910036
        // DD4 = 0.85328607
        // DD5 = 0.40315988
        // DD6 = 0.56975041
        // PO1 = 1.55877040
        // PO2 = 0.83752232
        // PO3 = 1.22831702
        // PO4 = 0.71617642
        // PO5 = 1.12877683
        // PO6 = 0.54394491
        // PO7 = 1.55877040
        // PO8 = 0.84359471
        // PO9 = 1.55877040
        assert!(
            (d_params.ddp[2] - 1.03469697).abs() < 1e-6,
            "DD2 mismatch: {}",
            d_params.ddp[2]
        );
        assert!(
            (d_params.ddp[3] - 1.30910036).abs() < 1e-6,
            "DD3 mismatch: {}",
            d_params.ddp[3]
        );
        assert!(
            (d_params.ddp[4] - 0.85328607).abs() < 1e-6,
            "DD4 mismatch: {}",
            d_params.ddp[4]
        );
        assert!(
            (d_params.ddp[5] - 0.40315988).abs() < 1e-6,
            "DD5 mismatch: {}",
            d_params.ddp[5]
        );
        assert!(
            (d_params.ddp[6] - 0.56975041).abs() < 1e-6,
            "DD6 mismatch: {}",
            d_params.ddp[6]
        );

        assert!(
            (d_params.po[1] - 1.55877040).abs() < 1e-6,
            "PO1 mismatch: {}",
            d_params.po[1]
        );
        assert!(
            (d_params.po[2] - 0.83752232).abs() < 1e-6,
            "PO2 mismatch: {}",
            d_params.po[2]
        );
        assert!(
            (d_params.po[3] - 1.22831702).abs() < 1e-6,
            "PO3 mismatch: {}",
            d_params.po[3]
        );
        assert!(
            (d_params.po[4] - 0.71617642).abs() < 1e-6,
            "PO4 mismatch: {}",
            d_params.po[4]
        );
        assert!(
            (d_params.po[5] - 1.12877683).abs() < 1e-6,
            "PO5 mismatch: {}",
            d_params.po[5]
        );
        assert!(
            (d_params.po[6] - 0.54394491).abs() < 1e-6,
            "PO6 mismatch: {}",
            d_params.po[6]
        );
        assert!(
            (d_params.po[7] - 1.55877040).abs() < 1e-6,
            "PO7 mismatch: {}",
            d_params.po[7]
        );
        assert!(
            (d_params.po[8] - 0.84359471).abs() < 1e-6,
            "PO8 mismatch: {}",
            d_params.po[8]
        );
        assert!(
            (d_params.po[9] - 1.55877040).abs() < 1e-6,
            "PO9 mismatch: {}",
            d_params.po[9]
        );
    }
}