Python vibrations.Vibrations类代码示例

    def __init__(self, atoms, Excitations,
                 indices=None,
                 gsname='rraman',  # name for ground state calculations
                 exname=None,      # name for excited state calculations
                 delta=0.01,
                 nfree=2,
                 directions=None,
                 approximation='Profeta',
                 observation={'geometry': '-Z(XX)Z'},
                 exkwargs={},      # kwargs to be passed to Excitations
                 exext='.ex.gz',   # extension for Excitation names
                 txt='-',
                 verbose=False,):
        assert(nfree == 2)
        Vibrations.__init__(self, atoms, indices, gsname, delta, nfree)
        self.name = gsname + '-d%.3f' % delta
        if exname is None:
            exname = gsname
        self.exname = exname + '-d%.3f' % delta
        self.exext = exext

        if directions is None:
            self.directions = np.array([0, 1, 2])
        else:
            self.directions = np.array(directions)

        self.approximation = approximation
        self.observation = observation
        self.exobj = Excitations
        self.exkwargs = exkwargs

        self.timer = Timer()
        self.txt = convert_string_to_fd(txt)

        self.verbose = verbose

def run(job, atoms):
    kwargs = {'mode': 0,
          'label': 'vib-{0}'.format(job, job),
          'xc': 'PBE', 

          'scf_guess': 'atomic',
          'max_scf': 500,
          'EPS_SCF': 5.0E-7,
          'added_mos': 500,

          'sme/method': 'fermi_dirac',
          'ELECTRONIC_TEMPERATURE': 300,
          
          'DIA/ALGORITHM': 'STANDARD',

          'mix/METHOD': 'BROYDEN_MIXING',
          'ALPHA': 0.1,
          'BETA': 1.5,
          'NBUFFER': 8,

          'cpu': 36,
          'cutoff': 300,
          'run_type': 'ENERGY_FORCE',  # ENERGY_FORCE, GEO_OPT, CELL_OPT, MD
          'atoms': atoms,
        }
    calc = CP2K(**kwargs)
    atoms.set_calculator(calc)
    vib = Vibrations(atoms, indices = [65, 69])
    vib.run()
    vib.summary()

    def __init__(self, atoms, Excitations,
                 indices=None,
                 gsname='rraman',  # name for ground state calculations
                 exname=None,      # name for excited state calculations
                 delta=0.01,
                 nfree=2,
                 directions=None,
                 exkwargs={},      # kwargs to be passed to Excitations
                 txt='-'):
        assert(nfree == 2)
        Vibrations.__init__(self, atoms, indices, gsname, delta, nfree)
        self.name = gsname + '-d%.3f' % delta
        if exname is None:
            exname = gsname
        self.exname = exname + '-d%.3f' % delta

        if directions is None:
            self.directions = np.array([0, 1, 2])
        else:
            self.directions = np.array(directions)

        self.exobj = Excitations
        self.exkwargs = exkwargs

        self.timer = Timer()
        self.txt = get_txt(txt, rank)

 def read(self, method='standard', direction='central'):
     """Read data from a pre-performed calculation."""
     if not hasattr(self, 'modes'):
         self.timer.start('read vibrations')
         Vibrations.read(self, method, direction)
         # we now have:
         # self.H     : Hessian matrix
         # self.im    : 1./sqrt(masses)
         # self.modes : Eigenmodes of the mass weighted H
         self.om_r = self.hnu.real    # energies in eV
         self.timer.stop('read vibrations')
     if not hasattr(self, 'ex0E_p'):
         self.read_excitations()

    def __init__(self, atoms, vibname, minfreq=None, maxfreq=None):
        """Input is a atoms object and the corresponding vibrations. 
        With minfreq and maxfreq frequencies can    
        be excluded from the calculation"""

        self.atoms = atoms
        # V = a * v is the combined atom and xyz-index
        self.mm05_V = np.repeat(1.0 / np.sqrt(atoms.get_masses()), 3)
        self.minfreq = minfreq
        self.maxfreq = maxfreq
        self.shape = (len(self.atoms), 3)

        vib = Vibrations(atoms, name=vibname)
        self.energies = np.real(vib.get_energies(method="frederiksen"))  # [eV]
        self.frequencies = np.real(vib.get_frequencies(method="frederiksen"))  # [cm^-1]
        self.modes = vib.modes
        self.H = vib.H

def run(label, atoms):
    #calc.mode = 1
    calc.directory =  'vib/pt/{0}'.format(label)
    calc.prefix = 'al2o3-pt-{0}'.format(label)
    calc.results = {}
    calc.CP2K_INPUT.FORCE_EVAL_list[0].DFT.Wfn_restart_file_name = 'al2o3-pt-{0}-RESTART.wfn'.format(label)
    calc.CP2K_INPUT.MOTION.CONSTRAINT.FIXED_ATOMS_list = []
    #===============================================================================
    atoms.set_calculator(calc)
    ###calc.write_input_file()
    #e = atoms.get_potential_energy()
    #t = calc.get_time()
    print(' {0} '.format(label))
    vib = Vibrations(atoms, indices = [120])
    vib.run()
    vib.summary()
    import os, shutil
    for file in os.listdir('.'):
        if "pckl" in file:
            shutil.move(file,calc.directory) 

               grid='nodisk',
               tolerances='tight',
               maxiter=777,
               convergence={'lshift': 0.0},
               basis=basis,
               basispar='spherical',
               direct='noio',
               raw='set int:txs:limxmem 134217728\nmemory total 8000 Mb noverify\n',
               label=label)
 atoms.set_calculator(calc)
 t = time.time()
 if len(atoms) > 1:
     # need a copy of atoms for calculation of vibrations
     vibatoms = atoms.copy()
     vibatoms.set_calculator(calc)
     vib = Vibrations(vibatoms, name=label + '_fixed')
     vib.run()
     f = vib.get_frequencies()[-1].real
     # clean nwchem restart
     if os.path.exists(name + '_' + code + '.db'):
         os.remove(name + '_' + code + '.db')
     atoms.get_potential_energy()
     c.write(atoms, name=name, relaxed=False, basis=basis,
             frequency=f,
             time=time.time()-t)
 else:
     atoms.get_potential_energy()
     c.write(atoms, name=name, relaxed=False, basis=basis,
             time=time.time()-t)
 if len(atoms) > 1:
     opt = BFGS(atoms,

# name of output file for free energies
output_name = 'out.energy'

### At 300K and 101325 Pa
### change for your operating conditions
T = 300     # K
P = 101325  # Pa

#########################################################################################################
#####                                     END                                                       #####
#########################################################################################################

energy = atoms.get_potential_energy()                 # caclulate the energy, to be used to determine G
vibrateatoms = [atom.index for atom in atoms if atom.symbol in ['H','N']]   # calculate the vibrational modes for all N and H atoms
# Calculate vibrations
vib = Vibrations(atoms,indices=vibrateatoms,delta=0.03)    # define a vibration calculation
vib.run()                                                  # run the vibration calculation
vib.summary(method='standard')                             # summarize the calculated results

for mode in range(len(vibrateatoms)*3):                    # Make trajectory files to visualize the modes.
    vib.write_mode(mode)

vibenergies=vib.get_energies()
vibenergies=[vib for vib in vibenergies if not isinstance(vib,complex)]  # only take the real modes
gibbs = HarmonicThermo(vib_energies = vibenergies, electronicenergy = energy)

freeenergy = gibbs.get_gibbs_energy(T,P)

f=open(output_name,'w')
f.write('Potential energy: '+str(energy)+'\n'+'Free energy: '+str(freeenergy)+'\n')
f.close

"""Calculate the vibrational modes of a H2O molecule."""

from ase.vibrations import Vibrations
from gpaw import GPAW

h2o = GPAW('h2o.gpw', txt=None).get_atoms()

# Create vibration calculator
vib = Vibrations(h2o)
vib.run()
vib.summary(method='frederiksen')

# Make trajectory files to visualize normal modes:
for mode in range(9):
    vib.write_mode(mode)

from __future__ import print_function
from ase import Atoms
from ase.calculators.emt import EMT
from ase.optimize import QuasiNewton
from ase.vibrations import Vibrations
from ase.thermochemistry import IdealGasThermo

n2 = Atoms('N2',
           positions=[(0, 0, 0), (0, 0, 1.1)],
           calculator=EMT())
QuasiNewton(n2).run(fmax=0.01)
vib = Vibrations(n2)
vib.run()
print(vib.get_frequencies())
vib.summary()
print(vib.get_mode(-1))
vib.write_mode(n=None, nimages=20)
vib_energies = vib.get_energies()

thermo = IdealGasThermo(vib_energies=vib_energies, geometry='linear',
                        atoms=n2, symmetrynumber=2, spin=0)
thermo.get_gibbs_energy(temperature=298.15, pressure=2 * 101325.)

assert vib.clean(empty_files=True) == 0
assert vib.clean() == 13

#!/usr/bin/env python
from ase.db import connect
from ase.calculators.aims import Aims
from ase.lattice.surface import fcc111, add_adsorbate
from ase.constraints import FixAtoms
from ase import Atoms, Atom
from ase.io.aims import read_aims
from ase.optimize import BFGS
from ase.vibrations import Vibrations

mydb = connect("mydb.db")
atoms = mydb.get_atoms(name='pt-co-relax')

calc = Aims(label='cluster/pt-co-vib-cons',
          xc='pbe',
          spin='none',
          relativistic = 'atomic_zora scalar',
          sc_accuracy_etot=1e-7,
          sc_accuracy_eev=1e-3,
          sc_accuracy_rho=1e-4,
          sc_accuracy_forces=1e-3)

atoms.set_calculator(calc)

vib = Vibrations(atoms, indices=[1, 2], name='cluster/pt-co-vib-cons/vib')
vib.run()
vib.summary()

def calculate(element, ref_data, p):
    values_dict = {}
    values_dict[p['xc']] = {}
    for XY, data in ref_data[p['xc']].items():
        X = XY.split('-')[0]
        Y = XY.split('-')[1]
        if (X == Y and X == element) or (X != Y and (X == element or Y == element)):
            # compound contains the requested element
            re_ref = data['re']
            we_ref = data['we']
            m0_ref = data.get('m0', 0.0)
            #
            compound = Atoms(X+Y,
                             [
                (0,       0,     0.5  ),
                (0,       0,     0.5+re_ref/a),
                ],
                             pbc=0)
            compound.set_cell([a, b, c], scale_atoms=1)
            compound.center()

            # calculation on the reference geometry
            calc = Calculator(**p)
            compound.set_calculator(calc)
            e_compound = compound.get_potential_energy()
            finegd = calc.density.finegd
            dip = finegd.calculate_dipole_moment(calc.density.rhot_g)*calc.a0
            vib = Vibrations(compound)
            vib.run()
            vib_compound = vib.get_frequencies(method='frederiksen').real[-1]
            world.barrier()
            vib_pckl = glob('vib.*.pckl')
            if rank == 0:
                for file in vib_pckl: remove(file)

            # calculation on the relaxed geometry
            qn = QuasiNewton(compound)
            #qn.attach(PickleTrajectory('compound.traj', 'w', compound).write)
            qn.run(fmax=0.05)
            e_compound_r = compound.get_potential_energy()
            dist_compound_r = compound.get_distance(0,1)
            dip_r = finegd.calculate_dipole_moment(calc.density.rhot_g)*calc.a0
            vib = Vibrations(compound)
            vib.run()
            vib_compound_r = vib.get_frequencies(method='frederiksen').real[-1]
            world.barrier()
            vib_pckl = glob('vib.*.pckl')
            if rank == 0:
                for file in vib_pckl: remove(file)

            del compound
            e = e_compound
            we = vib_compound
            m0 = dip
            e_r = e_compound_r
            we_r = vib_compound_r
            re_r = dist_compound_r
            m0_r = dip_r
            #
            values_dict[p['xc']][XY] = {'re': re_r, 'we': (we_r, we), 'm0': (m0_r, m0)}
    #
    return values_dict

H       5.320549047980879      3.220584852467720      3.974551561510350
H       7.723359150977955      3.224855971783890      4.574146712279462
H       7.580803493981530      5.034479218283977      4.877211530909463
"""

h = 0.3
atoms = Cluster(read_xyz(StringIO.StringIO(butadiene)))
atoms.minimal_box(3.0, h)
atoms.set_calculator(GPAW(h=h))
if 0:
    dyn = FIRE(atoms)
    dyn.run(fmax=0.05)
    atoms.write("butadiene.xyz")

vibname = "fcvib"
vib = Vibrations(atoms, name=vibname)
vib.run()

# Modul
a = FranckCondon(atoms, vibname, minfreq=250)

# excited state forces
F = np.array(
    [
        [-2.11413, 0.07317, -0.91682],
        [3.23569, -0.74520, 0.76758],
        [-3.44847, 0.63846, -0.81080],
        [2.77345, 0.01272, 0.74811],
        [-0.06544, -0.01078, -0.03209],
        [-0.01245, -0.01123, -0.00040],
        [0.00186, -0.05864, -0.00371],

#!/usr/bin/env python
from ase import Atom, Atoms
from ase.structure import molecule
from ase.calculators.cp2k import CP2K
from ase.optimize import BFGS
from ase.vibrations import Vibrations
from multiprocessing import Pool
import os

atoms = molecule('CO')
atoms.center(vacuum=2.0)


calc = CP2K(label = 'molecules/co/vib',
      xc='PBE')
atoms.set_calculator(calc)
gopt = BFGS(atoms, logfile=None)
gopt.run(fmax=1e-2)
e = atoms.get_potential_energy()
pos = atoms.get_positions()
d = ((pos[0] - pos[1])**2).sum()**0.5
print('{0:1.4f}  {1:1.4f} '.format( e, d))

vib = Vibrations(n2)
vib.run()
vib.summary()
vib.write_mode(-1)

                                      # the valence electrons
                      sigma=0.1,
                      psppath='/home/vossj/suncat/psp/gbrv1.5pbe',  # pseudopotential
                      convergence={'energy': 1e-5,
                                   'mixing': 0.1,
                                   'nmix': 10,
                                   'mix': 4,
                                   'maxsteps': 500,
                                   'diag': 'david'
                                   },  # convergence parameters
                      outdirprefix='calcdirv')  # output directory for Quantum Espresso files

atoms.set_calculator(calcvib)


vib = Vibrations(atoms, indices=vibrateatoms, delta=0.03)
vib.run()
vib.summary(method='standard')

# Make trajectory files to visualize the modes.
for mode in range(len(vibrateatoms) * 3):
    vib.write_mode(mode)

# Calculate free energy
vib_energies = vib.get_energies()
thermo = IdealGasThermo(vib_energies=vib_energies,
                        electronicenergy=energy,
                        atoms=atoms,
                        geometry='linear',
                        symmetrynumber=2, spin=0)