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 """ 

17 

18 from ase import Atoms, units 

19 from ase.geometry import cellpar_to_cell 

20 import re 

21 

22 # Read comment: 

23 fd.readline() 

24 

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

26 box = line.split() 

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

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

29 

30 keywords = [] 

31 positions = [] 

32 symbols = [] 

33 unit = 1.0 

34 

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

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

37 

38 while True: 

39 line = fd.readline() 

40 

41 if line == '': 

42 break # EOF 

43 

44 p = re_comment.match(line) 

45 if p is not None: 

46 # remove comment character at the beginning of line 

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

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

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

50 

51 elif re_node.match(line): 

52 unit = 1.0 

53 if not ("reduced" in keywords): 

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

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

56 unit = units.Bohr 

57 

58 fields = line.split() 

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

60 unit * float(fields[1]), 

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

62 symbols.append(fields[3]) 

63 

64 if ("surface" in keywords) or ("freeBC" in keywords): 

65 raise NotImplementedError 

66 

67 # create atoms object based on the information 

68 if "angdeg" in keywords: 

69 cell = cellpar_to_cell(box) 

70 else: 

71 unit = 1.0 

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

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

74 unit = units.Bohr 

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

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

77 cell *= unit 

78 

79 if "reduced" in keywords: 

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

81 else: 

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

83 

84 atoms.set_chemical_symbols(symbols) 

85 return atoms 

86 

87 

88@writer 

89def write_v_sim(fd, atoms): 

90 """Write V_Sim input file. 

91 

92 Writes the atom positions and unit cell. 

93 """ 

94 from ase.geometry import cellpar_to_cell, cell_to_cellpar 

95 

96 # Convert the lattice vectors to triangular matrix by converting 

97 # to and from a set of lengths and angles 

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

99 dxx = cell[0, 0] 

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

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

102 

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

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

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

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

107 

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

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

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

111 if np.alltrue(atoms.pbc): 

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

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

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

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

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

117 else: 

118 raise Exception( 

119 'Only supported boundary conditions are full PBC,' 

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

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

122 

123 # Add atoms (scaled positions) 

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

125 atoms.get_chemical_symbols()): 

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

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