Python Vibrations.run方法代码示例

# 需要导入模块: from ase.vibrations import Vibrations [as 别名]
# 或者: from ase.vibrations.Vibrations import run [as 别名]
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()

# 需要导入模块: from ase.vibrations import Vibrations [as 别名]
# 或者: from ase.vibrations.Vibrations import run [as 别名]
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) 

# 需要导入模块: from ase.vibrations import Vibrations [as 别名]
# 或者: from ase.vibrations.Vibrations import run [as 别名]
from __init__ import AnharmonicModes

slab = fcc111('Al', size=(2, 2, 2), vacuum=3.0)
CH3 = molecule('CH3')
add_adsorbate(slab, CH3, 2.5, 'ontop')

constraint = FixAtoms(mask=[a.symbol == 'Al' for a in slab])
slab.set_constraint(constraint)
slab.set_calculator(EMT())

dyn = QuasiNewton(slab, logfile='/dev/null')
dyn.run(fmax=0.05)

vib = Vibrations(slab, indices=[8, 9, 10, 11])
vib.run()
vib.summary(log='/dev/null')
vib.clean()

AM = AnharmonicModes(vibrations_object=vib)
rot_mode = AM.define_rotation(
    basepos=[0., 0., -1.],
    branch=[9, 10, 11],
    symnumber=3)

AM.run()
AM.summary(log='/dev/null')
AM.clean()

# print(AM.get_ZPE(), AM.get_entropic_energy())
assert abs(AM.get_ZPE() - 0.388) < 1e-3, AM.get_ZPE()

# 需要导入模块: from ase.vibrations import Vibrations [as 别名]
# 或者: from ase.vibrations.Vibrations import run [as 别名]
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

# 需要导入模块: from ase.vibrations import Vibrations [as 别名]
# 或者: from ase.vibrations.Vibrations import run [as 别名]
# 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