spherical_harmonics.cpp buffer overflow?

[aromanro]

Hi,

I'm trying to compile the project on windows and I noticed something:

In double SH::legendre (int n, double x) you have at line 58:

double v[n];

That would allow using indices for v from 0 to n-1. Despite that, v[n] is written and returned.

The simple solution would be: double v[n + 1]; (as you already have in SH::legendre_p).

Another issue is that the code as it is creates problems when compiling with VisualC++ (error because n is not a constant), The solution would be:

double *v = new double[n+1];

and just before return:

    const double val = v[n];

    delete[] v;

    return val;

and similarily in SH::legendre_p

PS Nice work. How long did it take you to implement it? I'm thinking of adding DFT to my Hartree-Fock program for molecules, I did implement some DFT code before, but with plane-waves basis: https://compphys.go.ro/dft-for-a-quantum-dot/

Adrian

[ifilot]

Hi Adrian,

Thank you for bringing this error to my attention! It is surprising that the code was working at all; I would have expected that it would have given a segmentation fault or something similar. Anyway, I have fixed it. In the end, I went for a more C++-style solution using std::vector<double> (which honestly I should have done in the first place), but your approach is of course equally valid.

Considering your other question. I guess I have spent about 1-2 years working on this code, but never full-time. I actually spent most time figuring out how to put the algorithm in the two papers of Becke into code. Cool to see that you wrote a PW-base DFT code. That is something that is also on my bucket list. :-)

I am closing this issue as it is fixed.

Thanks again! Ivo

[aromanro]

Hi,

Thank you for the info!

Probably I should take the simpler route first: compute for a single atom (that takes advantage of the spherical symmetry and reduces 3D to 1D), this is useful for computing pseudopotentials or going the APW/LAPW route with it (muffin tin).

Just a comment regarding the plane-waves project. It's a very good method to be used in crystals, but for many electrons/atom you need to use pseudopotentials. I was able for example to compute the H2O molecule without pseudopotentials, but it needs a lot of plane waves and it takes a lot of computing time. The solution is to use pseudopotentials, but that rises other issues. I implemented it with local pseudopotentials and I'm not happy with the results... using norm-conserving or ultrasoft pseudopotentials is out of the question right now for me, because I couldn't find an explicit enough description of files structure to use them and I don't have patience to figure out myself :)