## Python programing help!

This is the place for queries that don't fit in any of the other categories.

### Python programing help!

Helow every one.. I am new here.... Sory if crate the topic in a wrong direktory.

i Need some hlep for my project.

I work with molekylerdynamik simulations. .

Here is some eksamples; if eny one can this kind of stuff write me.

EKSEMPEL 1:
Code: Select all
`import numpy as npimport randomimport matplotlib.pyplot as plt import copyimport videodef distance(xi, yi, xj, yj):""" Calculate the distance between particle i and particle j """dx = xj-xidy = yj-yireturn np.sqrt(dx**2 + dy**2)def lennard_jones(x_positions, y_positions, n_particles):""" Calculate the force and energy based on the positions of the particles,in a Lennard-Jones potential"""energy = 0.0x_forces = [0.0 for _ in range(n_particles)]y_forces = [0.0 for _ in range(n_particles)]for i in range(n_particles):for j in range(n_particles):if i>j:# Distancerij = distance(x_positions[i],y_positions[i],x_positions[j],y_positions[j])# Energyenergy += 4*(1.0/rij**12 - 1.0/rij**6)# Forcefx = -48* (x_positions[j]-x_positions[i])/rij**2 * (1/rij**12 - 0.5 *1/rij**6)fy = -48* (y_positions[j]-y_positions[i])/rij**2 * (1/rij**12 - 0.5 *1/rij**6)x_forces[i] = x_forces[i] + fxy_forces[i] = y_forces[i] + fyx_forces[j] = x_forces[j] - fxy_forces[j] = y_forces[j] - fyreturn x_forces, y_forces, energydef initialize_particles(n_particles, box_width):""" Initialize particles in a grid position """sqrt_npart = int(np.ceil(np.sqrt(n_particles)))X = []Y = []for j in range(sqrt_npart):X += [i for i in range(sqrt_npart)]Y += [j for i in range(sqrt_npart)]# Remove excess particles.X = X[:n_particles]Y = Y[:n_particles]# Rescale particle positions to fit box.for i in range(n_particles):X[i] = (X[i] - 0.5*(sqrt_npart-1))*1.0/sqrt_npart*box_width*1.8#Y[i] = (Y[i] - 0.5*(sqrt_npart-1))*1.0/sqrt_npart*box_width*1.8Y[i] = (Y[i]) * 1.0/sqrt_npart*box_width*0.8# Initialize particle velocitiesVx = [2*(random.random() - 0.5) for i in range(n_particles)]Vy = [2*(random.random() - 0.5) for i in range(n_particles)]# Initialize particle forcesFx, Fy, energy = lennard_jones(X, Y, n_particles)return X, Y, Vx, Vy, Fx, Fy, energydef velo_verlet(x_positions, y_positions, x_velocities, y_velocities,x_forces, y_forces, box_width, n_particles, dt):""" Simulate particles movement for their positions and velocities in asingle time - step with time dt"""x_forces_old = copy.copy(x_forces)y_forces_old = copy.copy(y_forces)# Step 1:# Update positionsfor i in range(n_particles):x_positions[i] = x_positions[i] + dt*x_velocities[i] + 0.5*dt*dt*x_forces[i]y_positions[i] = y_positions[i] + dt*y_velocities[i] + 0.5*dt*dt*y_forces[i]if abs(x_positions[i]) > box_width:x_positions[i] -= dt*x_velocities[i]x_velocities[i] = -x_velocities[i]if abs(y_positions[i]) > box_width:y_positions[i] -= dt*y_velocities[i]y_velocities[i] = -y_velocities[i]# Step 2:# Update forcesx_forces, y_forces, energy = lennard_jones(x_positions, y_positions, n_particles)# Step 3:# Update velocitiesfor i in range(n_particles):x_velocities[i] = x_velocities[i] + 0.5*dt*(x_forces_old[i] + x_forces[i])y_velocities[i] = y_velocities[i] + 0.5*dt*(y_forces_old[i] + y_forces[i])return x_positions, y_positions, x_velocities, y_velocities, x_forces, y_forces, energydef plot_box(X, Y, Vx, Vy, box_width, filename):plt.plot(X, Y, 'ro')plt.quiver(X, Y, Vx, Vy)plt.axis((-box_width, box_width, -box_width, box_width))plt.savefig(filename)plt.clf()## Simulation Initializationbox_width = 10.0n_particles = 7*7n_particles = 42n_steps = 5000dt = 0.001# Empty energy listsenergy_list = []# x_test = [0.0, 0.0]# y_test = [0.0, 1.4]# print lennard_jones(x_test, y_test)## Initialize particlesX, Y, Vx, Vy, Fx, Fy, energy = initialize_particles(n_particles, box_width)## Take simulation stepsfor n in range(n_steps):X, Y, Vx, Vy, Fx, Fy, energy = velo_verlet(X, Y, Vx, Vy, Fx, Fy, box_width, n_particles, dt)if n == 0:plot_box(X, Y, Vx, Vy, box_width, 'first_step.png')plot_box(X, Y, Fx, Fy, box_width, 'first_step_f.png')if n % 1 == 0:energy_list.append(energy)if n % 10 == 0:video.add_frame(X, Y)if n % 100 == 0:print "step", n, "of", n_steps## Plot Results# Energy plot# plt.plot(energy_list)# plt.savefig('energy_time.png')# Save videovideo.save("week3video", box_width)`

EKSEMPEL 2:
Code: Select all
`import numpyimport timeimport copyfrom forces import lennardjones as fortran_ljdef initialize_positions(natms, rho):""" Initialize particle positions, based on number of atoms and density rho"""n = natmsp = int(numpy.floor(n**(1.0/3.0))) + 1box = (n/rho)**(1.0/3.0)du = box/pX = numpy.multiply.outer(numpy.ones(p), numpy.arange(p))X = numpy.multiply.outer(numpy.ones(p), X)Y = numpy.multiply.outer(numpy.arange(p), numpy.ones(p))Y = numpy.multiply.outer(numpy.ones(p), Y)Z = numpy.multiply.outer(numpy.arange(p), numpy.ones(p))Z = numpy.multiply.outer(Z, numpy.ones(p))U = numpy.reshape(numpy.ravel((X, Y, Z)), (3, p**3))U = U[:, :n] * du + du / 2return Udef initialize_box(n_atoms, rho):""" Initialize the box"""n = n_atomsp = int(numpy.floor(n**(1.0/3.0))) + 1box = (n/rho)**(1.0/3.0)du = box/pbox = numpy.array([box]*3)return boxdef initialize_velocities(U, T):""" Initialize random velocities"""# Calculate random velocities# Make the sum of velocities zerondim, n = numpy.shape(U)V = numpy.random.rand(ndim, n) - 0.5Vm = numpy.sum(V,1)/nV -= numpy.reshape(Vm, (ndim,1))# scale according to temperaturefs = numpy.dot( numpy.ravel(V), numpy.ravel(V) ) / n# the dot product of the velocities divided by the number of particlesfs = numpy.sqrt(ndim * T/fs)V *= fsreturn Vdef initialize_particles(n_atoms, temperature, rho):""" Initialize particles"""U = initialize_positions(n_atoms, rho)box = initialize_box(n_atoms, rho)V = initialize_velocities(U, temperature)epot, F, Vir = fortran_lj(U,box)return U, V, F, box[0]def lennard_jones(U, box):return fortran_lj(U, numpy.array([box,box,box]))`
Last edited by Yoriz on Tue May 13, 2014 12:02 pm, edited 1 time in total.
Reason: First post lock, Added code tags.
Thomas

Posts: 1
Joined: Tue May 13, 2014 11:48 am

### Re: Python programing help!

Hi welcome to the forum.
Did you have something specific you wanted help with?
Join the #python-forum IRC channel on irc.freenode.net!

Yoriz

Posts: 1465
Joined: Fri Feb 08, 2013 1:35 am
Location: UK