Python cubic.FaceCenteredCubic类代码示例

def equilShape(element,params,size=(10,10,10),distance=25.0,corrections=0,structure='fcc'):
    """
    this is to use the ratio of energies to calculate the equilibrium crystal shape, cycle through a bunch of (h,k,l) indices
    """
    slab = FaceCenteredCubic(element,directions=([[1,0,0],[0,1,0],[0,0,1]]),size=(10,10,10))    
    energy100 = fitFunction([1,0,0],params,corrections,structure)
    h100 = distance
    orig_positions = slab.get_positions()
    kept_positions = list(orig_positions)
    center = slab.get_center_of_mass()
    for h in range(-12,12):
        for k in range(0,9):
            for l in range(0,9):
                nvector=list([h,k,l]/numpy.sqrt(h**2+k**2+l**2))
                energyhkl = fitFunction(nvector,params,corrections,structure)
                distancehkl = energyhkl/energy100*h100
                for i in range(0,len(kept_positions)):
                    list_to_pop = []
                    if numpy.dot(kept_positions[i],nvector) > distancehkl:
                        list_to_pop.append(i)
                for i in list_to_pop:
                    kept_positions.pop(i)
    
    # set up new slab with new positions
    number_of_atoms = len(kept_positions)
    elstring = str(number_of_atoms)+'Pt' 
    new_slab = Atoms(elstring,positions=kept_positions)

    return new_slab

def relax(input_atoms, ref_db):
    atoms_string = input_atoms.get_chemical_symbols()

    # Open connection to the database with reference data
    db = connect(ref_db)

    # Load our model structure which is just FCC
    atoms = FaceCenteredCubic('X', latticeconstant=1.)
    atoms.set_chemical_symbols(atoms_string)

    # Compute the average lattice constant of the metals in this individual
    # and the sum of energies of the constituent metals in the fcc lattice
    # we will need this for calculating the heat of formation
    a = 0
    ei = 0
    for m in set(atoms_string):
        dct = db.get(metal=m)
        count = atoms_string.count(m)
        a += count * dct.latticeconstant
        ei += count * dct.energy_per_atom
    a /= len(atoms_string)
    atoms.set_cell([a, a, a], scale_atoms=True)

    # Since calculations are extremely fast with EMT we can also do a volume
    # relaxation
    atoms.set_calculator(EMT())
    eps = 0.05
    volumes = (a * np.linspace(1 - eps, 1 + eps, 9))**3
    energies = []
    for v in volumes:
        atoms.set_cell([v**(1. / 3)] * 3, scale_atoms=True)
        energies.append(atoms.get_potential_energy())

    eos = EquationOfState(volumes, energies)
    v1, ef, B = eos.fit()
    latticeconstant = v1**(1. / 3)

    # Calculate the heat of formation by subtracting ef with ei
    hof = (ef - ei) / len(atoms)

    # Place the calculated parameters in the info dictionary of the
    # input_atoms object
    input_atoms.info['key_value_pairs']['hof'] = hof
    
    # Raw score must always be set
    # Use one of the following two; they are equivalent
    input_atoms.info['key_value_pairs']['raw_score'] = -hof
    # set_raw_score(input_atoms, -hof)
    
    input_atoms.info['key_value_pairs']['latticeconstant'] = latticeconstant

    # Setting the atoms_string directly for easier analysis
    atoms_string = ''.join(input_atoms.get_chemical_symbols())
    input_atoms.info['key_value_pairs']['atoms_string'] = atoms_string

 def test_Grochola(self):
     a = FaceCenteredCubic('Au', size=[2,2,2])
     calc = EAM('Au-Grochola-JCP05.eam.alloy')
     a.set_calculator(calc)
     FIRE(StrainFilter(a, mask=[1,1,1,0,0,0]), logfile=None).run(fmax=0.001)
     a0 = a.cell.diagonal().mean()/2
     self.assertTrue(abs(a0-4.0701)<2e-5)
     self.assertTrue(abs(a.get_potential_energy()/len(a)+3.924)<0.0003)
     C, C_err = fit_elastic_constants(a, symmetry='cubic', verbose=False)
     C11, C12, C44 = Voigt_6x6_to_cubic(C)
     self.assertTrue(abs((C11-C12)/GPa-32.07)<0.7)
     self.assertTrue(abs(C44/GPa-45.94)<0.5)

def MakeCu(T=300, size=(29,29,30)):
    print "Preparing", T, "K Copper system."
    atoms = FaceCenteredCubic(directions=[[1,0,0],[0,1,0],[0,0,1]],
                              symbol='Cu', size=size)
    atoms.set_calculator(EMT())
    MaxwellBoltzmannDistribution(atoms, 2*T * units.kB)
    #dyn = VelocityVerlet(atoms, 5*units.fs)
    dyn = Langevin(atoms, 5*units.fs, T*units.kB, 0.05)
    dyn.run(50)
    print "Done.  Temperature =", temperature(atoms), \
          "K.  Number of atoms: ", len(atoms)
    return atoms

    def test_stress(self):
        a = FaceCenteredCubic('Au', size=[2,2,2])
        calc = EAM('Au-Grochola-JCP05.eam.alloy')
        a.set_calculator(calc)
        self.assertArrayAlmostEqual(a.get_stress(), calc.calculate_numerical_stress(a), tol=self.tol)

        sx, sy, sz = a.cell.diagonal()
        a.set_cell([sx, 0.9*sy, 1.2*sz], scale_atoms=True)
        self.assertArrayAlmostEqual(a.get_stress(), calc.calculate_numerical_stress(a), tol=self.tol)

        a.set_cell([[sx, 0.1*sx, 0], [0, 0.9*sy, 0], [0, -0.1*sy, 1.2*sz]], scale_atoms=True)
        self.assertArrayAlmostEqual(a.get_stress(), calc.calculate_numerical_stress(a), tol=self.tol)

    def makePerturbedLattice(self, shape=(2,2,2)):
        """
        to make perturbed bulk structure lattice positions
        """
        a = self.getFCCLattice()
        slab = FaceCenteredCubic(size = shape, symbol=self.element, latticeconstant=a)
        positions=slab.get_positions()

        perturbations = numpy.random.standard_normal(scipy.shape(positions))*a*0.05

        positions += perturbations

        return positions

def fcc211(symbol, size, a=None, vacuum=None, orthogonal=True):
    """FCC(211) surface.

    Does not currently support special adsorption sites.

    Currently only implemented for *orthogonal=True* with size specified
    as (i, j, k), where i, j, and k are number of atoms in each direction.
    i must be divisible by 3 to accommodate the step width.
    """
    if not orthogonal:
        raise NotImplementedError('Only implemented for orthogonal '
                                  'unit cells.')
    if size[0] % 3 != 0:
        raise NotImplementedError('First dimension of size must be '
                                  'divisible by 3.')
    atoms = FaceCenteredCubic(symbol,
                              directions=[[1, -1, -1],
                                          [0, 2, -2],
                                          [2, 1, 1]],
                              miller=(None, None, (2, 1, 1)),
                              latticeconstant=a,
                              size=(1, 1, 1),
                              pbc=True)
    z = (size[2] + 1) // 2
    atoms = atoms.repeat((size[0] // 3, size[1], z))
    if size[2] % 2:  # Odd: remove bottom layer and shrink cell.
        remove_list = [atom.index for atom in atoms
                       if atom.z < atoms[1].z]
        del atoms[remove_list]
        dz = atoms[0].z
        atoms.translate((0., 0., -dz))
        atoms.cell[2][2] -= dz

    atoms.cell[2] = 0.0
    atoms.pbc[1] = False
    if vacuum:
        atoms.center(vacuum, axis=2)

    # Renumber systematically from top down.
    orders = [(atom.index, round(atom.x, 3), round(atom.y, 3),
               -round(atom.z, 3), atom.index) for atom in atoms]
    orders.sort(key=itemgetter(3, 1, 2))
    newatoms = atoms.copy()
    for index, order in enumerate(orders):
        newatoms[index].position = atoms[order[0]].position.copy()
    return newatoms

def MakeAtoms(elem1, elem2=None):
    if elem2 is None:
        elem2 = elem1
    a1 = reference_states[elem1]['a']
    a2 = reference_states[elem2]['a']
    a0 = (0.5 * a1**3 + 0.5 * a2**3)**(1.0/3.0) * 1.03
    if ismaster:
        # 50*50*50 would be big enough, but some vacancies are nice.
        print "Z1 = %i,  Z2 = %i,  a0 = %.5f" % (elem1, elem2, a0)
        atoms = FaceCenteredCubic(symbol='Cu', size=(51,51,51))
        nremove = len(atoms) - 500000
        assert nremove > 0
        remove = np.random.choice(len(atoms), nremove, replace=False)
        del atoms[remove]
        if isparallel:
            atoms = atoms.repeat(cpuLayout)
        if elem1 != elem2:
            z = atoms.get_atomic_numbers()
            z[np.random.choice(len(atoms), len(atoms)/2, replace=False)] = elem2
            atoms.set_atomic_numbers(z)
    else:
        atoms = None
    if isparallel:
        atoms = MakeParallelAtoms(atoms, cpuLayout)
    MaxwellBoltzmannDistribution(atoms, T * units.kB)
    return atoms

def myfcc(symbol, pbc):
    "Make an fcc lattice with standard lattice constant"
    if ismaster:
        a = FaceCenteredCubic(symbol=symbol, size=(10,10,10), pbc=pbc)    
        dx = 0.01 * np.sin(0.1 * np.arange(len(a) * 3))
        dx.shape = (len(a),3)
        a.set_positions(a.get_positions() + dx)
        a.set_momenta(np.zeros((len(a),3)))
        #view(a)
    else:
        a = None
    if isparallel:
        a = MakeParallelAtoms(a, cpulayout)
    return a

if re.match("^n\d\d\d.dcsc.fysik.dtu.dk$", host):
    print "    This is a d512 node on Niflheim."
    fullhost = "niflheim-d512/%s" % (host.split(".")[0])
    host = "niflheim-d512"
elif re.match("^[stu]\d\d\d.dcsc.fysik.dtu.dk$", host):
    print "    This is an s50 node on Niflheim."
    fullhost = "niflheim-s50/%s" % (host.split(".")[0])
    host = "niflheim-s50"
else:
    fullhost = host
print "Current time is "+when
print ""

print "Preparing system"
initial = FaceCenteredCubic(directions=[[1,0,0],[0,1,0],[0,0,1]],
                            size=(30, 30, 30),
                            symbol="Pt")
ReportTest("Number of atoms", len(initial), 108000, 0)
r = initial.get_positions()
r.flat[:] += 0.14 * sin(arange(3*len(initial)))
initial.set_positions(r)

print "Running self-test."
atoms = Atoms(initial)
atoms.set_calculator(EMT2013(PtY_parameters))
e = atoms.get_potential_energies()
f = atoms.get_forces()
if os.access(selfcheckfilename, os.F_OK):
    olde, oldf = cPickle.load(open(selfcheckfilename))
    de = max(fabs(e - olde))
    df = max(fabs(f.flat[:] - oldf.flat[:]))

#!/usr/bin/env python
from ase.lattice.cubic import FaceCenteredCubic
from ase.calculators.aims import Aims
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np

atoms = FaceCenteredCubic('Pt')
kpts = [2, 3, 5, 7, 9, 11, 15, 19, 25] 
energies = []
ready = True
for k in kpts:
    calc = Aims(label='bulk/pt-kpts-{0}'.format(k),
              xc='pbe',
              spin='none',
              relativistic = 'atomic_zora scalar',
              kpts=(k, k, k),  # specifies the Monkhorst-Pack grid
              sc_accuracy_etot=1e-5,
              sc_accuracy_eev=1e-2,
              sc_accuracy_rho=1e-4,)
    atoms.set_calculator(calc)
    try:
        energies.append(atoms.get_potential_energy())
    except:
        print('aims failed when k = {0}'.format(k))
        ready = False
if not ready:
    import sys; sys.exit()

print '# pt-fcc-latt'

"""Test force calculations alone with the purpose of optimizing threading."""

from asap3 import *
from asap3.Timing import report_timing
from ase.lattice.cubic import FaceCenteredCubic
import time

nsteps = 10

start = time.time()
if len(sys.argv) >= 2 and sys.argv[1] == '-t':
    AsapThreads()   
if len(sys.argv) >= 2 and sys.argv[1] == '-T':
    AsapThreads(6)   

atoms = FaceCenteredCubic(size=(100,100,50), symbol='Cu')
atoms.set_calculator(EMT())
print "Number of atoms:", len(atoms)

d = 0.1
pos = atoms.arrays['positions']  # Nasty!

for i in range(nsteps):
    pos[50][0] += d
    d = -d
    f = atoms.get_forces()
    
report_timing()

print "Wall time elapsed:", time.time() - start

            j = int(jj)
            ReportTest.BoolTest("Atom %d on list %d (forward)" % (j, i),
                                j in fnb[i], silent=True)
            ReportTest.BoolTest("Atom %d on list %d (reverse)" % (i, j),
                                i in fnb[j], silent=True)
        if ReportTest.GetNumberOfErrors() > 10:
            print "*** Too many errors - giving up! ***"
            break
    

print_version(1)

element = "Cu"
latconst = ase.data.reference_states[ase.data.atomic_numbers[element]]['a']

atoms = FaceCenteredCubic(directions=[[1,0,0],[0,1,1],[0,0,1]], size=(9,7,5),
                          symbol=element, debug=0)
atoms.set_calculator(EMT(minimum_image=True))
epot = atoms.get_potential_energy()

nblist = atoms.get_calculator().get_neighborlist()
count = {}
for lst in nblist:
    n = len(lst)
    try:
        count[n] += 1
    except KeyError:
        count[n] = 1
# print "Histogram:"
numbers = count.keys()
numbers.sort()
sum = 0

conv = {'eigenstates' : 1e-4, 'density' : 1e-2, 'energy' : 1e-3}

# output benchmark parameters
if rank == 0:
    print("#"*60)
    print("GPAW benchmark: Copper Sheet")
    print("  dimensions: x=%d, y=%d, z=%d" % (x, y, z))
    print("  grid spacing: h=%f" % h)
    print("  Brillouin-zone sampling: kpts=" + str(kpts))
    print("  MPI task: %d out of %d" % (rank, size))
    print("  using MICs: " + str(use_mic))
    print("#"*60)
    print("")

# setup the system
atoms = FaceCenteredCubic(directions=[[1,-1,0], [1,1,-2], [1,1,1]],
        size=(x,y,z), symbol='Cu', pbc=(1,1,0))
#add_vacuum(atoms, 10.0)
atoms.center(vacuum=6.0, axis=2)
calc = GPAW(h=h, nbands=-20, width=0.2,
            kpts=kpts, xc='PBE',
            maxiter=maxiter,
            txt=txt, eigensolver=RMM_DIIS(niter=2),
            parallel={'sl_auto': True},
            mixer=Mixer(0.1, 5, 100),
           )
atoms.set_calculator(calc)

# execute the run
try:
    atoms.get_potential_energy()
except ConvergenceError:

from __future__ import print_function, division
from ase.lattice.cubic import FaceCenteredCubic
from ase.lattice.hexagonal import HexagonalClosedPacked
from ase.test import must_raise

with must_raise(ValueError):
    # The Miller indices of the surfaces are linearly dependent
    atoms = FaceCenteredCubic(symbol='Cu',
                              miller=[[1, 1, 0], [1, 1, 0], [0, 0, 1]])

# This one should be OK:
atoms = FaceCenteredCubic(symbol='Cu',
                          miller=[[1, 1, 0], [0, 1, 0], [0, 0, 1]])
print(atoms.get_cell())


with must_raise(ValueError):
    # The directions spanning the unit cell are linearly dependent
    atoms = FaceCenteredCubic(symbol='Cu',
                              directions=[[1, 1, 0], [1, 1, 0], [0, 0, 1]])

with must_raise(ValueError):
    # The directions spanning the unit cell are linearly dependent
    atoms = FaceCenteredCubic(symbol='Cu',
                              directions=[[1, 1, 0], [1, 0, 0], [0, 1, 0]])

# This one should be OK:
atoms = FaceCenteredCubic(symbol='Cu',
                          directions=[[1, 1, 0], [0, 1, 0], [0, 0, 1]])
print(atoms.get_cell())