[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!
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()))
#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
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)
More information about the SciPy-User
mailing list