1#!/usr/bin/env python 2# Copyright 2014-2019 The PySCF Developers. All Rights Reserved. 3# 4# Licensed under the Apache License, Version 2.0 (the "License"); 5# you may not use this file except in compliance with the License. 6# You may obtain a copy of the License at 7# 8# http://www.apache.org/licenses/LICENSE-2.0 9# 10# Unless required by applicable law or agreed to in writing, software 11# distributed under the License is distributed on an "AS IS" BASIS, 12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13# See the License for the specific language governing permissions and 14# limitations under the License. 15# 16# Author: Qiming Sun <osirpt.sun@gmail.com> 17# 18# 19# The symmetry detection method implemented here is not strictly follow the 20# point group detection flowchart. The detection is based on the degeneracy 21# of cartesian basis of multipole momentum, eg 22# http://symmetry.jacobs-university.de/cgi-bin/group.cgi?group=604&option=4 23# see the column of "linear functions, quadratic functions and cubic functions". 24# 25# Different point groups have different combinations of degeneracy for each 26# type of cartesian functions. Based on the degeneracy of cartesian function 27# basis, one can quickly filter out a few candidates of point groups for the 28# given molecule. Regular operations (rotation, mirror etc) can be applied 29# then to identify the symmetry. Current implementation only checks the 30# rotation functions and it's roughly enough for D2h and subgroups. 31# 32# There are special cases this detection method may break down, eg two H8 cube 33# molecules sitting on the same center but with random orientation. The 34# system is in C1 while this detection method gives O group because the 35# 3 rotation bases are degenerated. In this case, the code use the regular 36# method (point group detection flowchart) to detect the point group. 37# 38 39import sys 40import re 41import numpy 42import scipy.linalg 43from pyscf.gto import mole 44from pyscf.lib import norm 45from pyscf.lib import logger 46from pyscf.lib.exceptions import PointGroupSymmetryError 47from pyscf.symm.param import OPERATOR_TABLE 48from pyscf import __config__ 49 50TOLERANCE = getattr(__config__, 'symm_geom_tol', 1e-5) 51 52# For code compatiblity in python-2 and python-3 53if sys.version_info >= (3,): 54 unicode = str 55 56def parallel_vectors(v1, v2, tol=TOLERANCE): 57 if numpy.allclose(v1, 0, atol=tol) or numpy.allclose(v2, 0, atol=tol): 58 return True 59 else: 60 cos = numpy.dot(_normalize(v1), _normalize(v2)) 61 return (abs(cos-1) < TOLERANCE) | (abs(cos+1) < TOLERANCE) 62 63def argsort_coords(coords, decimals=None): 64 if decimals is None: 65 decimals = int(-numpy.log10(TOLERANCE)) - 1 66 coords = numpy.around(coords, decimals=decimals) 67 idx = numpy.lexsort((coords[:,2], coords[:,1], coords[:,0])) 68 return idx 69 70def sort_coords(coords, decimals=None): 71 if decimals is None: 72 decimals = int(-numpy.log10(TOLERANCE)) - 1 73 coords = numpy.asarray(coords) 74 idx = argsort_coords(coords, decimals=decimals) 75 return coords[idx] 76 77# ref. http://en.wikipedia.org/wiki/Rotation_matrix 78def rotation_mat(vec, theta): 79 '''rotate angle theta along vec 80 new(x,y,z) = R * old(x,y,z)''' 81 vec = _normalize(vec) 82 uu = vec.reshape(-1,1) * vec.reshape(1,-1) 83 ux = numpy.array(( 84 ( 0 ,-vec[2], vec[1]), 85 ( vec[2], 0 ,-vec[0]), 86 (-vec[1], vec[0], 0 ))) 87 c = numpy.cos(theta) 88 s = numpy.sin(theta) 89 r = c * numpy.eye(3) + s * ux + (1-c) * uu 90 return r 91 92# reflection operation with householder 93def householder(vec): 94 vec = _normalize(vec) 95 return numpy.eye(3) - vec[:,None]*vec*2 96 97def closest_axes(axes, ref): 98 xcomp, ycomp, zcomp = numpy.einsum('ix,jx->ji', axes, ref) 99 zmax = numpy.amax(abs(zcomp)) 100 zmax_idx = numpy.where(abs(abs(zcomp)-zmax)<TOLERANCE)[0] 101 z_id = numpy.amax(zmax_idx) 102 #z_id = numpy.argmax(abs(zcomp)) 103 xcomp[z_id] = ycomp[z_id] = 0 # remove z 104 xmax = numpy.amax(abs(xcomp)) 105 xmax_idx = numpy.where(abs(abs(xcomp)-xmax)<TOLERANCE)[0] 106 x_id = numpy.amax(xmax_idx) 107 #x_id = numpy.argmax(abs(xcomp)) 108 ycomp[x_id] = 0 # remove x 109 y_id = numpy.argmax(abs(ycomp)) 110 return x_id, y_id, z_id 111 112def alias_axes(axes, ref): 113 '''Rename axes, make it as close as possible to the ref axes 114 ''' 115 x_id, y_id, z_id = closest_axes(axes, ref) 116 new_axes = axes[[x_id,y_id,z_id]] 117 if numpy.linalg.det(new_axes) < 0: 118 new_axes = axes[[y_id,x_id,z_id]] 119 return new_axes 120 121 122def detect_symm(atoms, basis=None, verbose=logger.WARN): 123 '''Detect the point group symmetry for given molecule. 124 125 Return group name, charge center, and nex_axis (three rows for x,y,z) 126 ''' 127 if isinstance(verbose, logger.Logger): 128 log = verbose 129 else: 130 log = logger.Logger(sys.stdout, verbose) 131 132 tol = TOLERANCE / numpy.sqrt(1+len(atoms)) 133 decimals = int(-numpy.log10(tol)) 134 log.debug('geometry tol = %g', tol) 135 136 rawsys = SymmSys(atoms, basis) 137 w1, u1 = rawsys.cartesian_tensor(1) 138 axes = u1.T 139 log.debug('principal inertia moments %s', w1) 140 charge_center = rawsys.charge_center 141 142 if numpy.allclose(w1, 0, atol=tol): 143 gpname = 'SO3' 144 return gpname, charge_center, numpy.eye(3) 145 146 elif numpy.allclose(w1[:2], 0, atol=tol): # linear molecule 147 if rawsys.has_icenter(): 148 gpname = 'Dooh' 149 else: 150 gpname = 'Coov' 151 return gpname, charge_center, axes 152 153 else: 154 w1_degeneracy, w1_degen_values = _degeneracy(w1, decimals) 155 w2, u2 = rawsys.cartesian_tensor(2) 156 w2_degeneracy, w2_degen_values = _degeneracy(w2, decimals) 157 log.debug('2d tensor %s', w2) 158 159 n = None 160 c2x = None 161 mirrorx = None 162 if 3 in w1_degeneracy: # T, O, I 163 # Because rotation vectors Rx Ry Rz are 3-degenerated representation 164 # See http://www.webqc.org/symmetrypointgroup-td.html 165 166 w3, u3 = rawsys.cartesian_tensor(3) 167 w3_degeneracy, w3_degen_values = _degeneracy(w3, decimals) 168 log.debug('3d tensor %s', w3) 169 if (5 in w2_degeneracy and 170 4 in w3_degeneracy and len(w3_degeneracy) == 3): # I group 171 gpname, new_axes = _search_i_group(rawsys) 172 if gpname is not None: 173 return gpname, charge_center, _refine(new_axes) 174 175 elif 3 in w2_degeneracy and len(w2_degeneracy) <= 3: # T/O group 176 gpname, new_axes = _search_ot_group(rawsys) 177 if gpname is not None: 178 return gpname, charge_center, _refine(new_axes) 179 180 elif (2 in w1_degeneracy and 181 numpy.any(w2_degeneracy[w2_degen_values>0] >= 2)): 182 if numpy.allclose(w1[1], w1[2], atol=tol): 183 axes = axes[[1,2,0]] 184 n = rawsys.search_c_highest(axes[2])[1] 185 if n == 1: 186 n = None 187 else: 188 c2x = rawsys.search_c2x(axes[2], n) 189 mirrorx = rawsys.search_mirrorx(axes[2], n) 190 191 else: 192 n = -1 # tag as D2h and subgroup 193 194# They must not be I/O/T group, at most one C3 or higher rotation axis 195 if n is None: 196 zaxis, n = rawsys.search_c_highest() 197 if n > 1: 198 c2x = rawsys.search_c2x(zaxis, n) 199 mirrorx = rawsys.search_mirrorx(zaxis, n) 200 if c2x is not None: 201 axes = _make_axes(zaxis, c2x) 202 elif mirrorx is not None: 203 axes = _make_axes(zaxis, mirrorx) 204 else: 205 for axis in numpy.eye(3): 206 if not parallel_vectors(axis, zaxis): 207 axes = _make_axes(zaxis, axis) 208 break 209 else: # Ci or Cs or C1 with degenerated w1 210 mirror = rawsys.search_mirrorx(None, 1) 211 if mirror is not None: 212 xaxis = numpy.array((1.,0.,0.)) 213 axes = _make_axes(mirror, xaxis) 214 else: 215 axes = numpy.eye(3) 216 217 log.debug('Highest C_n = C%d', n) 218 if n >= 2: 219 if c2x is not None: 220 if rawsys.has_mirror(axes[2]): 221 gpname = 'D%dh' % n 222 elif rawsys.has_improper_rotation(axes[2], n): 223 gpname = 'D%dd' % n 224 else: 225 gpname = 'D%d' % n 226 # yaxis = numpy.cross(axes[2], c2x) 227 axes = _make_axes(axes[2], c2x) 228 elif mirrorx is not None: 229 gpname = 'C%dv' % n 230 axes = _make_axes(axes[2], mirrorx) 231 elif rawsys.has_mirror(axes[2]): 232 gpname = 'C%dh' % n 233 elif rawsys.has_improper_rotation(axes[2], n): 234 gpname = 'S%d' % (n*2) 235 else: 236 gpname = 'C%d' % n 237 return gpname, charge_center, _refine(axes) 238 239 else: 240 is_c2x = rawsys.has_rotation(axes[0], 2) 241 is_c2y = rawsys.has_rotation(axes[1], 2) 242 is_c2z = rawsys.has_rotation(axes[2], 2) 243# rotate to old axes, as close as possible? 244 if is_c2z and is_c2x and is_c2y: 245 if rawsys.has_icenter(): 246 gpname = 'D2h' 247 else: 248 gpname = 'D2' 249 axes = alias_axes(axes, numpy.eye(3)) 250 elif is_c2z or is_c2x or is_c2y: 251 if is_c2x: 252 axes = axes[[1,2,0]] 253 if is_c2y: 254 axes = axes[[2,0,1]] 255 if rawsys.has_mirror(axes[2]): 256 gpname = 'C2h' 257 elif rawsys.has_mirror(axes[0]): 258 gpname = 'C2v' 259 else: 260 gpname = 'C2' 261 else: 262 if rawsys.has_icenter(): 263 gpname = 'Ci' 264 elif rawsys.has_mirror(axes[0]): 265 gpname = 'Cs' 266 axes = axes[[1,2,0]] 267 elif rawsys.has_mirror(axes[1]): 268 gpname = 'Cs' 269 axes = axes[[2,0,1]] 270 elif rawsys.has_mirror(axes[2]): 271 gpname = 'Cs' 272 else: 273 gpname = 'C1' 274 axes = numpy.eye(3) 275 charge_center = numpy.zeros(3) 276 return gpname, charge_center, axes 277 278 279# reduce to D2h and its subgroups 280# FIXME, CPL, 209, 506 281def get_subgroup(gpname, axes): 282 if gpname in ('D2h', 'D2' , 'C2h', 'C2v', 'C2' , 'Ci' , 'Cs' , 'C1'): 283 return gpname, axes 284 elif gpname in ('SO3',): 285 #return 'D2h', alias_axes(axes, numpy.eye(3)) 286 return 'SO3', axes 287 elif gpname in ('Dooh',): 288 #return 'D2h', alias_axes(axes, numpy.eye(3)) 289 return 'Dooh', axes 290 elif gpname in ('Coov',): 291 #return 'C2v', axes 292 return 'Coov', axes 293 elif gpname in ('Oh',): 294 return 'D2h', alias_axes(axes, numpy.eye(3)) 295 elif gpname in ('O',): 296 return 'D2', alias_axes(axes, numpy.eye(3)) 297 elif gpname in ('Ih',): 298 return 'Ci', alias_axes(axes, numpy.eye(3)) 299 elif gpname in ('I',): 300 return 'C1', axes 301 elif gpname in ('Td', 'T', 'Th'): 302 return 'D2', alias_axes(axes, numpy.eye(3)) 303 elif re.search(r'S\d+', gpname): 304 n = int(re.search(r'\d+', gpname).group(0)) 305 if n % 2 == 0: 306 return 'C%d'%(n//2), axes 307 else: 308 return 'Ci', axes 309 else: 310 n = int(re.search(r'\d+', gpname).group(0)) 311 if n % 2 == 0: 312 if re.search(r'D\d+d', gpname): 313 subname = 'D2' 314 elif re.search(r'D\d+h', gpname): 315 subname = 'D2h' 316 elif re.search(r'D\d+', gpname): 317 subname = 'D2' 318 elif re.search(r'C\d+h', gpname): 319 subname = 'C2h' 320 elif re.search(r'C\d+v', gpname): 321 subname = 'C2v' 322 else: 323 subname = 'C2' 324 else: 325 # rotate axes and 326 # Dnh -> C2v 327 # Dn -> C2 328 # Dnd -> Ci 329 # Cnh -> Cs 330 # Cnv -> Cs 331 if re.search(r'D\d+h', gpname): 332 subname = 'C2v' 333 axes = axes[[1,2,0]] 334 elif re.search(r'D\d+d', gpname): 335 subname = 'C2h' 336 axes = axes[[1,2,0]] 337 elif re.search(r'D\d+', gpname): 338 subname = 'C2' 339 axes = axes[[1,2,0]] 340 elif re.search(r'C\d+h', gpname): 341 subname = 'Cs' 342 elif re.search(r'C\d+v', gpname): 343 subname = 'Cs' 344 axes = axes[[1,2,0]] 345 else: 346 subname = 'C1' 347 return subname, axes 348subgroup = get_subgroup 349 350def as_subgroup(topgroup, axes, subgroup=None): 351 from pyscf.symm import std_symb 352 from pyscf.symm.param import SUBGROUP 353 354 groupname, axes = get_subgroup(topgroup, axes) 355 356 if isinstance(subgroup, (str, unicode)): 357 subgroup = std_symb(subgroup) 358 if groupname == 'C2v' and subgroup == 'Cs': 359 axes = numpy.einsum('ij,kj->ki', rotation_mat(axes[1], numpy.pi/2), axes) 360 361 elif (groupname == 'D2' and re.search(r'D\d+d', topgroup) and 362 subgroup in ('C2v', 'Cs')): 363 # Special treatment for D2d, D4d, .... get_subgroup gives D2 by 364 # default while C2v is also D2d's subgroup. 365 groupname = 'C2v' 366 axes = numpy.einsum('ij,kj->ki', rotation_mat(axes[2], numpy.pi/4), axes) 367 368 elif topgroup in ('Td', 'T', 'Th') and subgroup == 'C2v': 369 x, y, z = axes 370 x = _normalize(x+y) 371 y = numpy.cross(z, x) 372 axes = numpy.array((x,y,z)) 373 374 elif subgroup not in SUBGROUP[groupname]: 375 raise PointGroupSymmetryError('%s not in Ablien subgroup of %s' % 376 (subgroup, topgroup)) 377 378 groupname = subgroup 379 return groupname, axes 380 381def symm_ops(gpname, axes=None): 382 if axes is not None: 383 raise PointGroupSymmetryError('TODO: non-standard orientation') 384 op1 = numpy.eye(3) 385 opi = -1 386 387 opc2z = -numpy.eye(3) 388 opc2z[2,2] = 1 389 opc2x = -numpy.eye(3) 390 opc2x[0,0] = 1 391 opc2y = -numpy.eye(3) 392 opc2y[1,1] = 1 393 394 opcsz = numpy.dot(opc2z, opi) 395 opcsx = numpy.dot(opc2x, opi) 396 opcsy = numpy.dot(opc2y, opi) 397 opdic = {'E' : op1, 398 'C2z': opc2z, 399 'C2x': opc2x, 400 'C2y': opc2y, 401 'i' : opi, 402 'sz' : opcsz, 403 'sx' : opcsx, 404 'sy' : opcsy,} 405 return opdic 406 407def symm_identical_atoms(gpname, atoms): 408 ''' Requires ''' 409 from pyscf import gto 410 # Dooh Coov for linear molecule 411 if gpname == 'Dooh': 412 coords = numpy.array([a[1] for a in atoms], dtype=float) 413 idx0 = argsort_coords(coords) 414 coords0 = coords[idx0] 415 opdic = symm_ops(gpname) 416 newc = numpy.dot(coords, opdic['sz']) 417 idx1 = argsort_coords(newc) 418 dup_atom_ids = numpy.sort((idx0,idx1), axis=0).T 419 uniq_idx = numpy.unique(dup_atom_ids[:,0], return_index=True)[1] 420 eql_atom_ids = dup_atom_ids[uniq_idx] 421 eql_atom_ids = [list(sorted(set(i))) for i in eql_atom_ids] 422 return eql_atom_ids 423 elif gpname == 'Coov': 424 eql_atom_ids = [[i] for i,a in enumerate(atoms)] 425 return eql_atom_ids 426 427 coords = numpy.array([a[1] for a in atoms]) 428 429# charges = numpy.array([gto.charge(a[0]) for a in atoms]) 430# center = numpy.einsum('z,zr->r', charges, coords)/charges.sum() 431# if not numpy.allclose(center, 0, atol=TOLERANCE): 432# sys.stderr.write('WARN: Molecular charge center %s is not on (0,0,0)\n' 433# % center) 434 opdic = symm_ops(gpname) 435 ops = [opdic[op] for op in OPERATOR_TABLE[gpname]] 436 idx = argsort_coords(coords) 437 coords0 = coords[idx] 438 439 dup_atom_ids = [] 440 for op in ops: 441 newc = numpy.dot(coords, op) 442 idx = argsort_coords(newc) 443 if not numpy.allclose(coords0, newc[idx], atol=TOLERANCE): 444 raise PointGroupSymmetryError('Symmetry identical atoms not found') 445 dup_atom_ids.append(idx) 446 447 dup_atom_ids = numpy.sort(dup_atom_ids, axis=0).T 448 uniq_idx = numpy.unique(dup_atom_ids[:,0], return_index=True)[1] 449 eql_atom_ids = dup_atom_ids[uniq_idx] 450 eql_atom_ids = [list(sorted(set(i))) for i in eql_atom_ids] 451 return eql_atom_ids 452 453def check_given_symm(gpname, atoms, basis=None): 454 # more strict than symm_identical_atoms, we required not only the coordinates 455 # match, but also the symbols and basis functions 456 457 #FIXME: compare the basis set when basis is given 458 if gpname == 'Dooh': 459 coords = numpy.array([a[1] for a in atoms], dtype=float) 460 if numpy.allclose(coords[:,:2], 0, atol=TOLERANCE): 461 opdic = symm_ops(gpname) 462 rawsys = SymmSys(atoms, basis) 463 return rawsys.has_icenter() and numpy.allclose(rawsys.charge_center, 0) 464 else: 465 return False 466 elif gpname == 'Coov': 467 coords = numpy.array([a[1] for a in atoms], dtype=float) 468 return numpy.allclose(coords[:,:2], 0, atol=TOLERANCE) 469 470 opdic = symm_ops(gpname) 471 ops = [opdic[op] for op in OPERATOR_TABLE[gpname]] 472 rawsys = SymmSys(atoms, basis) 473 for lst in rawsys.atomtypes.values(): 474 coords = rawsys.atoms[lst,1:] 475 idx = argsort_coords(coords) 476 coords0 = coords[idx] 477 478 for op in ops: 479 newc = numpy.dot(coords, op) 480 idx = argsort_coords(newc) 481 if not numpy.allclose(coords0, newc[idx], atol=TOLERANCE): 482 return False 483 return True 484 485def shift_atom(atoms, orig, axis): 486 c = numpy.array([a[1] for a in atoms]) 487 c = numpy.dot(c - orig, numpy.array(axis).T) 488 return [[atoms[i][0], c[i]] for i in range(len(atoms))] 489 490class RotationAxisNotFound(PointGroupSymmetryError): 491 pass 492 493class SymmSys(object): 494 def __init__(self, atoms, basis=None): 495 self.atomtypes = mole.atom_types(atoms, basis) 496 # fake systems, which treates the atoms of different basis as different atoms. 497 # the fake systems do not have the same symmetry as the potential 498 # it's only used to determine the main (Z-)axis 499 chg1 = numpy.pi - 2 500 coords = [] 501 fake_chgs = [] 502 idx = [] 503 for k, lst in self.atomtypes.items(): 504 idx.append(lst) 505 coords.append([atoms[i][1] for i in lst]) 506 ksymb = mole._rm_digit(k) 507 if ksymb != k: 508 # Put random charges on the decorated atoms 509 fake_chgs.append([chg1] * len(lst)) 510 chg1 *= numpy.pi-2 511 elif mole.is_ghost_atom(k): 512 if ksymb == 'X' or ksymb.upper() == 'GHOST': 513 fake_chgs.append([.3] * len(lst)) 514 elif ksymb[0] == 'X': 515 fake_chgs.append([mole.charge(ksymb[1:])+.3] * len(lst)) 516 elif ksymb[:5] == 'GHOST': 517 fake_chgs.append([mole.charge(ksymb[5:])+.3] * len(lst)) 518 else: 519 fake_chgs.append([mole.charge(ksymb)] * len(lst)) 520 coords = numpy.array(numpy.vstack(coords), dtype=float) 521 fake_chgs = numpy.hstack(fake_chgs) 522 self.charge_center = numpy.einsum('i,ij->j', fake_chgs, coords)/fake_chgs.sum() 523 coords = coords - self.charge_center 524 525 idx = numpy.argsort(numpy.hstack(idx)) 526 self.atoms = numpy.hstack((fake_chgs.reshape(-1,1), coords))[idx] 527 528 self.group_atoms_by_distance = [] 529 decimals = int(-numpy.log10(TOLERANCE)) - 1 530 for index in self.atomtypes.values(): 531 index = numpy.asarray(index) 532 c = self.atoms[index,1:] 533 dists = numpy.around(norm(c, axis=1), decimals) 534 u, idx = numpy.unique(dists, return_inverse=True) 535 for i, s in enumerate(u): 536 self.group_atoms_by_distance.append(index[idx == i]) 537 538 def cartesian_tensor(self, n): 539 z = self.atoms[:,0] 540 r = self.atoms[:,1:] 541 ncart = (n+1)*(n+2)//2 542 natm = len(z) 543 tensor = numpy.sqrt(numpy.copy(z).reshape(natm,-1) / z.sum()) 544 for i in range(n): 545 tensor = numpy.einsum('zi,zj->zij', tensor, r).reshape(natm,-1) 546 e, c = scipy.linalg.eigh(numpy.dot(tensor.T,tensor)) 547 return e[-ncart:], c[:,-ncart:] 548 549 def symmetric_for(self, op): 550 for lst in self.group_atoms_by_distance: 551 r0 = self.atoms[lst,1:] 552 r1 = numpy.dot(r0, op) 553# FIXME: compare whehter two sets of coordinates are identical 554 yield all((_vec_in_vecs(x, r0) for x in r1)) 555 556 def has_icenter(self): 557 return all(self.symmetric_for(-1)) 558 559 def has_rotation(self, axis, n): 560 op = rotation_mat(axis, numpy.pi*2/n).T 561 return all(self.symmetric_for(op)) 562 563 def has_mirror(self, perp_vec): 564 return all(self.symmetric_for(householder(perp_vec).T)) 565 566 def has_improper_rotation(self, axis, n): 567 s_op = numpy.dot(householder(axis), rotation_mat(axis, numpy.pi/n)).T 568 return all(self.symmetric_for(s_op)) 569 570 def search_possible_rotations(self, zaxis=None): 571 '''If zaxis is given, the rotation axis is parallel to zaxis''' 572 maybe_cn = [] 573 for lst in self.group_atoms_by_distance: 574 natm = len(lst) 575 if natm > 1: 576 coords = self.atoms[lst,1:] 577# possible C2 axis 578 for i in range(1, natm): 579 if abs(coords[0]+coords[i]).sum() > TOLERANCE: 580 maybe_cn.append((coords[0]+coords[i], 2)) 581 else: # abs(coords[0]-coords[i]).sum() > TOLERANCE: 582 maybe_cn.append((coords[0]-coords[i], 2)) 583 584# atoms of equal distances may be associated with rotation axis > C2. 585 r0 = coords - coords[0] 586 distance = norm(r0, axis=1) 587 eq_distance = abs(distance[:,None] - distance) < TOLERANCE 588 for i in range(2, natm): 589 for j in numpy.where(eq_distance[i,:i])[0]: 590 cos = numpy.dot(r0[i],r0[j]) / (distance[i]*distance[j]) 591 ang = numpy.arccos(cos) 592 nfrac = numpy.pi*2 / (numpy.pi-ang) 593 n = int(numpy.around(nfrac)) 594 if abs(nfrac-n) < TOLERANCE: 595 maybe_cn.append((numpy.cross(r0[i],r0[j]),n)) 596 597 # remove zero-vectors and duplicated vectors 598 vecs = numpy.vstack([x[0] for x in maybe_cn]) 599 idx = norm(vecs, axis=1) > TOLERANCE 600 ns = numpy.hstack([x[1] for x in maybe_cn]) 601 vecs = _normalize(vecs[idx]) 602 ns = ns[idx] 603 604 if zaxis is not None: # Keep parallel rotation axes 605 cos = numpy.dot(vecs, _normalize(zaxis)) 606 vecs = vecs[(abs(cos-1) < TOLERANCE) | (abs(cos+1) < TOLERANCE)] 607 ns = ns[(abs(cos-1) < TOLERANCE) | (abs(cos+1) < TOLERANCE)] 608 609 possible_cn = [] 610 seen = numpy.zeros(len(vecs), dtype=bool) 611 for k, v in enumerate(vecs): 612 if not seen[k]: 613 where1 = numpy.einsum('ix->i', abs(vecs[k:] - v)) < TOLERANCE 614 where1 = numpy.where(where1)[0] + k 615 where2 = numpy.einsum('ix->i', abs(vecs[k:] + v)) < TOLERANCE 616 where2 = numpy.where(where2)[0] + k 617 seen[where1] = True 618 seen[where2] = True 619 620 vk = _normalize((numpy.einsum('ix->x', vecs[where1]) - 621 numpy.einsum('ix->x', vecs[where2]))) 622 for n in (set(ns[where1]) | set(ns[where2])): 623 possible_cn.append((vk,n)) 624 return possible_cn 625 626 def search_c2x(self, zaxis, n): 627 '''C2 axis which is perpendicular to z-axis''' 628 decimals = int(-numpy.log10(TOLERANCE)) - 1 629 for lst in self.group_atoms_by_distance: 630 if len(lst) > 1: 631 r0 = self.atoms[lst,1:] 632 zcos = numpy.around(numpy.einsum('ij,j->i', r0, zaxis), 633 decimals=decimals) 634 uniq_zcos = numpy.unique(zcos) 635 maybe_c2x = [] 636 for d in uniq_zcos: 637 if d > TOLERANCE: 638 mirrord = abs(zcos+d)<TOLERANCE 639 if mirrord.sum() == (zcos==d).sum(): 640 above = r0[zcos==d] 641 below = r0[mirrord] 642 nelem = len(below) 643 maybe_c2x.extend([above[0] + below[i] 644 for i in range(nelem)]) 645 elif abs(d) < TOLERANCE: # plane which crosses the orig 646 r1 = r0[zcos==d][0] 647 maybe_c2x.append(r1) 648 r2 = numpy.dot(rotation_mat(zaxis, numpy.pi*2/n), r1) 649 if abs(r1+r2).sum() > TOLERANCE: 650 maybe_c2x.append(r1+r2) 651 else: 652 maybe_c2x.append(r2-r1) 653 654 if len(maybe_c2x) > 0: 655 idx = norm(maybe_c2x, axis=1) > TOLERANCE 656 maybe_c2x = _normalize(maybe_c2x)[idx] 657 maybe_c2x = _remove_dupvec(maybe_c2x) 658 for c2x in maybe_c2x: 659 if (not parallel_vectors(c2x, zaxis) and 660 self.has_rotation(c2x, 2)): 661 return c2x 662 663 def search_mirrorx(self, zaxis, n): 664 '''mirror which is parallel to z-axis''' 665 if n > 1: 666 for lst in self.group_atoms_by_distance: 667 natm = len(lst) 668 r0 = self.atoms[lst[0],1:] 669 if natm > 1 and not parallel_vectors(r0, zaxis): 670 r1 = numpy.dot(rotation_mat(zaxis, numpy.pi*2/n), r0) 671 mirrorx = _normalize(r1-r0) 672 if self.has_mirror(mirrorx): 673 return mirrorx 674 else: 675 for lst in self.group_atoms_by_distance: 676 natm = len(lst) 677 r0 = self.atoms[lst,1:] 678 if natm > 1: 679 maybe_mirror = [r0[i]-r0[0] for i in range(1, natm)] 680 for mirror in _normalize(maybe_mirror): 681 if self.has_mirror(mirror): 682 return mirror 683 684 def search_c_highest(self, zaxis=None): 685 possible_cn = self.search_possible_rotations(zaxis) 686 nmax = 1 687 cmax = numpy.array([0.,0.,1.]) 688 for cn, n in possible_cn: 689 if n > nmax and self.has_rotation(cn, n): 690 nmax = n 691 cmax = cn 692 return cmax, nmax 693 694 695def _normalize(vecs): 696 vecs = numpy.asarray(vecs) 697 if vecs.ndim == 1: 698 return vecs / (numpy.linalg.norm(vecs) + 1e-200) 699 else: 700 return vecs / (norm(vecs, axis=1).reshape(-1,1) + 1e-200) 701 702def _vec_in_vecs(vec, vecs): 703 norm = numpy.sqrt(len(vecs)) 704 return min(numpy.einsum('ix->i', abs(vecs-vec))/norm) < TOLERANCE 705 706def _search_i_group(rawsys): 707 possible_cn = rawsys.search_possible_rotations() 708 c5_axes = [c5 for c5, n in possible_cn 709 if n == 5 and rawsys.has_rotation(c5, 5)] 710 if len(c5_axes) <= 1: 711 return None,None 712 713 zaxis = c5_axes[0] 714 cos = numpy.dot(c5_axes, zaxis) 715 assert(numpy.all((abs(cos[1:]+1/numpy.sqrt(5)) < TOLERANCE) | 716 (abs(cos[1:]-1/numpy.sqrt(5)) < TOLERANCE))) 717 718 if rawsys.has_icenter(): 719 gpname = 'Ih' 720 else: 721 gpname = 'I' 722 723 c5 = c5_axes[1] 724 if numpy.dot(c5, zaxis) < 0: 725 c5 = -c5 726 c5a = numpy.dot(rotation_mat(zaxis, numpy.pi*6/5), c5) 727 xaxis = c5a + c5 728 return gpname, _make_axes(zaxis, xaxis) 729 730def _search_ot_group(rawsys): 731 possible_cn = rawsys.search_possible_rotations() 732 c4_axes = [c4 for c4, n in possible_cn 733 if n == 4 and rawsys.has_rotation(c4, 4)] 734 735 if len(c4_axes) > 0: # T group 736 assert(len(c4_axes) > 1) 737 if rawsys.has_icenter(): 738 gpname = 'Oh' 739 else: 740 gpname = 'O' 741 return gpname, _make_axes(c4_axes[0], c4_axes[1]) 742 743 else: # T group 744 c3_axes = [c3 for c3, n in possible_cn 745 if n == 3 and rawsys.has_rotation(c3, 3)] 746 if len(c3_axes) <= 1: 747 return None, None 748 749 cos = numpy.dot(c3_axes, c3_axes[0]) 750 assert(numpy.all((abs(cos[1:]+1./3) < TOLERANCE) | 751 (abs(cos[1:]-1./3) < TOLERANCE))) 752 753 if rawsys.has_icenter(): 754 gpname = 'Th' 755# Because C3 axes are on the mirror of Td, two C3 can determine a mirror. 756 elif rawsys.has_mirror(numpy.cross(c3_axes[0], c3_axes[1])): 757 gpname = 'Td' 758 else: 759 gpname = 'T' 760 761 c3a = c3_axes[0] 762 if numpy.dot(c3a, c3_axes[1]) > 0: 763 c3a = -c3a 764 c3b = numpy.dot(rotation_mat(c3a,-numpy.pi*2/3), c3_axes[1]) 765 c3c = numpy.dot(rotation_mat(c3a, numpy.pi*2/3), c3_axes[1]) 766 zaxis, xaxis = c3a+c3b, c3a+c3c 767 return gpname, _make_axes(zaxis, xaxis) 768 769def _degeneracy(e, decimals): 770 e = numpy.around(e, decimals) 771 u, idx = numpy.unique(e, return_inverse=True) 772 degen = numpy.array([numpy.count_nonzero(idx==i) for i in range(len(u))]) 773 return degen, u 774 775def _pseudo_vectors(vs): 776 idy0 = abs(vs[:,1])<TOLERANCE 777 idz0 = abs(vs[:,2])<TOLERANCE 778 vs = vs.copy() 779 # ensure z component > 0 780 vs[vs[:,2]<0] *= -1 781 # if z component == 0, ensure y component > 0 782 vs[(vs[:,1]<0) & idz0] *= -1 783 # if y and z component == 0, ensure x component > 0 784 vs[(vs[:,0]<0) & idy0 & idz0] *= -1 785 return vs 786 787def _remove_dupvec(vs): 788 def rm_iter(vs): 789 if len(vs) <= 1: 790 return vs 791 else: 792 x = numpy.sum(abs(vs[1:]-vs[0]), axis=1) 793 rest = rm_iter(vs[1:][x>TOLERANCE]) 794 return numpy.vstack((vs[0], rest)) 795 return rm_iter(_pseudo_vectors(vs)) 796 797def _make_axes(z, x): 798 y = numpy.cross(z, x) 799 x = numpy.cross(y, z) # because x might not perp to z 800 return _normalize(numpy.array((x,y,z))) 801 802def _refine(axes): 803 # Make sure the axes can be rotated from continuous unitary transformation 804 if axes[2,2] < 0: 805 axes[2] *= -1 806 if abs(axes[0,0]) > abs(axes[1,0]): 807 x_id, y_id = 0, 1 808 else: 809 x_id, y_id = 1, 0 810 if axes[x_id,0] < 0: 811 axes[x_id] *= -1 812 if numpy.linalg.det(axes) < 0: 813 axes[y_id] *= -1 814 return axes 815 816 817if __name__ == "__main__": 818 atom = [["O" , (1. , 0. , 0. ,)], 819 ['H' , (0. , -.757 , 0.587,)], 820 ['H' , (0. , 0.757 , 0.587,)] ] 821 gpname, orig, axes = detect_symm(atom) 822 atom = shift_atom(atom, orig, axes) 823 print(gpname, symm_identical_atoms(gpname, atom)) 824 825 atom = [['H', (0,0,0)], ['H', (0,0,-1)], ['H', (0,0,1)]] 826 gpname, orig, axes = detect_symm(atom) 827 print(gpname, orig, axes) 828 atom = shift_atom(atom, orig, axes) 829 print(gpname, symm_identical_atoms(gpname, atom)) 830 831 atom = [['H', (0., 0., 0.)], 832 ['H', (0., 0., 1.)], 833 ['H', (0., 1., 0.)], 834 ['H', (1., 0., 0.)], 835 ['H', (-1, 0., 0.)], 836 ['H', (0.,-1., 0.)], 837 ['H', (0., 0.,-1.)]] 838 gpname, orig, axes = detect_symm(atom) 839 print(gpname, orig, axes) 840 atom = shift_atom(atom, orig, axes) 841 print(gpname, symm_identical_atoms(subgroup(gpname, axes)[0], atom)) 842