Source code for ase.calculators.nwchem

"""This module defines an ASE interface to NWchem

http://www.nwchem-sw.org/
"""
import os
import re
import numpy as np

from warnings import warn
from ase.atoms import Atoms
from ase.units import Hartree, Bohr
from ase.io.nwchem import write_nwchem
from ase.calculators.calculator import (FileIOCalculator, Parameters,
                                        ReadError, EigenvalOccupationMixin)


class KPoint:
    def __init__(self, s):
        self.s = s
        self.eps_n = []
        self.f_n = []


[docs]class NWChem(FileIOCalculator, EigenvalOccupationMixin): implemented_properties = ['energy', 'forces', 'dipole', 'magmom', 'orbitals'] command = 'nwchem PREFIX.nw > PREFIX.out' default_parameters = dict( xc='LDA', smearing=None, charge=None, task='gradient', # Warning: nwchem centers atoms by default # see ase-developers/2012-March/001356.html geometry='nocenter noautosym', convergence={'energy': None, 'density': None, 'gradient': None, 'lshift': None, # set lshift to 0.0 for nolevelshifting 'damp': None, }, basis='3-21G', basispar=None, ecp=None, so=None, spinorbit=False, tddft=False, odft=False, raw='') # additional outside of dft block control string def __init__(self, restart=None, ignore_bad_restart_file=False, label='nwchem', atoms=None, **kwargs): """Construct NWchem-calculator object.""" FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, label, atoms, **kwargs) def set(self, **kwargs): changed_parameters = FileIOCalculator.set(self, **kwargs) if changed_parameters: self.reset() def check_state(self, atoms): system_changes = FileIOCalculator.check_state(self, atoms) # Ignore unit cell and boundary conditions: if 'cell' in system_changes: system_changes.remove('cell') if 'pbc' in system_changes: system_changes.remove('pbc') return system_changes def write_input(self, atoms, properties=None, system_changes=None): FileIOCalculator.write_input(self, atoms, properties, system_changes) p = self.parameters p.initial_magmoms = atoms.get_initial_magnetic_moments().tolist() p.write(self.label + '.ase') del p['initial_magmoms'] f = open(self.label + '.nw', 'w') if p.charge is not None: f.write('charge %s\n' % p.charge) write_nwchem(f, atoms, p.geometry) f.write('start\n') if p.basispar is not None: basispar = 'basis ' + p.basispar else: basispar = 'basis' def format_basis_set(string, tag=basispar): formatted = tag + '\n' lines = string.split('\n') if len(lines) > 1: formatted += string else: formatted += ' * library ' + string return formatted + '\nend\n' basis = format_basis_set(p.basis) if p.ecp is not None: basis += format_basis_set(p.ecp, 'ecp') if p.so is not None: basis += format_basis_set(p.so, 'so') f.write(basis) if p.xc == 'RHF': task = 'scf' elif p.xc == 'MP2': task = 'mp2' else: if p.spinorbit: task = 'sodft' elif p.tddft: task = 'tddft' else: task = 'dft' xc = {'LDA': 'slater pw91lda', 'PBE': 'xpbe96 cpbe96', 'revPBE': 'revpbe cpbe96', 'RPBE': 'rpbe cpbe96', 'CAM-B3LYP': 'xcamb88 1.00 lyp 0.81 vwn_5 0.19 hfexch 1.00'}.get(p.xc, p.xc) if p.tddft: f.write('\n' + 'dft' + '\n') else: f.write('\n' + task + '\n') f.write(' xc ' + xc + '\n') if p.xc == 'CAM-B3LYP': f.write(' cam 0.33 cam_alpha 0.19 cam_beta 0.46\n') f.write(' direct\n') for key in p.convergence: if p.convergence[key] is not None: if key == 'lshift': if p.convergence[key] <= 0.0: f.write(' convergence nolevelshifting\n') else: f.write(' convergence %s %s\n' % (key, p.convergence[key] / Hartree)) else: f.write(' convergence %s %s\n' % (key, p.convergence[key])) if p.smearing is not None: assert p.smearing[0].lower() == 'gaussian', p.smearing f.write(' smear %s\n' % (p.smearing[1] / Hartree)) if 'mult' not in p: # Obtain multiplicity from magnetic momenta: tot_magmom = atoms.get_initial_magnetic_moments().sum() if tot_magmom < 0: mult = tot_magmom - 1 # fill minority bands else: mult = tot_magmom + 1 else: mult = p.mult if mult != int(mult): raise RuntimeError('Noninteger multiplicity not possible. ' + 'Check initial magnetic moments.') f.write(' mult %d\n' % mult) if p.odft: f.write(' odft\n') # open shell aka spin polarized dft for key in sorted(p.keys()): if key in ['charge', 'geometry', 'basis', 'basispar', 'ecp', 'so', 'xc', 'spinorbit', 'convergence', 'smearing', 'raw', 'mult', 'task', 'odft', 'tddft']: continue f.write(u" {0} {1}\n".format(key, p[key])) f.write('end\n') if p.raw: f.write(p.raw + '\n') f.write('\ntask ' + task + ' ' + p.task + '\n') f.close() def read(self, label): FileIOCalculator.read(self, label) if not os.path.isfile(self.label + '.out'): raise ReadError f = open(self.label + '.nw') for line in f: if line.startswith('geometry'): break symbols = [] positions = [] for line in f: if line.startswith('end'): break words = line.split() symbols.append(words[0]) positions.append([float(word) for word in words[1:]]) self.parameters = Parameters.read(self.label + '.ase') self.atoms = Atoms(symbols, positions, magmoms=self.parameters.pop('initial_magmoms')) self.read_results() def read_results(self): self.read_energy() if self.parameters.task.find('gradient') > -1: self.read_forces() if self.parameters.task.find('optimize') > -1: self.read_coordinates() self.read_forces() if self.parameters.task.find('qmd') > -1: self.read_coordinates() self.read_forces() self.niter = self.read_number_of_iterations() self.nelect = self.read_number_of_electrons() self.nvector = self.read_number_of_bands() self.results['magmom'] = self.read_magnetic_moment() dipole = self.read_dipole_moment() if dipole is not None: self.results['dipole'] = dipole self.read_mos() def get_ibz_k_points(self): return np.array([0., 0., 0.]) def get_number_of_bands(self): return self.nvector def read_number_of_bands(self): nvector = 0 for line in open(self.label + '.out'): if line.find('Vector ') != -1: # count all printed vectors nvector += 1 if not nvector: nvector = None return nvector def get_number_of_electrons(self): return self.nelect def read_number_of_electrons(self): nelect = None for line in open(self.label + '.out'): # find last one if line.find('of electrons') != -1: nelect = float(line.split(':')[1].strip()) return nelect def get_number_of_iterations(self): return self.niter def read_number_of_iterations(self): niter = 0 for line in open(self.label + '.out'): if line.find('d= ') != -1: # count all iterations niter += 1 if not niter: niter = None return niter def read_magnetic_moment(self): magmom = None for line in open(self.label + '.out'): if line.find('Spin multiplicity') != -1: # last one magmom = float(line.split(':')[-1].strip()) if magmom < 0: magmom += 1 else: magmom -= 1 return magmom def read_dipole_moment(self): dipolemoment = [] for line in open(self.label + '.out'): for component in ['1 1 0 0', '1 0 1 0', '1 0 0 1']: if line.find(component) != -1: value = float(line.split(component)[1].split()[0]) value = value * Bohr dipolemoment.append(value) if len(dipolemoment) == 0: if len(self.atoms) == 1: dipolemoment = [0.0, 0.0, 0.0] else: return None return np.array(dipolemoment) def read_energy(self): """Read Energy from nwchem output file.""" text = open(self.label + '.out', 'r').read() lines = iter(text.split('\n')) # Energy: estring = 'Total ' if self.parameters.xc == 'RHF': estring += 'SCF' elif self.parameters.xc == 'MP2': estring += 'MP2' else: estring += 'DFT' estring += ' energy' for line in lines: if line.find(estring) >= 0: energy = float(line.split()[-1]) # Excited state energy replaces ground state, if present if (self.parameters.tddft): lines = text.split('\n') estring = 'Excited state energy' for line in lines: if line.find(estring) >= 0: energy = float(line.split()[-1]) self.results['energy'] = energy * Hartree # All lines have been 'eaten' while iterating; loop from scratch. # (We could insert a break, but I (askhl) don't know if multiple # energies might be listed; then we would not get the last one. lines = text.split('\n') # Eigenstates spin = -1 kpts = [] for line in lines: if line.find('Molecular Orbital Analysis') >= 0: last_eps = -99999.0 spin += 1 kpts.append(KPoint(spin)) if spin >= 0: if line.find('Vector') >= 0: line = line.lower().replace('d', 'e') line = line.replace('=', ' ') word = line.split() this_occ = float(word[3]) this_eps = float(word[5]) kpts[spin].f_n.append(this_occ) kpts[spin].eps_n.append(this_eps) if this_occ < 0.1 and this_eps < last_eps: warn('HOMO above LUMO - if this is not an exicted ' + 'state - this might be introduced by levelshift.', RuntimeWarning) last_eps = this_eps self.kpts = kpts def read_forces(self): """Read Forces from nwchem output file.""" file = open(self.label + '.out', 'r') lines = file.readlines() file.close() for i, line in enumerate(lines): if line.find('ENERGY GRADIENTS') >= 0: gradients = [] for j in range(i + 4, i + 4 + len(self.atoms)): word = lines[j].split() gradients.append([float(word[k]) for k in range(5, 8)]) self.results['forces'] = -np.array(gradients) * Hartree / Bohr def read_coordinates(self): """Read updated coordinates from nwchem output file.""" file = open(self.label + '.out', 'r') lines = file.readlines() file.close() for i, line in enumerate(lines): if line.find('ENERGY GRADIENTS') >= 0: positions = [] for j in range(i + 4, i + 4 + len(self.atoms)): word = lines[j].split() positions.append([float(word[k]) for k in range(2, 5)]) self.atoms.set_positions(np.array(positions) * Bohr) def read_mos(self): """ Read the molecular orbitals and molecular orbital energies The output can be requested by specification of the print keyword {'print': '\"final vectors analysis\" \"final vectors\"'} when constructing the calculator object. """ orbitals = [] orbital_coefficients = [] orbital_coefficients_line = [] filename = self.label + '.out' orb_number = None mo_read = False mo_init = False ev_number = 0 for line in open(filename, 'r'): if 'Final MO vectors' in line: mo_read = True mo_init = False orb_number = None orbital_coefficients = [] orbital_coefficients_line = [] continue if mo_read: regex = ( r'^\s*(\d+)(\s+[+-]?\d+\.\d+)(\s+[+-]?\d+\.\d+)?' r'(\s+[+-]?\d+\.\d+)?(\s+[+-]?\d+\.\d+)?' r'(\s+[+-]?\d+\.\d+)?(\s+[+-]?\d+\.\d+)?' ) match = re.search(regex, line) if match: orb_number = int(match.group(1)) for i in range(2, 8): if match.group(i): coeff = float(match.group(i)) orbital_coefficients_line.append(coeff) if mo_init: orbital_coefficients[orb_number-1].extend( orbital_coefficients_line ) else: orbital_coefficients.append(orbital_coefficients_line) orbital_coefficients_line = [] elif orb_number: if orb_number == len(orbital_coefficients): mo_init = True if orb_number == len(orbital_coefficients[0]): mo_read = False continue if 'Vector' in line: regex = ( r'Vector\s+(\d+)\s+Occ=\s*([\w\.+-]+)\s+E=\s*([+-\.\w]+)' ) match = re.search(regex, line) if match: ev_number += 1 # for the case of geometry optimization output if int(match.group(1)) < ev_number: orbitals = [] ev_number = int(match.group(1)) orbitals.append({ 'number': int(match.group(1)), 'occupation': float(re.sub('[dD]', 'E', match.group(2))), 'energy': float(re.sub('[dD]', 'E', match.group(3))) }) if len(orbital_coefficients) > 0: # this is to transpose the nested list mos = [[x[i] for x in orbital_coefficients] for i in range(len(orbital_coefficients[0]))] for orbital in orbitals: orbital['coefficients'] = mos.pop(0) self.results['orbitals'] = orbitals def get_eigenvalues(self, kpt=0, spin=0): """Return eigenvalue array.""" return np.array(self.kpts[spin].eps_n) * Hartree def get_occupation_numbers(self, kpt=0, spin=0): """Return occupation number array.""" return self.kpts[spin].f_n def get_number_of_spins(self): """Return the number of spins in the calculation. Spin-paired calculations: 1, spin-polarized calculation: 2.""" return len(self.kpts) def get_spin_polarized(self): """Is it a spin-polarized calculation?""" return len(self.kpts) == 2 def get_mos(self, atoms): """ return eigenvectors and eigenvalues as matrices (numpy arrays) """ orbital_coefficients = [] eigen_energies = [] for orbital in self.get_property('orbitals'): orbital_coefficients.append(orbital['coefficients']) eigen_energies.append(orbital['energy']) return [ np.array(orbital_coefficients), np.diag(np.array(eigen_energies)*Hartree) ]