root / ase / calculators / abinit.py @ 1
Historique | Voir | Annoter | Télécharger (21,37 ko)
| 1 |
"""This module defines an ASE interface to ABINIT.
|
|---|---|
| 2 |
|
| 3 |
http://www.abinit.org/
|
| 4 |
"""
|
| 5 |
|
| 6 |
import os |
| 7 |
from glob import glob |
| 8 |
from os.path import join, isfile, islink |
| 9 |
|
| 10 |
import numpy as np |
| 11 |
|
| 12 |
from ase.data import chemical_symbols |
| 13 |
from ase.data import atomic_numbers |
| 14 |
from ase.units import Bohr, Hartree |
| 15 |
|
| 16 |
|
| 17 |
class Abinit: |
| 18 |
"""Class for doing ABINIT calculations.
|
| 19 |
|
| 20 |
The default parameters are very close to those that the ABINIT
|
| 21 |
Fortran code would use. These are the exceptions::
|
| 22 |
|
| 23 |
calc = Abinit(label='abinit', xc='LDA', pulay=5, mix=0.1)
|
| 24 |
|
| 25 |
Use the set_inp method to set extra INPUT parameters::
|
| 26 |
|
| 27 |
calc.set_inp('nstep', 30)
|
| 28 |
|
| 29 |
"""
|
| 30 |
def __init__(self, label='abinit', xc='LDA', kpts=None, nbands=0, |
| 31 |
width=0.04*Hartree, ecut=None, charge=0, |
| 32 |
pulay=5, mix=0.1, pps='fhi', toldfe=1.0e-6 |
| 33 |
): |
| 34 |
"""Construct ABINIT-calculator object.
|
| 35 |
|
| 36 |
Parameters
|
| 37 |
==========
|
| 38 |
label: str
|
| 39 |
Prefix to use for filenames (label.in, label.txt, ...).
|
| 40 |
Default is 'abinit'.
|
| 41 |
xc: str
|
| 42 |
Exchange-correlation functional. Must be one of LDA, PBE,
|
| 43 |
revPBE, RPBE.
|
| 44 |
kpts: list of three int
|
| 45 |
Monkhost-Pack sampling.
|
| 46 |
nbands: int
|
| 47 |
Number of bands.
|
| 48 |
Default is 0.
|
| 49 |
width: float
|
| 50 |
Fermi-distribution width in eV.
|
| 51 |
Default is 0.04 Hartree.
|
| 52 |
ecut: float
|
| 53 |
Planewave cutoff energy in eV.
|
| 54 |
No default.
|
| 55 |
charge: float
|
| 56 |
Total charge of the system.
|
| 57 |
Default is 0.
|
| 58 |
pulay: int
|
| 59 |
Number of old densities to use for Pulay mixing.
|
| 60 |
mix: float
|
| 61 |
Mixing parameter between zero and one for density mixing.
|
| 62 |
|
| 63 |
Examples
|
| 64 |
========
|
| 65 |
Use default values:
|
| 66 |
|
| 67 |
>>> h = Atoms('H', calculator=Abinit())
|
| 68 |
>>> h.center(vacuum=3.0)
|
| 69 |
>>> e = h.get_potential_energy()
|
| 70 |
|
| 71 |
"""
|
| 72 |
|
| 73 |
if not nbands > 0: |
| 74 |
raise ValueError('Number of bands (nbands) not set') |
| 75 |
|
| 76 |
if ecut is None: |
| 77 |
raise ValueError('Planewave cutoff energy in eV (ecut) not set') |
| 78 |
|
| 79 |
self.label = label#################### != out |
| 80 |
self.xc = xc
|
| 81 |
self.kpts = kpts
|
| 82 |
self.nbands = nbands
|
| 83 |
self.width = width
|
| 84 |
self.ecut = ecut
|
| 85 |
self.charge = charge
|
| 86 |
self.pulay = pulay
|
| 87 |
self.mix = mix
|
| 88 |
self.pps = pps
|
| 89 |
self.toldfe = toldfe
|
| 90 |
if not pps in ['fhi', 'hgh', 'hgh.sc']: |
| 91 |
raise ValueError('Unexpected PP identifier %s' % pps) |
| 92 |
|
| 93 |
self.converged = False |
| 94 |
self.inp = {}
|
| 95 |
self.n_entries_int = 20 # integer entries per line |
| 96 |
self.n_entries_float = 8 # float entries per line |
| 97 |
|
| 98 |
def update(self, atoms): |
| 99 |
if (not self.converged or |
| 100 |
len(self.numbers) != len(atoms) or |
| 101 |
(self.numbers != atoms.get_atomic_numbers()).any()):
|
| 102 |
self.initialize(atoms)
|
| 103 |
self.calculate(atoms)
|
| 104 |
elif ((self.positions != atoms.get_positions()).any() or |
| 105 |
(self.pbc != atoms.get_pbc()).any() or |
| 106 |
(self.cell != atoms.get_cell()).any()):
|
| 107 |
self.calculate(atoms)
|
| 108 |
|
| 109 |
def initialize(self, atoms): |
| 110 |
self.numbers = atoms.get_atomic_numbers().copy()
|
| 111 |
self.species = []
|
| 112 |
for a, Z in enumerate(self.numbers): |
| 113 |
if Z not in self.species: |
| 114 |
self.species.append(Z)
|
| 115 |
|
| 116 |
if 'ABINIT_PP_PATH' in os.environ: |
| 117 |
pppaths = os.environ['ABINIT_PP_PATH'].split(':') |
| 118 |
else:
|
| 119 |
pppaths = [] |
| 120 |
|
| 121 |
self.ppp_list = []
|
| 122 |
if self.xc != 'LDA': |
| 123 |
xcname = 'GGA'
|
| 124 |
else:
|
| 125 |
xcname = 'LDA'
|
| 126 |
|
| 127 |
for Z in self.species: |
| 128 |
symbol = chemical_symbols[abs(Z)]
|
| 129 |
number = atomic_numbers[symbol] |
| 130 |
|
| 131 |
pps = self.pps
|
| 132 |
if pps == 'fhi': |
| 133 |
name = '%02d-%s.%s.fhi' % (number, symbol, xcname)
|
| 134 |
elif pps in ('hgh', 'hgh.sc'): |
| 135 |
hghtemplate = '%d%s.%s.hgh' # E.g. "42mo.6.hgh" |
| 136 |
# There might be multiple files with different valence
|
| 137 |
# electron counts, so we must choose between
|
| 138 |
# the ordinary and the semicore versions for some elements.
|
| 139 |
#
|
| 140 |
# Therefore we first use glob to get all relevant files,
|
| 141 |
# then pick the correct one afterwards.
|
| 142 |
name = hghtemplate % (number, symbol.lower(), '*')
|
| 143 |
|
| 144 |
found = False
|
| 145 |
for path in pppaths: |
| 146 |
if pps.startswith('hgh'): |
| 147 |
filenames = glob(join(path, name)) |
| 148 |
if not filenames: |
| 149 |
continue
|
| 150 |
assert len(filenames) in [0, 1, 2] |
| 151 |
if pps == 'hgh': |
| 152 |
selector = min # Lowest possible valence electron count |
| 153 |
else:
|
| 154 |
assert pps == 'hgh.sc' |
| 155 |
selector = max # Semicore - highest electron count |
| 156 |
Z = selector([int(os.path.split(name)[1].split('.')[1]) |
| 157 |
for name in filenames]) |
| 158 |
name = hghtemplate % (number, symbol.lower(), str(Z))
|
| 159 |
filename = join(path, name) |
| 160 |
if isfile(filename) or islink(filename): |
| 161 |
found = True
|
| 162 |
self.ppp_list.append(filename)
|
| 163 |
break
|
| 164 |
if not found: |
| 165 |
raise RuntimeError('No pseudopotential for %s!' % symbol) |
| 166 |
|
| 167 |
self.converged = False |
| 168 |
|
| 169 |
def get_potential_energy(self, atoms, force_consistent=False): |
| 170 |
self.update(atoms)
|
| 171 |
|
| 172 |
if force_consistent:
|
| 173 |
return self.efree |
| 174 |
else:
|
| 175 |
# Energy extrapolated to zero Kelvin:
|
| 176 |
return (self.etotal + self.efree) / 2 |
| 177 |
|
| 178 |
def get_number_of_bands(self): |
| 179 |
return self.nbands |
| 180 |
|
| 181 |
def get_kpts_info(self, kpt=0, spin=0, mode='eigenvalues'): |
| 182 |
return self.read_kpts_info(kpt, spin, mode) |
| 183 |
|
| 184 |
def get_k_point_weights(self): |
| 185 |
return self.get_kpts_info(kpt=0, spin=0, mode='k_point_weights') |
| 186 |
|
| 187 |
def get_bz_k_points(self): |
| 188 |
raise NotImplementedError |
| 189 |
|
| 190 |
def get_ibz_k_points(self): |
| 191 |
return self.get_kpts_info(kpt=0, spin=0, mode='ibz_k_points') |
| 192 |
|
| 193 |
def get_fermi_level(self): |
| 194 |
return self.read_fermi() |
| 195 |
|
| 196 |
def get_eigenvalues(self, kpt=0, spin=0): |
| 197 |
return self.get_kpts_info(kpt, spin, 'eigenvalues') |
| 198 |
|
| 199 |
def get_occupations(self, kpt=0, spin=0): |
| 200 |
return self.get_kpts_info(kpt, spin, 'occupations') |
| 201 |
|
| 202 |
def get_forces(self, atoms): |
| 203 |
self.update(atoms)
|
| 204 |
return self.forces.copy() |
| 205 |
|
| 206 |
def get_stress(self, atoms): |
| 207 |
self.update(atoms)
|
| 208 |
return self.stress.copy() |
| 209 |
|
| 210 |
def calculate(self, atoms): |
| 211 |
self.positions = atoms.get_positions().copy()
|
| 212 |
self.cell = atoms.get_cell().copy()
|
| 213 |
self.pbc = atoms.get_pbc().copy()
|
| 214 |
|
| 215 |
self.write_files()
|
| 216 |
|
| 217 |
self.write_inp(atoms)
|
| 218 |
|
| 219 |
abinit = os.environ['ABINIT_SCRIPT']
|
| 220 |
locals = {'label': self.label}
|
| 221 |
|
| 222 |
# Now, because (stupidly) abinit when it finds a name it uses nameA
|
| 223 |
# and when nameA exists it uses nameB, etc.
|
| 224 |
# we need to rename our *.txt file to *.txt.bak
|
| 225 |
filename = self.label + '.txt' |
| 226 |
if islink(filename) or isfile(filename): |
| 227 |
os.rename(filename, filename+'.bak')
|
| 228 |
|
| 229 |
execfile(abinit, {}, locals) |
| 230 |
exitcode = locals['exitcode'] |
| 231 |
if exitcode != 0: |
| 232 |
raise RuntimeError(('Abinit exited with exit code: %d. ' + |
| 233 |
'Check %s.log for more information.') %
|
| 234 |
(exitcode, self.label))
|
| 235 |
|
| 236 |
self.read()
|
| 237 |
|
| 238 |
self.converged = True |
| 239 |
|
| 240 |
def write_files(self): |
| 241 |
"""Write input parameters to files-file."""
|
| 242 |
fh = open(self.label + '.files', 'w') |
| 243 |
|
| 244 |
import getpass |
| 245 |
#find a suitable default scratchdir (should be writeable!)
|
| 246 |
username=getpass.getuser() |
| 247 |
|
| 248 |
if os.access("/scratch/"+username,os.W_OK): |
| 249 |
scratch = "/scratch/"+username
|
| 250 |
elif os.access("/scratch/",os.W_OK): |
| 251 |
scratch = "/scratch/"
|
| 252 |
else:
|
| 253 |
if os.access(os.curdir,os.W_OK):
|
| 254 |
scratch = os.curdir #if no /scratch use curdir
|
| 255 |
else:
|
| 256 |
raise IOError,"No suitable scratch directory and no write access to current dir" |
| 257 |
|
| 258 |
fh.write('%s\n' % (self.label+'.in')) # input |
| 259 |
fh.write('%s\n' % (self.label+'.txt')) # output |
| 260 |
fh.write('%s\n' % (self.label+'i')) # input |
| 261 |
fh.write('%s\n' % (self.label+'o')) # output |
| 262 |
# scratch files
|
| 263 |
fh.write('%s\n' % (os.path.join(scratch, self.label+'.abinit'))) |
| 264 |
# Provide the psp files
|
| 265 |
for ppp in self.ppp_list: |
| 266 |
fh.write('%s\n' % (ppp)) # psp file path |
| 267 |
|
| 268 |
fh.close() |
| 269 |
|
| 270 |
def set_inp(self, key, value): |
| 271 |
"""Set INPUT parameter."""
|
| 272 |
self.inp[key] = value
|
| 273 |
|
| 274 |
def write_inp(self, atoms): |
| 275 |
"""Write input parameters to in-file."""
|
| 276 |
fh = open(self.label + '.in', 'w') |
| 277 |
|
| 278 |
inp = {
|
| 279 |
#'SystemLabel': self.label,
|
| 280 |
#'LatticeConstant': 1.0,
|
| 281 |
'natom': len(atoms), |
| 282 |
'charge': self.charge, |
| 283 |
'nband': self.nbands, |
| 284 |
#'DM.UseSaveDM': self.converged,
|
| 285 |
#'SolutionMethod': 'diagon',
|
| 286 |
'npulayit': self.pulay, # default 7 |
| 287 |
'diemix': self.mix |
| 288 |
} |
| 289 |
|
| 290 |
if self.ecut is not None: |
| 291 |
inp['ecut'] = str(self.ecut)+' eV' # default Ha |
| 292 |
|
| 293 |
if self.width is not None: |
| 294 |
inp['tsmear'] = str(self.width)+' eV' # default Ha |
| 295 |
fh.write('occopt 3 # Fermi-Dirac smearing\n')
|
| 296 |
|
| 297 |
inp['ixc'] = { # default 1 |
| 298 |
'LDA': 7, |
| 299 |
'PBE': 11, |
| 300 |
'revPBE': 14, |
| 301 |
'RPBE': 15, |
| 302 |
'WC': 23 |
| 303 |
}[self.xc]
|
| 304 |
|
| 305 |
magmoms = atoms.get_initial_magnetic_moments() |
| 306 |
if magmoms.any():
|
| 307 |
inp['nsppol'] = 2 |
| 308 |
fh.write('spinat\n')
|
| 309 |
for n, M in enumerate(magmoms): |
| 310 |
if M != 0: |
| 311 |
fh.write('%.14f %.14f %.14f\n' % (0, 0, M)) |
| 312 |
else:
|
| 313 |
inp['nsppol'] = 1 |
| 314 |
#fh.write('spinat\n')
|
| 315 |
#for n, M in enumerate(magmoms):
|
| 316 |
# if M != 0:
|
| 317 |
# fh.write('%.14f\n' % (M))
|
| 318 |
|
| 319 |
inp.update(self.inp)
|
| 320 |
|
| 321 |
for key, value in inp.items(): |
| 322 |
if value is None: |
| 323 |
continue
|
| 324 |
|
| 325 |
if isinstance(value, list): |
| 326 |
fh.write('%block %s\n' % key)
|
| 327 |
for line in value: |
| 328 |
fh.write(' '.join(['%s' % x for x in line]) + '\n') |
| 329 |
fh.write('%endblock %s\n' % key)
|
| 330 |
|
| 331 |
unit = keys_with_units.get(inpify(key)) |
| 332 |
if unit is None: |
| 333 |
fh.write('%s %s\n' % (key, value))
|
| 334 |
else:
|
| 335 |
if 'fs**2' in unit: |
| 336 |
value /= fs**2
|
| 337 |
elif 'fs' in unit: |
| 338 |
value /= fs |
| 339 |
fh.write('%s %f %s\n' % (key, value, unit))
|
| 340 |
|
| 341 |
fh.write('#Definition of the unit cell\n')
|
| 342 |
fh.write('acell\n')
|
| 343 |
fh.write('%.14f %.14f %.14f Angstrom\n' % (1.0, 1.0, 1.0)) |
| 344 |
fh.write('rprim\n')
|
| 345 |
for v in self.cell: |
| 346 |
fh.write('%.14f %.14f %.14f\n' % tuple(v)) |
| 347 |
|
| 348 |
fh.write('chkprim 0 # Allow non-primitive cells\n')
|
| 349 |
|
| 350 |
fh.write('#Definition of the atom types\n')
|
| 351 |
fh.write('ntypat %d\n' % (len(self.species))) |
| 352 |
fh.write('znucl')
|
| 353 |
for n, Z in enumerate(self.species): |
| 354 |
fh.write(' %d' % (Z))
|
| 355 |
fh.write('\n')
|
| 356 |
fh.write('#Enumerate different atomic species\n')
|
| 357 |
fh.write('typat')
|
| 358 |
fh.write('\n')
|
| 359 |
self.types = []
|
| 360 |
for Z in self.numbers: |
| 361 |
for n, Zs in enumerate(self.species): |
| 362 |
if Z == Zs:
|
| 363 |
self.types.append(n+1) |
| 364 |
for n, type in enumerate(self.types): |
| 365 |
fh.write(' %d' % (type)) |
| 366 |
if n > 1 and ((n % self.n_entries_int) == 1): |
| 367 |
fh.write('\n')
|
| 368 |
fh.write('\n')
|
| 369 |
|
| 370 |
fh.write('#Definition of the atoms\n')
|
| 371 |
fh.write('xangst\n')
|
| 372 |
a = 0
|
| 373 |
for pos, Z in zip(self.positions, self.numbers): |
| 374 |
a += 1
|
| 375 |
fh.write('%.14f %.14f %.14f\n' % tuple(pos)) |
| 376 |
|
| 377 |
if self.kpts is not None: |
| 378 |
fh.write('kptopt %d\n' % (1)) |
| 379 |
fh.write('ngkpt ')
|
| 380 |
fh.write('%d %d %d\n' % tuple(self.kpts)) |
| 381 |
|
| 382 |
fh.write('#Definition of the SCF procedure\n')
|
| 383 |
fh.write('toldfe %.1g\n' % self.toldfe) |
| 384 |
fh.write('chkexit 1 # abinit.exit file in the running directory terminates after the current SCF\n')
|
| 385 |
|
| 386 |
fh.close() |
| 387 |
|
| 388 |
def read_fermi(self): |
| 389 |
"""Method that reads Fermi energy in Hartree from the output file
|
| 390 |
and returns it in eV"""
|
| 391 |
E_f=None
|
| 392 |
filename = self.label + '.txt' |
| 393 |
text = open(filename).read().lower()
|
| 394 |
assert 'error' not in text |
| 395 |
for line in iter(text.split('\n')): |
| 396 |
if line.rfind('fermi (or homo) energy (hartree) =') > -1: |
| 397 |
E_f = float(line.split('=')[1].strip().split()[0]) |
| 398 |
return E_f*Hartree
|
| 399 |
|
| 400 |
def read_kpts_info(self, kpt=0, spin=0, mode='eigenvalues'): |
| 401 |
""" Returns list of eigenvalues, occupations, kpts weights, or
|
| 402 |
kpts coordinates for given kpt and spin"""
|
| 403 |
# output may look like this (or without occupation entries); 8 entries per line:
|
| 404 |
#
|
| 405 |
# Eigenvalues (hartree) for nkpt= 20 k points:
|
| 406 |
# kpt# 1, nband= 3, wtk= 0.01563, kpt= 0.0625 0.0625 0.0625 (reduced coord)
|
| 407 |
# -0.09911 0.15393 0.15393
|
| 408 |
# occupation numbers for kpt# 1
|
| 409 |
# 2.00000 0.00000 0.00000
|
| 410 |
# kpt# 2, nband= 3, wtk= 0.04688, kpt= 0.1875 0.0625 0.0625 (reduced coord)
|
| 411 |
# ...
|
| 412 |
#
|
| 413 |
assert mode in ['eigenvalues' , 'occupations', 'ibz_k_points', 'k_point_weights'], 'mode not in [\'eigenvalues\' , \'occupations\', \'ibz_k_points\', \'k_point_weights\']' |
| 414 |
# number of lines of eigenvalues/occupations for a kpt
|
| 415 |
n_entry_lines = max(1, int(self.nbands/self.n_entries_float)) |
| 416 |
#
|
| 417 |
filename = self.label + '.txt' |
| 418 |
text = open(filename).read().lower()
|
| 419 |
assert 'error' not in text |
| 420 |
lines = iter(text.split('\n')) |
| 421 |
text_list = [] |
| 422 |
# find the begining line of eigenvalues
|
| 423 |
contains_eigenvalues = False
|
| 424 |
for line in lines: |
| 425 |
#if line.rfind('eigenvalues (hartree) for nkpt') > -1: #MDTMP
|
| 426 |
if line.rfind('eigenvalues ( ev ) for nkpt') > -1: |
| 427 |
n_kpts = int(line.split('nkpt=')[1].strip().split()[0]) |
| 428 |
contains_eigenvalues = True
|
| 429 |
break
|
| 430 |
# find the end line of eigenvalues starting from linenum
|
| 431 |
for line in lines: |
| 432 |
text_list.append(line) |
| 433 |
if not line.strip(): # find a blank line |
| 434 |
break
|
| 435 |
# remove last (blank) line
|
| 436 |
text_list = text_list[:-1]
|
| 437 |
|
| 438 |
assert contains_eigenvalues, 'No eigenvalues found in the output' |
| 439 |
|
| 440 |
# join text eigenvalues description with eigenvalues
|
| 441 |
# or occupation numbers for kpt# with occupations
|
| 442 |
contains_occupations = False
|
| 443 |
for line in text_list: |
| 444 |
if line.rfind('occupation numbers') > -1: |
| 445 |
contains_occupations = True
|
| 446 |
break
|
| 447 |
if mode == 'occupations': |
| 448 |
assert contains_occupations, 'No occupations found in the output' |
| 449 |
|
| 450 |
if contains_occupations:
|
| 451 |
range_kpts = 2*n_kpts
|
| 452 |
else:
|
| 453 |
range_kpts = n_kpts |
| 454 |
#
|
| 455 |
values_list = [] |
| 456 |
offset = 0
|
| 457 |
for kpt_entry in range(range_kpts): |
| 458 |
full_line = ''
|
| 459 |
for entry_line in range(n_entry_lines+1): |
| 460 |
full_line = full_line+str(text_list[offset+entry_line])
|
| 461 |
first_line = text_list[offset] |
| 462 |
if mode == 'occupations': |
| 463 |
if first_line.rfind('occupation numbers') > -1: |
| 464 |
# extract numbers
|
| 465 |
full_line = [float(v) for v in full_line.split('#')[1].strip().split()[1:]] |
| 466 |
values_list.append(full_line) |
| 467 |
elif mode in ['eigenvalues', 'ibz_k_points', 'k_point_weights']: |
| 468 |
if first_line.rfind('reduced coord') > -1: |
| 469 |
# extract numbers
|
| 470 |
if mode == 'eigenvalues': |
| 471 |
#full_line = [Hartree*float(v) for v in full_line.split(')')[1].strip().split()[:]] # MDTMP
|
| 472 |
full_line = [float(v) for v in full_line.split(')')[1].strip().split()[:]] |
| 473 |
elif mode == 'ibz_k_points': |
| 474 |
full_line = [float(v) for v in full_line.split('kpt=')[1].strip().split('(')[0].split()] |
| 475 |
else:
|
| 476 |
full_line = float(full_line.split('wtk=')[1].strip().split(',')[0].split()[0]) |
| 477 |
values_list.append(full_line) |
| 478 |
offset = offset+n_entry_lines+1
|
| 479 |
#
|
| 480 |
if mode in ['occupations', 'eigenvalues']: |
| 481 |
return np.array(values_list[kpt])
|
| 482 |
else:
|
| 483 |
return np.array(values_list)
|
| 484 |
|
| 485 |
def read(self): |
| 486 |
"""Read results from ABINIT's text-output file."""
|
| 487 |
filename = self.label + '.txt' |
| 488 |
text = open(filename).read().lower()
|
| 489 |
assert 'error' not in text |
| 490 |
assert 'was not enough scf cycles to converge' not in text |
| 491 |
# some consistency ckecks
|
| 492 |
for line in iter(text.split('\n')): |
| 493 |
if line.rfind('natom ') > -1: |
| 494 |
natom = int(line.split()[-1]) |
| 495 |
assert natom == len(self.numbers) |
| 496 |
for line in iter(text.split('\n')): |
| 497 |
if line.rfind('znucl ') > -1: |
| 498 |
znucl = [float(Z) for Z in line.split()[1:]] |
| 499 |
for n, Z in enumerate(self.species): |
| 500 |
assert Z == znucl[n]
|
| 501 |
lines = text.split('\n')
|
| 502 |
for n, line in enumerate(lines): |
| 503 |
if line.rfind(' typat ') > -1: |
| 504 |
nlines = len(self.numbers) / self.n_entries_int |
| 505 |
typat = [int(t) for t in line.split()[1:]] # first line |
| 506 |
for nline in range(nlines): # remaining lines |
| 507 |
for t in lines[1 + n + nline].split()[:]: |
| 508 |
typat.append(int(t))
|
| 509 |
for n, t in enumerate(self.types): |
| 510 |
assert t == typat[n]
|
| 511 |
|
| 512 |
lines = iter(text.split('\n')) |
| 513 |
# Stress:
|
| 514 |
# Printed in the output in the following format [Hartree/Bohr^3]:
|
| 515 |
# sigma(1 1)= 4.02063464E-04 sigma(3 2)= 0.00000000E+00
|
| 516 |
# sigma(2 2)= 4.02063464E-04 sigma(3 1)= 0.00000000E+00
|
| 517 |
# sigma(3 3)= 4.02063464E-04 sigma(2 1)= 0.00000000E+00
|
| 518 |
for line in lines: |
| 519 |
if line.rfind('cartesian components of stress tensor (hartree/bohr^3)') > -1: |
| 520 |
self.stress = np.empty((3, 3)) |
| 521 |
for i in range(3): |
| 522 |
entries = lines.next().split() |
| 523 |
self.stress[i,i] = float(entries[2]) |
| 524 |
self.stress[min(3, 4-i)-1, max(1, 2-i)-1] = float(entries[5]) |
| 525 |
self.stress[max(1, 2-i)-1, min(3, 4-i)-1] = float(entries[5]) |
| 526 |
self.stress = self.stress*Hartree/Bohr**3 |
| 527 |
break
|
| 528 |
else:
|
| 529 |
raise RuntimeError |
| 530 |
|
| 531 |
# Energy [Hartree]:
|
| 532 |
# Warning: Etotal could mean both electronic energy and free energy!
|
| 533 |
for line in iter(text.split('\n')): |
| 534 |
if line.rfind('>>>>> internal e=') > -1: |
| 535 |
self.etotal = float(line.split('=')[-1])*Hartree |
| 536 |
for line1 in iter(text.split('\n')): |
| 537 |
if line1.rfind('>>>>>>>>> etotal=') > -1: |
| 538 |
self.efree = float(line1.split('=')[-1])*Hartree |
| 539 |
break
|
| 540 |
else:
|
| 541 |
raise RuntimeError |
| 542 |
break
|
| 543 |
else:
|
| 544 |
for line2 in iter(text.split('\n')): |
| 545 |
if line2.rfind('>>>>>>>>> etotal=') > -1: |
| 546 |
self.etotal = float(line2.split('=')[-1])*Hartree |
| 547 |
self.efree = self.etotal |
| 548 |
break
|
| 549 |
else:
|
| 550 |
raise RuntimeError |
| 551 |
|
| 552 |
# Forces:
|
| 553 |
for line in lines: |
| 554 |
if line.rfind('cartesian forces (ev/angstrom) at end:') > -1: |
| 555 |
forces = [] |
| 556 |
for i in range(len(self.numbers)): |
| 557 |
forces.append(np.array([float(f) for f in lines.next().split()[1:]])) |
| 558 |
self.forces = np.array(forces)
|
| 559 |
break
|
| 560 |
else:
|
| 561 |
raise RuntimeError |
| 562 |
|
| 563 |
def inpify(key): |
| 564 |
return key.lower().replace('_', '').replace('.', '').replace('-', '') |
| 565 |
|
| 566 |
|
| 567 |
keys_with_units = {
|
| 568 |
} |
| 569 |
#keys_with_units = {
|
| 570 |
# 'paoenergyshift': 'eV',
|
| 571 |
# 'zmunitslength': 'Bohr',
|
| 572 |
# 'zmunitsangle': 'rad',
|
| 573 |
# 'zmforcetollength': 'eV/Ang',
|
| 574 |
# 'zmforcetolangle': 'eV/rad',
|
| 575 |
# 'zmmaxdispllength': 'Ang',
|
| 576 |
# 'zmmaxdisplangle': 'rad',
|
| 577 |
# 'ecut': 'eV',
|
| 578 |
# 'dmenergytolerance': 'eV',
|
| 579 |
# 'electronictemperature': 'eV',
|
| 580 |
# 'oneta': 'eV',
|
| 581 |
# 'onetaalpha': 'eV',
|
| 582 |
# 'onetabeta': 'eV',
|
| 583 |
# 'onrclwf': 'Ang',
|
| 584 |
# 'onchemicalpotentialrc': 'Ang',
|
| 585 |
# 'onchemicalpotentialtemperature': 'eV',
|
| 586 |
# 'mdmaxcgdispl': 'Ang',
|
| 587 |
# 'mdmaxforcetol': 'eV/Ang',
|
| 588 |
# 'mdmaxstresstol': 'eV/Ang**3',
|
| 589 |
# 'mdlengthtimestep': 'fs',
|
| 590 |
# 'mdinitialtemperature': 'eV',
|
| 591 |
# 'mdtargettemperature': 'eV',
|
| 592 |
# 'mdtargetpressure': 'eV/Ang**3',
|
| 593 |
# 'mdnosemass': 'eV*fs**2',
|
| 594 |
# 'mdparrinellorahmanmass': 'eV*fs**2',
|
| 595 |
# 'mdtaurelax': 'fs',
|
| 596 |
# 'mdbulkmodulus': 'eV/Ang**3',
|
| 597 |
# 'mdfcdispl': 'Ang',
|
| 598 |
# 'warningminimumatomicdistance': 'Ang',
|
| 599 |
# 'rcspatial': 'Ang',
|
| 600 |
# 'kgridcutoff': 'Ang',
|
| 601 |
# 'latticeconstant': 'Ang'}
|