Coverage for /builds/kinetik161/ase/ase/io/v_sim.py: 86.84%

1""" 

2This module contains functionality for reading and writing an ASE 

3Atoms object in V_Sim 3.5+ ascii format. 

4 

5""" 

6 

7import numpy as np 

8 

9from ase.utils import reader, writer 

10 

11 

12@reader 

13def read_v_sim(fd): 

14 """Import V_Sim input file. 

15 

16 Reads cell, atom positions, etc. from v_sim ascii file. 

17 V_sim format is specified here: 

18 https://l_sim.gitlab.io/v_sim/sample.html#sample_ascii 

19 """ 

20 

21 import re 

22 

23 from ase import Atoms, units 

24 from ase.geometry import cellpar_to_cell 

25 

26 # Read comment: 

27 fd.readline() 

28 

29 line = fd.readline() + ' ' + fd.readline() 

30 box = line.split() 

31 for i in range(len(box)): 

32 box[i] = float(box[i]) 

33 

34 keywords = [] 

35 positions = [] 

36 symbols = [] 

37 unit = 1.0 

38 

39 re_comment = re.compile(r'^\s*[#!]') 

40 re_node = re.compile(r'^\s*\S+\s+\S+\s+\S+\s+\S+') 

41 

42 while True: 

43 line = fd.readline() 

44 

45 if line == '': 

46 break # EOF 

47 

48 p = re_comment.match(line) 

49 if p is not None: 

50 # remove comment character at the beginning of line 

51 line = line[p.end():].replace(',', ' ').lower() 

52 if line[:8] == "keyword:": 

53 keywords.extend(line[8:].split()) 

54 

55 elif re_node.match(line): 

56 unit = 1.0 

57 if "reduced" not in keywords: 

58 if (("bohr" in keywords) or ("bohrd0" in keywords) or 

59 ("atomic" in keywords) or ("atomicd0" in keywords)): 

60 unit = units.Bohr 

61 

62 fields = line.split() 

63 positions.append([unit * float(fields[0]), 

64 unit * float(fields[1]), 

65 unit * float(fields[2])]) 

66 symbols.append(fields[3]) 

67 

68 # create atoms object based on the information 

69 if "angdeg" in keywords: 

70 cell = cellpar_to_cell(box) 

71 else: 

72 unit = 1.0 

73 if (("bohr" in keywords) or ("bohrd0" in keywords) or 

74 ("atomic" in keywords) or ("atomicd0" in keywords)): 

75 unit = units.Bohr 

76 cell = np.zeros((3, 3)) 

77 cell.flat[[0, 3, 4, 6, 7, 8]] = box[:6] 

78 cell *= unit 

79 

80 if "reduced" in keywords: 

81 atoms = Atoms(cell=cell, scaled_positions=positions) 

82 else: 

83 atoms = Atoms(cell=cell, positions=positions) 

84 

85 if "periodic" in keywords: 

86 atoms.pbc = [True, True, True] 

87 elif "freebc" in keywords: 

88 atoms.pbc = [False, False, False] 

89 elif "surface" in keywords: 

90 atoms.pbc = [True, False, True] 

91 else: # default is periodic boundary conditions 

92 atoms.pbc = [True, True, True] 

93 

94 atoms.set_chemical_symbols(symbols) 

95 return atoms 

96 

97 

98@writer 

99def write_v_sim(fd, atoms): 

100 """Write V_Sim input file. 

101 

102 Writes the atom positions and unit cell. 

103 V_sim format is specified here: 

104 https://l_sim.gitlab.io/v_sim/sample.html#sample_ascii 

105 """ 

106 from ase.geometry import cell_to_cellpar, cellpar_to_cell 

107 

108 # Convert the lattice vectors to triangular matrix by converting 

109 # to and from a set of lengths and angles 

110 cell = cellpar_to_cell(cell_to_cellpar(atoms.cell)) 

111 dxx = cell[0, 0] 

112 dyx, dyy = cell[1, 0:2] 

113 dzx, dzy, dzz = cell[2, 0:3] 

114 

115 fd.write('===== v_sim input file created using the' 

116 ' Atomic Simulation Environment (ASE) ====\n') 

117 fd.write(f'{dxx} {dyx} {dyy}\n') 

118 fd.write(f'{dzx} {dzy} {dzz}\n') 

119 

120 # Use v_sim 3.5 keywords to indicate scaled positions, etc. 

121 fd.write('#keyword: reduced\n') 

122 fd.write('#keyword: angstroem\n') 

123 if np.all(atoms.pbc): 

124 fd.write('#keyword: periodic\n') 

125 elif not np.any(atoms.pbc): 

126 fd.write('#keyword: freeBC\n') 

127 elif np.array_equiv(atoms.pbc, [True, False, True]): 

128 fd.write('#keyword: surface\n') 

129 else: 

130 raise Exception( 

131 'Only supported boundary conditions are full PBC,' 

132 ' no periodic boundary, and surface which is free in y direction' 

133 ' (i.e. Atoms.pbc = [True, False, True]).') 

134 

135 # Add atoms (scaled positions) 

136 for position, symbol in zip(atoms.get_scaled_positions(), 

137 atoms.get_chemical_symbols()): 

138 fd.write('{} {} {} {}\n'.format( 

139 position[0], position[1], position[2], symbol))