Coverage for /builds/hweiske/ase/ase/io/gpw.py: 11.86%

1"""Read gpw-file from GPAW.""" 

2import ase.io.ulm as ulm 

3from ase import Atoms 

4from ase.calculators.singlepoint import (SinglePointDFTCalculator, 

5 SinglePointKPoint, all_properties) 

6from ase.io.trajectory import read_atoms 

7from ase.units import Bohr, Hartree 

8 

9 

10def read_gpw(filename): 

11 try: 

12 reader = ulm.open(filename) 

13 except ulm.InvalidULMFileError: 

14 return read_old_gpw(filename) 

15 

16 atoms = read_atoms(reader.atoms, _try_except=False) 

17 

18 wfs = reader.wave_functions 

19 kpts = wfs.get('kpts') 

20 if kpts is None: 

21 ibzkpts = None 

22 bzkpts = None 

23 bz2ibz = None 

24 else: 

25 ibzkpts = kpts.ibzkpts 

26 bzkpts = kpts.get('bzkpts') 

27 bz2ibz = kpts.get('bz2ibz') 

28 

29 if reader.version >= 3: 

30 efermi = reader.wave_functions.fermi_levels.mean() 

31 else: 

32 efermi = reader.occupations.fermilevel 

33 

34 atoms.calc = SinglePointDFTCalculator( 

35 atoms, 

36 efermi=efermi, 

37 ibzkpts=ibzkpts, 

38 bzkpts=bzkpts, 

39 bz2ibz=bz2ibz, 

40 # New gpw-files may have "non_collinear_magmom(s)" which ASE 

41 # doesn't know: 

42 **{property: value 

43 for property, value in reader.results.asdict().items() 

44 if property in all_properties}) 

45 

46 if kpts is not None: 

47 atoms.calc.kpts = [] 

48 for spin, (eps_kn, f_kn) in enumerate(zip(wfs.eigenvalues, 

49 wfs.occupations)): 

50 for kpt, (weight, eps_n, f_n) in enumerate(zip(kpts.weights, 

51 eps_kn, f_kn)): 

52 atoms.calc.kpts.append( 

53 SinglePointKPoint(weight, spin, kpt, eps_n, f_n)) 

54 reader.close() 

55 

56 return atoms 

57 

58 

59def read_old_gpw(filename): 

60 from gpaw.io.tar import Reader 

61 r = Reader(filename) 

62 positions = r.get('CartesianPositions') * Bohr 

63 numbers = r.get('AtomicNumbers') 

64 cell = r.get('UnitCell') * Bohr 

65 pbc = r.get('BoundaryConditions') 

66 tags = r.get('Tags') 

67 magmoms = r.get('MagneticMoments') 

68 energy = r.get('PotentialEnergy') * Hartree 

69 

70 if r.has_array('CartesianForces'): 

71 forces = r.get('CartesianForces') * Hartree / Bohr 

72 else: 

73 forces = None 

74 

75 atoms = Atoms(positions=positions, 

76 numbers=numbers, 

77 cell=cell, 

78 pbc=pbc) 

79 if tags.any(): 

80 atoms.set_tags(tags) 

81 

82 if magmoms.any(): 

83 atoms.set_initial_magnetic_moments(magmoms) 

84 magmom = magmoms.sum() 

85 else: 

86 magmoms = None 

87 magmom = None 

88 

89 atoms.calc = SinglePointDFTCalculator(atoms, energy=energy, 

90 forces=forces, 

91 magmoms=magmoms, 

92 magmom=magmom) 

93 kpts = [] 

94 if r.has_array('IBZKPoints'): 

95 for w, kpt, eps_n, f_n in zip(r.get('IBZKPointWeights'), 

96 r.get('IBZKPoints'), 

97 r.get('Eigenvalues'), 

98 r.get('OccupationNumbers')): 

99 kpts.append(SinglePointKPoint(w, kpt[0], kpt[1], 

100 eps_n[0], f_n[0])) 

101 atoms.calc.kpts = kpts 

102 

103 return atoms