#Simple md. Force is soft-sphere: V(r)=e*(s/r)**12 #Units: m=1, e=1, s=1 from __future__ import division, print_function from visual import * from visual.graph import * def forces(ux,uy,N): #compute forces on all the particles fx=zeros((N,N)) fy=zeros((N,N)) rx=ux-ux[:,newaxis] #This obscure command gives a matrix rx[i,j] = ux[i]-ux[j] ry=uy-uy[:,newaxis] d=rx*rx+ry*ry + eye(N) #Add and subtract indentity matrix to avoid d[i,i]=0. dinv=1./d -eye(N) fx=12*rx*dinv**7 fy=12.*ry*dinv**7 fx= sum(fx,axis=0) #sum rows of force matrix -- total force on i due to all j fy = sum(fy,axis=0) return fx,fy maxsteps=100000 #Set up quantities needed dt=0.002 N=25 #number of particles L=12 #Box size t=0. Nsqrt=ceil(sqrt(N)) ux=zeros(N) uy=zeros(N) vx=zeros(N) vy=zeros(N) atoms=[] for ij in arange(N): #populate a lattice i,j=divmod(ij,Nsqrt) ux[ij]=1.5*i -L/2.+1 uy[ij]=1.5*j -L/2.+1 atoms=atoms+[sphere(pos=(ux[ij],uy[ij],0), color=color.yellow,radius=0.5)] vx[0]=1. #initial conditions: one particle moves vy[0]=3. ax,ay=forces(ux,uy,N) for m in arange(maxsteps): #main loop: all the computation is done here rate(10000) t=t+dt #verlet update. ux=ux+vx*dt+ax*dt*dt/2. uy=uy+vy*dt+ay*dt*dt/2. vx=vx+ax*dt/2. vy=vy+ay*dt/2. ax,ay=forces(ux,uy,N) vx=vx+ax*dt/2. vy=vy+ay*dt/2. for i in arange(N): atoms[i].pos=(ux[i],uy[i],0.) if abs(ux[i]) > L/2.-.5: #bounce off walls vx[i]=-vx[i] if abs(uy[i]) >L/2.-.5: vy[i]=-vy[i]