1from ase.atoms import Atoms 

2from ase.utils import reader, writer 

3 

4 

5@writer 

6def write_crystal(fd, atoms): 

7 """Method to write atom structure in crystal format 

8 (fort.34 format) 

9 """ 

10 

11 ispbc = atoms.get_pbc() 

12 box = atoms.get_cell() 

13 

14 # here it is assumed that the non-periodic direction are z 

15 # in 2D case, z and y in the 1D case. 

16 

17 if ispbc[2]: 

18 fd.write('%2s %2s %2s %23s \n' % 

19 ('3', '1', '1', 'E -0.0E+0 DE 0.0E+0( 1)')) 

20 elif ispbc[1]: 

21 fd.write('%2s %2s %2s %23s \n' % 

22 ('2', '1', '1', 'E -0.0E+0 DE 0.0E+0( 1)')) 

23 box[2, 2] = 500. 

24 elif ispbc[0]: 

25 fd.write('%2s %2s %2s %23s \n' % 

26 ('1', '1', '1', 'E -0.0E+0 DE 0.0E+0( 1)')) 

27 box[2, 2] = 500. 

28 box[1, 1] = 500. 

29 else: 

30 fd.write('%2s %2s %2s %23s \n' % 

31 ('0', '1', '1', 'E -0.0E+0 DE 0.0E+0( 1)')) 

32 box[2, 2] = 500. 

33 box[1, 1] = 500. 

34 box[0, 0] = 500. 

35 

36 # write box 

37 # crystal dummy 

38 fd.write(' %.17E %.17E %.17E \n' 

39 % (box[0][0], box[0][1], box[0][2])) 

40 fd.write(' %.17E %.17E %.17E \n' 

41 % (box[1][0], box[1][1], box[1][2])) 

42 fd.write(' %.17E %.17E %.17E \n' 

43 % (box[2][0], box[2][1], box[2][2])) 

44 

45 # write symmetry operations (not implemented yet for 

46 # higher symmetries than C1) 

47 fd.write(' %2s \n' % (1)) 

48 fd.write(' %.17E %.17E %.17E \n' % (1, 0, 0)) 

49 fd.write(' %.17E %.17E %.17E \n' % (0, 1, 0)) 

50 fd.write(' %.17E %.17E %.17E \n' % (0, 0, 1)) 

51 fd.write(' %.17E %.17E %.17E \n' % (0, 0, 0)) 

52 

53 # write coordinates 

54 fd.write(' %8s \n' % (len(atoms))) 

55 coords = atoms.get_positions() 

56 tags = atoms.get_tags() 

57 atomnum = atoms.get_atomic_numbers() 

58 for iatom, coord in enumerate(coords): 

59 fd.write('%5i %19.16f %19.16f %19.16f \n' 

60 % (atomnum[iatom] + tags[iatom], 

61 coords[iatom][0], coords[iatom][1], coords[iatom][2])) 

62 

63 

64@reader 

65def read_crystal(fd): 

66 """Method to read coordinates form 'fort.34' files 

67 additionally read information about 

68 periodic boundary condition 

69 """ 

70 lines = fd.readlines() 

71 

72 atoms_pos = [] 

73 anumber_list = [] 

74 my_pbc = [False, False, False] 

75 mycell = [] 

76 

77 if float(lines[4]) != 1: 

78 raise ValueError('High symmetry geometry is not allowed.') 

79 

80 if float(lines[1].split()[0]) < 500.0: 

81 cell = [float(c) for c in lines[1].split()] 

82 mycell.append(cell) 

83 my_pbc[0] = True 

84 else: 

85 mycell.append([1, 0, 0]) 

86 

87 if float(lines[2].split()[1]) < 500.0: 

88 cell = [float(c) for c in lines[2].split()] 

89 mycell.append(cell) 

90 my_pbc[1] = True 

91 else: 

92 mycell.append([0, 1, 0]) 

93 

94 if float(lines[3].split()[2]) < 500.0: 

95 cell = [float(c) for c in lines[3].split()] 

96 mycell.append(cell) 

97 my_pbc[2] = True 

98 else: 

99 mycell.append([0, 0, 1]) 

100 

101 natoms = int(lines[9].split()[0]) 

102 for i in range(natoms): 

103 index = 10 + i 

104 anum = int(lines[index].split()[0]) % 100 

105 anumber_list.append(anum) 

106 

107 position = [float(p) for p in lines[index].split()[1:]] 

108 atoms_pos.append(position) 

109 

110 atoms = Atoms(positions=atoms_pos, numbers=anumber_list, 

111 cell=mycell, pbc=my_pbc) 

112 

113 return atoms