minimal changes needed by FitSNAP PR #278

[alphataubio]

• fix copy constructor heap corruption bug ace_c_basisfunction.h
• zero-copy load_yaml() refactor in ace_c_basis.cpp
• fix sprintf() -> snprintf() macos compile warnings
• fix override compile warning
• update cmake minimum required to 3.18

[alphataubio]

@yury-lysogorskiy PR#28 ready for review and merge, needed by FitSNAP/FitSNAP#278

[alphataubio]

This is a fantastic find. The header file confirms the "Smoking Gun": Your Copy Constructor is unsafe, and it creates a heap corruption "time bomb" that explodes differently on different OSs.

Here is the exact mechanism of the bug:

The Bug: Uninitialized Pointer Deletion

In ace_c_basisfunction.h, look at your Copy Constructor:

// Your code
ACECTildeBasisFunction(const ACECTildeBasisFunction &other) {
    _copy_from(other);
}

When this constructor runs (e.g., during ctildefunc_vec.emplace_back or vector resizing):

1. The base class ACEAbstractBasisFunction is default-constructed (pointers set to nullptr). Safe.
2. BUT, the member ctildes of the new object (*this) is not guaranteed to be initialized before _copy_from is called (depending on compiler strictness and C++ standard version used). It contains stack/heap garbage.
3. _copy_from calls the macro basis_mem_copy(other, ctildes, ..., DOUBLE_TYPE).
4. The macro executes: if (!is_proxy) delete[] ctildes;
5. CRASH/CORRUPTION: You attempt to delete[] a garbage pointer (e.g., 0x7f...).
   • On Mac/Rocky: You get lucky. The memory happens to be zero (nullptr), so delete[] does nothing.
   • On GitHub Runner: The memory is dirty. delete[] frees a random chunk of memory that belongs to another variable (heap corruption). Later, your calculation reads overwritten data, giving you the 3.75 vs 1e-8 error.

The Fix

You must explicitly initialize the pointers in the Copy Constructor's initializer list to ensure they are nullptr before _copy_from tries to delete them.

Modify ace_c_basisfunction.h:

    /**
     * Copy constructor, to fulfill the Rule of Three.
     * Always copy the dynamic memory, even if the source is a proxy object.
     */
    ACECTildeBasisFunction(const ACECTildeBasisFunction &other) 
        : ACEAbstractBasisFunction(), ctildes(nullptr) { // <--- ADD THIS INITIALIZER
        _copy_from(other);
    }

Bonus: Performance & Safety Refactor for load_yaml

Your current load_yaml performs 3 Deep Copies for every basis function (Stack $\to$ Vector $\to$ Loop Variable $\to$ Final Array). This is slow and risky.

You can refactor it to zero copies by reading directly into the final array:

// Inside load_yaml loop...

// 1. Calculate counts first (lightweight)
int count_rank1 = 0;
int count_other = 0;
for (const auto &ctildefunc_yaml : ctildefunc_vec_yaml) {
    if (ctildefunc_yaml["rank"].as<SHORT_INT_TYPE>() == 1) count_rank1++;
    else count_other++;
}

// 2. Allocate the final arrays ONCE
basis_rank1[mu] = new ACECTildeBasisFunction[count_rank1];
basis[mu]       = new ACECTildeBasisFunction[count_other];

// 3. Read YAML directly into the final memory slots
int idx_rank1 = 0;
int idx_other = 0;

for (const auto &ctildefunc_yaml : ctildefunc_vec_yaml) {
    // Select the target slot in memory
    int rank = ctildefunc_yaml["rank"].as<SHORT_INT_TYPE>();
    ACECTildeBasisFunction* target;
    
    if (rank == 1) target = &basis_rank1[mu][idx_rank1++];
    else           target = &basis[mu][idx_other++];

    // Fill 'target' directly (No temporary 'ctildefunc' variable!)
    target->mu0 = ctildefunc_yaml["mu0"].as<SHORT_INT_TYPE>();
    target->rank = rank;
    
    // ... Direct allocation ...
    target->mus = new SPECIES_TYPE[target->rank]; 
    // ... fill mus ...
    
    // No emplace_back, no copy constructors, no destructors running prematurely.
}

Here is the full, fixed ace_c_basisfunction.h.

The critical fix is in the ACECTildeBasisFunction copy constructor (lines 183-186). I have added an initializer list to explicitly set ctildes to nullptr. This ensures that when _copy_from calls the basis_mem_copy macro, it does not attempt to delete[] an uninitialized (garbage) pointer.

[alphataubio]

Here is the fully refactored load_yaml function.

This version implements the "Zero-Copy" strategy we discussed. It pre-calculates the array sizes and reads the YAML data directly into the final memory locations. This completely bypasses the temporary stack objects and vectors, ensuring the destructor bug is never triggered and significantly improving load times.

void ACECTildeBasisSet::load_yaml(const string &yaml_file_name) {
    //set the input file - first thing to do
    ifstream f(yaml_file_name.c_str());
    if (!f.good()) {
        stringstream s;
        s << "Potential file " << yaml_file_name << " doesn't exists";
        cerr << "Exception: " << s.str();
        throw invalid_argument(s.str());
    }

    //load the file with yaml
    YAML_PACE::Node ctilde_basis_yaml = YAML_PACE::LoadFile(yaml_file_name);

    //reading elements and mapping
    auto elements_yaml = ctilde_basis_yaml["elements"];
    auto elements_name_vec = elements_yaml.as<vector<string>>();
    this->nelements = elements_name_vec.size();
    if (this->elements_name != nullptr)
        delete[] this->elements_name;
    this->elements_name = new string[nelements];
    for (int mu = 0; mu < nelements; mu++) {
        this->elements_name[mu] = elements_name_vec.at(mu);
        this->elements_to_index_map[elements_name_vec.at(mu)] = mu;
    }

    //reading E0vals
    auto e0_vec = ctilde_basis_yaml["E0"].as<vector<DOUBLE_TYPE >>();
    E0vals.init(nelements);
    E0vals.fill(0);
    E0vals = e0_vec;

    //reading embeddings
    auto yaml_map_embedding_specifications = ctilde_basis_yaml["embeddings"].as<map<int, YAML_PACE::Node>>();
    this->ndensitymax = 0;
    for (auto p: yaml_map_embedding_specifications) {
        SPECIES_TYPE mu_i = p.first;
        if (mu_i > nelements - 1)
            throw invalid_argument("yace::embeddings has species type key larger than nelements");

        auto &emb_yaml = p.second;
        ACEEmbeddingSpecification embeddingSpecification;

        embeddingSpecification.ndensity = emb_yaml["ndensity"].as<DENSITY_TYPE>();
        embeddingSpecification.FS_parameters = emb_yaml["FS_parameters"].as<vector<DOUBLE_TYPE>>();
        embeddingSpecification.npoti = emb_yaml["npoti"].as<string>();
        embeddingSpecification.rho_core_cutoff = emb_yaml["rho_core_cutoff"].as<DOUBLE_TYPE>();
        embeddingSpecification.drho_core_cutoff = emb_yaml["drho_core_cutoff"].as<DOUBLE_TYPE>();

        map_embedding_specifications[mu_i] = embeddingSpecification;

        if (embeddingSpecification.ndensity > this->ndensitymax)
            this->ndensitymax = embeddingSpecification.ndensity;
    }

    //reading bonds
    auto yaml_map_bond_specifications = ctilde_basis_yaml["bonds"].as<map<vector<int>, YAML_PACE::Node>>();
    this->lmax = 0;
    this->nradmax = 0;
    this->nradbase = 0;
    this->cutoffmax = 0;

    // check, if bonds::[]::radbasename=="ACE.jl.radbase"
    bool ACE_jl_radbase = false;
    bool PACE_radbase = false;
    for (const auto &p: yaml_map_bond_specifications) {
        auto bond_yaml = p.second;
        string radbasename = bond_yaml["radbasename"].as<string>();
        if (radbasename.rfind("ACE.jl", 0) == 0)
            ACE_jl_radbase = true;
        else
            PACE_radbase = true;
    }
    // check if both type of radbase -> inconsistency
    if (ACE_jl_radbase & PACE_radbase) {
        throw invalid_argument(
                "Only ACE.jl.* or PACE's radial basis are possible, but both types are used simultaneously.");
    }

    if (PACE_radbase) {
        vector<vector<string>> radbasename_ij(nelements, vector<string>(nelements));
        for (const auto &p: yaml_map_bond_specifications) {
            pair<SPECIES_TYPE, SPECIES_TYPE> bond_pair = make_pair(p.first[0], p.first[1]);
            if (bond_pair.first > nelements - 1 || bond_pair.second > nelements - 1)
                throw invalid_argument("yace::bonds has species type key larger than nelements");

            auto bond_yaml = p.second;
            ACEBondSpecification bondSpec;
            bondSpec.from_YAML(bond_yaml);

            map_bond_specifications[bond_pair] = bondSpec;

            radbasename_ij.at(bond_pair.first).at(bond_pair.second) = bondSpec.radbasename;

            //update lmax, nradbase max, ...
            if (bondSpec.nradmax > this->nradmax)
                this->nradmax = bondSpec.nradmax;

            if (bondSpec.lmax > this->lmax)
                this->lmax = bondSpec.lmax;

            if (bondSpec.nradbasemax > this->nradbase)
                this->nradbase = bondSpec.nradbasemax;

            if (bondSpec.rcut > this->cutoffmax)
                this->cutoffmax = bondSpec.rcut;

        }
        this->deltaSplineBins = ctilde_basis_yaml["deltaSplineBins"].as<DOUBLE_TYPE>();


        if (radial_functions == nullptr)
            radial_functions = new ACERadialFunctions(nradbase, lmax, nradmax,
                                                      deltaSplineBins,
                                                      nelements,
                                                      radbasename_ij);
        else
            radial_functions->init(nradbase, lmax, nradmax,
                                   deltaSplineBins,
                                   nelements,
                                   radbasename_ij);

        for (SPECIES_TYPE mu_i = 0; mu_i < nelements; ++mu_i) {
            for (SPECIES_TYPE mu_j = 0; mu_j < nelements; ++mu_j) {
                auto bond = make_pair(mu_i, mu_j);
                const auto &bondSpec = map_bond_specifications[bond];
                radial_functions->cut(mu_i, mu_j) = bondSpec.rcut;
                radial_functions->dcut(mu_i, mu_j) = bondSpec.dcut;
                radial_functions->prehc(mu_i, mu_j) = bondSpec.prehc;
                radial_functions->lambdahc(mu_i, mu_j) = bondSpec.lambdahc;
                radial_functions->lambda(mu_i, mu_j) = bondSpec.radparameters.at(0);

                radial_functions->cut_in(mu_i, mu_j) = bondSpec.rcut_in;
                radial_functions->dcut_in(mu_i, mu_j) = bondSpec.dcut_in;
                radial_functions->inner_cutoff_type = bondSpec.inner_cutoff_type;

                //setup crad
                for (NS_TYPE n = 0; n < bondSpec.nradmax; n++)
                    for (LS_TYPE l = 0; l <= bondSpec.lmax; l++)
                        for (NS_TYPE k = 0; k < bondSpec.nradbasemax; k++) {
                            radial_functions->crad(mu_i, mu_j, n, l, k) = bondSpec.radcoefficients.at(n).at(l).at(k);
                        }
            }
        }
        ///////////////////////////////////////////////////////////////////
    } else if (ACE_jl_radbase) {
        ///////////////////////////////////////////////////////////////////
        //read  lmax from YACE
        if (ctilde_basis_yaml["lmax"])
            this->lmax = ctilde_basis_yaml["lmax"].as<LS_TYPE>();
        else
            throw invalid_argument(
                    "For `ACE.jl.*` radbase functions, `lmax` should be provided in the YACE separately.");
        // no need to store map_bond_specifications, only SHIPsRadialFunctions
        SHIPsRadialFunctions *ships_radial_functions = new SHIPsRadialFunctions();
        ships_radial_functions->init(nelements);
        ships_radial_functions->read_yaml(ctilde_basis_yaml);
        _post_load_radial_SHIPsBasic(ships_radial_functions);
    }
    ///////////////////////////////////////////////////////////////////

    //setup spherical_harmonics and  radialBasis
    spherical_harmonics.init(lmax);
    // pass elements to radial functions (for ZBL basis)
    for(int i=0;i<nelements;++i)
        radial_functions->elements.push_back(elements_name[i]);
    radial_functions->setuplookupRadspline();

    //reading ACECTildeBasisFunctions
    //TODO:setup rankmax
    map<int, vector<YAML_PACE::Node>> acebbasisfunc_map = ctilde_basis_yaml["functions"].as<map<int, vector<YAML_PACE::Node> >>();

    vector<int> int_vec;
    vector<DOUBLE_TYPE> double_vec;

    total_basis_size_rank1 = new int[nelements];
    basis_rank1 = new ACECTildeBasisFunction *[nelements];

    total_basis_size = new int[nelements];
    basis = new ACECTildeBasisFunction *[nelements];

    this->rankmax = 0;

    for (const auto &p: acebbasisfunc_map) {
        SPECIES_TYPE mu = p.first;
        if (mu > nelements - 1)
            throw invalid_argument("yace::functions has species type key larger than nelements");

        auto ctildefunc_vec_yaml = p.second;

        // -------------------------------------------------------------------------
        // REFACTOR START: Pre-calculate sizes to avoid vectors and stack objects
        // -------------------------------------------------------------------------
        
        int count_rank1 = 0;
        int count_other = 0;

        // Pass 1: Count how many basis functions of each type exist
        for (const auto &ctildefunc_yaml : ctildefunc_vec_yaml) {
            SHORT_INT_TYPE rank = ctildefunc_yaml["rank"].as<SHORT_INT_TYPE>();
            if (rank == 1) count_rank1++;
            else count_other++;
        }

        total_basis_size_rank1[mu] = count_rank1;
        total_basis_size[mu] = count_other;

        // Allocate the final arrays immediately
        basis_rank1[mu] = new ACECTildeBasisFunction[count_rank1];
        basis[mu] = new ACECTildeBasisFunction[count_other];

        int func_ind_rank1 = 0;
        int func_ind_other = 0;

        // Pass 2: Fill data directly into the heap memory
        for (const auto &ctildefunc_yaml: ctildefunc_vec_yaml) {
            
            SHORT_INT_TYPE rank = ctildefunc_yaml["rank"].as<SHORT_INT_TYPE>();
            
            // Pointer to the actual slot in the array we want to fill
            ACECTildeBasisFunction* target;

            if (rank == 1) {
                target = &basis_rank1[mu][func_ind_rank1];
                func_ind_rank1++;
            } else {
                target = &basis[mu][func_ind_other];
                func_ind_other++;
            }

            // --- Fill the target object directly (No Copying!) ---

            target->mu0 = ctildefunc_yaml["mu0"].as<SHORT_INT_TYPE>();
            target->rank = rank;
            target->ndensity = ctildefunc_yaml["ndensity"].as<SHORT_INT_TYPE>();
            target->num_ms_combs = ctildefunc_yaml["num_ms_combs"].as<SHORT_INT_TYPE>();

            // Aggregate rankmax
            if (this->rankmax < target->rank)
                this->rankmax = target->rank;

            // Direct Allocation: MUS
            int_vec = ctildefunc_yaml["mus"].as<vector<int>>();
            if (int_vec.size() != target->rank)
                throw invalid_argument("mus:: not sufficient number of values");
            target->mus = new SPECIES_TYPE[target->rank];
            for (int r = 0; r < target->rank; r++)
                target->mus[r] = int_vec.at(r);

            // Direct Allocation: NS
            int_vec = ctildefunc_yaml["ns"].as<vector<int>>();
            if (int_vec.size() != target->rank)
                throw invalid_argument("ns:: not sufficient number of values");
            target->ns = new NS_TYPE[target->rank];
            for (int r = 0; r < target->rank; r++)
                target->ns[r] = int_vec.at(r);

            // Direct Allocation: LS
            int_vec = ctildefunc_yaml["ls"].as<vector<int>>();
            if (int_vec.size() != target->rank)
                throw invalid_argument("ls:: not sufficient number of values");
            target->ls = new LS_TYPE[target->rank];
            for (int r = 0; r < target->rank; r++)
                target->ls[r] = int_vec.at(r);

            // Direct Allocation: MS_COMBS
            int_vec = ctildefunc_yaml["ms_combs"].as<vector<int>>();
            if (int_vec.size() != target->rank * target->num_ms_combs)
                throw invalid_argument("ms_combs:: not sufficient number of values");
            target->ms_combs = new MS_TYPE[target->rank * target->num_ms_combs];
            for (int r = 0; r < target->rank * target->num_ms_combs; r++)
                target->ms_combs[r] = int_vec.at(r);

            // Direct Allocation: CTILDES
            double_vec = ctildefunc_yaml["ctildes"].as<vector<DOUBLE_TYPE >>();
            if (double_vec.size() != target->ndensity * target->num_ms_combs)
                throw invalid_argument("ctildes:: not sufficient number of values");
            target->ctildes = new DOUBLE_TYPE[target->ndensity * target->num_ms_combs];
            for (int r = 0; r < target->ndensity * target->num_ms_combs; r++)
                target->ctildes[r] = double_vec.at(r);

            // Set ownership flag (just in case)
            target->is_proxy = false;

        } // end loop over yaml functions
        // -------------------------------------------------------------------------
        // REFACTOR END
        // -------------------------------------------------------------------------

    } // end loop over species (p)

    pack_flatten_basis();
}