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

1# fmt: off 

2 

3""" 

4This module contains functionality for reading and writing an ASE 

5Atoms object in V_Sim 3.5+ ascii format. 

6 

7""" 

8 

9import numpy as np 

10 

11from ase.utils import reader, writer 

12 

13 

14@reader 

15def read_v_sim(fd): 

16 """Import V_Sim input file. 

17 

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

19 V_sim format is specified here: 

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

21 """ 

22 

23 import re 

24 

25 from ase import Atoms, units 

26 from ase.geometry import cellpar_to_cell 

27 

28 # Read comment: 

29 fd.readline() 

30 

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

32 box = line.split() 

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

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

35 

36 keywords = [] 

37 positions = [] 

38 symbols = [] 

39 unit = 1.0 

40 

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

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

43 

44 while True: 

45 line = fd.readline() 

46 

47 if line == '': 

48 break # EOF 

49 

50 p = re_comment.match(line) 

51 if p is not None: 

52 # remove comment character at the beginning of line 

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

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

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

56 

57 elif re_node.match(line): 

58 unit = 1.0 

59 if "reduced" not in keywords: 

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

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

62 unit = units.Bohr 

63 

64 fields = line.split() 

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

66 unit * float(fields[1]), 

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

68 symbols.append(fields[3]) 

69 

70 # create atoms object based on the information 

71 if "angdeg" in keywords: 

72 cell = cellpar_to_cell(box) 

73 else: 

74 unit = 1.0 

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

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

77 unit = units.Bohr 

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

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

80 cell *= unit 

81 

82 if "reduced" in keywords: 

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

84 else: 

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

86 

87 if "periodic" in keywords: 

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

89 elif "freebc" in keywords: 

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

91 elif "surface" in keywords: 

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

93 else: # default is periodic boundary conditions 

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

95 

96 atoms.set_chemical_symbols(symbols) 

97 return atoms 

98 

99 

100@writer 

101def write_v_sim(fd, atoms): 

102 """Write V_Sim input file. 

103 

104 Writes the atom positions and unit cell. 

105 V_sim format is specified here: 

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

107 """ 

108 from ase.geometry import cell_to_cellpar, cellpar_to_cell 

109 

110 # Convert the lattice vectors to triangular matrix by converting 

111 # to and from a set of lengths and angles 

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

113 dxx = cell[0, 0] 

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

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

116 

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

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

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

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

121 

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

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

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

125 if np.all(atoms.pbc): 

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

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

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

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

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

131 else: 

132 raise Exception( 

133 'Only supported boundary conditions are full PBC,' 

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

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

136 

137 # Add atoms (scaled positions) 

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

139 atoms.get_chemical_symbols()): 

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

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