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 

8from ase.utils import reader, writer 

9 

10 

11@reader 

12def read_v_sim(fd): 

13 """Import V_Sim input file. 

14 

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

16 V_sim format is specified here: 

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

18 """ 

19 

20 from ase import Atoms, units 

21 from ase.geometry import cellpar_to_cell 

22 import re 

23 

24 # Read comment: 

25 fd.readline() 

26 

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

28 box = line.split() 

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

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

31 

32 keywords = [] 

33 positions = [] 

34 symbols = [] 

35 unit = 1.0 

36 

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

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

39 

40 while True: 

41 line = fd.readline() 

42 

43 if line == '': 

44 break # EOF 

45 

46 p = re_comment.match(line) 

47 if p is not None: 

48 # remove comment character at the beginning of line 

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

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

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

52 

53 elif re_node.match(line): 

54 unit = 1.0 

55 if not ("reduced" in keywords): 

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

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

58 unit = units.Bohr 

59 

60 fields = line.split() 

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

62 unit * float(fields[1]), 

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

64 symbols.append(fields[3]) 

65 

66 # create atoms object based on the information 

67 if "angdeg" in keywords: 

68 cell = cellpar_to_cell(box) 

69 else: 

70 unit = 1.0 

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

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

73 unit = units.Bohr 

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

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

76 cell *= unit 

77 

78 if "reduced" in keywords: 

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

80 else: 

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

82 

83 if "periodic" in keywords: 

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

85 elif "freebc" in keywords: 

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

87 elif "surface" in keywords: 

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

89 else: # default is periodic boundary conditions 

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

91 

92 atoms.set_chemical_symbols(symbols) 

93 return atoms 

94 

95 

96@writer 

97def write_v_sim(fd, atoms): 

98 """Write V_Sim input file. 

99 

100 Writes the atom positions and unit cell. 

101 V_sim format is specified here: 

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

103 """ 

104 from ase.geometry import cellpar_to_cell, cell_to_cellpar 

105 

106 # Convert the lattice vectors to triangular matrix by converting 

107 # to and from a set of lengths and angles 

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

109 dxx = cell[0, 0] 

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

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

112 

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

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

115 fd.write('{0} {1} {2}\n'.format(dxx, dyx, dyy)) 

116 fd.write('{0} {1} {2}\n'.format(dzx, dzy, dzz)) 

117 

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

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

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

121 if np.alltrue(atoms.pbc): 

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

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

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

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

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

127 else: 

128 raise Exception( 

129 'Only supported boundary conditions are full PBC,' 

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

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

132 

133 # Add atoms (scaled positions) 

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

135 atoms.get_chemical_symbols()): 

136 fd.write('{0} {1} {2} {3}\n'.format( 

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