git.psu.edu will be upgraded to critical security release 13.7.4 Monday, 11/18 between 9 and 10pm. Please email support@git.psu.edu if you have trouble with anything Gitlab-related. Please see the git.psu.edu Yammer group for more information.

Commit 95f6cc7d authored by Dundar Yilmaz's avatar Dundar Yilmaz

MPI VERSION

parent 1861e7db
# ------------------------------------------------------
#
# DipoleAnalyzer: Reax_Class
#
# ------------------------------------------------------
import reax_class as reax
import numpy as np
from timeit import default_timer as timer
import mmap
import sys
def new_read_data_adf(posfile_adf, bondfile_adf,
initial_frame, final_frame):
frames = []
f = open(posfile_adf)
start = timer()
s = mmap.mmap(f.fileno(), 0, access=mmap.ACCESS_READ)
br = timer()
q = 0
cnt = 0
cursor = 0
parts = []
timesteps = []
while cnt < final_frame:
line0 = s.readline()
try:
natoms = int(line0)
except:
print ('file error')
if cnt >= initial_frame and cnt < final_frame:
part = [None] * (natoms + 2)
part[0] = line0
part[1] = s.readline()
words = part[1].split()
timestep = int(words[1])
timesteps.append(timestep)
for i in range(natoms):
part[i + 2] = s.readline()
parts.append(part)
del part
else:
s.readline()
for i in range(natoms):
s.readline()
cnt += 1
ft = timer()
s.close()
del s
print ('Position file mapped in %6.4f seconds'
'parted in %6.4f seconds' % (ft-br, br-start))
f = open(bondfile_adf)
start = timer()
s = mmap.mmap(f.fileno(), 0, access=mmap.ACCESS_READ)
br = timer()
q = 0
cnt = 0
cursor = 0
bondparts = [None] * len(parts)
tq = 0
while cnt < final_frame:
shiftcursor = 50
initial = s.find('Iteration:', cursor)
if initial < shiftcursor:
shiftcursor = initial
if initial == -1:
break
cursor = initial + 10
final = s.find('Iteration:', cursor)
back = s[initial - shiftcursor:initial]
kk = back.find('\n') + 1
shiftcursor = shiftcursor - kk
part = s[initial - shiftcursor:final - shiftcursor]
cursor = final
if cnt >= initial_frame and cnt < final_frame:
lines = part.split('\n')
words = lines[0].split()
timestep = int(words[3])
if timestep < timesteps[tq]:
cnt -= 1
else:
bondparts[tq] = lines
tq += 1
del part
del lines
else:
del part
cnt += 1
ft = timer()
print ('Bonds file mapped in %6.4f '
'seconds parted in %6.4f seconds' % (ft-br, br-start))
s.close()
del s
for part, bondpart in zip(parts, bondparts):
natoms = int(part[0])
words = part[1].split()
xlo = float(0.0)
xhi = float(words[3])
ylo = float(0.0)
yhi = float(words[4])
zlo = float(0.0)
zhi = float(words[5])
timestep = int(words[1])
frame = reax.Frame(timestep)
frame.set_box([xlo, ylo, zlo], [xhi, yhi, zhi])
frame.natoms = natoms
for j in range(natoms):
words = part[2+j].split()
atomindex = j
symbol = words[0]
t = np.array(words[1:4])
coord = t.astype(float)
frame.atoms.append(reax.Atom(atomindex, symbol, coord))
chargeinfo = []
line = bondpart[0]
words = line.split()
timestep = int(words[3])
for j in range(natoms):
words = bondpart[1 + j].split()
index = int(words[0])
ch = float(words[-1])
frame.atoms[index - 1].charge = ch
frames.append(frame)
del frame
del parts
del bondparts
return timesteps, frames
def new_data_parter_lammps(posfile_lmp, bondfile_lmp,
initial_frame, final_frame):
frames = []
f = open(posfile_lmp)
start = timer()
s = mmap.mmap(f.fileno(), 0, access=mmap.ACCESS_READ)
br = timer()
q = 0
cnt = 0
cursor = 0
parts = []
timesteps = []
while cnt <= final_frame:
initial = s.find('ITEM: TIMESTEP', cursor)
if initial == -1:
break
cursor = initial + 14
final = s.find('ITEM: TIMESTEP', cursor)
part = s[initial:final]
cursor = final
if cnt >= initial_frame and cnt <= final_frame:
lines = part.split('\n')
timestep = int(lines[1])
parts.append(lines)
timesteps.append(timestep)
del part
del lines
else:
del part
cnt += 1
s.close()
ft = timer()
print ('Position file mapped in %6.4f seconds '
'parted in %6.4f seconds' % (ft-br, br-start))
f = open(bondfile_lmp)
start = timer()
s = mmap.mmap(f.fileno(), 0, access=mmap.ACCESS_READ)
br = timer()
q = 0
cnt = 0
cursor = 0
bondparts = [None] * len(parts)
tq = 0
while tq < len(timesteps):
initial = s.find('# Timestep', cursor)
if initial == -1:
break
cnt += 1
header = s[cursor:cursor + 20]
h = header.split('\n')
q = h[0]
del h
c = q.split()
del q
timestep = int(c[2])
del c
cursor = initial+10
final = s.find('# Timestep', cursor)
cursor = final
if timestep == timesteps[tq]:
part = s[initial:final]
lines = part.split('\n')
words = lines[0].split()
bondparts[tq] = lines
tq += 1
del part
del lines
ft = timer()
print ('Bonds file mapped in %6.4f seconds '
'parted in %6.4f seconds' % (ft-br, br-start))
return parts, bondparts, timesteps
def new_data_parter_lammps_charges(posfile_lmp,
initial_frame, final_frame):
frames = []
f = open(posfile_lmp)
start = timer()
s = mmap.mmap(f.fileno(), 0, access=mmap.ACCESS_READ)
br = timer()
q = 0
cnt = 0
cursor = 0
parts = []
timesteps = []
while cnt <= final_frame:
initial = s.find(b'ITEM: TIMESTEP', cursor)
if initial == -1:
break
cursor = initial + 14
final = s.find(b'ITEM: TIMESTEP', cursor)
part = s[initial:final]
cursor = final
if cnt >= initial_frame and cnt <= final_frame:
lines = part.split(b'\n')
timestep = int(lines[1])
parts.append(lines)
timesteps.append(timestep)
del part
del lines
else:
del part
cnt += 1
s.close()
ft = timer()
print ('Position file mapped in %6.4f seconds '
'parted in %6.4f seconds' % (ft-br, br-start))
return parts, timesteps
def new_data_reader_lammps_charges(datagroup):
st=timer()
frames=[]
for i,data in enumerate(datagroup):
part, timestep = data
natoms = int(part[3])
words = part[5].split()
xlo = float(words[0])
xhi = float(words[1])
words = part[6].split()
ylo = float(words[0])
yhi = float(words[1])
words = part[7].split()
zlo = float(words[0])
zhi = float(words[1])
timestep = int(part[1])
frame = reax.Frame(timestep)
frame.set_box([xlo, ylo, zlo], [xhi, yhi, zhi])
frame.natoms = natoms
for j in range(natoms):
words = part[9 + j].split()
atomindex = int(words[0])
symbol = words[1]
t = np.array(words[3:6])
coord = t.astype(float)
charge = float(words[2])
frame.atoms.append(reax.Atom(atomindex, symbol, coord, charge))
frames.append(frame)
return frames
def new_data_reader_lammps(datagroup):
st = timer()
frames = []
for i, data in enumerate(datagroup):
part, bondpart, timestep = data
natoms = int(part[3])
words = part[5].split()
xlo = float(words[0])
xhi = float(words[1])
words = part[6].split()
ylo = float(words[0])
yhi = float(words[1])
words = part[7].split()
zlo = float(words[0])
zhi = float(words[1])
timestep = int(part[1])
frame = reax.Frame(timestep)
frame.set_box([xlo, ylo, zlo], [xhi, yhi, zhi])
frame.natoms = natoms
#print 'FRAME : ',frame.boxlo[:],frame.boxhi[:]
for j in range(natoms):
words = part[9 + j].split()
atomindex = int(words[0])
symbol = words[1]
t = np.array(words[2:5])
coord = t.astype(float)
frame.atoms.append(reax.Atom(atomindex, symbol, coord))
for line in bondpart:
words = line.split()
if len(words) <= 1:
continue
if words[0] == '#':
if words[1] == 'Timestep':
timestep = int(words[2])
#print 'bond time ',timestep
if words[1] == 'Number':
natoms = int(words[4])
charges=np.zeros(natoms)
continue
#print i,timestep, words
index=int(words[0])
ch=float(words[-1])
frame.atoms[index-1].charge=ch
frames.append(frame)
laptime=timer()
if (i+1)%50 == 0:
print_progress(i+1,len(datagroup),(laptime-st))
#if i%50 == 0 :
# print 'Elapsed time %6.2f for %6d of %d frames. Remaining time: %6.2f '%\
# ((laptime-st),i+1,len(datagroup),float(len(datagroup)-(i+1))*(laptime-st)/float(i+1))
del frame
#for timestep,frame in zip(timesteps,frames):
# print timestep, frame.timestep
ft=timer()
print( '\n')
print( 'Positions and charges read ', ft-st)
return frames
def read_positions_lammps(posfile_lmp,initial_frame,final_frame):
frames=[]
# Parting Trajectory File
with open(posfile_lmp, 'r') as trajectory:
pstart=timer()
g=trajectory.read()
parts=g.split("ITEM: TIMESTEP")
pfinish=timer()
print ('There are %d frames '% len(parts))
print ('File read and splitted in %6.4f seconds' % (pfinish-pstart))
for i,part in enumerate(g.split("ITEM: TIMESTEP")):
if len(part) == 0 :
continue
if i > final_frame : continue
if i < initial_frame : continue
lines=part.split('\n')
lines[0]="ITEM: TIMESTEP\n"
timestep=int(lines[1])
natoms=int(lines[3])
#print timestep, natoms
words=lines[5].split()
xlo=float(words[0]);xhi=float(words[1])
words=lines[6].split()
ylo=float(words[0]);yhi=float(words[1])
words=lines[7].split()
zlo=float(words[0]);zhi=float(words[1])
frame=reax.Frame(timestep)
frame.set_box([xlo,ylo,zlo],[xhi,yhi,zhi])
frame.natoms=natoms
for j in range(natoms):
words=lines[9+j].split()
atomindex=int(words[0])
symbol=words[1]
t=np.array(words[2:5])
coord=t.astype(float)
frame.atoms.append(reax.Atom(atomindex,symbol,coord))
frames.append(frame)
del frame
rfinish=timer()
print ('Coordinates read in %6.4f seconds' % (rfinish-pfinish))
return frames
def read_bondinfo_lammps(bondsfile_lmp,initial_frame,final_frame):
with open(bondsfile_lmp, 'r') as bonds:
pstart=timer()
g=bonds.read()
chargeinfo=[]
timesteps=[]
parts=g.split("# Timestep")
pfinish=timer()
print ('There are %d frames in bonds file'% len(parts))
print ('File read and splitted in %6.4f seconds' % (pfinish-pstart))
for i,part in enumerate(g.split("# Timestep")):
if len(part) == 0 :
continue
if i > final_frame : continue
if i < initial_frame : continue
part="# Timestep"+part
lines=part.split('\n')
#timestep=int(lines[0])
#print lines[0]
for line in lines:
words=line.split()
#print words
if len(words) <= 1 : continue
if words[0] == '#' :
if words[1] == 'Timestep': timestep=int(words[2])
if words[1] == 'Number':
natoms=int(words[4])
charges=np.zeros(natoms)
continue
index=int(words[0])
ch=float(words[-1])
charges[index-1]=ch
#print timestep,natoms
chargeinfo.append(charges)
timesteps.append(timestep)
rfinish=timer()
print ('Bonds read in %6.4f seconds' % (rfinish-pfinish))
return timesteps,chargeinfo
def merge_charges_with_frames(frames,timesteps,charges):
nframes=len(frames)
#print nframes,len(timesteps)
lag=0
for i,frame in enumerate(frames):
if timesteps[i+lag] < frame.timestep : lag+=1
frame.set_charges(charges[i+lag])
#frame.calculate_dipole_moments()
#frame.calculate_total_polarization()
#print i,lag,frame.timestep,timesteps[i+lag]
#frame.print_atoms()
def print_progress(iteration, total,elapsed_time, prefix='', suffix='', decimals=1, bar_length=100):
"""
Call in a loop to create terminal progress bar
@params:
iteration - Required : current iteration (Int)
total - Required : total iterations (Int)
prefix - Optional : prefix string (Str)
suffix - Optional : suffix string (Str)
decimals - Optional : positive number of decimals in percent complete (Int)
bar_length - Optional : character length of bar (Int)
"""
str_format = "{0:." + str(decimals) + "f}"
percents = str_format.format(100 * (iteration / float(total)))
remaining_time=elapsed_time/float(iteration)*float(total-iteration)
filled_length = int(round(bar_length * iteration / float(total)))
bar = 'X' * filled_length + '-' * (bar_length - filled_length)
#print iteration,total,elapsed_time,remaining_time
sys.stdout.write('\r%s |%s| %s%s %s Remaining t: %6.2f ' % (prefix, bar, percents, '%', suffix,remaining_time)),
if iteration == total:
sys.stdout.write('\n')
sys.stdout.flush()
from timeit import default_timer as timer
import numpy as np
from periodic_kdtree import PeriodicCKDTree
import multiprocessing as mp
import time
from multiprocessing import Pool
from functools import partial
from numba import cuda
def batchquery(data):
T,vectors= data
neighbors=[]
for vec in vectors:
d,i = T.query(vec,k=8)
#print 'p', i
neighbors.append(i)
return neighbors
bits = 64
np_type = np.float64
@cuda.jit("void(float{}[:, :],float{}[:,:], float{}[:, :])".format(bits, bits, bits))
def distance_matrixAB(matA, matB, out):
m = matA.shape[0]
n = matB.shape[0]
print(n)
i, j = cuda.grid(2)
d = 0
if i < m and j < n:
for k in range(3):
tmp = matA[i, k] + matB[j, k]
d += tmp * tmp
out[i, j] = d
@cuda.jit("void(float{}[:, :], float{}[:, :])".format(bits, bits))
def distance_matrix(mat, out):
m = mat.shape[0]
n = mat.shape[1]
i, j = cuda.grid(2)
d = 0
if i < m and j < m:
for k in range(n):
tmp = mat[i, k] - mat[j, k]
d += tmp * tmp
out[i, j] = d
def gpu_dist_matrix(mat):
rows = mat.shape[0]
block_dim = (16, 16)
grid_dim = (int(rows/block_dim[0] + 1), int(rows/block_dim[1] + 1))
stream = cuda.stream()
mat2 = cuda.to_device(np.asarray(mat, dtype=np_type), stream=stream)
out2 = cuda.device_array((rows, rows))
distance_matrix[grid_dim, block_dim](mat2, out2)
out = out2.copy_to_host(stream=stream)
return out
def gpu_dist_matrixAB(matA,matB):
rows = matB.shape[0]
n = matB.shape[0]
m = matA.shape[0]
block_dim = (16, 16)
grid_dim = (int(rows/block_dim[0] + 1), int(rows/block_dim[1] + 1))
stream = cuda.stream()
mat2A = cuda.to_device(np.asarray(matA, dtype=np_type), stream=stream)
mat2B = cuda.to_device(np.asarray(matB, dtype=np_type), stream=stream)
out2 = cuda.device_array((n, m))
distance_matrix[grid_dim, block_dim](mat2A,mat2A, out2)
out = out2.copy_to_host(stream=stream)
return out
def find_neighbors(frames):
for frame in frames:
print (frame.boxlo, frame.boxhi)
bounds = frame.boxhi
frame = frames[0]
bounds = frame.boxhi
start = timer()
neighborlist = [None] * frame.natoms
amap = []
XA = []
XB = []
for i in range(frame.natoms):
XB.append(frame.atoms[i].coord)
if frame.atoms[i].symbol == b'Ti' :
XA.append(frame.atoms[i].coord)
amap.append(i)
XA = np.asarray(XA, dtype=np.float32)
XB = np.asarray(XB, dtype=np.float32)
print (XA.shape, XB.shape)
XChunk = XA[0:4,:]
XCB = XB[0:4,:]
start = timer()
#dists = -2*np.dot(XB,XChunk.T)+np.sum(XChunk**2,axis=1)+np.sum(XB**2,axis=1)[:,np.newaxis]
print (timer()-start)
start = timer()
#dist = gpu_dist_matrix(XChunk)
print (timer()-start)
dist2 = gpu_dist_matrixAB(XChunk,XCB)
#for d,d2 in zip(dist,dist2):
# print(d,d2)
print(dist2)
"""
nXA = XA.shape[0]
nXB = XB.shape[0]
start = timer()
T = PeriodicCKDTree(bounds,XB)
print( 'Tree built in %8.4f seconds'% (timer()-start))
start = timer()
neighborlist=[]
for i in range(40):
vec = XA[i,:]
d, i = T.query(vec,k=8)
neighborlist.append(i)
#print i
#print i
#print '----'
neighborlist = np.asarray(neighborlist, dtype=int)
print( 'Query completed %8.4f seconds'% (timer()-start))
start = timer()
nump=4
sliced_XA=[XA[i:i+nump,:] for i in range(0,40,nump)]
TS = [T for i in range(nump)]
data=[]
for A,B in zip(TS,sliced_XA):
d= A, B
data.append(d)
q=Pool(4)
s=q.map_async(batchquery,data)
q.close()
q.join()
data = s.get()
data = np.asarray(data,dtype=int)
data = np.reshape(data,(-1,8))
#print(data)
print( 'Query completed parallel %8.4f seconds'% (timer()-start))
for n,d in zip(neighborlist,data):
print 'Serial : ', n
print 'Paralel: ', d
# for n in neighborlist:
# veci = XA[n[0],:]
# vecj =
"""
#--------------------------------------------------------------------------------------------------
#
# DipoleAnalyzer: MAIN
#
#--------------------------------------------------------------------------------------------------
import numpy as np
import frame_reader as fr
import reax_class as reax
import sys
import os
import time
from timeit import default_timer as timer
from ktree_test import find_neighbors
def main():