Python neb.NEB类代码示例

    def write(self):
        method = self.method()
        images = self.getOriginalImages()
        new_images = []
        for i in range(len(images)-1):
            initial, final = images[i:i+2]
            diff = final.positions - initial.positions
            diff = vectorLength(diff, 1)
            n_images = int(np.ceil(diff.max()/0.3)) - 1
            print 'n_images=%d' % n_images
            t_images = [initial]
            for j in range(n_images):
                temp = t_images[0].copy()
                t_images.append(temp)
            t_images.append(final)
            neb = NEB(t_images)
            neb.interpolate()
            new_images.extend(neb.images[:-1])

        new_images.append(images[-1])
        neb = NEB(new_images)
        kpts = self.baseTask().parameters().kpts()
        self.method().writeNEB(
                runfile=self.runFile(),
                neb=neb,
                kpts=kpts,
                )

def run_neb_calculation(cpu):
    images = [PickleTrajectory('H.traj')[-1]]
    for i in range(nimages):
        images.append(images[0].copy())
    images[-1].positions[6, 1] = 2 - images[0].positions[6, 1]
    neb = NEB(images, parallel=True, world=cpu)
    neb.interpolate()

    images[cpu.rank + 1].set_calculator(Calculator())

    dyn = BFGS(neb)
    dyn.run(fmax=fmax)

    if cpu.rank == 1:
        results.append(images[2].get_potential_energy())

def get_atoms():
    # 2x2-Al(001) surface with 3 layers and an
    # Au atom adsorbed in a hollow site:
    slab = fcc100('Al', size=(2, 2, 3))
    add_adsorbate(slab, 'Au', 1.7, 'hollow')
    slab.center(axis=2, vacuum=4.0)

    # Fix second and third layers:
    mask = [atom.tag > 1 for atom in slab]
    slab.set_constraint(FixAtoms(mask=mask))

    # Use EMT potential:
    slab.set_calculator(EMT())

    # Initial state:
    qn = QuasiNewton(slab, logfile=None)
    qn.run(fmax=0.05)
    initial = slab.copy()

    # Final state:
    slab[-1].x += slab.get_cell()[0, 0] / 2
    qn = QuasiNewton(slab, logfile=None)
    qn.run(fmax=0.05)
    final = slab.copy()

    # Setup a NEB calculation
    constraint = FixAtoms(mask=[atom.tag > 1 for atom in initial])

    images = [initial]
    for i in range(3):
        image = initial.copy()
        image.set_constraint(constraint)
        images.append(image)

    images.append(final)

    neb = NEB(images, parallel=mpi.parallel)
    neb.interpolate()

    def set_calculator(calc):
        i = 0
        for image in neb.images[1:-1]:
            if not mpi.parallel or mpi.rank // (mpi.size // 3) == i:
                image.set_calculator(calc)
            i += 1
    neb.set_calculator = set_calculator

    return neb

def interp_by_neb(initial, final, n_images, interp='idpp',
                  calculator=None, constraint=None, fmax=0.5, steps=10,
                  **kwargs):
    """Interpolate between the initial and final points by NEB method.

    :param initial: The initial image.
    :param final: The final image.
    :param int n_images: The number of image between the initial and final
        images.
    :param string interp: The interpolation method.

    :param Callable calculator: The callable to generate calculators for the
        images.  It can be set to None if no real optimization is intended.
    :param constraint: The constraint for the images.
    :param fmax: The maximal force to stop the optimization.
    :param steps: The number of optimization steps allowed.

    :param kwargs:  All other keyword arguments are forwarded to the
        initializer of the ASE NEB class.

    :return: The list of images between the points.
    """

    # To circumvent the error from interpolation when we have no middle images.
    if n_images == 0:
        return [initial, final]

    images = [initial]
    for _ in range(n_images):
        image = initial.copy()
        images.append(image)
        if calculator is not None:
            image.set_calculator(calculator())
        if constraint is not None:
            image.set_constraint(constraint)
        continue
    images.append(final)

    neb = NEB(images, **kwargs)
    neb.interpolate(method=interp)

    if calculator is not None:
        dyn = MDMin(neb)
        dyn.run(fmax=fmax, steps=steps)

    return neb.images

def nebmake(directory, start, final, images, tolerance=0, ci=False, poscar_override=[], linear=False, write=True):

    if type(start) == str:
        start_POSCAR = os.path.join(start, 'CONTCAR') if os.path.exists(os.path.join(start, 'CONTCAR')) and os.path.getsize(os.path.join(start, 'CONTCAR')) > 0 else os.path.join(start, 'POSCAR')
        final_POSCAR = os.path.join(final, 'CONTCAR') if os.path.exists(os.path.join(final, 'CONTCAR')) and os.path.getsize(os.path.join(final, 'CONTCAR')) > 0 else os.path.join(final, 'POSCAR')
        p1 = Poscar.from_file(start_POSCAR)
        p2 = Poscar.from_file(final_POSCAR)
        s1 = p1.structure
        s2 = p2.structure
    else:
        s1 = start
        s2 = final
    # s1.sort()
    # s2.sort()
    atoms = []
    if poscar_override:
        for i in range(int(len(poscar_override)/2)):
            atoms.append( (poscar_override[i*2], poscar_override[i*2+1]) )
        (s1, s2) = reorganize_structures(s1, s2, atoms=atoms, autosort_tol=tolerance)
        tolerance=0
    try:
        structures = s1.interpolate(s2, images, autosort_tol=tolerance)
    except Exception as e:
        a=input('Failed.  Type y to sort --> ')
        if a=='y':
            s1.sort()
            s2.sort()
        else:
            raise e
        structures = s1.interpolate(s2, images, autosort_tol=tolerance)


    if not linear:
        from pymatgen.io.ase import AseAtomsAdaptor
        from ase.neb import NEB
        structures_ase = [ AseAtomsAdaptor.get_atoms(struc) for struc in structures ]
        neb = NEB(structures_ase)
        neb.interpolate('idpp') # type: NEB
        structures = [ AseAtomsAdaptor.get_structure(atoms) for atoms in neb.images ]

    if write:
        start_OUTCAR = os.path.join(start, 'OUTCAR')
        final_OUTCAR = os.path.join(final, 'OUTCAR')
        incar = Incar.from_file(os.path.join(start, 'INCAR'))
        kpoints = Kpoints.from_file(os.path.join(start, 'KPOINTS'))
        potcar = Potcar.from_file(os.path.join(start, 'POTCAR'))
        incar['ICHAIN'] = 0
        incar['IMAGES'] = images-1
        incar['LCLIMB'] = ci

        for i, s in enumerate(structures):
            folder = os.path.join(directory, str(i).zfill(2))
            os.mkdir(folder)
            Poscar(s, selective_dynamics=p1.selective_dynamics).write_file(os.path.join(folder, 'POSCAR'))
            if i == 0:
                shutil.copy(start_OUTCAR, os.path.join(folder, 'OUTCAR'))
            if i == images:
                shutil.copy(final_OUTCAR, os.path.join(folder, 'OUTCAR'))
            i += 1

        incar.write_file(os.path.join(directory, 'INCAR'))
        kpoints.write_file(os.path.join(directory, 'KPOINTS'))
        potcar.write_file(os.path.join(directory, 'POTCAR'))
    return structures

    def __init__(self, images, control, k=1.0, climb=False, parallel=False, minmodes=None, decouple_modes=False):
        self.control = control
        NEB.__init__(self, images, k, climb, parallel)

        self.spring_force = 'full'

        # Set up MinModeAtoms objects for each image and make individual logfiles for each
        # NB: Shouldn't there be a ERM_Control class that takes care of this crap?
        self.images = []
        for i in range(self.nimages):
            min_control = control.copy()

            i_num = ('%0' + str(len(str(self.nimages))) + 'i') % i
            d_logfile_old = self.control.get_logfile()
            m_logfile_old = self.control.get_eigenmode_logfile()
            if d_logfile_old not in ['-', None]:
                if type(d_logfile_old) == str:
                    d_logfile_old = d_logfile_old.split('.')
                else:
                    d_logfile_old = d_logfile_old.name.split('.')
                d_logfile_old.insert(-1, i_num)
                d_logfile_new = '-'.join(['.'.join(d_logfile_old[:-2]), '.'.join(d_logfile_old[-2:])])
            else:
                d_logfile_new = d_logfile_old
            if m_logfile_old not in ['-', None]:
                if type(m_logfile_old) == str:
                    m_logfile_old = m_logfile_old.split('.')
                else:
                    m_logfile_old = m_logfile_old.name.split('.')
                m_logfile_old.insert(-1, i_num)
                m_logfile_new = '-'.join(['.'.join(m_logfile_old[:-2]), '.'.join(m_logfile_old[-2:])])
            else:
                m_logfile_new = m_logfile_old

            if i in [0, self.nimages - 1]:
                write_rank = 0
            else:
                write_rank = (i - 1) * size // (self.nimages - 2)

            min_control.set_write_rank(write_rank)
            min_control.initialize_logfiles(logfile = d_logfile_new, eigenmode_logfile = m_logfile_new)
            if minmodes is None:
                minmode = None
            else:
                minmodes = np.array(minmodes)
                if minmodes.shape == (self.nimages, self.natoms, 3):
                    # Assume one minmode for each image
                    raise NotImplementedError()
                elif minmodes.shape == (2, self.natoms, 3):
                    # Assume end images minmodes and interpolate
                    raise NotImplementedError()
                elif minmodes.shape == (self.natoms, 3):
                    minmode = [minmodes.copy()]
                else:
                    raise ValueError('ERM did not understand the minmodes given to it.')

            image = MinModeAtoms(images[i], min_control, eigenmodes = minmode)
            self.images.append(image)

        self.forces['dimer'] = np.zeros((self.nimages, self.natoms, 3))

        # Populate the tangents
        for i in range(1, self.nimages - 1):
            p_m = self.images[i - 1].get_positions()
            p_p = self.images[i + 1].get_positions()
            t = (p_p - p_m) / 2.0
            if 0.0 in t:
                # Assume a linear interpolation
                # HACK/BUG: Currently the last or first "free" image will yield p[-1] - p[0]
                t = self.images[-1].get_positions() - self.images[0].get_positions()
                t /= (self.nimages - 1.0)
            self.tangents[i] = t
        self.tangents[0] = t
        self.tangents[-1] = -t

        # Save user variables
        self.decouple_modes = decouple_modes # Release the orthogonality constraint of the minmode and tanget.

        # Development stuff
        self.plot_devplot = False
        self.plot_subplot = False
        self.plot_animate = 0
        self.plot_x = None
        self.plot_y = None
        self.plot_e = None
        self.xrange = None
        self.yrange = None

import os

from ase.io import read
from ase.neb import NEB
from ase.calculators.turbomole import Turbomole
from ase.optimize import BFGS

initial = read('initial.coord')
final = read('final.coord')
os.system('rm -f coord; cp initial.coord coord')

# Make a band consisting of 5 configs:
configs = [initial]
configs += [initial.copy() for i in range(3)]
configs += [final]

band = NEB(configs, climb=True)
# Interpolate linearly the positions of the not-endpoint-configs:
band.interpolate()

#Set calculators
for config in configs:
    config.set_calculator(Turbomole())

# Optimize the Path:
relax = BFGS(band, trajectory='neb.traj')
relax.run(fmax=0.05)

    initial.set_calculator(LennardJones())

    images = [initial]

    # Set calculator
    for i in range(nimages):
        image = initial.copy()
        image.set_calculator(LennardJones())
        images.append(image)

    images.append(final)

    # Define the NEB and make a linear interpolation
    # with removing translational
    # and rotational degrees of freedom
    neb = NEB(images,
              remove_rotation_and_translation=remove_rotation_and_translation)
    neb.interpolate()
    # Test used these old defaults which are not optimial, but work
    # in this particular system
    neb.idpp_interpolate(fmax=0.1, optimizer=BFGS)

    qn = FIRE(neb, dt=0.005, maxmove=0.05, dtmax=0.1)
    qn.run(steps=20)

    # Switch to CI-NEB, still removing the external degrees of freedom
    # Also spesify the linearly varying spring constants
    neb = NEB(images, climb=True,
              remove_rotation_and_translation=remove_rotation_and_translation)
    qn = FIRE(neb, dt=0.005, maxmove=0.05, dtmax=0.1)
    qn.run(fmax=fmax)

images = [slab]
for i in range(6):
    image = slab.copy()
    # Set constraints and calculator:
    image.set_constraint(constraint)
    image.calc = EMT()
    images.append(image)

# Displace last image:
image[-2].position = image[-1].position
image[-1].x = d
image[-1].y = d / sqrt(3)

dyn = QuasiNewton(images[-1])
dyn.run(fmax=0.05)
neb = NEB(images, climb=not True)

# Interpolate positions between initial and final states:
neb.interpolate(method='idpp')

for image in images:
    print(image.positions[-1], image.get_potential_energy())

dyn = BFGS(neb, maxstep=0.04, trajectory='mep.traj')
dyn.run(fmax=0.05)

for image in images:
    print(image.positions[-1], image.get_potential_energy())

if locals().get('display'):
    import os

from ase.build import molecule
from ase.neb import NEB

initial = molecule('C2H6')
final = initial.copy()
final.positions[2:5] = initial.positions[[3, 4, 2]]

images = [initial]
for i in range(5):
    images.append(initial.copy())
images.append(final)

neb = NEB(images)
d0 = images[3].get_distance(2, 3)
neb.interpolate()
d1 = images[3].get_distance(2, 3)
neb.idpp_interpolate(fmax=0.005)
d2 = images[3].get_distance(2, 3)
print(d0, d1, d2)
assert abs(d2 - 1.74) < 0.01

#optimise the initial state
# Atom above step
slab[-1].position = (x3,y2+1,z2+3.5)
final = slab.copy()
final.set_calculator(EMT())
relax = QuasiNewton(final)
relax.run(fmax=0.05)
view(final)

#create a list of images for interpolation
images = [initial]
for i in range(nimages):
    images.append(initial.copy())

for image in images:
    image.set_calculator(EMT())

images.append(final)
view(images)

#carry out idpp interpolation
neb = NEB(images)
neb.interpolate('idpp')

#Run NEB calculation
qn = QuasiNewton(neb, trajectory='N_diffusion.traj',  logfile='N_diffusion.log')
qn.run(fmax=0.05)

# Approximate height of Ag atom on Cu(100) surfece:
h0 = 2.2373
initial += Atom('Pt', (10.96, 11.074, h0))
initial += Atom('Pt', (13.7, 11.074, h0))
initial += Atom('Pt', (9.59, 8.701, h0))
initial += Atom('Pt', (12.33, 8.701, h0))
initial += Atom('Pt', (15.07, 8.701, h0))
initial += Atom('Pt', (10.96, 6.328, h0))
initial += Atom('Pt', (13.7, 6.328, h0))

if 0:
    view(initial)

# Make band:
images = [initial.copy() for i in range(7)]
neb = NEB(images)

# Set constraints and calculator:
indices = np.compress(initial.positions[:, 2] < -5.0, range(len(initial)))
constraint = FixAtoms(indices)
for image in images:
    image.set_calculator(ASAP())
    image.constraints.append(constraint)

# Displace last image:
for i in xrange(1,8,1):
    images[-1].positions[-i] += (d/2, -h1/3, 0)

write('initial.traj', images[0])
# Relax height of Ag atom for initial and final states:
for image in [images[0], images[-1]]:

    def run(self):
        '''Run the AutoNEB optimization algorithm.'''
        n_cur = self.__initialize__()
        while len(self.all_images) < self.n_simul + 2:
            if isinstance(self.k, (float, int)):
                self.k = [self.k] * (len(self.all_images) - 1)
            if self.world.rank == 0:
                print('Now adding images for initial run')
            # Insert a new image where the distance between two images is
            # the largest
            spring_lengths = []
            for j in range(n_cur - 1):
                spring_vec = self.all_images[j + 1].get_positions() - \
                    self.all_images[j].get_positions()
                spring_lengths.append(np.linalg.norm(spring_vec))
            jmax = np.argmax(spring_lengths)

            if self.world.rank == 0:
                print('Max length between images is at ', jmax)

            # The interpolation used to make initial guesses
            # If only start and end images supplied make all img at ones
            if len(self.all_images) == 2:
                n_between = self.n_simul
            else:
                n_between = 1

            toInterpolate = [self.all_images[jmax]]
            for i in range(n_between):
                toInterpolate += [toInterpolate[0].copy()]
            toInterpolate += [self.all_images[jmax + 1]]

            neb = NEB(toInterpolate)
            neb.interpolate(method=self.interpolate_method)

            tmp = self.all_images[:jmax + 1]
            tmp += toInterpolate[1:-1]
            tmp.extend(self.all_images[jmax + 1:])

            self.all_images = tmp

            # Expect springs to be in equilibrium
            k_tmp = self.k[:jmax]
            k_tmp += [self.k[jmax] * (n_between + 1)] * (n_between + 1)
            k_tmp.extend(self.k[jmax + 1:])
            self.k = k_tmp

            # Run the NEB calculation with the new image included
            n_cur += n_between

        # Determine if any images do not have a valid energy yet
        energies = self.get_energies()

        n_non_valid_energies = len([e for e in energies if e != e])

        if self.world.rank == 0:
            print('Start of evaluation of the initial images')

        while n_non_valid_energies != 0:
            if isinstance(self.k, (float, int)):
                self.k = [self.k] * (len(self.all_images) - 1)

            # First do one run since some energie are non-determined
            to_run, climb_safe = self.which_images_to_run_on()
            self.execute_one_neb(n_cur, to_run, climb=False)

            energies = self.get_energies()
            n_non_valid_energies = len([e for e in energies if e != e])

        if self.world.rank == 0:
            print('Finished initialisation phase.')

        # Then add one image at a time until we have n_max images
        while n_cur < self.n_max:
            if isinstance(self.k, (float, int)):
                self.k = [self.k] * (len(self.all_images) - 1)
            # Insert a new image where the distance between two images
            # is the largest OR where a higher energy reselution is needed
            if self.world.rank == 0:
                print('****Now adding another image until n_max is reached',
                      '({0}/{1})****'.format(n_cur, self.n_max))
            spring_lengths = []
            for j in range(n_cur - 1):
                spring_vec = self.all_images[j + 1].get_positions() - \
                    self.all_images[j].get_positions()
                spring_lengths.append(np.linalg.norm(spring_vec))

            total_vec = self.all_images[0].get_positions() - \
                self.all_images[-1].get_positions()
            tl = np.linalg.norm(total_vec)

            fR = max(spring_lengths) / tl

            e = self.get_energies()
            ed = []
            emin = min(e)
            enorm = max(e) - emin
            for j in range(n_cur - 1):
                delta_E = (e[j + 1] - e[j]) * (e[j + 1] + e[j] - 2 *
                                               emin) / 2 / enorm
#.........这里部分代码省略.........

    def execute_one_neb(self, n_cur, to_run, climb=False, many_steps=False):
        '''Internal method which executes one NEB optimization.'''
        self.iteration += 1
        # First we copy around all the images we are not using in this
        # neb (for reproducability purposes)
        if self.world.rank == 0:
            for i in range(n_cur):
                if i not in to_run[1: -1]:
                    filename = '%s%03d.traj' % (self.prefix, i)
                    self.all_images[i].write(filename)
                    filename_ref = self.iter_folder + \
                        '/%s%03diter%03d.traj' % (self.prefix, i,
                                                  self.iteration)
                    if os.path.isfile(filename):
                        shutil.copy2(filename, filename_ref)
        if self.world.rank == 0:
            print('Now starting iteration %d on ' % self.iteration, to_run)
        # Attach calculators to all the images we will include in the NEB
        self.attach_calculators([self.all_images[i] for i in to_run[1: -1]])
        neb = NEB([self.all_images[i] for i in to_run],
                  k=[self.k[i] for i in to_run[0:-1]],
                  method=self.method,
                  parallel=self.parallel,
                  remove_rotation_and_translation=self
                  .remove_rotation_and_translation,
                  climb=climb)

        # Do the actual NEB calculation
        qn = self.optimizer(neb,
                            logfile=self.iter_folder +
                            '/%s_log_iter%03d.log' % (self.prefix,
                                                      self.iteration))

        # Find the ranks which are masters for each their calculation
        if self.parallel:
            nneb = to_run[0]
            nim = len(to_run) - 2
            n = self.world.size // nim      # number of cpu's per image
            j = 1 + self.world.rank // n    # my image number
            assert nim * n == self.world.size
            traj = Trajectory('%s%03d.traj' % (self.prefix, j + nneb), 'w',
                              self.all_images[j + nneb],
                              master=(self.world.rank % n == 0))
            filename_ref = self.iter_folder + \
                '/%s%03diter%03d.traj' % (self.prefix,
                                          j + nneb, self.iteration)
            trajhist = Trajectory(filename_ref, 'w',
                                  self.all_images[j + nneb],
                                  master=(self.world.rank % n == 0))
            qn.attach(traj)
            qn.attach(trajhist)
        else:
            num = 1
            for i, j in enumerate(to_run[1: -1]):
                filename_ref = self.iter_folder + \
                    '/%s%03diter%03d.traj' % (self.prefix, j, self.iteration)
                trajhist = Trajectory(filename_ref, 'w', self.all_images[j])
                qn.attach(seriel_writer(trajhist, i, num).write)

                traj = Trajectory('%s%03d.traj' % (self.prefix, j), 'w',
                                  self.all_images[j])
                qn.attach(seriel_writer(traj, i, num).write)
                num += 1

        if isinstance(self.maxsteps, (list, tuple)) and many_steps:
            steps = self.maxsteps[1]
        elif isinstance(self.maxsteps, (list, tuple)) and not many_steps:
            steps = self.maxsteps[0]
        else:
            steps = self.maxsteps

        if isinstance(self.fmax, (list, tuple)) and many_steps:
            fmax = self.fmax[1]
        elif isinstance(self.fmax, (list, tuple)) and not many_steps:
            fmax = self.fmax[0]
        else:
            fmax = self.fmax
        qn.run(fmax=fmax, steps=steps)

        # Remove the calculators and replace them with single
        # point calculators and update all the nodes for
        # preperration for next iteration
        neb.distribute = types.MethodType(store_E_and_F_in_spc, neb)
        neb.distribute()

del atoms[5]
print(atoms)
assert len(atoms.constraints[0].index) == 5




fmax = 0.05
nimages = 3

print([a.get_potential_energy() for a in Trajectory('H.traj')])
images = [Trajectory('H.traj')[-1]]
for i in range(nimages):
    images.append(images[0].copy())
images[-1].positions[6, 1] = 2 - images[0].positions[6, 1]
neb = NEB(images)
neb.interpolate()
if 0:  # verify that initial images make sense
    from ase.visualize import view
    view(neb.images)

for image in images:
    image.set_calculator(MorsePotential())

dyn = BFGS(neb, trajectory='mep.traj')  # , logfile='mep.log')

dyn.run(fmax=fmax)

for a in neb.images:
    print(a.positions[-1], a.get_potential_energy())