[SciPy-User] Get array of separation vectors from an array of vectors and use it to compute the force in a MD code

[josef.pktd at gmail.com]

On Thu, Jan 14, 2010 at 7:26 PM, davide_fiocco at yahoo.it
<davide_fiocco at yahoo.it> wrote:
> Thanks Josef!
> I post here the code i wrote to compute the matrix ff of the forces
> between all the pairs of particles in a given set interacting under
> the Lennard-Jones potential.
> Note that:
> - coordinates of the i-th particle is stored in self.txyz[i,1:].
> - the returned matrix ff contains at f[i,j,:] the three components of
> the force due to the interaction between i and j.
> - the for loop is the way I used to rebuild a triangular matrix from
> its reduced representation

When you are rebuilding the triu, or the full symmetric distance
matrix ff from the vectorized version then you can use again the
intial triu indices I,J, and inplace add the transpose.
might require a bit of thinking to get the 3rd axis right, but
something like this:

ff[I,J,:] = f      # unless numpy switches axis
ff += np.swapaxis(ff,2,1)  # diagonal is zero so not duplicate to worry about

You might want to try on a simple example, but I'm pretty sure
something like this should work

Josef

>
> I guess it can't be considered good code...and it'd be cool if someone
> could point out its major flaws!
> Thanks a lot again!
>
> Davide
>
>        def get_forces(self):
>                if self.pair_style == 'lj/cut':
>                        #Josef suggestion to get the reduced array of separation vectors R
>                        I, J = numpy.nonzero(numpy.triu(numpy.ones((self.natoms,
> self.natoms)), k=1))
>                        R = self.atoms.txyz[I,1:] - self.atoms.txyz[J,1:]
>                        #invoking a vectorized function to apply the
> minimum image convention to the separation vectors
>                        R = minimum_image(R, self.boxes[-1].bounds)
>                        #compute the array of inverse distances
>                        S = 1/numpy.sqrt(numpy.add.reduce((R*R).transpose()))

isn't the transpose here just choosing the axis ? 1/numpy.sqrt(((R*R).sum(0)))
it won't make much difference but I find it easier to read

>                        #in f I will store the information about the upper triangular part
> of the matrix of forces
>                        f = numpy.zeros((S.size, 3))
>                        invcut = 1./2.5
>                        #compute Lennard Jones force for distances
> below a given cutoff
>                        f[S > invcut, :] = (R[S > invcut,:])*((24.*(-2.*S[S > invcut]**13 +
> S[S > invcut]**7))*S[S > invcut]).reshape(-1,1)

you might want to replace  the repeated comparison with a temp
variable:  mask = S > invcut

Josef

>                        ff = numpy.zeros((self.natoms, self.natoms, 3))
>                        #convert reduced array of forces into an
> antisymmetric matrix ff (f contains all the information about its
> triu)
>                        for i in range(self.natoms):
>                                ff[i,i+1:,:] =  f[self.natoms*i - i*(i+1)/2:self.natoms*(i+1) - (i
> + 1)*(i + 2)/2,:]
>                                ff[i+1:,i,:] = -f[self.natoms*i - i*(i+1)/2:self.natoms*(i+1) - (i
> + 1)*(i + 2)/2,:]
>
>                        return ff
>        #apply the minimum image convention
>        def minimum_image_scalar(dx, box):
>                dx = dx - int(round(dx/box))*box
>                return dx
>        minimum_image = numpy.vectorize(minimum_image_scalar)
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>