Python NeighborList.update方法代码示例

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

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

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

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import update [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 update [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 update [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 update [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 update [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 update [as 别名]
class OPLSff:
    def __init__(self, fileobj=None, warnings=0):
        self.warnings = warnings
        self.data = {}
        if fileobj is not None:
            self.read(fileobj)

    def read(self, fileobj, comments='#'):
        if isinstance(fileobj, str):
            fileobj = open(fileobj)

        def read_block(name, 
                       symlen, # length of the symbol
                       nvalues # of values expected
                       ):
            if name not in self.data:
                self.data[name] = {}
            data = self.data[name]

            def add_line():
                line = fileobj.readline()
                if len(line) <= 1: # end of the block
                    return False
                line = line.split('#')[0] # get rid of comments
                if len(line) > symlen:
                    symbol = line[:symlen]
                    words = line[symlen:].split()
                    if len(words) >=  nvalues:
                        if nvalues == 1:
                            data[symbol] = float(words[0])
                        else:
                            data[symbol] = [float(word) 
                                            for word in words[:nvalues]]
                return True

            while add_line():
                pass
 
        read_block('one',      2, 3)
        read_block('bonds',      5, 2)
        read_block('angles',     8, 2)
        read_block('dihedrals', 11, 3)

        read_block('cutoffs',      5, 1)

        self.bonds = BondData(self.data['bonds'])
        self.angles = AnglesData(self.data['angles'])
        self.cutoffs = CutoffList(self.data['cutoffs'])

    def write_lammps(self, atoms):
        btypes, atypes = self.write_lammps_atoms(atoms)
        self.write_lammps_definitions(atoms, btypes, atypes)

    def write_lammps_atoms(self, atoms):
        """Write atoms infor for LAMMPS"""
        
        fileobj = 'lammps_atoms'
        if isinstance(fileobj, str):
            fileobj = open(fileobj, 'w')
        write_lammps_data(fileobj, atoms, 
                          specorder=atoms.types,
                          speclist=(atoms.get_tags() + 1),
                          )

        # masses
        fileobj.write('\nMasses\n\n')
        for i, typ in enumerate(atoms.types):
            cs = atoms.split_symbol(typ)[0]
            fileobj.write('%6d %g # %s -> %s\n' % 
                          (i + 1, 
                           atomic_masses[chemical_symbols.index(cs)],
                           typ, cs))
  
        # bonds
        btypes, blist = self.get_bonds(atoms)
        fileobj.write('\n' + str(len(btypes)) + ' bond types\n')
        fileobj.write(str(len(blist)) + ' bonds\n')
        fileobj.write('\nBonds\n\n')
        
        for ib, bvals in enumerate(blist):
            fileobj.write('%8d %6d %6d %6d\n' %
                          (ib + 1, bvals[0] + 1, bvals[1] + 1, bvals[2] + 1))

        # 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)
#.........这里部分代码省略.........

# 需要导入模块: from ase.calculators.neighborlist import NeighborList [as 别名]
# 或者: from ase.calculators.neighborlist.NeighborList import update [as 别名]
class LennardJones(Calculator):
    implemented_properties = ['energy', 'forces', 'stress']
    default_parameters = {'epsilon': 1.0,
                          'sigma': 1.0,
                          'rc': None}
    nolabel = True

    def __init__(self, **kwargs):
        Calculator.__init__(self, **kwargs)

    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]]

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

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

The breakdown of energy contribution from the indvidual components are
stored in the calculator instance ``.results['energy_components']``

Arguments
=========

=========================  ====================================================
Keyword                    Description
=========================  ====================================================
``potential``              file of potential in ``.alloy`` or ``.adp`` format
                           (This is generally all you need to supply)

``elements[N]``            array of N element abbreviations

``embedded_energy[N]``     arrays of embedded energy functions

``electron_density[N]``    arrays of electron density functions

``phi[N,N]``               arrays of pair potential functions

``d_embedded_energy[N]``   arrays of derivative embedded energy functions

``d_electron_density[N]``  arrays of derivative electron density functions

``d_phi[N,N]``             arrays of derivative pair potentials functions

``d[N,N], q[N,N]``         ADP dipole and quadrupole function

``d_d[N,N], d_q[N,N]``     ADP dipole and quadrupole derivative functions

``skin``                   skin distance passed to NeighborList(). If no atom
                           has moved more than the skin-distance since the last
                           call to the ``update()`` method then the neighbor
                           list can be reused. Defaults to 1.0.

``form``                   the form of the potential ``alloy`` or ``adp``. This
                           will be determined from the file suffix or must be
                           set if using equations

=========================  ====================================================


Additional parameters for writing potential files
=================================================

The following parameters are only required for writing a potential in
``.alloy`` or ``.adp`` format file.

=========================  ====================================================
Keyword                    Description
=========================  ====================================================
``header``                 Three line text header. Default is standard message.

``Z[N]``                   Array of atomic number of each element

``mass[N]``                Atomic mass of each element

``a[N]``                   Array of lattice parameters for each element

``lattice[N]``             Lattice type

``nrho``                   No. of rho samples along embedded energy curve

``drho``                   Increment for sampling density