Skip to content

Commit

Permalink
Merge pull request #1190 from nisse3000/temp_q_sq_face_cap_trig_pris
Browse files Browse the repository at this point in the history 
Added a new order paramter (q_sq_face_cap_trig_pris) and unit test
  • Loading branch information
@shyuep
shyuep committed Jul 7, 2018
2 parents c86cc6a + 9cbebb8 commit 74392fa
Show file tree
Hide file tree
Showing 3 changed files with 87 additions and 13 deletions.
59 changes: 47 additions & 12 deletions pymatgen/analysis/local_env.py
View file
Expand Up @@ -33,7 +33,7 @@
__status__ = "Production"
__date__ = "August 17, 2017"

from math import pow, pi, asin, atan, sqrt, exp, sin, cos, acos, fabs
from math import pow, pi, asin, atan, sqrt, exp, sin, cos, acos, fabs, atan2
import numpy as np

try:
Expand Down Expand Up @@ -1573,7 +1573,8 @@ class LocalStructOrderParams(object):
"sq_plan_max", "pent_plan", "pent_plan_max", "sq", "tet", "tet_max", "tri_pyr", \
"sq_pyr", "sq_pyr_legacy", "tri_bipyr", "sq_bipyr", "oct", \
"oct_legacy", "pent_pyr", "hex_pyr", "pent_bipyr", "hex_bipyr", \
"T", "cuboct", "cuboct_max", "see_saw_rect", "bcc", "q2", "q4", "q6", "oct_max", "hex_plan_max")
"T", "cuboct", "cuboct_max", "see_saw_rect", "bcc", "q2", "q4", "q6", "oct_max", \
"hex_plan_max", "sq_face_cap_trig_pris")

def __init__(self, types, parameters=None, cutoff=-10.0):
"""
Expand Down Expand Up @@ -1684,6 +1685,7 @@ def __init__(self, types, parameters=None, cutoff=-10.0):
t + ")!")
self._types = tuple(types)

self._comp_azi = False
self._params = []
for i, t in enumerate(self._types):
d = deepcopy(default_op_params[t]) if default_op_params[t] is not None \
Expand All @@ -1708,8 +1710,11 @@ def __init__(self, types, parameters=None, cutoff=-10.0):
"sq_plan", "pent_plan", "tri_pyr", "pent_pyr", "hex_pyr",
"pent_bipyr", "hex_bipyr", "T", "cuboct", "oct_max", "tet_max",
"tri_plan_max", "sq_plan_max", "pent_plan_max", "cuboct_max",
"bent", "see_saw_rect", "hex_plan_max"]):
"bent", "see_saw_rect", "hex_plan_max",
"sq_face_cap_trig_pris"]):
self._computerijs = self._geomops = True
if "sq_face_cap_trig_pris" in self._types:
self._comp_azi = True
if not set(self._types).isdisjoint(["reg_tri", "sq"]):
self._computerijs = self._computerjks = self._geomops2 = True
if not set(self._types).isdisjoint(["q2", "q4", "q6"]):
Expand Down Expand Up @@ -2387,12 +2392,18 @@ def get_order_parameters(self, structure, n, indices_neighs=None, \
tmp = max(
-1.0, min(np.inner(zaxis, rijnorm[k]), 1.0))
thetak = acos(tmp)
xaxistmp = gramschmidt(rijnorm[k], zaxis)
if np.linalg.norm(xaxistmp) < very_small:
xaxis = gramschmidt(rijnorm[k], zaxis)
if np.linalg.norm(xaxis) < very_small:
flag_xaxis = True
else:
xaxis = xaxistmp / np.linalg.norm(xaxistmp)
xaxis = xaxis / np.linalg.norm(xaxis)
flag_xaxis = False
if self._comp_azi:
flag_yaxis = True
yaxis = np.cross(zaxis, xaxis)
if np.linalg.norm(yaxis) > very_small:
yaxis = yaxis / np.linalg.norm(yaxis)
flag_yaxis = False

# Contributions of j-i-k angles, where i represents the
# central atom and j and k two of the neighbors.
Expand Down Expand Up @@ -2449,6 +2460,12 @@ def get_order_parameters(self, structure, n, indices_neighs=None, \
qsptheta[i][j][kc] += (self._params[i]['w_SPP'] *
exp(-0.5 * tmp * tmp))
norms[i][j][kc] += self._params[i]['w_SPP']
elif t == "sq_face_cap_trig_pris":
if thetak < self._params[i]['TA3']:
tmp = self._params[i]['IGW_TA1'] * (
thetak * ipi - self._params[i]['TA1'])
qsptheta[i][j][kc] += exp(-0.5 * tmp * tmp)
norms[i][j][kc] += 1

for m in range(nneigh):
if (m != j) and (m != k) and (not flag_xaxis):
Expand All @@ -2465,7 +2482,11 @@ def get_order_parameters(self, structure, n, indices_neighs=None, \
-1.0,
min(np.inner(xtwoaxis, xaxis), 1.0)))
flag_xtwoaxis = False

if self._comp_azi:
phi2 = atan2(
np.dot(xtwoaxis, yaxis),
np.dot(xtwoaxis, xaxis))
#print('{} {}'.format(180*phi/pi, 180*phi2/pi))
# South pole contributions of m.
if t in ["tri_bipyr", "sq_bipyr", "pent_bipyr",
"hex_bipyr", "oct_max", "sq_plan_max",
Expand Down Expand Up @@ -2587,13 +2608,27 @@ def get_order_parameters(self, structure, n, indices_neighs=None, \
qsptheta[i][j][kc] += exp(-0.5 * tmp * tmp) * \
exp(-0.5 * tmp2 * tmp2)
norms[i][j][kc] += 1.0
elif t == "sq_face_cap_trig_pris" and not flag_yaxis:
if thetak < self._params[i]['TA3']:
if thetam < self._params[i]['TA3']:
tmp = cos(self._params[i]['fac_AA1'] * \
phi2) ** self._params[i]['exp_cos_AA1']
tmp2 = self._params[i]['IGW_TA1'] * (
thetam * ipi - self._params[i]['TA1'])
else:
tmp = cos(self._params[i]['fac_AA2'] * \
(phi2 + self._params[i]['shift_AA2'])) ** \
self._params[i]['exp_cos_AA2']
tmp2 = self._params[i]['IGW_TA2'] * (
thetam * ipi - self._params[i]['TA2'])
#print("phi2 {} phi2+shift {} tmp {} ".format(phi2*180/pi, 180*(phi2 + self._params[i]['shift_AA2'])/pi, tmp))
qsptheta[i][j][kc] += tmp * exp(-0.5 * tmp2 * tmp2)
norms[i][j][kc] += 1

kc += 1

# Normalize Peters-style OPs.
for i, t in enumerate(self._types):
#if t == "pent_plan":
# ops[i] = ops[i] / sum(norms[i]) \
# if sum(norms[i]) > 1.0e-12 else None
if t in ["tri_plan", "tet", "bent", "sq_plan",
"oct", "oct_legacy", "cuboct", "pent_plan"]:
ops[i] = tmp_norm = 0.0
Expand All @@ -2604,7 +2639,8 @@ def get_order_parameters(self, structure, n, indices_neighs=None, \
elif t in ["T", "tri_pyr", "see_saw_rect", "sq_pyr", "tri_bipyr",
"sq_bipyr", "pent_pyr", "hex_pyr", "pent_bipyr",
"hex_bipyr", "oct_max", "tri_plan_max", "tet_max",
"sq_plan_max", "pent_plan_max", "cuboct_max", "hex_plan_max"]:
"sq_plan_max", "pent_plan_max", "cuboct_max", "hex_plan_max",
"sq_face_cap_trig_pris"]:
ops[i] = None
if nneigh > 1:
for j in range(nneigh):
Expand All @@ -2613,7 +2649,6 @@ def get_order_parameters(self, structure, n, indices_neighs=None, \
if norms[i][j][k] > 1.0e-12 else 0.0
ops[i] = max(qsptheta[i][j]) if j == 0 \
else max(ops[i], max(qsptheta[i][j]))
#ops[i] = max(qsptheta[i]) if len(qsptheta[i]) > 0 else None
elif t == "bcc":
ops[i] = 0.0
for j in range(nneigh):
Expand Down
21 changes: 21 additions & 0 deletions pymatgen/analysis/op_params.yaml
View file
Expand Up @@ -287,3 +287,24 @@ hex_plan_max:
IGW_TA:
25
TA: 0.16666667
sq_face_cap_trig_pris:
TA1:
0.37662414278557504
TA2:
0.7728144741714951
TA3:
1.8055339573923717
IGW_TA1:
15.5
IGW_TA2:
13
fac_AA1:
2
exp_cos_AA1:
2
fac_AA2:
1
shift_AA2:
-0.7853981633974483
exp_cos_AA2:
2
20 changes: 19 additions & 1 deletion pymatgen/analysis/tests/test_local_env.py
View file
Expand Up @@ -696,6 +696,18 @@ def setUp(self):
[0.0, 0.0, -1.0], [-1.0, 0.0, 0.0]],
validate_proximity=False, to_unit_cell=False,
coords_are_cartesian=True, site_properties=None)
self.sq_face_capped_trig_pris = Structure(
Lattice.from_lengths_and_angles(
[30, 30, 30], [90, 90, 90]),
["H", "H", "H", "H", "H", "H", "H", "H"],
[[0, 0, 0], [-0.6546536707079771, -0.37796447300922725, 0.6546536707079771],
[0.6546536707079771, -0.37796447300922725, 0.6546536707079771],
[0.0, 0.7559289460184545, 0.6546536707079771],
[-0.6546536707079771, -0.37796447300922725, -0.6546536707079771],
[0.6546536707079771, -0.37796447300922725, -0.6546536707079771],
[0.0, 0.7559289460184545, -0.6546536707079771], [0.0, -1.0, 0.0]],
validate_proximity=False, to_unit_cell=False,
coords_are_cartesian=True, site_properties=None)

def test_init(self):
self.assertIsNotNone(
Expand All @@ -714,7 +726,7 @@ def test_get_order_parameters(self):
"tri_plan", "sq_plan", "pent_plan", "sq_pyr", "tri_pyr", \
"pent_pyr", "hex_pyr", "pent_bipyr", "hex_bipyr", "T", "cuboct", \
"see_saw_rect", "hex_plan_max", "tet_max", "oct_max", "tri_plan_max", "sq_plan_max", \
"pent_plan_max", "cuboct_max", "tet_max"]
"pent_plan_max", "cuboct_max", "tet_max", "sq_face_cap_trig_pris"]
op_params = [None for i in range(len(op_types))]
op_params[1] = {'TA': 1, 'IGW_TA': 1./0.0667}
op_params[2] = {'TA': 45./180, 'IGW_TA': 1./0.0667}
Expand Down Expand Up @@ -893,6 +905,12 @@ def test_get_order_parameters(self):
self.hexagonal_planar, 0, indices_neighs=[1,2,3,4,5,6])
self.assertAlmostEqual(int(op_vals[26] * 1000 + 0.5), 1000)

# Square face capped trigonal prism.
op_vals = ops_101.get_order_parameters(
self.sq_face_capped_trig_pris, 0,
indices_neighs=[i for i in range(1, 8)])
self.assertAlmostEqual(int(op_vals[34] * 1000 + 0.5), 1000)

# Test providing explicit neighbor lists.
op_vals = ops_101.get_order_parameters(self.bcc, 0, indices_neighs=[1])
self.assertIsNotNone(op_vals[0])
Expand Down

0 comments on commit 74392fa

Please sign in to comment.