Skip to content

Commit

Permalink
submit align code
Browse files Browse the repository at this point in the history
  • Loading branch information
@ensemble-learning
ensemble-learning committed Mar 10, 2017
1 parent 121383d commit a2d7c43
Show file tree
Hide file tree
Showing 24 changed files with 3,049 additions and 0 deletions.
133 changes: 133 additions & 0 deletions examples/04_align_axis/align.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,133 @@
"""
"""

import math
import numpy as np
import scipy.linalg

DEBUG = 0

def read_coords(fname):
f = open(fname, "r")
lines = f.readlines()
coords = []
for i in range(2, len(lines)):
tokens = lines[i].strip().split()
if len(tokens) > 0:
coords.append([float(j) for j in tokens[1:]])
return coords

def get_surface_atoms(coords):
sur = []
for i in range(32):
sur.append(coords[i])
return sur

def get_com(coords):
com = [0.0, 0.0, 0.0]
for i in range(len(coords)):
for j in range(3):
com[j] += coords[i][j]
for i in range(3):
com[i] = com[i]/len(coords)

return com

def center_coords(coords, s):
coords_new = []
for i in range(len(coords)):
coords_new.append([])
for j in range(3):
coords_new[i].append(coords[i][j] - s[j])
return coords_new

def fit_plannar(coords):
# regular grid covering the domain of the data
data = np.array(coords)
X,Y = np.meshgrid(np.arange(-10.0, 10.0, 1), np.arange(-10.0, 10.0, 1))
# best-fit linear plane
A = np.c_[data[:,0], data[:,1], np.ones(data.shape[0])]
C,_,_,_ = scipy.linalg.lstsq(A, data[:,2]) # coefficients
if DEBUG:
o = open("plannar.xyz", "w")
o.write("%d\n\n"%(len(coords)*2))
for i in range(len(coords)):
x = coords[i][0]
y = coords[i][1]
z0 = coords[i][2]
z1 = C[0]*x + C[1]*y + C[2]
o.write("Au %12.4f%12.4f%12.4f\n"%(x, y, z0))
o.write("Fe %12.4f%12.4f%12.4f\n"%(x, y, z1))
o.close()
return C

def align_axis(v1, v2):
w = np.cross(v1,v2)
w = w/np.linalg.norm(w)
w_hat = np.matrix([[ 0, -w[2], w[1]],
[ w[2], 0, -w[0]],
[-w[1], w[0], 0]])
cos_tht = np.dot(v1, v2)/np.linalg.norm(v1)/np.linalg.norm(v2)
tht = np.arccos(cos_tht)
R = np.identity(3) + w_hat*np.sin(tht) + w_hat*w_hat*(1-np.cos(tht))
return R

def rotate_coords(coords, R):
coords_new = []
for i in range(len(coords)):
coords_new.append(np.dot(R, coords[i]).A1)
return coords_new

def output_coords(coords, fname):
o = open(fname, "w")
o.write("%d\n\n"%(len(coords)+20))
for i in range(len(coords)):
o.write("Au ")
for j in range(3):
o.write("%12.4f"%coords[i][j])
o.write("\n")
for i in range(20):
o.write("O 0 0 %12.4f\n"%(i-10))
o.close()

def output_coords_debug(coords, fname, pc):
o = open(fname, "w")
o.write("%d\n\n"%(len(coords)+20))
for i in range(len(coords)):
o.write("Au ")
for j in range(3):
o.write("%12.4f"%coords[i][j])
o.write("\n")
for i in range(20):
x = pc[0]*i
y = pc[1]*i
z = -i
o.write("O %12.4f%12.4f%12.4f\n"%(x, y, z))
o.close()

def main():
# read the xyz file
coords = read_coords("input.xyz")
# get the surface atoms
sur = get_surface_atoms(coords)
# calculate the center of mass (com)
com = get_com(sur)
# move the com to zero
sur_zero = center_coords(sur, com)
# determine the surface and surface vector
pc = fit_plannar(sur_zero)
v1 = [pc[0], pc[1], -1]
v0 = [0, 0, 1] # rotate to z axis
coords_zero = center_coords(coords, com)
if DEBUG:
output_coords_debug(coords_zero, "vector.xyz", pc)
# calculate the rotation matrix
R = align_axis(v1, v0)
# apply the rotation
coords_new = rotate_coords(coords_zero, R)
# output the data
output_coords(coords_zero, "t0.xyz")
output_coords(coords_new, "t1.xyz")

if __name__ == "__main__":
main()
84 changes: 84 additions & 0 deletions examples/04_align_axis/input.xyz
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,84 @@
82
Nanocluster
Pt 4.303200 10.058001 5.519400
Pt 5.532200 9.051200 3.205000
Pt 1.845400 12.071800 10.148200
Pt 3.074400 11.065001 7.833800
Pt 5.597800 11.509000 8.980600
Pt 6.826800 10.502000 6.666200
Pt 8.055600 9.495200 4.351800
Pt 9.284600 8.488199 2.037400
Pt 11.808000 8.932200 3.184200
Pt 9.350201 10.945999 7.813000
Pt 10.579200 9.939200 5.498600
Pt 13.102799 10.383200 6.645400
Pt 14.331599 9.376200 4.331000
Pt 16.855200 9.820200 5.477800
Pt 3.140000 13.522600 13.609400
Pt 4.369000 12.515800 11.295000
Pt 6.892400 12.959800 12.441799
Pt 8.121400 11.953000 10.127400
Pt 5.663600 13.966599 14.756201
Pt 8.187200 14.410600 15.903000
Pt 9.416000 13.403799 13.588600
Pt 10.710600 14.854601 17.049801
Pt 10.645000 12.396801 11.274199
Pt 11.873799 11.390000 8.959800
Pt 14.397400 11.834001 10.106601
Pt 15.626201 10.827200 7.792200
Pt 18.149799 11.271000 8.939000
Pt 11.939600 13.847800 14.735399
Pt 13.168401 12.840800 12.421000
Pt 15.692000 13.284800 13.567801
Pt 16.920799 12.278000 11.253400
Pt 14.463000 14.291800 15.882200
Pt 6.147400 7.950400 5.713000
Pt 3.689600 9.964200 10.341801
Pt 4.918600 8.957399 8.027400
Pt 7.442000 9.401400 9.174200
Pt 7.376400 6.943600 3.398600
Pt 9.899799 7.387600 4.545400
Pt 8.671000 8.394400 6.859800
Pt 11.194400 8.838400 8.006600
Pt 12.423400 7.831600 5.692200
Pt 14.947000 8.275600 6.839000
Pt 2.460800 10.971201 12.656200
Pt 4.984200 11.415200 13.803000
Pt 6.213200 10.408200 11.488600
Pt 8.736600 10.852200 12.635401
Pt 7.507800 11.859200 14.949800
Pt 10.031400 12.302999 16.096600
Pt 9.965600 9.845400 10.320999
Pt 12.489200 10.289400 11.467800
Pt 13.717999 9.282400 9.153400
Pt 16.241600 9.726399 10.300200
Pt 17.470400 8.719600 7.985800
Pt 11.260200 11.296200 13.782200
Pt 13.783799 11.740200 14.929000
Pt 15.012600 10.733400 12.614599
Pt 12.554800 12.747001 17.243401
Pt 5.533800 7.856600 10.535601
Pt 6.762800 6.849800 8.221200
Pt 7.991600 5.842800 5.906800
Pt 9.286200 7.293800 9.368000
Pt 10.515200 6.286800 7.053600
Pt 13.038599 6.730800 8.200399
Pt 4.305000 8.863600 12.849999
Pt 6.828400 9.307600 13.996600
Pt 8.057400 8.300600 11.682400
Pt 9.351999 9.751600 15.143400
Pt 10.580799 8.744600 12.829000
Pt 11.809799 7.737800 10.514801
Pt 14.333400 8.181800 11.661400
Pt 15.562201 7.174800 9.347000
Pt 11.875401 10.195600 16.290199
Pt 13.104400 9.188601 13.975801
Pt 7.378000 5.749000 10.729199
Pt 8.607000 4.742200 8.414801
Pt 11.130400 5.186200 9.561600
Pt 6.149200 6.756000 13.043600
Pt 8.672600 7.200000 14.190399
Pt 9.901600 6.193000 11.876000
Pt 12.424999 6.637000 13.022799
Pt 13.653999 5.630200 10.708400
Pt 11.196199 7.644000 15.337200
104 changes: 104 additions & 0 deletions examples/04_align_axis/output/t0.xyz
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,104 @@
102

Au -5.6944 -1.6134 -4.0242
Au -4.4654 -2.6203 -6.3386
Au -8.1522 0.4003 0.6046
Au -6.9232 -0.6064 -1.7098
Au -4.3998 -0.1625 -0.5630
Au -3.1708 -1.1695 -2.8774
Au -1.9420 -2.1763 -5.1918
Au -0.7130 -3.1833 -7.5062
Au 1.8104 -2.7393 -6.3594
Au -0.6474 -0.7255 -1.7306
Au 0.5816 -1.7323 -4.0450
Au 3.1052 -1.2883 -2.8982
Au 4.3340 -2.2953 -5.2126
Au 6.8576 -1.8513 -4.0658
Au -6.8576 1.8511 4.0658
Au -5.6286 0.8443 1.7514
Au -3.1052 1.2883 2.8982
Au -1.8762 0.2815 0.5838
Au -4.3340 2.2951 5.2126
Au -1.8104 2.7391 6.3594
Au -0.5816 1.7323 4.0450
Au 0.7130 3.1832 7.5062
Au 0.6474 0.7254 1.7306
Au 1.8762 -0.2815 -0.5838
Au 4.3998 0.1626 0.5630
Au 5.6286 -0.8443 -1.7514
Au 8.1522 -0.4005 -0.6046
Au 1.9420 2.1763 5.1918
Au 3.1708 1.1693 2.8774
Au 5.6944 1.6133 4.0242
Au 6.9232 0.6065 1.7098
Au 4.4654 2.6203 6.3386
Au -3.8502 -3.7211 -3.8306
Au -6.3080 -1.7073 0.7982
Au -5.0790 -2.7141 -1.5162
Au -2.5556 -2.2701 -0.3694
Au -2.6212 -4.7279 -6.1450
Au -0.0978 -4.2839 -4.9982
Au -1.3266 -3.2771 -2.6838
Au 1.1968 -2.8331 -1.5370
Au 2.4258 -3.8399 -3.8514
Au 4.9494 -3.3959 -2.7046
Au -7.5368 -0.7002 3.1126
Au -5.0134 -0.2563 4.2594
Au -3.7844 -1.2633 1.9450
Au -1.2610 -0.8193 3.0918
Au -2.4898 0.1877 5.4062
Au 0.0338 0.6315 6.5530
Au -0.0320 -1.8261 0.7774
Au 2.4916 -1.3821 1.9242
Au 3.7204 -2.3891 -0.3902
Au 6.2440 -1.9451 0.7566
Au 7.4728 -2.9519 -1.5578
Au 1.2626 -0.3753 4.2386
Au 3.7862 0.0687 5.3854
Au 5.0150 -0.9381 3.0710
Au 2.5572 1.0756 7.6998
Au -4.4638 -3.8149 0.9920
Au -3.2348 -4.8217 -1.3224
Au -2.0060 -5.8287 -3.6368
Au -0.7114 -4.3777 -0.1756
Au 0.5176 -5.3847 -2.4900
Au 3.0410 -4.9407 -1.3432
Au -5.6926 -2.8079 3.3064
Au -3.1692 -2.3639 4.4530
Au -1.9402 -3.3709 2.1388
Au -0.6456 -1.9199 5.5998
Au 0.5832 -2.9269 3.2854
Au 1.8122 -3.9337 0.9712
Au 4.3358 -3.4897 2.1178
Au 5.5646 -4.4967 -0.1966
Au 1.8778 -1.4759 6.7466
Au 3.1068 -2.4828 4.4322
Au -2.6196 -5.9225 1.1856
Au -1.3906 -6.9293 -1.1288
Au 1.1328 -6.4853 0.0180
Au -3.8484 -4.9155 3.5000
Au -1.3250 -4.4715 4.6468
Au -0.0960 -5.4785 2.3324
Au 2.4274 -5.0345 3.4792
Au 3.6564 -6.0413 1.1648
Au 1.1986 -4.0275 5.7936
O 0 0 -10.0000
O 0 0 -9.0000
O 0 0 -8.0000
O 0 0 -7.0000
O 0 0 -6.0000
O 0 0 -5.0000
O 0 0 -4.0000
O 0 0 -3.0000
O 0 0 -2.0000
O 0 0 -1.0000
O 0 0 0.0000
O 0 0 1.0000
O 0 0 2.0000
O 0 0 3.0000
O 0 0 4.0000
O 0 0 5.0000
O 0 0 6.0000
O 0 0 7.0000
O 0 0 8.0000
O 0 0 9.0000
Loading

0 comments on commit a2d7c43

Please sign in to comment.