1""" write gromos96 geometry files 

2(the exact file format is copied from the freely available 

3gromacs package, http://www.gromacs.org 

4its procedure src/gmxlib/confio.c (write_g96_conf) 

5""" 

6 

7import numpy as np 

8 

9from ase import Atoms 

10from ase.data import chemical_symbols 

11from ase.utils import reader, writer 

12 

13 

14@reader 

15def read_gromos(fileobj): 

16 """Read gromos geometry files (.g96). 

17 Reads: 

18 atom positions, 

19 and simulation cell (if present) 

20 tries to set atom types 

21 """ 

22 

23 lines = fileobj.readlines() 

24 read_pos = False 

25 read_box = False 

26 tmp_pos = [] 

27 symbols = [] 

28 mycell = None 

29 for line in lines: 

30 if (read_pos and ('END' in line)): 

31 read_pos = False 

32 if (read_box and ('END' in line)): 

33 read_box = False 

34 if read_pos: 

35 symbol, dummy, x, y, z = line.split()[2:7] 

36 tmp_pos.append((10 * float(x), 10 * float(y), 10 * float(z))) 

37 if (len(symbol) != 2): 

38 symbols.append(symbol[0].lower().capitalize()) 

39 else: 

40 symbol2 = symbol[0].lower().capitalize() + \ 

41 symbol[1] 

42 if symbol2 in chemical_symbols: 

43 symbols.append(symbol2) 

44 else: 

45 symbols.append(symbol[0].lower().capitalize()) 

46 if symbols[-1] not in chemical_symbols: 

47 raise RuntimeError("Symbol '{}' not in chemical symbols" 

48 .format(symbols[-1])) 

49 if read_box: 

50 try: 

51 grocell = list(map(float, line.split())) 

52 except ValueError: 

53 pass 

54 else: 

55 mycell = np.diag(grocell[:3]) 

56 if len(grocell) >= 9: 

57 mycell.flat[[1, 2, 3, 5, 6, 7]] = grocell[3:9] 

58 mycell *= 10. 

59 if ('POSITION' in line): 

60 read_pos = True 

61 if ('BOX' in line): 

62 read_box = True 

63 

64 gmx_system = Atoms(symbols=symbols, positions=tmp_pos, cell=mycell) 

65 if mycell is not None: 

66 gmx_system.pbc = True 

67 return gmx_system 

68 

69 

70@writer 

71def write_gromos(fileobj, atoms): 

72 """Write gromos geometry files (.g96). 

73 Writes: 

74 atom positions, 

75 and simulation cell (if present) 

76 """ 

77 

78 from ase import units 

79 

80 natoms = len(atoms) 

81 try: 

82 gromos_residuenames = atoms.get_array('residuenames') 

83 except KeyError: 

84 gromos_residuenames = [] 

85 for idum in range(natoms): 

86 gromos_residuenames.append('1DUM') 

87 try: 

88 gromos_atomtypes = atoms.get_array('atomtypes') 

89 except KeyError: 

90 gromos_atomtypes = atoms.get_chemical_symbols() 

91 

92 pos = atoms.get_positions() 

93 pos = pos / 10.0 

94 

95 vel = atoms.get_velocities() 

96 if vel is None: 

97 vel = pos * 0.0 

98 else: 

99 vel *= 1000.0 * units.fs / units.nm 

100 

101 fileobj.write('TITLE\n') 

102 fileobj.write('Gromos96 structure file written by ASE \n') 

103 fileobj.write('END\n') 

104 fileobj.write('POSITION\n') 

105 count = 1 

106 rescount = 0 

107 oldresname = '' 

108 for resname, atomtype, xyz in zip(gromos_residuenames, 

109 gromos_atomtypes, 

110 pos): 

111 if resname != oldresname: 

112 oldresname = resname 

113 rescount = rescount + 1 

114 okresname = resname.lstrip('0123456789 ') 

115 fileobj.write('%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n' % 

116 (rescount, okresname, atomtype, count, 

117 xyz[0], xyz[1], xyz[2])) 

118 count = count + 1 

119 fileobj.write('END\n') 

120 

121 if atoms.get_pbc().any(): 

122 fileobj.write('BOX\n') 

123 mycell = atoms.get_cell() 

124 grocell = mycell.flat[[0, 4, 8, 1, 2, 3, 5, 6, 7]] * 0.1 

125 fileobj.write(''.join(['{:15.9f}'.format(x) for x in grocell])) 

126 fileobj.write('\nEND\n')