Implement orbital-based calculation

[JaeHyeokLee98]

copy cioverlap.pyx to cioverlap_orb.pyx
copy tdnac.c to tdnac_orb.c

[skmin-lab]

@gemini review

[gemini-code-assist]

Code Review

This pull request introduces an orbital-based calculation by adding cioverlap_orb.pyx and tdnac_orb.c, which appear to be copies of existing files. However, the new C code in tdnac_orb.c contains several critical issues, including incorrect function signatures, calls to undefined functions, and inconsistencies between function prototypes, implementations, and calls. These issues will prevent the code from compiling. The Cython wrapper cioverlap_orb.pyx also needs to be updated to correctly interface with the new C code once the C code is fixed. I have left detailed comments on these critical issues. Additionally, there's a potential for memory leaks in the Cython file and a critical bug combined with a performance issue in a numerical loop in the C code.

[gemini-code-assist]

This extern block has several issues:

1. It refers to "tdnac.c", but it should probably be "tdnac_orb.c" to use the new orbital-based implementation.
2. The function declared is TD_NAC, but the function in tdnac_orb.c is named TD_NAC_ORB.
3. The signature of TD_NAC declared here does not match the intended signature of TD_NAC_ORB in tdnac_orb.c (which itself has issues).

These issues need to be fixed in conjunction with the issues in tdnac_orb.c. After fixing tdnac_orb.c, this block should be updated to correctly declare TD_NAC_ORB.

[gemini-code-assist]

The function call is to TD_NAC, but it should be TD_NAC_ORB to match the new C implementation in tdnac_orb.c. This needs to be updated after fixing the cdef extern block and the C function itself.

    TD_NAC_ORB(istep, nst, nbasis, norb, nocc, nvirt, dt, orb_ini, orb_final, nacme, \
        ao_overlap, mo_coef_old, mo_coef_new, ci_coef_old, ci_coef_new)

[gemini-code-assist]

The function prototype for calc_MO_over does not match its implementation (line 278) or its usage (line 41). The prototype has incorrect parameter types and names (nst instead of norb, ph_fac instead of permut_mat) and an extra parameter l_spin_restricted that is not in the implementation. This will cause compilation errors.

Suggested change

	static void calc_MO_over(int nbasis, int nst, double **mo_overlap, double *ph_fac,
	double **ao_overlap, double **mo_coef_old, double **mo_coef_new, int l_spin_restricted);
	static void calc_MO_over(int nbasis, int norb, double **mo_overlap, double **permut_mat,
	double **ao_overlap, double **mo_coef_old, double **mo_coef_new);

[gemini-code-assist]

The function prototype for MO_phase_order is inconsistent with its usage on line 85. The prototype specifies (int nbasis, int nst, double **mo_coef_new, double *ph_fac), but the call passes norb for nst and a double** (permut_mat) for ph_fac (double*). This will cause a compilation error. The prototype should be updated to match the intended usage and the implementation.

Suggested change

	static void MO_phase_order(int nbasis, int nst, double **mo_coef_new, double *ph_fac);
	static void MO_phase_order(int nbasis, int norb, double **mo_coef_new, double **permut_mat);

[gemini-code-assist]

The function TD_NAC_ORB is missing several parameters in its signature which are used within its body. Variables such as norb, nocc, nvirt, orb_ini, orb_final, ci_coef_old, and ci_coef_new are used but not declared or passed as arguments. This will lead to compilation errors. The signature should be updated to include these missing parameters. Also, the parameter l_spin_restricted is declared but not used within the function.

Suggested change

	static void TD_NAC_ORB(int istep, int nst, int nbasis, double dt,
	double **nacme, double **ao_overlap, double **mo_coef_old, double **mo_coef_new, int l_spin_restricted){
	static void TD_NAC_ORB(int istep, int nst, int nbasis, int norb, int nocc, int nvirt, double dt,
	int *orb_ini, int *orb_final, double **nacme, double **ao_overlap, double **mo_coef_old, double **mo_coef_new, double ***ci_coef_old, double ***ci_coef_new, int l_spin_restricted){

[gemini-code-assist]

This function calls CI_phase_order, but this function is not defined in this file (tdnac_orb.c). This will cause linker errors. The same issue exists for state_phase (line 133) and norm_CI_coef (lines 148, 150). It seems that parts of tdnac.c were not copied over to tdnac_orb.c. Please ensure all necessary functions are defined or included.

[gemini-code-assist]

The calculation of tmp_sign is incorrect. It is being determined based on the partial sum of mo_overlap[iorb][jorb] inside the innermost loops. The sign should be determined only after the full value of mo_overlap[iorb][jorb] has been computed, i.e., after the loops over ibasis and jbasis are complete.

Additionally, the four nested loops result in O(N⁴) complexity, which is inefficient. This can be optimized to O(N³) by using intermediate matrices, similar to what the LAPACK/MKL implementation does with dgemm.

    double **tmp_mat = malloc(norb * sizeof(double*));
    for(iorb = 0; iorb < norb; iorb++){
        tmp_mat[iorb] = malloc(nbasis * sizeof(double));
    }

    // tmp_mat = mo_coef_old * ao_overlap
    for(iorb = 0; iorb < norb; iorb++){
        for(jbasis = 0; jbasis < nbasis; jbasis++){
            tmp_mat[iorb][jbasis] = 0.0;
            for(ibasis = 0; ibasis < nbasis; ibasis++){
                tmp_mat[iorb][jbasis] += mo_coef_old[iorb][ibasis] * ao_overlap[ibasis][jbasis];
            }
        }
    }

    // mo_overlap = tmp_mat * mo_coef_new^T
    for(iorb = 0; iorb < norb; iorb++){
        for(jorb = 0; jorb < norb; jorb++){
            mo_overlap[iorb][jorb] = 0.0;
            for(jbasis = 0; jbasis < nbasis; jbasis++){
                mo_overlap[iorb][jorb] += tmp_mat[iorb][jbasis] * mo_coef_new[jorb][jbasis];
            }
            // Correctly determine sign after full calculation
            if(mo_overlap[iorb][jorb] < 0.0){
                tmp_sign[iorb][jorb] = -1.0;
            }
            else{
                tmp_sign[iorb][jorb] = 1.0;
            }
        }
    }

    for(iorb = 0; iorb < norb; iorb++){
        free(tmp_mat[iorb]);
    }
    free(tmp_mat);

[gemini-code-assist]

The function signature of MO_phase_order is incorrect. It uses variables norb and permut_mat in its body, but they are not passed as parameters. This will cause compilation errors. The signature should be updated to accept these variables to match the prototype and usage.

Suggested change

	static void MO_phase_order(int nbasis, int nst, double **mo_coef_new, double *ph_fac){
	static void MO_phase_order(int nbasis, int norb, double **mo_coef_new, double **permut_mat){

[gemini-code-assist]

The manual memory management with PyMem_Malloc and PyMem_Free is prone to memory leaks. If any operation between allocation and deallocation raises an exception (for example, in the TD_NAC call), the PyMem_Free calls at the end of the function will be skipped, leading to a memory leak. It is highly recommended to wrap the code that can raise exceptions in a try...finally block to ensure that memory is always deallocated.

Example:

# ... allocations ...
try:
    # Assign NACME variables from python to C
    # ...
    # Calculate TDNAC term for CIoverlap
    TD_NAC(...)
    # Assign NACME variables from C to python
    # ...
finally:
    # Deallocate NACME variables
    PyMem_Free(orb_ini)
    # ... all other PyMem_Free calls

This would make the code more robust against memory leaks.