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

davide_fiocco@yah... davide_fiocco@yah...
Thu Jan 14 18:26:33 CST 2010

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

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!


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

More information about the SciPy-User mailing list