Stack Overflow Asked by FuzzyFiso on December 20, 2021
I am trying to write some codes for my research project. I have 3 different arrays which are x
, y
, and z
. These together represent the positions of the atoms that belong to my system. When I start the simulation, the code prints the positions for every atom. I was able to do this with this code for a 6-atom system via this:
for i in range(number_of_steps):
x, z = RunMD(0.1, x, z, v_x, v_z, m)
print(num_particles)
print "Ni"
print "Ni", x[0], y[0], z[0]
print "Ni", x[1], y[1], z[1]
print "Ni", x[2], y[2], z[2]
print "Ni", x[3], y[3], z[3]
print "Ni", x[4], y[4], z[4]
print "Ni", x[5], y[5], z[5]
However, now I need to study with much more number of atoms. How can I handle this printing problem for 100-atom system? All the code is below. Thanks
import numpy as np
import math
'''
units: eV / A / ps
'''
def cart_to_polar(x,z):
rho = math.sqrt(x**2 + z**2)
phi = math.atan2(z, x)
return rho, phi
def polar_to_cart(rho, theta):
x = rho * math.cos(theta)
z = rho * math.sin(theta)
return x, z
def GetLJForce(r, epsilon, sigma):
return 48 * epsilon * np.power(
sigma, 12) / np.power(
r, 13) - 24 * epsilon * np.power(
sigma, 6) / np.power(r, 7)
def GetSpringForce(k, radius):
return -1 * k * radius
def GetAcc(x_positions, z_positions, m):
xAcc = [0]*len(x_positions)
zAcc = [0]*len(z_positions)
for i in range(0, len(x_positions)-1):
for j in range(i+1, len(x_positions)):
delta_x = x_positions[j] - x_positions[i]
delta_z = z_positions[j] - z_positions[i]
radius, theta = cart_to_polar(delta_x, delta_z)
if(radius == 0):
radius = 1
force_mag = GetLJForce(radius, 0.84, 2.56) + GetSpringForce(0, radius)
force_x, force_z = polar_to_cart(force_mag, theta)
xAcc[i] += force_x / m[i]
zAcc[i] += force_z / m[i]
force_x, force_z = polar_to_cart(-force_mag, theta)
xAcc[j] = force_x / m[j]
zAcc[j] = force_z / m[j]
return xAcc, zAcc
def UpdatePos(x, v, a, dt):
return x + v*dt + 0.5*a*dt**2
def UpdateVel(v, a, dt):
return v + a*dt
def RunMD(dt, x, z, v_x, v_z, m):
num_particles = len(x)
a_x, a_z = GetAcc(x, z, m)
for i in range(num_particles):
v_x[i] = UpdateVel(v_x[i], a_x[i], dt)
v_z[i] = UpdateVel(v_z[i], a_z[i], dt)
for i in range(num_particles):
x[i] = UpdatePos(x[i], v_x[i], a_x[i], dt)
z[i] = UpdatePos(z[i], v_z[i], a_z[i], dt)
if(z[1], z[2], z[3], z[4], z[5] != 0):
z[1] = 0
z[2] = 0
z[3] = 0
z[4] = 0
z[5] = 0
return x, z
num_particles = 6
m = [1] * num_particles
x = [4, 5, 10, 15, 20, 25]
y = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
z = [5, 0.0, 0.0, 0.0, 0.0, 0.0]
v_x = [0] * num_particles
v_z = [0] * num_particles
v_x[0] = 5
number_of_steps = 100
for i in range(number_of_steps):
x, z = RunMD(0.1, x, z, v_x, v_z, m)
print(num_particles)
print "Ni"
print "Ni", x[0], y[0], z[0]
print "Ni", x[1], y[1], z[1]
print "Ni", x[2], y[2], z[2]
print "Ni", x[3], y[3], z[3]
print "Ni", x[4], y[4], z[4]
print "Ni", x[5], y[5], z[5]
Here is a simple way to do it:
for i in range(number_of_steps):
x, z = RunMD(0.1, x, z, v_x, v_z, m)
print(num_particles)
print "Ni"
for j in range(len(x)):
print("Ni",x[j],y[j],z[j])
I hope this is helpful ;)
Answered by Hakim CHERIF on December 20, 2021
You can try the zip
operator
for z in zip(arr_1,arr_2,arr_3):
print(*z)
Answered by Spandan Brahmbhatt on December 20, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP