Coverage for /builds/ase/ase/ase/io/gromos.py: 95.35%

1# fmt: off 

2 

3""" write gromos96 geometry files 

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

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

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

7""" 

8 

9import numpy as np 

10 

11from ase import Atoms 

12from ase.data import chemical_symbols 

13from ase.utils import reader, writer 

14 

15 

16@reader 

17def read_gromos(fileobj): 

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

19 Reads: 

20 atom positions, 

21 and simulation cell (if present) 

22 tries to set atom types 

23 """ 

24 

25 lines = fileobj.readlines() 

26 read_pos = False 

27 read_box = False 

28 tmp_pos = [] 

29 symbols = [] 

30 mycell = None 

31 for line in lines: 

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

33 read_pos = False 

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

35 read_box = False 

36 if read_pos: 

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

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

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

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

41 else: 

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

43 symbol[1] 

44 if symbol2 in chemical_symbols: 

45 symbols.append(symbol2) 

46 else: 

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

48 if symbols[-1] not in chemical_symbols: 

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

50 .format(symbols[-1])) 

51 if read_box: 

52 try: 

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

54 except ValueError: 

55 pass 

56 else: 

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

58 if len(grocell) >= 9: 

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

60 mycell *= 10. 

61 if ('POSITION' in line): 

62 read_pos = True 

63 if ('BOX' in line): 

64 read_box = True 

65 

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

67 if mycell is not None: 

68 gmx_system.pbc = True 

69 return gmx_system 

70 

71 

72@writer 

73def write_gromos(fileobj, atoms): 

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

75 Writes: 

76 atom positions, 

77 and simulation cell (if present) 

78 """ 

79 

80 from ase import units 

81 

82 natoms = len(atoms) 

83 try: 

84 gromos_residuenames = atoms.get_array('residuenames') 

85 except KeyError: 

86 gromos_residuenames = [] 

87 for _ in range(natoms): 

88 gromos_residuenames.append('1DUM') 

89 try: 

90 gromos_atomtypes = atoms.get_array('atomtypes') 

91 except KeyError: 

92 gromos_atomtypes = atoms.get_chemical_symbols() 

93 

94 pos = atoms.get_positions() 

95 pos = pos / 10.0 

96 

97 vel = atoms.get_velocities() 

98 if vel is None: 

99 vel = pos * 0.0 

100 else: 

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

102 

103 fileobj.write('TITLE\n') 

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

105 fileobj.write('END\n') 

106 fileobj.write('POSITION\n') 

107 count = 1 

108 rescount = 0 

109 oldresname = '' 

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

111 gromos_atomtypes, 

112 pos): 

113 if resname != oldresname: 

114 oldresname = resname 

115 rescount = rescount + 1 

116 okresname = resname.lstrip('0123456789 ') 

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

118 (rescount, okresname, atomtype, count, 

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

120 count = count + 1 

121 fileobj.write('END\n') 

122 

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

124 fileobj.write('BOX\n') 

125 mycell = atoms.get_cell() 

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

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

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