Python programing help!

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

Python programing help!

Postby Thomas » Tue May 13, 2014 11:55 am

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 np
import random
import matplotlib.pyplot as plt import copy

import video


def distance(xi, yi, xj, yj):
""" Calculate the distance between particle i and particle j """
dx = xj-xi
dy = yj-yi
return 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.0

x_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:

# Distance
rij = distance(x_positions[i],
y_positions[i],
x_positions[j],
y_positions[j])

# Energy
energy += 4*(1.0/rij**12 - 1.0/rij**6)

# Force
fx = -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] + fx
y_forces[i] = y_forces[i] + fy

x_forces[j] = x_forces[j] - fx
y_forces[j] = y_forces[j] - fy


return x_forces, y_forces, energy


def 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.8
Y[i] = (Y[i]) * 1.0/sqrt_npart*box_width*0.8


# Initialize particle velocities
Vx = [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 forces
Fx, Fy, energy = lennard_jones(X, Y, n_particles)

return X, Y, Vx, Vy, Fx, Fy, energy


def 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 a
single time - step with time dt
"""

x_forces_old = copy.copy(x_forces)
y_forces_old = copy.copy(y_forces)

# Step 1:
# Update positions
for 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 forces
x_forces, y_forces, energy = lennard_jones(x_positions, y_positions, n_particles)

# Step 3:
# Update velocities
for 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, energy


def 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 Initialization
box_width = 10.0
n_particles = 7*7
n_particles = 42
n_steps = 5000
dt = 0.001

# Empty energy lists
energy_list = []

# x_test = [0.0, 0.0]
# y_test = [0.0, 1.4]
# print lennard_jones(x_test, y_test)

## Initialize particles
X, Y, Vx, Vy, Fx, Fy, energy = initialize_particles(n_particles, box_width)

## Take simulation steps
for 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 video
video.save("week3video", box_width)

EKSEMPEL 2:
Code: Select all
import numpy
import time
import copy


from forces import lennardjones as fortran_lj


def initialize_positions(natms, rho):
""" Initialize particle positions, based on number of atoms and density rho
"""
n = natms
p = int(numpy.floor(n**(1.0/3.0))) + 1
box = (n/rho)**(1.0/3.0)
du = box/p
X = 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 / 2
return U


def initialize_box(n_atoms, rho):
""" Initialize the box
"""
n = n_atoms
p = int(numpy.floor(n**(1.0/3.0))) + 1
box = (n/rho)**(1.0/3.0)
du = box/p
box = numpy.array([box]*3)
return box


def initialize_velocities(U, T):
""" Initialize random velocities
"""
# Calculate random velocities
# Make the sum of velocities zero
ndim, n = numpy.shape(U)
V = numpy.random.rand(ndim, n) - 0.5
Vm = numpy.sum(V,1)/n
V -= numpy.reshape(Vm, (ndim,1))

# scale according to temperature
fs = numpy.dot( numpy.ravel(V), numpy.ravel(V) ) / n

# the dot product of the velocities divided by the number of particles
fs = numpy.sqrt(ndim * T/fs)
V *= fs

return V


def 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!

Postby Yoriz » Tue May 13, 2014 12:03 pm

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.
User avatar
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: Google [Bot] and 7 guests