Is there a better way to optimize the lennard jones potential function?

Question

In actual fact, its the derivative of the Lennard Jones potential. The reason for is that I am writing a Molecular Dynamics program and at least 80% of the time is spent in the following function, even with the most aggressive compiler options (gcc ** -O3).

double ljd(double r) /* Derivative of Lennard Jones Potential for Argon with 
                        respect to distance (r) */ 
{  
    double temp;  
    temp = Si/r;  
    temp = temp*temp;
    temp = temp*temp*temp;  
    return ( (24*Ep/r)*(temp-(2 * pow(temp,2))) );  
}  

This code is from a file "functs.h", which I import into my main file. I thought that using temporary variables in this way would make the function faster, but I am worried that creating them is too wasteful. Should I use static? Also the code is written in parallel using openmp, so I can't really declare temp as a global variable?

The variables Ep and Si are defined (using #define). I have only been using C for about 1 month. I tried to look at the assembler code generated by gcc, but I was completely lost.\

Answer 1

I would get rid of the call to pow() for a start:

double ljd(double r) /* Derivative of Lennard Jones Potential for Argon with 
                        respect to distance (r) */ 
{  
    double temp;  
    temp = Si / r;  
    temp = temp * temp;
    temp = temp * temp * temp;  
    return ( (24.0 * Ep / r) * (temp - (2.0 * temp * temp)) );  
}  

Answer 2

On my architecture (intel Centrino Duo, MinGW-gcc 4.5.2 on Windows XP), non-optimized code using pow()

static inline double ljd(double r)
{
    return 24 * Ep / Si * (pow(Si / r, 7) - 2 * pow(Si / r, 13));
}

actually outperforms your version if -ffast-math is provided.

The generated assembly (using some arbitrary values for Ep and Si) looks like this:

fldl    LC0
fdivl   8(%ebp)
fld     %st(0)
fmul    %st(1), %st
fmul    %st, %st(1)
fld     %st(0)
fmul    %st(1), %st
fmul    %st(2), %st
fxch    %st(1)
fmul    %st(2), %st
fmul    %st(0), %st
fmulp   %st, %st(2)
fxch    %st(1)
fadd    %st(0), %st
fsubrp  %st, %st(1)
fmull   LC1

Answer 3

Well, as I've said before, compilers suck at optimising floating point code for many reasons. So, here's an Intel assembly version that should be faster (compiled using DevStudio 2005):

const double Si6 = /*whatever pow(Si,6) is*/;
const double Si_value = /*whatever Si is*/; /* need _value as Si is a register name! */
const double Ep24 = /*whatever 24.Ep is*/;

double ljd (double r)
{
  double result;
  __asm
  {
    fld qword ptr [r]
    fld st(0)
    fmul st(0),st(0)
    fld st(0)
    fmul st(0),st(0)
    fmulp st(1),st(0)
    fld qword ptr [Si6]
    fdivrp st(1),st(0)
    fld st(0)
    fld1
    fsub st(0),st(1)
    fsubrp st(1),st(0)
    fmulp st(1),st(0)
    fld qword ptr [Ep24]
    fmulp st(1),st(0)
    fdivrp st(1),st(0)
    fstp qword ptr [result]
  }

  return result;
}

This version will produce slightly different results to the version posted. The compiler will probably be writing intermediate results to RAM in the original code. This will lose precision since the (Intel) FPU operates at 80bits internally whereas the double type is only 64bits. The above assembler will not lose precision in the intermediate results, it is all done at 80bits. Only the final result is rounded to 64bits.

Answer 4

Ah, that brings me back some memories... I've done MD with Lennard Jones potential years ago.

In my scenario (not huge systems) it was enough to replace the pow() with several multiplications, as suggested by another answer. I also restricted the range of neighbours, effectively truncating the potential at about r ~ 3.5 and applying some standard thermodynamic correction afterwards.

But if all this is not enough for you, I suggest to precompute the function for closely spaced values of r and simply interpolate (linear or quadratic, I'd say).

Answer 5

Is your application structured in such a way that you could profitably vectorise this function, calculating several independent values in parallel? This would allow you to utilise hardware vector units, such as SSE.

It also seems like you would be better off keeping 1/r values around, rather than r itself.

---

This is an example explicitly using SSE2 instructions to implement the function. ljd() calculates two values at once.

static __m128d ljd(__m128d r)
{
    static const __m128d two = { 2.0, 2.0 };
    static const __m128d si = { Si, Si };
    static const __m128d ep24 = { 24 * Ep, 24 * Ep };

    __m128d temp2, temp3;
    __m128d temp = _mm_div_pd(si, r);
    __m128d ep24_r = _mm_div_pd(ep24, r);

    temp = _mm_mul_pd(temp, temp);
    temp2 = _mm_mul_pd(temp, temp);
    temp2 = _mm_mul_pd(temp2, temp);
    temp3 = _mm_mul_pd(temp2, temp2);
    temp3 = _mm_mul_pd(temp3, two);

    return _mm_mul_pd(ep24_r, _mm_sub_pd(temp2, temp3));
}

/* Requires `out` and `in` to be 16-byte aligned */
void ljd_array(double out[], const double in[], int n)
{
    int i;

    for (i = 0; i < n; i += 2)
    {
        _mm_store_pd(out + i, ljd(_mm_load_pd(in + i)));
    }
}

However, it is important to note that recent versions of GCC are often able to vectorise functions like this automatically, as long as you're targetting the right architecture and have optimisation enabled. If you're targetting 32-bit x86, try compiling with -msse2 -O3, and adjust things such that the input and output arrays are 16-byte aligned.

Alignment for static and automatic arrays can be achieved under gcc with the type attribute __attribute__ ((aligned (16))), and for dynamic arrays using the posix_memalign() function.

Answer 6

The local variable is just fine. It doesn't cost anything. Leave it alone.

As others said, get rid of the pow call. It can't be any faster than simply squaring the number, and it could be a lot slower.

That said, just because the function is active 80+% of the time does not mean it's a problem. It only means if there is something you can optimize, it's either in there, or in something it calls (like pow) or in something that calls it.

If you try random pausing, which is a method of stack-sampling, you will see that routine on 80+% of samples, plus the lines within it that are responsible for the time, plus its callers that are responsible for the time, and their callers, and so on. All the lines of code on the stack are jointly responsible for the time.

Optimality is not when nothing take a large percent of time, it is when nothing you can fix takes a large percent of time.