Python optimize.QuasiNewton类代码示例

def aual100(site, height, calc=None):

    slab = fcc100('Al', size=(2, 2, 1))
    add_adsorbate(slab, 'Au', height, site)
    slab.center(axis=2, vacuum=3.0)
    mask = [atom.symbol == 'Al' for atom in slab]
    fixlayer = FixAtoms(mask=mask)
    slab.set_constraint(fixlayer)

    if calc is None:
        calc = GPAW(h=0.25, kpts=(2, 2, 1), xc='PBE', txt=site + '.txt')

    slab.set_calculator(calc)
    qn = QuasiNewton(slab, trajectory=site + '.traj')
    qn.run(fmax=0.05)

    if isinstance(calc, GPAW):
        calc.write(site + '.gpw')

    return slab.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 evaluate(self, imember):
        entry = self.population.get_entry(imember)
        pcm_structure = pychemia.Structure.from_dict(entry['structure'])

        ase_structure = pychemia.external.ase.pychemia2ase(pcm_structure)
        ase_structure.set_calculator(LennardJones())

        dyn = QuasiNewton(ase_structure)
        dyn.run()

        ase_structure.set_constraint(FixAtoms(mask=[True for atom in ase_structure]))
        ucf = UnitCellFilter(ase_structure)
        qn = QuasiNewton(ucf)
        qn.run()

        new_structure = pychemia.external.ase.ase2pychemia(ase_structure)
        energy = ase_structure.get_potential_energy()
        forces = ase_structure.get_forces()
        stress = ase_structure.get_stress()
        new_properties = {'energy': float(energy), 'forces': generic_serializer(forces),
                          'stress': generic_serializer(stress)}

        self.population.db.update(imember, structure=new_structure, properties=new_properties)

from ase.optimize import QuasiNewton


def set_work_cell(molecule,h,vacuum):
    molecule.center(vacuum=vacuum)
    cell = molecule.get_cell()
    for i in [0,1,2]:
        cell[i,i] = int(cell[i,i]/4./h)*4.*h
    molecule.set_cell(cell)
    return

h=0.2
vacuum=4.5
atoms = read('init.traj')
#atoms = read('testMolecule1nm.xyz')
#atoms.set_pbc(True)
#atoms.set_cell([10.0,10.0,10.0])
#set_work_cell(atoms,h,vacuum)

## gpaw calculator:
calc = GPAW(h=h, nbands=-20, xc='PBE', txt='relax.txt',
            spinpol=True, charge=+1)
atoms.set_calculator(calc)

relax = QuasiNewton(atoms, logfile='relax.log',trajectory='relax.traj')
relax.run(fmax=0.05)

calc.write('relax.gpw',mode='all')

import numpy as np
from ase import Atoms
from ase.units import fs
from ase.calculators.test import TestPotential
from ase.calculators.emt import EMT
from ase.md import VelocityVerlet
from ase.io import PickleTrajectory, read
from ase.optimize import QuasiNewton

np.seterr(all='raise')
a = Atoms('4X',
          masses=[1, 2, 3, 4],
          positions=[(0, 0, 0),
                     (1, 0, 0),
                     (0, 1, 0),
                     (0.1, 0.2, 0.7)],
          calculator=TestPotential())
print a.get_forces()
md = VelocityVerlet(a, dt=0.5 * fs, logfile='-', loginterval=500)
traj = PickleTrajectory('4N.traj', 'w', a)
md.attach(traj.write, 100)
e0 = a.get_total_energy()
md.run(steps=10000)
del traj
assert abs(read('4N.traj').get_total_energy() - e0) < 0.0001

qn = QuasiNewton(a)
qn.run(0.001)
assert abs(a.get_potential_energy() - 1.0) < 0.000002

    def add_displacement_energy(self, displacement):
        """Add the groundstate energy for a displacements along the
        translational path, and adds it to an_mode['displacement_energies'].

        Args:
            displacement (float): How much to follow translational path.
        """

        # Will otherwise do a groundstate calculation at initial positions
        if displacement:
            if displacement != self.an_mode['transition_path_length']:
                self.atoms.set_positions(
                    self.get_translation_positions(displacement))

                # Do 1D optimization
                fix_environment = FixAtoms(mask=[
                    i not in self.an_mode['indices']
                    for i in range(len(self.atoms))])

                axis_relax = self.an_mode.get('relax_axis')
                if axis_relax:
                    if self.use_forces:
                        warnings.warn(' '.join([
                            "relax along axis and force_consistent",
                            "should only be used with ase releases after",
                            "Jan 2017. See",
                            "https://gitlab.com/ase/ase/merge_requests/354"
                        ]))
                    c = []
                    for i in self.an_mode['indices']:
                        c.append(FixedLine(i, axis_relax))
                    # Fixing everything that is not the vibrating part
                    c.append(fix_environment)
                    self.atoms.set_constraint(c)

                    # Optimization
                    dyn = QuasiNewton(self.atoms, logfile='/dev/null')
                    dyn.run(fmax=self.settings.get('fmax', 0.05))

                    self.atoms.set_constraint(fix_environment)

        if not self.an_mode.get('displacement_energies'):
            self.an_mode['displacement_energies'] = list()

        if self.use_forces:
            e = self.atoms.get_potential_energy(force_consistent=True)

            # For the forces, we need the projection of the forces
            # on the normal mode of the rotation at the current angle
            v_force = self.atoms.get_forces()[
                self.an_mode['indices']].reshape(-1)

            f = float(np.dot(
                v_force, self.an_mode['mode_tangent']))

            if not self.an_mode.get('displacement_forces'):
                self.an_mode['displacement_forces'] = [f]
            else:
                self.an_mode['displacement_forces'].append(f)
        else:
            e = self.atoms.get_potential_energy()

        if self.traj is not None:
            self.traj.write(self.atoms)

        self.an_mode['displacement_energies'].append(e)

        # adding to trajectory:
        if self.traj is not None:
            self.traj.write(self.atoms)

        self.atoms.set_positions(self.groundstate_positions)

        # save to backup file:
        if self.an_filename:
            self.save_to_backup()

"""Test that the filter and trajectories are playing well together."""

from ase.build import molecule
from ase.constraints import Filter
from ase.optimize import QuasiNewton
from ase.calculators.emt import EMT

atoms = molecule('CO2')
atoms.set_calculator(EMT())
filter = Filter(atoms, indices=[1, 2])

opt = QuasiNewton(filter, trajectory='filter-test.traj', logfile='filter-test.log')
opt.run()

# Initial state:
# 2x2-Al(001) surface with 1 layer and an
# Au atom adsorbed in a hollow site:
slab = fcc100('Al', size=(2, 2, 2))
slab.center(axis=2, vacuum=3.0)
add_adsorbate(slab, 'Au', 1.6, 'hollow')

# Make sure the structure is correct:
view(slab)

# Fix the Al atoms:
mask = [atom.symbol == 'Al' for atom in slab]
print(mask)
fixlayer = FixAtoms(mask=mask)
slab.set_constraint(fixlayer)

# Use GPAW:
calc = GPAW(mode=PW(200), kpts=(2, 2, 1), xc='PBE', txt='hollow.txt')
slab.set_calculator(calc)

qn = QuasiNewton(slab, trajectory='hollow.traj')

# Find optimal height.  The stopping criterion is: the force on the
# Au atom should be less than 0.05 eV/Ang
qn.run(fmax=0.05)

calc.write('hollow.gpw')  # Write gpw output after the minimization

print('energy:', slab.get_potential_energy())
print('height:', slab.positions[-1, 2] - slab.positions[0, 2])

#!/usr/bin/env python
#

import os

from ase import Atoms
from ase.calculators.dftb import Dftb
from ase.optimize import QuasiNewton
from ase.io import write, read

from ase.structure import molecule
test = molecule('H2O')
test.set_calculator(Dftb(label='h2o',atoms=test,
run_manyDftb_steps = True,
Driver_='ConjugateGradient',
Driver_MaxForceComponent='1E-4',
Driver_MaxSteps=1000,
Hamiltonian_MaxAngularMomentum_ = '',
Hamiltonian_MaxAngularMomentum_O = '"p"',
Hamiltonian_MaxAngularMomentum_H = '"s"',
))

dyn = QuasiNewton(test, trajectory='test.traj')
dyn.run(fmax=100, steps=0)
test = read('geo_end.gen')
write('test.final.xyz', test)

        cell=[10, 10, 10])

c_basis = """2 nodes 1.00
0 1 S 0.20 P 1 0.20 6.00
5.00
1.00
1 2 S 0.20 P 1 E 0.20 6.00
6.00 5.00
1.00 0.95"""

specie = Specie(symbol='C', basis_set=PAOBasisBlock(c_basis))
calc = Siesta(
    label='ch4',
    basis_set='SZ',
    xc='LYP',
    mesh_cutoff=300 * Ry,
    species=[specie],
    restart='ch4.XV',
    ignore_bad_restart_file=True,
    fdf_arguments={'DM.Tolerance': 1E-5,
                   'DM.MixingWeight': 0.15,
                   'DM.NumberPulay': 3,
                   'MaxSCFIterations': 200,
                   'ElectronicTemperature': 0.02585 * eV,  # 300 K
                   'SaveElectrostaticPotential': True})

bud.set_calculator(calc)
dyn = QuasiNewton(bud, trajectory=traj)
dyn.run(fmax=0.02)
e = bud.get_potential_energy()

#		OP_f = [OP[i:i+2] for i in xrange(0,len(OP),2)]
#		OP_f_text = np.savetxt('out.tmp', OP_f)
#		OP_f_infile = open('out.tmp','r')
#		OP_f_text = OP_f_infile.read()
#		text_file.write(OP_f_text)
#		text_file.close()
#		exit()

	x_dist = a.get_distance(a0 = 0, a1 = 1)
	OP.append(x_dist)	
	epot = a.get_potential_energy() / len(a)
	ekin = a.get_kinetic_energy() / len(a)
	ETOT = epot+ekin
	OP.append(ETOT)

qn = QuasiNewton(atoms, trajectory='qn.traj')
qn.attach(print_distance)
qn.run(fmax=0.001)
qn.attach(print_distance)
write('qn.final.xyz', atoms)

# Set the momenta corresponding to T=300K
MaxwellBoltzmannDistribution(atoms, 300*units.kB)
print 'Removing linear momentum and angular momentum'
Stationary(atoms) # zero linear momentum
ZeroRotation(atoms) # zero angular momentum

# We want to run MD using the VelocityVerlet algorithm.
dyn = VelocityVerlet(atoms, 0.1*units.fs, trajectory='moldyn4.traj') # save trajectory.

#Function to print the potential, kinetic and total energy.

from ase.data.molecules import molecule
from ase.optimize import QuasiNewton
from gpaw import GPAW

for name in ['H2', 'N2', 'O2', 'NO']:
    mol = molecule(name)
    mol.center(vacuum=5.0)
    if name == 'NO':
        mol.translate((0, 0.1, 0))
    calc = GPAW(xc='PBE',
                h=0.2,
                stencils=(3, 3),
                txt=name + '.txt')
    mol.set_calculator(calc)
  
    opt = QuasiNewton(mol, logfile=name + '.log', trajectory=name + '.traj')
    opt.run(fmax=0.05)
    calc.write(name)

                         Hamiltonian_MaxAngularMomentum_='',
                         Hamiltonian_MaxAngularMomentum_H='"s"',
                         ))

c = Hookean_Always(a1=0, a2=1, k=0, rt=0.6)
atoms.set_constraint(c)
OP = []

def print_distance(a=atoms):
        distance = a.get_distance(a0=0,a1=1)
	OP.append(distance)
        epot = a.get_potential_energy() / len(a)
        ekin = a.get_kinetic_energy() / len(a)
        ETOT = epot+ekin
        OP.append(ETOT)

dyn = QuasiNewton(atoms, trajectory='atoms.traj')
dyn.attach(print_distance)
dyn.run(steps=10)
print atoms.get_distance(a0=0, a1=1)
write('test.final.xyz', atoms)

OP_f = [OP[i:i+2] for i in xrange(0,len(OP),2)]
OP_f_text = np.savetxt('out.tmp', OP_f)         # reads array and saves into a file as a string
OP_f_infile = open('out.tmp','r')       # open file to write string array onto
OP_f_text = OP_f_infile.read()

text_file.write(OP_f_text)
text_file.close()
                     

# see the module for the required format of reactions definition
from ase.test.cmr.reactions import reactions

# assure that all reactions define a reaction_id
for r in reactions:
    assert r[-1][0] == 'reaction_id'

optimize = True

calculator = EMT()

# find names of compounds
# (in one of the most obscure ways - python list flattening)
compounds = [c[0] for c in sum([r[:-1] for r in reactions], [])]
# unique
compounds = list(set(compounds))

for formula in compounds:
    m = molecule(formula)
    m.set_calculator(calculator)
    if optimize:
        dyn = QuasiNewton(m,
                          logfile=('%s.log' % formula),
                          trajectory=('%s.traj' % formula),
                          )
        dyn.run()
    else:
        e = m.get_potential_energy()
        write(filename=('%s.traj' % formula), images=m, format='traj')

d = a / sqrt(2)
y = d * sqrt(3) / 2
fcc111 = Atoms('Cu',
               cell=[(d, 0, 0),
                     (d / 2, y, 0),
                     (d / 2, y / 3, -a / sqrt(3))],
               pbc=True)
slab = fcc111 * (2, 2, 4)
slab.set_cell([2 * d, 2 * y, 1])
slab.set_pbc((1, 1, 0))
slab.set_calculator(EMT())
Z = slab.get_positions()[:, 2]
indices = [i for i, z in enumerate(Z) if z < Z.mean()]
constraint = FixAtoms(indices=indices)
slab.set_constraint(constraint)
dyn = QuasiNewton(slab)
dyn.run(fmax=0.05)
Z = slab.get_positions()[:, 2]
print Z[0] - Z[1]
print Z[1] - Z[2]
print Z[2] - Z[3]

b = 1.2
h = 2.0
slab += Atom('C', (d, 2 * y / 3, h))
slab += Atom('O', (3 * d / 2, y / 3, h))
traj = PickleTrajectory('initial.traj', 'w', slab)
dyn = QuasiNewton(slab)
dyn.attach(traj.write)
dyn.run(fmax=0.05)
#view(slab)