0

I have a while loop that continues as long as energy variable (type double) has not converged to below a certain threshold. One of the variables needed to calculate this energy is an Armadillo matrix of doubles, named f_mo. In the while loop, this f_mo updates iteratively, so I calculate f_mo at the beginning of each loop as:

    arma::mat f_mo = h_core_mo;//h_core_mo is an Armadillo matrix of doubles
    for(size_t p = 0; p < n_mo; p++) {//n_mo is of type size_t
    for(size_t q = 0; q < n_mo; q++) {
        double sum = 0.0;
            for(size_t i = 0; i < n_occ; i++) {//n_occ is of type size_t
                //f_mo(p,q) += 2.0*g_mat_full_qp1_qp1_mo(p*n_mo + q, i*n_mo + i)-g_mat_full_qp1_qp1_mo(p*n_mo+i,i*n_mo+q); //all g_mat_ are Armadillo matrices of doubles
                sum += 2.0*g_mat_full_qp1_qp1_mo(p*n_mo + q, i*n_mo + i)-g_mat_full_qp1_qp1_mo(p*n_mo+i,i*n_mo+q);
            }
            for(size_t i2 = 0; i2 < n_occ2; i2++) {//n_occ2 is of type size_t
                //f_mo(p,q) -= 1.0*g_mat_full_qp1_qp2_mo(p*n_mo + q, i2*n_mo2 + i2);
                sum -= 1.0*g_mat_full_qp1_qp2_mo(p*n_mo + q, i2*n_mo2 + i2);
            }
        
        f_mo(p,q) +=sum;
    }}

But say I replace the sum (which I add at the end to f_mo(p,q)) with addition to f_mo(p,q) directly (the commented out code). The output f_mo matrices are identical to machine precision. Nothing about the code should change. The only variables affected in the loop are sum and f_mo. And YET, the code converges to a different energy and in vastly different number of while loop iterations. I am at a loss as to the cause of the difference. When I run the same code 2,3,4,5 times, I get the same result every time. When I recompile with no optimization, I get the same issue. When I run on a different computer (controlling for environment), I yet again get a discrepancy in # of while loops despite identical f_mo, but the total number of iterations for each method (sum += and f_mo(p,q) += ) differ.

It is worth noting that the point at which the code outputs differ is always g_mat_full_qp1_qp2_mo, which is recalculated later in the while loop. HOWEVER, every variable going into the calculation of g_mat_full_qp1_qp2_mo is identical between the two codes. This leads me to think there is something more profound about C++ that I do not understand. I welcome any ideas as to how you would proceed in debugging this behavior (I am all but certain it is not a typical bug, and I've controlled for environment and optimization)

user213305
  • 402
  • 2
  • 10
jezzo
  • 193
  • 10
  • This could be a precision error. `sum` starts out as zero, so adding small magnitude numbers to 0 will keep the full precision. Because `f_mo(p,q)` has a value in it, adding small magnitude numbers can result in a loss of precision because the lower precision bits get lost. See [this question](https://stackoverflow.com/questions/6699066/in-which-order-should-floats-be-added-to-get-the-most-precise-result). – 1201ProgramAlarm Apr 18 '21 at 21:02
  • @1201ProgramAlarm Thank you for the tip. I should add that the first iteration where the energy values differ, they differ on the order of 10^-10 (and in the previous iteration they match to 10^-50). Please correct me if my interpretation is wrong, but it therefore seems unlikely to me that this is caused by a precision error. – jezzo Apr 18 '21 at 21:21
  • 2
    Without your code and a [mre] it is only possible to speculate on the cause. – 1201ProgramAlarm Apr 18 '21 at 21:27

1 Answers1

1

I'm going to assume this a Hartree-Fock, or some other kind of electronic structure calculation where you adding the two-electron integrals to the core Hamiltonian, and apply some domain knowledge.

Part of that assume is the individual elements of the two-electron integrals are very small, in particular compared to the core Hamiltonian. Hence as 1201ProgramAlarm mentions in their comment, the order of addition will matter. You will get a more accurate result if you add smaller numbers together first to avoid loosing precision when adding two numbers many orders of magintude apart.. Because you iterate this processes until the Fock matrix f_mo has tightly converged you eventually converge to the same value.

In order to add up the numbers in a more accurate order, and hopefully converge faster, most electronic structure programs have a seperate routine to calculate the two-electron integrals, and then add them to the core Hamiltonian, which is what you are doing, element by element, in your example code.

Presentation on numerical computing.

user213305
  • 402
  • 2
  • 10