How to write it in the most effective way

Question

I have an array (N = 10 ^ 4), and I need to find the difference between each of the two entries (calculating the potential taking into account the coordinates of the atoms) Here is the code that I write in pure python, but its really inefficient, can anyone tell me how to speed it up? (using numpy or weaving). Here x, y are the atomic coordinates arrays (just a simple 1D array)

def potential(r):
   U = 4.*(np.power(r,-12) - np.power(r,-6))
   return U
def total_energy(x):
   E = 0.
   #need to speed up this part
   for i in range(N-1):
     for j in range(i):
         E += potential(np.sqrt((x[i]-x[j])**2))  
return E

Answer 1

you can use array arithmetic first:

def potential(r):
    return 4.*(r**(-12) - r**(-6))

def total_energy(x):
    E = 0.
    for i in range(N-1):
        E += potential(np.sqrt((x[i]-x[:i])**2)).sum()
    return E

or you can test the fully vectorized version:

def total_energy(x):
    b=np.diag(x).cumsum(1)-x
    return potential(abs(b[np.triu_indices_from(b,1)])).sum()

Answer 2

scipy.spatial.distance. pdist, , .

, (Nx3), :

def potential(r):
       U = 4.*(np.power(r,-12) - np.power(r,-6))
       return U
def total_energy(x):
   E = 0.
   #need to speed up this part
   for i in range(N):                                    #To N here
     for j in range(i):
         E += potential(np.sqrt(np.sum((x[i]-x[j])**2))) #Add sum here
   return E

:

import scipy.spatial.distance as sd


def scipy_LJ(arr, sigma=None):
    """
    Computes the Lennard-Jones potential for an array (M x N) of M points
    in N dimensional space. Usage of a sigma parameter is optional.
    """

    if len(arr.shape)==1:
        arr = arr[:,None]

    r = sd.pdist(arr)

    if sigma==None:
        np.power(r, -6, out=r)
        return np.sum(r**2 - r)*4

    else:
       r *= sigma
       np.power(r, -6, out=r)
       return np.sum(r**2 - r)*4

:

N = 1000
points = np.random.rand(N,3)+0.1

np.allclose(total_energy(points), scipy_LJ(points))
Out[43]: True

%timeit total_energy(points)
1 loops, best of 3: 13.6 s per loop

%timeit scipy_LJ(points)
10 loops, best of 3: 24.3 ms per loop

~ 500 !

N = 10000
points = np.random.rand(N,3)+0.1

%timeit scipy_LJ(points)
1 loops, best of 3: 3.05 s per loop

~ 2 .

Answer 3

0) ( )

In [16]: %timeit total_energy(points)
1 loops, best of 3: 14.9 s per loop

1) SciPy

In [9]: %timeit scipy_LJ(points)
10 loops, best of 3: 44 ms per loop

1-2)

 %timeit sum( potential(np.sqrt((x[i]-x[:i])**2 + (y[i]-y[:i])**2 + (z[i] - z[:i])**2)).sum() for i in range(N-1))
10 loops, best of 3: 126 ms per loop

2) Fortran (! - )

    subroutine EnergyForces(Pos, PEnergy, Dim, NAtom)
    implicit none
    integer, intent(in) :: Dim, NAtom
    real(8), intent(in), dimension(0:NAtom-1, 0:Dim-1) :: Pos
!    real(8), intent(in) :: L
    real(8), intent(out) :: PEnergy
    real(8), dimension(Dim) :: rij, Posi
    real(8) :: d2, id2, id6, id12
    real(8) :: rc2, Shift
    integer :: i, j
    PEnergy = 0.
    do i = 0, NAtom - 1
        !store Pos(i,:) in a temporary array for faster access in j loop
        Posi = Pos(i,:)
        do j = i + 1, NAtom - 1
            rij = Pos(j,:) - Posi
!            rij = rij - L * dnint(rij / L)
            !compute only the squared distance and compare to squared cut
            d2 = sum(rij * rij)
            id2 = 1. / d2            !inverse squared distance
            id6 = id2 * id2 * id2    !inverse sixth distance
            id12 = id6 * id6         !inverse twelvth distance
            PEnergy = PEnergy + 4. * (id12 - id6)
      enddo
    enddo
end subroutine

In [14]: %timeit ljlib.energyforces(points.transpose(), 3, N)
10000 loops, best of 3: 61 us per loop

3) Fortran 1000 , scipy 3000 , numpy, , python. , Scipy , , Fortran .

Answer 4

. .

• return sum( potential(np.sqrt((x[i]-x[:i])**2)).sum() for i in range(N-1))
  

• .

• , , - f2py, Fortran ( ), python , : program_lj.f90

$gfortran -c program_lj.f90

fortran, .

$f2py -c -m program_lj program_lj.f90

python. python:

import program_lj
result = program_lj.subroutine_in_program(parameters)

, -.