[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 8:01 PM,  <josef.pktd at gmail.com> wrote:
> 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
typo, I think it's axis=1 in sum

And thanks for posting, it's nice to see whether my answers are helpful or not

Josef

>
>>                        #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
>>
>