Python neighborlist.NeighborList类代码示例

 def __init__(self, filename):
     atoms = ase.io.read(filename, format='vasp')
     cell = atoms.get_cell()
     symb = atoms.get_atomic_numbers()
     
     cutoffs = np.ones((len(atoms))) * 1.35  # radius around each atom (half the max bondlength)
     nl = NeighborList(cutoffs, 0., self_interaction = False, bothways = False)
     nl.update(atoms)
     
     std_len = 2.35805097505
     
     # Loop over atoms
     bond_no = 0
     bonds = np.zeros((nl.nneighbors,4))
     
     for atom1idx, (indices, offsets, atom1pos) in enumerate(zip(nl.neighbors, nl.displacements, nl.positions)):
         # Loop over bonds
         for atom2idx, offset in zip(indices, offsets):
             #if symb[atom1idx] == symb[atom2idx] == 'Si': 
             atom2pos = nl.positions[atom2idx] + np.dot(offset, cell)
             #print "pos[%i] = %.2f %.2f %.2f" % (atom1idx, atom1pos[0],atom1pos[1],atom1pos[2])
             #print "pos[%i] = %.2f %.2f %.2f" % (atom2idx, atom2pos[0],atom2pos[1],atom2pos[2])
             bond = atom2pos - atom1pos
             bond_center = atom1pos + bond/2.
             bond_z = bond_center[2]
             bondlength = np.sqrt(np.dot(bond,bond))    
             bonds[bond_no] = [bond_z, bondlength, symb[atom1idx], symb[atom2idx]]
             bond_no += 1
     self.bonds = bonds
     print bonds

def get_bondpairs(atoms, radius=1.1):
    """Get all pairs of bonding atoms

    Return all pairs of atoms which are closer than radius times the
    sum of their respective covalent radii.  The pairs are returned as
    tuples::

      (a, b, (i1, i2, i3))

    so that atoms a bonds to atom b displaced by the vector::

        _     _     _
      i c + i c + i c ,
       1 1   2 2   3 3

    where c1, c2 and c3 are the unit cell vectors and i1, i2, i3 are
    integers."""
    
    from ase.data import covalent_radii
    from ase.calculators.neighborlist import NeighborList
    cutoffs = radius * covalent_radii[atoms.numbers]
    nl = NeighborList(cutoffs=cutoffs, self_interaction=False)
    nl.update(atoms)
    bondpairs = []
    for a in range(len(atoms)):
        indices, offsets = nl.get_neighbors(a)
        bondpairs.extend([(a, a2, offset)
                          for a2, offset in zip(indices, offsets)])
    return bondpairs

def find_tip_coordination(a, bondlength=2.6, bulk_nn=4):
    """
    Find position of tip in crack cluster from coordination
    """
    nl = NeighborList([bondlength/2.0]*len(a),
                      skin=0.0,
                      self_interaction=False,
                      bothways=True)
    nl.update(a)
    nn = np.array([len(nl.get_neighbors(i)[0]) for i in range(len(a))])
    a.set_array('n_neighb', nn)
    g = a.get_array('groups')

    y = a.positions[:, 1]
    above = (nn < bulk_nn) & (g != 0) & (y > a.cell[1,1]/2.0)
    below = (nn < bulk_nn) & (g != 0) & (y < a.cell[1,1]/2.0)

    a.set_array('above', above)
    a.set_array('below', below)

    bond1 = np.asscalar(above.nonzero()[0][a.positions[above, 0].argmax()])
    bond2 = np.asscalar(below.nonzero()[0][a.positions[below, 0].argmax()])

    # These need to be ints, otherwise they are no JSON serializable.
    a.info['bond1'] = bond1
    a.info['bond2'] = bond2

    return bond1, bond2

    def bind(self, frame):
        if not self.ui.get_widget('/MenuBar/ViewMenu/ShowBonds'
                                  ).get_active():
            self.bonds = np.empty((0, 5), int)
            return
        
        from ase.atoms import Atoms
        from ase.calculators.neighborlist import NeighborList
        nl = NeighborList(self.images.r * 1.5, skin=0, self_interaction=False)
        nl.update(Atoms(positions=self.images.P[frame],
                        cell=(self.images.repeat[:, np.newaxis] *
                              self.images.A[frame]),
                        pbc=self.images.pbc))
        nb = nl.nneighbors + nl.npbcneighbors
        self.bonds = np.empty((nb, 5), int)
        if nb == 0:
            return
        
        n1 = 0
        for a in range(self.images.natoms):
            indices, offsets = nl.get_neighbors(a)
            n2 = n1 + len(indices)
            self.bonds[n1:n2, 0] = a
            self.bonds[n1:n2, 1] = indices
            self.bonds[n1:n2, 2:] = offsets
            n1 = n2

        i = self.bonds[:n2, 2:].any(1)
        self.bonds[n2:, 0] = self.bonds[i, 1]
        self.bonds[n2:, 1] = self.bonds[i, 0]
        self.bonds[n2:, 2:] = -self.bonds[i, 2:]

 def calculate(self, image, key):
     cutoff = self.globals.cutoff
     n = NeighborList(cutoffs=[cutoff / 2.] * len(image),
                      self_interaction=False,
                      bothways=True,
                      skin=0.)
     n.update(image)
     return [n.get_neighbors(index) for index in xrange(len(image))]

 def buildNeighborList(self):
     atoms = self.atoms
     nl = NeighborList([0.8 for atom in atoms],
                       self_interaction=False, bothways=True)
     nl.update(atoms)
     self.nl = nl
     neigh = []
     for i in range(self.natom):
         neigh.append(self.sort_nei(i))
         # print neigh[i]
     self.neigh = neigh

	def testNeighborlist(self):
		atoms=read('range',format='lammps')
		atoms. set_pbc((1, 1, 1))
		nl = NeighborList([0.8 for atom in atoms],self_interaction=False,bothways=True)
		nl.update(atoms)
		ang=[]		
		for i in xrange(3):
			indices, offsets = nl. get_neighbors(i)
			angs=[]
			for j, offset in zip(indices, offsets):
				pos= atoms. positions[j] + dot(offset, atoms. get_cell())-atoms.positions[i]
				ang1=atan2(pos[1],pos[0])+pi
				angs.append((j,ang1))
			newangs=sorted(angs,key=lambda d:d[1])
			print newangs

    def update(self, atoms):
        # check all the elements are available in the potential
        elements = np.unique(atoms.get_chemical_symbols())
        unavailable = np.logical_not(
            np.array([item in self.elements for item in elements]))
        if np.any(unavailable):
            raise RuntimeError('These elements are not in the potential: %s' %
                               elements[unavailable])
                    
        # check if anything has changed to require re-calculation
        if (self.positions is None or
            len(self.positions) != len(atoms.get_positions()) or
            (self.positions != atoms.get_positions()).any() or
            (self.cell != atoms.get_cell()).any() or
            (self.pbc != atoms.get_pbc()).any()):

            # cutoffs need to be a vector for NeighborList
            cutoffs = self.cutoff * np.ones(len(atoms))

            # convert the elements to an index of the position
            # in the eam format
            self.index = np.array([atoms.calc.elements.index(el)
                                   for el in atoms.get_chemical_symbols()])
            self.pbc = atoms.get_pbc()

            # since we need the contribution of all neighbors to the
            # local electron density we cannot just calculate and use
            # one way neighbors
            self.neighbors = NeighborList(cutoffs, skin=self.skin,
                                          self_interaction=False,
                                          bothways=True)
            self.neighbors.update(atoms)
            self.calculate(atoms)

    def find_connected(self, index, dmax=None, scale=1.5):
        """Find the atoms connected to self[index] and return them.

        If dmax is not None:
        Atoms are defined to be connected if they are nearer than dmax
        to each other.

        If dmax is None:
        Atoms are defined to be connected if they are nearer than the
        sum of their covalent radii * scale to each other.

        """

        # set neighbor lists
        neighborlist = []
        if dmax is None:
            # define neighbors according to covalent radii
            radii = scale * covalent_radii[self.get_atomic_numbers()]
            for atom in self:
                positions = self.positions - atom.position
                distances = np.sqrt(np.sum(positions**2, axis=1))
                radius = scale * covalent_radii[atom.get_atomic_number()]
                neighborlist.append(np.where(distances < radii + radius)[0])
        else:
            # define neighbors according to distance
            nl = NeighborList([0.5 * dmax] * len(self), skin=0)
            nl.update(self)
            for i, atom in enumerate(self):
                neighborlist.append(list(nl.get_neighbors(i)[0]))

        connected = list(neighborlist[index])
        isolated = False
        while not isolated:
            isolated = True
            for i in connected:
                for j in neighborlist[i]:
                    if j in connected:
                        pass
                    else:
                        connected.append(j)
                        isolated = False 
 
        atoms = Cluster()
        for i in connected:
            atoms.append(self[i])

        return atoms

    def find_connected(self, index, dmax=None, scale=1.5):
        """Find the atoms connected to self[index] and return them.

        If dmax is not None:
        Atoms are defined to be connected if they are nearer than dmax
        to each other.

        If dmax is None:
        Atoms are defined to be connected if they are nearer than the
        sum of their covalent radii * scale to each other.

        """

        if index < 0:
            index = len(self) + index

        # set neighbor lists
        if dmax is None:
            # define neighbors according to covalent radii
            radii = scale * covalent_radii[self.get_atomic_numbers()]
        else:
            # define neighbors according to distance
            radii = [0.5 * dmax] * len(self)
        nl = NeighborList(radii, skin=0, self_interaction=False, bothways=True)
        nl.update(self)

        connected = [index] + list(nl.get_neighbors(index)[0])
        isolated = False
        while not isolated:
            isolated = True
            for i in connected:
                for j in nl.get_neighbors(i)[0]:
                    if j in connected:
                        pass
                    else:
                        connected.append(j)
                        isolated = False

        atoms = Cluster()
        for i in connected:
            atoms.append(self[i])

        return atoms

    def calculate(self, atoms=None,
                  properties=['energy'],
                  system_changes=all_changes):
        Calculator.calculate(self, atoms, properties, system_changes)

        natoms = len(self.atoms)

        sigma = self.parameters.sigma
        epsilon = self.parameters.epsilon
        rc = self.parameters.rc
        if rc is None:
            rc = 3 * sigma
        
        if 'numbers' in system_changes:
            self.nl = NeighborList([rc / 2] * natoms, self_interaction=False)

        self.nl.update(self.atoms)
        
        positions = self.atoms.positions
        cell = self.atoms.cell
        
        e0 = 4 * epsilon * ((sigma / rc)**12 - (sigma / rc)**6)
        
        energy = 0.0
        forces = np.zeros((natoms, 3))
        stress = np.zeros((3, 3))

        for a1 in range(natoms):
            neighbors, offsets = self.nl.get_neighbors(a1)
            cells = np.dot(offsets, cell)
            d = positions[neighbors] + cells - positions[a1]
            r2 = (d**2).sum(1)
            c6 = (sigma**2 / r2)**3
            c6[r2 > rc**2] = 0.0
            energy -= e0 * (c6 != 0.0).sum()
            c12 = c6**2
            energy += 4 * epsilon * (c12 - c6).sum()
            f = (24 * epsilon * (2 * c12 - c6) / r2)[:, np.newaxis] * d
            #print d
            #print r2**.5
            #print offsets
            #print f
            #print neighbors
            forces[a1] -= f.sum(axis=0)
            for a2, f2 in zip(neighbors, f):
                forces[a2] += f2
            stress += np.dot(f.T, d)
        
        #stress = np.dot(stress, cell)
        stress += stress.T.copy()
        stress *= -0.5 / self.atoms.get_volume()
        
        self.results['energy'] = energy
        self.results['forces'] = forces
        self.results['stress'] = stress.flat[[0, 4, 8, 5, 2, 1]]

class NeighborPairs:
    """Class for looping over pairs of atoms using a neighbor list."""
    def __init__(self, cutoff_a, cell_cv, pbc_c, self_interaction):
        self.neighbors = NeighborList(cutoff_a, skin=0, sorted=True, 
                                      self_interaction=self_interaction)
        self.atoms = Atoms('X%d' % len(cutoff_a), cell=cell_cv, pbc=pbc_c)
        # Warning: never use self.atoms.get_scaled_positions() for
        # anything.  Those positions suffer from roundoff errors!
        
    def set_positions(self, spos_ac):
        self.spos_ac = spos_ac
        self.atoms.set_scaled_positions(spos_ac)
        self.neighbors.update(self.atoms)

    def iter(self):
        cell_cv = self.atoms.cell
        for a1, spos1_c in enumerate(self.spos_ac):
            a2_a, offsets = self.neighbors.get_neighbors(a1)
            for a2, offset in zip(a2_a, offsets):
                spos2_c = self.spos_ac[a2] + offset
                R_c = np.dot(spos2_c - spos1_c, cell_cv)
                yield a1, a2, R_c, offset

    def initialize(self, atoms):
        self.par = {}
        self.rc = 0.0
        self.numbers = atoms.get_atomic_numbers()
        maxseq = max(par[1] for par in parameters.values()) * Bohr
        rc = self.rc = beta * maxseq * 0.5 * (sqrt(3) + sqrt(4))
        rr = rc * 2 * sqrt(4) / (sqrt(3) + sqrt(4))
        self.acut = np.log(9999.0) / (rr - rc)
        for Z in self.numbers:
            if Z not in self.par:
                p = parameters[chemical_symbols[Z]]
                s0 = p[1] * Bohr
                eta2 = p[3] / Bohr
                kappa = p[4] / Bohr
                x = eta2 * beta * s0
                gamma1 = 0.0
                gamma2 = 0.0
                for i, n in enumerate([12, 6, 24]):
                    r = s0 * beta * sqrt(i + 1)
                    x = n / (12 * (1.0 + exp(self.acut * (r - rc))))
                    gamma1 += x * exp(-eta2 * (r - beta * s0))
                    gamma2 += x * exp(-kappa / beta * (r - beta * s0))

                self.par[Z] = {'E0': p[0],
                               's0': s0,
                               'V0': p[2],
                               'eta2': eta2,
                               'kappa': kappa,
                               'lambda': p[5] / Bohr,
                               'n0': p[6] / Bohr**3,
                               'rc': rc,
                               'gamma1': gamma1,
                               'gamma2': gamma2}
                #if rc + 0.5 > self.rc:
                #    self.rc = rc + 0.5

        self.ksi = {}
        for s1, p1 in self.par.items():
            self.ksi[s1] = {}
            for s2, p2 in self.par.items():
                #self.ksi[s1][s2] = (p2['n0'] / p1['n0'] *
                #                    exp(eta1 * (p1['s0'] - p2['s0'])))
                self.ksi[s1][s2] = p2['n0'] / p1['n0']
                
        self.forces = np.empty((len(atoms), 3))
        self.sigma1 = np.empty(len(atoms))
        self.deds = np.empty(len(atoms))
                    
        self.nl = NeighborList([0.5 * self.rc + 0.25] * len(atoms),
                               self_interaction=False)

def check_min_dist(totalsol, type='Defect', nat=None, min_len=0.7, STR=''):
    if type=='Defect' or type=='Crystal' or type=='Surface':
        if nat==None:
            nat=len(totalsol)
        cutoffs=[2.0 for one in totalsol]
        nl=NeighborList(cutoffs,bothways=True,self_interaction=False)
        nl.update(totalsol)
        for one in totalsol[0:nat]:
            nbatoms=Atoms()
            nbatoms.append(one)
            indices, offsets=nl.get_neighbors(one.index)
            for index, d in zip(indices,offsets):
                index = int(index)
                sym=totalsol[index].symbol
                pos=totalsol[index].position + numpy.dot(d,totalsol.get_cell())
                at=Atom(symbol=sym,position=pos)
                nbatoms.append(at)
            while True:
                dflag=False
                for i in range(1,len(nbatoms)):
                    d=nbatoms.get_distance(0,i)
                    if d < min_len:
                        nbatoms.set_distance(0,i,min_len+.01,fix=0.5)
                        STR+='--- WARNING: Atoms too close (<0.7A) - Implement Move ---\n'
                        dflag=True
                if dflag==False:
                    break
            for i in range(len(indices)):
                totalsol[indices[i]].position=nbatoms[i+1].position
            totalsol[one.index].position=nbatoms[0].position
            nl.update(totalsol)
    elif type=='Cluster':
        for i in range(len(totalsol)):
            for j in range(len(totalsol)):
                if i != j:
                    d=totalsol.get_distance(i,j)
                    if d < min_len:
                        totalsol.set_distance(i,j,min_len,fix=0.5)
                        STR+='--- WARNING: Atoms too close (<0.7A) - Implement Move ---\n'
    else:
        print 'WARNING: In Check_Min_Dist in EvalEnergy: Structure Type not recognized'
    return totalsol, STR

def check_min_dist(Optimizer, totalsol, type='Defect', nat=None, min_len=0.7, STR=''):
        if type=='Defect' or type=='Crystal' or type=='Surface':
            if nat==None:
                nat=len(totalsol)
            cutoffs=[2.0 for one in totalsol]
            nl=NeighborList(cutoffs,bothways=True,self_interaction=False)
            nl.update(totalsol)
            for one in totalsol[0:nat]:
                nbatoms=Atoms()
                nbatoms.append(one)
                indices, offsets=nl.get_neighbors(one.index)
                for index, d in zip(indices,offsets):
                    index = int(index)
                    sym=totalsol[index].symbol
                    pos=totalsol[index].position + numpy.dot(d,totalsol.get_cell())
                    at=Atom(symbol=sym,position=pos)
                    nbatoms.append(at)
                while True:
                    dflag=False
                    for i in range(1,len(nbatoms)):
                        d=nbatoms.get_distance(0,i)
                        if d < min_len:
                            nbatoms.set_distance(0,i,min_len+.01,fix=0.5)
                            STR+='--- WARNING: Atoms too close (<0.7A) - Implement Move ---\n'
                            dflag=True
                    if dflag==False:
                        break
                for i in range(len(indices)):
                    totalsol[indices[i]].position=nbatoms[i+1].position
                totalsol[one.index].position=nbatoms[0].position
                nl.update(totalsol)
        elif type=='Cluster':
            if not 'LAMMPS' in Optimizer.modules:
                for i in range(len(totalsol)):
                    for j in range(len(totalsol)):
                        if i != j:
                            d=totalsol.get_distance(i,j)
                            if d < min_len:
                                totalsol.set_distance(i,j,min_len,fix=0.5)
                                STR+='--- WARNING: Atoms too close (<0.7A) - Implement Move ---\n'
            else:
                rank = MPI.COMM_WORLD.Get_rank()
                #logger = logging.getLogger(Optimizer.loggername)
                R = totalsol.arrays['positions']
                tol = 0.01
                epsilon = 0.05
                fix = 0.5
                
                closelist = numpy.arange(len(totalsol))
                iter = 0
                while len(closelist) > 0 and iter<2:
                    iter+=1
                    closelist = []
                    dist=spatial.distance.cdist(R,R)
                    numpy.fill_diagonal(dist,1.0)
                    smalldist = numpy.where(dist < min_len-tol)
   
                    for ind in range(len(smalldist[0])):
                        i = smalldist[0][ind]
                        j = smalldist[1][ind]
                        if i < j and dist[i][j] < min_len-tol:
                            closelist.append(i)
                            closelist.append(j)
                            if dist[i][j] > epsilon:
                                x = 1.0 - min_len / dist[i][j]
                                D = R[j]-R[i]
                                R[i] += (x * fix) * D
                                R[j] -= (x * (1.0 - fix)) * D
                            else:
                                R[i] += [0.2, 0.0, 0.0]
                                R[j] -= [0.2, 0.0, 0.0]
                            R2P = [R[i],R[j]]
                            dist2P=spatial.distance.cdist(R2P,R)
                            dist[i] = dist2P[0]
                            dist[j] = dist2P[1]
                            for k in range(len(R)):
                                dist[k][i] = dist[i][k]
                                dist[k][j] = dist[j][k]
                    closelist=list(set(closelist))
                    closelist.sort()
                    #if len(closelist) != 0:
                    #    logger.info('M:iter {0}, closelist size {1}'.format(iter,len(closelist)))
        else:
            print 'WARNING: In Check_Min_Dist in EvalEnergy: Structure Type not recognized'
        return totalsol, STR