""" This module is the EMBED module for ASE implemented by T. Kerber Torsten Kerber, ENS LYON: 2011, 07, 11 This work is supported by Award No. UK-C0017, made by King Abdullah University of Science and Technology (KAUST) """ import math from ase import Atom, Atoms from ase.data import covalent_radii, atomic_numbers import numpy as np class Embed(Atoms): #--- constructor of the Embed class --- def __init__(self, system=None, cluster=None, cell_cluster = "Auto"): super(Embed, self).__init__() # define the atom map self.atom_map_sys_cl = [] self.atom_map_cl_sys = [] self.linkatoms = [] # cluster dimensions self.xyz_cl_min = None self.xyz_cl_max = None # set the search radius for link atoms self.d = 10 # define the systems for calculations self.set_system(system) self.set_cluster(cluster) # set the cell of the system self.set_cell(system.get_cell()) self.cell_cluster = cell_cluster return #--- set the cluster --- def set_cluster(self, atoms): import copy # set the min/max cluster dimensions self.xyz_cl_min = atoms[0].get_position() self.xyz_cl_max = atoms[0].get_position() for atom in atoms: # assign the label "Cluster (10)" in atom.TAG atom.set_tag(10) xyz=atom.get_position() for i in xrange(3): # set the min/max cluster dimensions if xyz[i] < self.xyz_cl_min[i]: self.xyz_cl_min[i] = xyz[i] if xyz[i] > self.xyz_cl_max[i]: self.xyz_cl_max[i] = xyz[i] # add self.d around min/max cluster dimensions self.xyz_cl_min -= [self.d, self.d, self.d] self.xyz_cl_max += [self.d, self.d, self.d] # set the cluster for low and high level calculation self.atoms_cluster = atoms return #--- set the system --- def set_system(self, atoms): self.atoms_system = atoms # assign the label "Cluster (10)" in atom.TAG for atom in atoms: atom.set_tag(0) # update search radius for link atoms dx = 0 for atom in atoms: r = covalent_radii[atom.get_atomic_number()] if (r > dx): dx = r self.d = dx * 2.1 return #--- return cluster --- def get_cluster(self): return self.atoms_cluster def get_system(self): return self.atoms_system #--- Embedding --- def embed(self): # is the cluster and the host system definied ? if self.atoms_cluster is None or self.atoms_system is None: return self.find_cluster() self.set_linkatoms() print "link atoms found: ", len(self.linkatoms) if self.cell_cluster == "System": self.atoms_cluster.set_cell(self.atoms_system.get_cell()) elif self.cell_cluster == "Auto": positions = self.atoms_cluster.get_positions() #find the biggest dimensions of the cluster in x,y,z direction l = 0 for idir in xrange(3): l = max(l, positions[:, idir].max() - positions[:, idir].min()) # calculate the box parameters (cluster + min 5 Ang) l = (math.floor(l/2.5)+1)*2.5 + 5.0 # build cell cell = np.zeros((3, 3), float) # apply cell parameters for idir in xrange(3): cell[idir, idir] = l # set parameters to cluster self.atoms_cluster.set_cell(cell) # print information on the screen print "length of box surrounding the cluster: ", print l*10, print "pm" else: self.atoms_cluster.set_cell(self.cell_cluster) def find_cluster(self): # set tolerance d = 0.001 #atoms xyzs_cl=[] for atom_cl in self.atoms_cluster: xyzs_cl.append(atom_cl.get_position()) xyzs_sys=[] for atom_sys in self.atoms_system: xyzs_sys.append(atom_sys.get_position()) self.atom_map_sys_cl = np.zeros(len(self.atoms_system), int) self.atom_map_cl_sys = np.zeros(len(self.atoms_cluster), int) # loop over cluster atoms atom_sys for iat_sys in xrange(len(self.atoms_system)): # get the coordinates of the system atom atom_sys xyz_sys = xyzs_sys[iat_sys] # bSysOnly: no identical atom has been found bSysOnly = True # loop over system atoms atom_cl for iat_cl in xrange(len(self.atoms_cluster)): # difference vector between both atoms xyz_diff = np.abs(xyzs_sys[iat_sys]-xyzs_cl[iat_cl]) # identical atoms if xyz_diff[0] < d and xyz_diff[1] < d and xyz_diff[2] < d: # set tag (CLUSTER+HOST: 10) to atom_sys self.atoms_system[iat_sys].set_tag(10) # map the atom in the atom list self.atom_map_sys_cl[iat_sys] = iat_cl self.atom_map_cl_sys[iat_cl] = iat_sys # atom has been identified bSysOnly = False break if bSysOnly: self.atom_map_sys_cl[iat_sys] = -1 def set_linkatoms(self, tol=15., linkAtom=None, debug=False): # local copies of xyz coordinates to avoid massive copying of xyz objects xyzs_cl=[] for atom_cl in self.atoms_cluster: xyzs_cl.append(atom_cl.get_position()) xyzs_sys=[] for atom_sys in self.atoms_system: xyzs_sys.append(atom_sys.get_position()) # set the standard link atom if linkAtom is None: linkAtom ='H' # number of atoms in the cluster and the system nat_cl=len(self.atoms_cluster) nat_sys=len(self.atoms_system) # system has pbc? pbc = self.get_pbc() # set the bond table bonds = [] # set the 27 cell_vec, starting with the (0,0,0) vector for the unit cell cells_L = [(0.0, 0.0, 0.0)] # get the cell vectors of the host system and build up a 3 by 3 supercell to search for neighbors in the surrounding cell = self.atoms_system.get_cell() if self.atoms_system.get_pbc().any(): for ix in xrange(-1, 2): for iy in xrange(-1, 2): for iz in xrange(-1, 2): if ix == 0 and iy == 0 and iz == 0: continue cells_L.append(np.dot([ix, iy, iz], cell)) # save the radius of system atoms rs_sys = [] for atom in self.atoms_system: rs_sys.append(covalent_radii[atom.get_atomic_number()]) # sum over cluster atoms (iat_cl) for iat_cl in xrange(nat_cl): # get the cluster atom atom_cl=self.atoms_cluster[iat_cl] # ignore link atoms if atom_cl.get_tag() == 50: continue # xyz coordinates and covalent radius of the cluster atom iat_cl xyz_cl = xyzs_cl[iat_cl] r_cl = covalent_radii[atom_cl.get_atomic_number()] # sum over system atoms (iat_sys) for iat_sys in xrange(nat_sys): # avoid cluster atoms if self.atoms_system[iat_sys].get_tag()==10: continue # sum over all cell_L for cell_L in cells_L: # xyz coordinates and covalent radius of the system atom iat_sys xyz_sys = xyzs_sys[iat_sys]+cell_L # go only in distance self.d around the cluster lcont = True for i in xrange(3): if (xyz_sys[i] < self.xyz_cl_min[i] or xyz_sys[i] > self.xyz_cl_max[i]): lcont = False break if not lcont: continue # xyz coordinates and covalent radius of the system atom iat_sys r_sys = rs_sys[iat_sys] # diff vector xyz_diff = xyz_sys - xyz_cl # distance between the atoms r = np.sqrt(np.dot(xyz_diff, xyz_diff)) # ratio of the distance to the sum of covalent radius f = r / (r_cl + r_sys) if debug: print "Covalent radii = ",r_cl, r_sys print "Distance ", f print "tol = ",(1+tol/100.),(1-tol/100.),(1-2*tol/100.) if f <= (1+tol/100.) and f >= (1-2*tol/100.): s = cell_L, self.atom_map_cl_sys[iat_cl], iat_sys, r_cl bonds.append(s) break if f <= (1-2*tol/100.): raise RuntimeError("QMX: The cluster atom", iat_cl, " and the system atom", iat_sys, "came too close") r_h = covalent_radii[atomic_numbers[linkAtom]] for bond in bonds: cell_L, iat_cl_sys, iat_sys, r_cl = bond # assign the tags for the border atoms atom_sys.set_tag(1) atom_cl.set_tag(11) #difference vector for the link atom, scaling xyz_diff = xyzs_sys[iat_sys] + cell_L - xyzs_sys[iat_cl_sys] r = (r_cl + r_h) xyz_diff *= r / np.sqrt(np.dot(xyz_diff, xyz_diff)) # determine position of the link atom xyz_diff += xyzs_sys[iat_cl_sys] # create link atom atom = Atom(symbol=linkAtom, position=xyz_diff, tag=50) # add atom to cluster self.atoms_cluster.append(atom) # add atom to the linkatoms s = cell_L, iat_cl_sys, iat_sys, r, len(self.atoms_cluster)-1 self.linkatoms.append(s) return def set_positions(self, positions_new): # number of atoms nat_sys=len(self.atoms_system) # go over all pairs of atoms for iat_sys in xrange(nat_sys): xyz = positions_new[iat_sys] self.atoms_system[iat_sys].set_position(xyz) iat_cl = self.atom_map_sys_cl[iat_sys] if iat_cl > -1: self.atoms_cluster[iat_cl].set_position(xyz) for cell_L, iat_cl_sys, iat_sys, r, iat in self.linkatoms: # determine position of the link atom xyz_cl = positions_new[iat_cl_sys] xyz = positions_new[iat_sys] - xyz_cl + cell_L xyz *= r / np.sqrt(np.dot(xyz, xyz)) xyz += xyz_cl # update xyz coordinates of the cluster self.atoms_cluster[iat].set_position(xyz) def __getitem__(self, i): return self.atoms_system.__getitem__(i) def get_positions(self): return self.atoms_system.get_positions() def __add__(self, other): return self.atoms_system.__add__(other) def __delitem__(self, i): return self.atoms_system.__delitem__(i) def __len__(self): return self.atoms_system.__len__() def get_chemical_symbols(self, reduce=False): return self.atoms_system.get_chemical_symbols(reduce)