## 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.
Please read the new users read this link in my signature.
Did you have something specific you wanted help with?
Due to the reasons discussed here we will be moving to python-forum.io/ on October 1 2016
This forum will be locked down and no one will be able to post/edit/create threads, etc. here from thereafter. Please create an account at the new site to continue discussion.

Yoriz

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

Return to General Coding Help

### Who is online

Users browsing this forum: No registered users and 7 guests