Python NeighborList.get_neighbors方法代码示例

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
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

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
    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:]

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
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

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
 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))]

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
    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

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
    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

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
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

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
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

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
def find_defects(solid, bulko, rcutoff, atomlistcheck = False, trackvacs = False, trackswaps = False, debug = False, dcheck = 0.6):
    """Function to find interstitials, vacancies, and substitutional atoms (swaps) in a defected structure.
    Identifies species by comparison to perfect structure.
    Inputs:
        solid = ASE atoms class for defected structure
        bulko = ASE atoms class for perfect structure
        rcutoff = float value of distance to surrounding atoms to include
        atomlistcheck = False/list of atom types and concentrations according to atomlist format
        trackvacs = True/False whether or not to identify vacancies in defect
        trackswaps = True/False whether or not to identify substitutional defects
        debug = False/file object to write debug structures"""

    # Combine perfect and defect structures together
    b = bulko.copy()
    b.extend(solid)
    b.set_pbc(True)

    # Debug: Write solid and bulko to file
    if debug:
        print len(bulko)
        write_xyz(debug, b, 'Find Ints: Solid and Bulko')

    # Identify nearest neighbor atoms for each atom in perfect structure
    ntot = len(bulko)
    ctoff1 = [1.2 for one in b]
    nl = NeighborList(ctoff1, bothways = True, self_interaction = False)
    nl.update(b)
    slist = []
    blist = []
    wlist = []

    # Loop over each atom in perfect structure
    for one in range(ntot):
        indices, offsets = nl.get_neighbors(one)
        for index, d in zip(indices, offsets):
            index = int(index)
            if index >= ntot:
                pos = b[index].position + numpy.dot(d, bulko.get_cell())
                natm1 = Atom(position = pos)
                dist, dx, dy, dz = calc_dist(b[one], natm1)
                if dist <= dcheck:
                    # Assume atoms closer than 0.6 Angstroms to be equivalent
                    slist.append(index-ntot)
                    blist.append(one)
                    if b[one].symbol == b[index].symbol:
                        wlist.append(index-ntot)

    # Identify those atoms corresponding to interstitials, vacancies, and substitutions
    oslist = [atm.index for atm in solid if atm.index not in slist]
    vlist = [atm.index for atm in bulko if atm.index not in blist]
    swlist = [atm.index for atm in solid if atm.index not in wlist and atm.index not in oslist]

    # Create Atoms objects for each identified defect
    ntot = len(solid)
    cluster = Atoms()
    for one in oslist:
        cluster.append(solid[one])
    vacant = Atoms()
    if trackvacs == True:
        for one in vlist:
            vacant.append(Atom(symbol = bulko[one].symbol, position = bulko[one].position))
            solid.append(Atom(symbol = 'X', position = bulko[one].position))
            oslist.append(len(solid)-1)
        stro = 'Cluster Identified with length = {0}\nIdentified {1} vacancies\n'.format(len(cluster), len(vlist))
    swaps = Atoms()
    if trackswaps == True:
        for one in swlist:
            swaps.append(solid[one])
            oslist.append(one)
        stro = 'Cluster Identified with length = {0}\nIdentified {1} swaps\n'.format(len(cluster), len(swlist))
    else:
        stro = 'Cluster Identified with length = {0}\n'.format(len(cluster))

    # Debug: write cluster to file
    if debug:
        b = cluster.copy()
        write_xyz(debug, b, 'Find Ints: Cluster')
        debug.flush()
        print('Found cluster size = ', len(b))

    # Identify atoms surrounding the identified defects in the defected structure
    box = Atoms()
    bulki = Atoms()
    if rcutoff != 0:
        if rcutoff > 2.0:
            cutoffs = [rcutoff for one in solid]
        else:
            cutoffs = [2.0 for one in solid]
        solid.set_pbc(True)
        nl = NeighborList(cutoffs, bothways = True, self_interaction = False)
        nl.update(solid)
        nbatmsd = []
        repinds = []
        for one in oslist:
            if one not in repinds:
                if one < ntot:
                    nbatmsd.append((0, one))
                    repinds.append(one)
        for one in oslist:
            indices, offsets = nl.get_neighbors(one)
#.........这里部分代码省略.........

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]

#.........这里部分代码省略.........
        # check we have all the properties requested
        for property in properties:
            if property not in self.results:
                if property is "energy":
                    self.calculate_energy(self.atoms)

                if property is "forces":
                    self.calculate_forces(self.atoms)

        # we need to remember the previous state of parameters

    #        if 'potential' in parameter_changes and potential != None:
    #                self.read_potential(potential)

    def calculate_energy(self, atoms):
        """Calculate the energy
        the energy is made up of the ionic or pair interaction and
        the embedding energy of each atom into the electron cloud
        generated by its neighbors
        """

        pair_energy = 0.0
        embedding_energy = 0.0
        mu_energy = 0.0
        lam_energy = 0.0
        trace_energy = 0.0

        self.total_density = np.zeros(len(atoms))
        if self.form == "adp":
            self.mu = np.zeros([len(atoms), 3])
            self.lam = np.zeros([len(atoms), 3, 3])

        for i in range(len(atoms)):  # this is the atom to be embedded
            neighbors, offsets = self.neighbors.get_neighbors(i)
            offset = np.dot(offsets, atoms.get_cell())

            rvec = atoms.positions[neighbors] + offset - atoms.positions[i]

            ## calculate the distance to the nearest neighbors
            r = np.sqrt(np.sum(np.square(rvec), axis=1))  # fast
            #            r = np.apply_along_axis(np.linalg.norm, 1, rvec)  # sloow

            nearest = np.arange(len(r))[r <= self.cutoff]
            for j_index in range(self.Nelements):
                use = self.index[neighbors[nearest]] == j_index
                if not use.any():
                    continue
                pair_energy += np.sum(self.phi[self.index[i], j_index](r[nearest][use])) / 2.0

                density = np.sum(self.electron_density[j_index](r[nearest][use]))
                self.total_density[i] += density

                if self.form == "adp":
                    self.mu[i] += self.adp_dipole(r[nearest][use], rvec[nearest][use], self.d[self.index[i], j_index])

                    self.lam[i] += self.adp_quadrupole(
                        r[nearest][use], rvec[nearest][use], self.q[self.index[i], j_index]
                    )

            # add in the electron embedding energy
            embedding_energy += self.embedded_energy[self.index[i]](self.total_density[i])

        components = dict(pair=pair_energy, embedding=embedding_energy)

        if self.form == "adp":
            mu_energy += np.sum(self.mu ** 2) / 2.0

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]

#.........这里部分代码省略.........
                               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 get_forces(self, atoms):
        # calculate the forces based on derivatives of the three EAM functions
        self.update(atoms)
        self.forces = np.zeros((len(atoms), 3))

        for i in range(len(atoms)):  # this is the atom to be embedded
            neighbors, offsets = self.neighbors.get_neighbors(i)
            offset = np.dot(offsets, atoms.get_cell())
            # create a vector of relative positions of neighbors
            rvec = atoms.positions[neighbors] + offset - atoms.positions[i]
            r = np.sqrt(np.sum(np.square(rvec), axis=1))
            nearest = np.arange(len(r))[r < self.cutoff]

            d_embedded_energy_i = self.d_embedded_energy[
                self.index[i]](self.total_density[i])
            urvec = rvec.copy()  # unit directional vector

            for j in np.arange(len(neighbors)):
                urvec[j] = urvec[j] / r[j]
    
            for j_index in range(self.Nelements):
                use = self.index[neighbors[nearest]] == j_index
                if not use.any():
                    continue
                rnuse = r[nearest][use]
                density_j = self.total_density[neighbors[nearest][use]]
                scale = (self.d_phi[self.index[i], j_index](rnuse) +
                         (d_embedded_energy_i *
                          self.d_electron_density[j_index](rnuse)) +
                         (self.d_embedded_energy[j_index](density_j) *
                          self.d_electron_density[self.index[i]](rnuse)))
            
                self.forces[i] += np.dot(scale, urvec[nearest][use])
                
                if (self.form == 'adp'):
                    adp_forces = self.angular_forces(
                        self.mu[i],
                        self.mu[neighbors[nearest][use]],
                        self.lam[i],

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
class EMT(Calculator):
    implemented_properties = ['energy', 'forces']

    nolabel = True

    def __init__(self):
        Calculator.__init__(self)

    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 calculate(self, atoms=None, properties=['energy'],
                  system_changes=all_changes):
        Calculator.calculate(self, atoms, properties, system_changes)

        if 'numbers' in system_changes:
            self.initialize(self.atoms)

        positions = self.atoms.positions
        numbers = self.atoms.numbers
        cell = self.atoms.cell
        
        self.nl.update(self.atoms)
        
        self.energy = 0.0
        self.sigma1[:] = 0.0
        self.forces[:] = 0.0

        natoms = len(self.atoms)

        for a1 in range(natoms):
            Z1 = numbers[a1]
            p1 = self.par[Z1]
            ksi = self.ksi[Z1]
            neighbors, offsets = self.nl.get_neighbors(a1)
            offsets = np.dot(offsets, cell)
            for a2, offset in zip(neighbors, offsets):
                d = positions[a2] + offset - positions[a1]
                r = sqrt(np.dot(d, d))
                if r < self.rc + 0.5:
                    Z2 = numbers[a2]
                    p2 = self.par[Z2]
                    self.interact1(a1, a2, d, r, p1, p2, ksi[Z2])
                                
        for a in range(natoms):
            Z = numbers[a]
            p = self.par[Z]
            try:
                ds = -log(self.sigma1[a] / 12) / (beta * p['eta2'])
            except (OverflowError, ValueError):
                self.deds[a] = 0.0
                self.energy -= p['E0']
#.........这里部分代码省略.........

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
def swap_int_local(indiv, Optimizer):
    """Move function to perform n atom structure swap based on swaplist for cluster in local region
    Inputs:
        indiv = Individual class object to be altered
        Optimizer = Optimizer class object with needed parameters
    Outputs:
        indiv = Altered Individual class object
    """
    if 'MU' in Optimizer.debug:
        debug = True
    else:
        debug = False
    Optimizer.output.write('Local Cluster Swap Mutation performed on individual\n')
    Optimizer.output.write('Index = '+repr(indiv.index)+'\n')
    # Identify the local atoms available for swapping
    atmsst = indiv[0].copy()
    nat = len(atmsst)
    totalsol = atmsst.copy()
    if Optimizer.structure=='Defect':
        totalsol.extend(indiv.bulki.copy())
        ctoff = [Optimizer.sf for one in totalsol]
    else:
        ctoff = [1.5 for one in totalsol]
    nl = NeighborList(ctoff, bothways=True, self_interaction=False)
    nl.update(totalsol)
    cluster = Atoms()
    indices = [] #store indices of the atoms for later
    for i in range(nat):
        indices1, offsets = nl.get_neighbors(i)
        for index, d in zip(indices1,offsets):
            index = int(index)
            #Make sure to prevent repeats of the same atom
            if index not in indices:
                pos = totalsol[index].position + numpy.dot(d,totalsol.get_cell())
                cluster.append(Atom(symbol=totalsol[index].symbol,position=pos))
                indices.append(index)
    # Identify Number of atoms to swap
    swapmax=sum([c for sym,c in indiv.swaplist])
    if swapmax >1:
        natomsswap=random.randint(1,swapmax)
    elif swapmax==1:
        natomsswap=1
    else:
        Optimizer.output.write('Maximum number of swaps is '+repr(swapmax)+'. Unable to perform mutation')
        natomsswap=0
    if len(indiv.swaplist)==1:
        Optimizer.output.write('WARNING: Swap Mutation attempted on single atom structure system\n')
        natomsswap=0
    Optimizer.output.write('Number of swaps = '+repr(natomsswap)+'\n')
    # Identify symbols that can be swapped
    syms=[sym for sym,c in indiv.swaplist]
    slist = copy.deepcopy(indiv.swaplist)
    if debug: print 'Starting swaplist = '+repr(indiv.swaplist)
    #Sanity check
    sanch = [[sym,0] for sym in syms]
    for i in range(len(sanch)):
        nc = len([atm for atm in cluster if atm.symbol==sanch[i][0]])
        nc += [c for sym,c in slist if sym==sanch[i][0]][0]
        sanch[i][1]=nc
    # Swap Atoms
    for i in range(natomsswap):
        while True:
            at1 = random.choice(cluster)
            osym = at1.symbol
            nsyml=[sym for sym,c in slist if sym != osym and c > 0]
            if len(nsyml) > 0:
                break
        nsym = random.choice(nsyml)
        cluster[at1.index].symbol = nsym
        Optimizer.output.write('Swapped '+nsym+' atom with '+osym+'\n')
        for i in range(len(slist)):
            if slist[i][0]==nsym:
                slist[i][1]-=1
            elif slist[i][0]==osym:
                slist[i][1]+=1
    # Apply the new swap to the actual atoms in the structure
    for i in range(len(indices)):
        totalsol[indices[i]].symbol = cluster[i].symbol
    # Separate out bulk from defects
    indiv[0] = totalsol[0:nat]
    if Optimizer.structure=='Defect':
        indiv.bulki = totalsol[nat::]
    # Update new swaplist
    indiv.swaplist = slist
    #Sanity check
    sanchn = [[sym,0] for sym in syms]
    for i in range(len(sanchn)):
        nc = len([atm for atm in cluster if atm.symbol==sanchn[i][0]])
        nc += [c for sym,c in slist if sym==sanchn[i][0]][0]
        sanchn[i][1]=nc
    if debug:
        for i in range(len(sanch)):
            print 'INTSWAP: Starting Atom structure '+repr(sanch[i][0])+' = '+repr(sanch[i][1])
            print 'INTSWAP: After Mutation Atom structure '+repr(sanchn[i][0])+' = '+repr(sanchn[i][1])
    Optimizer.output.write(repr(indiv[0])+'\n')
    muttype='SCL'+repr(natomsswap)
    if indiv.energy==0:
        indiv.history_index=indiv.history_index+'m'+muttype
    else:
        indiv.history_index=repr(indiv.index)+'m'+muttype
#.........这里部分代码省略.........

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]

#.........这里部分代码省略.........

        # angles
        atypes, alist = self.get_angles()
        fileobj.write('\n' + str(len(atypes)) + ' angle types\n')
        fileobj.write(str(len(alist)) + ' angles\n')
        fileobj.write('\nAngles\n\n')
        
        for ia, avals in enumerate(alist):
            fileobj.write('%8d %6d %6d %6d %6d\n' %
                          (ia + 1, avals[0] + 1, 
                           avals[1] + 1, avals[2] + 1, avals[3] + 1))

        return btypes, atypes

    def update_neighbor_list(self, atoms):
        cut = 0.5 * max(self.data['cutoffs'].values())
        self.nl = NeighborList([cut] * len(atoms), skin=0, bothways=True)
        self.nl.update(atoms)
        self.atoms = atoms
    
    def get_bonds(self, atoms):
        """Find bonds and return them and their types"""
        cutoffs = CutoffList(self.data['cutoffs'])
        self.update_neighbor_list(atoms)

        types = atoms.get_types()
        tags = atoms.get_tags()
        cell = atoms.get_cell()
        positions = atoms.get_positions()
        bond_list = []
        bond_types = []
        for i, atom in enumerate(atoms):
            iname = types[tags[i]]
            indices, offsets = self.nl.get_neighbors(i)
            for j, offset in zip(indices, offsets):
                if j <= i:
                    continue # do not double count
                jname = types[tags[j]]
                cut = cutoffs.value(iname, jname)
                if cut is None:
                    if self.warnings > 1:
                        print ('Warning: cutoff %s-%s not found'
                               % (iname, jname))
                    continue # don't have it
                dist = np.linalg.norm(atom.position - atoms[j].position
                                      - np.dot(offset, cell))
                if dist > cut:
                    continue # too far away
                name, val = self.bonds.name_value(iname, jname)
                if name is None:
                    if self.warnings:
                        print ('Warning: potential %s-%s not found'
                               % (iname, jname))
                    continue # don't have it
                if name not in bond_types:
                    bond_types.append(name)
                bond_list.append([bond_types.index(name), i, j])
        return bond_types, bond_list
                
    def get_angles(self, atoms=None):
        cutoffs = CutoffList(self.data['cutoffs'])
        if atoms is not None:
            self.update_neighbor_list(atoms)
        else:
            atoms = self.atoms
         

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import get_neighbors [as 别名]
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':
        rank = MPI.COMM_WORLD.Get_rank()
        logger = logging.getLogger(Optimizer.loggername)
        R = totalsol.arrays['positions']
        tol = 0.01
        epsilon = 0.05
        fix = 0.5
        if Optimizer.forcing == 'EllipoidShape' or Optimizer.forcing == 'FreeNatom': 
          com = totalsol.get_center_of_mass()       
          cmax = numpy.maximum.reduce(R)
          cmin = numpy.minimum.reduce(R)
          rmax= (cmax-cmin)/2.0 
          if Optimizer.forcing == 'FreeNatom':
             rcutoff = 44.0
             cutoff = [44.0,44.0,20.0]        
          else:
             rcutoff = 11.0
             cutoff = [12.0,12.0,12.0]        
          rcutoff = 44.0
          cutoff = [44.0,44.0,20.0]        
          #check if atoms are isolated outside of cluster
          cutoffs=[3.0 for one in totalsol]
          nl=NeighborList(cutoffs,bothways=True,self_interaction=False)
          nl.update(totalsol)
          for i in range(len(totalsol)):
             indices, offsets=nl.get_neighbors(i)
             D = R[i]-com
             radius = (numpy.dot(D,D))**0.5 #numpy.linalg.norm(D)
             if len(indices) < 12 or radius > rcutoff :
                # logger.info('M:move atoms back when indice {0} or radius {1}'.format(len(indices),radius))
                # R[i] = [com[j] + D[j]/radius*rcutoff for j in range(3)]
                theta=math.radians(random.uniform(0,360))
                phi=math.radians(random.uniform(0,180))
                R[i][0] = com[0] + (rmax[0]+2.5)*math.sin(theta)*math.cos(phi) #allow atoms expend by 2.5 ang
                R[i][1] = com[1] + (rmax[1]+2.5)*math.sin(theta)*math.sin(phi)
                R[i][2] = com[2] + rmax[2]*math.cos(theta)
                # logger.info('M:move atoms back new pos {0} {1} {2}'.format(rmax[0]*math.sin(theta)*math.cos(phi),rmax[1]*math.sin(theta)*math.sin(phi),rmax[2]*math.cos(theta)))
               
          # check if atoms are within cluster region 
          for i in range(0,len(totalsol)):
            # D = R[i]-com
            for j in range(3):                 
               if D[j] > cutoff[j] or D[j] < -cutoff[j]:
         #         logger.info('M:before move R {0} & com {1}'.format(R[i][j],com[j]))
                  #if rank==0:
                  #   print "before:",j,R[i][j],com[j] 
                  R[i][j] = com[j]+numpy.sign(D[j])*cutoff[j]*random.random()
         #         logger.info('M:after move R {0} '.format(R[i][j]))
                  #if rank==0:
                  #   print "after:",R[i][j]
              #    STR+='--- WARNING: Atoms too far along x-y (>44A) - Implement Move ---\n'          
                  D = R[i]-com
        #    radius = (numpy.dot(D,D))**0.5 #numpy.linalg.norm(D)
             #  STR+='--- WARNING: Atoms too far (>56A) - Implement Move ---\n'          
              

        closelist = numpy.arange(len(totalsol))
        iter = 0
        while len(closelist) > 0 and iter<2:
          iter+=1 
         # checklist = numpy.copy(closelist)
          closelist = []  
          dist=scipy.spatial.distance.cdist(R,R)       
          numpy.fill_diagonal(dist,1.0)
          smalldist = numpy.where(dist < min_len-tol)
         # for i in checklist:
            # for j in range(i+1,len(totalsol)):
         #    if len(checklist) == len(totalsol):
         #       jstart = i+1
#.........这里部分代码省略.........