1"""This module defines an ASE interface to GROMACS. 

2 

3http://www.gromacs.org/ 

4It is VERY SLOW compared to standard Gromacs 

5(due to slow formatted io required here). 

6 

7Mainly intended to be the MM part in the ase QM/MM 

8 

9Markus.Kaukonen@iki.fi 

10 

11To be done: 

121) change the documentation for the new file-io-calculator (test works now) 

132) change gromacs program names 

14-now: hard coded 

15-future: set as dictionary in params_runs 

16 

17""" 

18 

19import os 

20import subprocess 

21from glob import glob 

22from shutil import which 

23 

24import numpy as np 

25 

26from ase import units 

27from ase.calculators.calculator import (EnvironmentError, 

28 FileIOCalculator, 

29 all_changes) 

30from ase.io.gromos import read_gromos, write_gromos 

31 

32 

33def parse_gromacs_version(output): 

34 import re 

35 match = re.search(r'GROMACS version\:\s*(\S+)', output, re.M) 

36 return match.group(1) 

37 

38 

39def get_gromacs_version(executable): 

40 output = subprocess.check_output([executable, '--version'], 

41 encoding='utf-8') 

42 return parse_gromacs_version(output) 

43 

44 

45def do_clean(name='#*'): 

46 """ remove files matching wildcards """ 

47 myfiles = glob(name) 

48 for myfile in myfiles: 

49 try: 

50 os.remove(myfile) 

51 except OSError: 

52 pass 

53 

54 

55class Gromacs(FileIOCalculator): 

56 """Class for doing GROMACS calculations. 

57 Before running a gromacs calculation you must prepare the input files 

58 separately (pdb2gmx and grompp for instance.) 

59 

60 Input parameters for gromacs runs (the .mdp file) 

61 are given in self.params and can be set when initializing the calculator 

62 or by method set_own. 

63 for example:: 

64 

65 CALC_MM_RELAX = Gromacs() 

66 CALC_MM_RELAX.set_own_params('integrator', 'steep', 

67 'use steepest descent') 

68 

69 Run command line arguments for gromacs related programs: 

70 pdb2gmx, grompp, mdrun, energy, traj. These can be given as:: 

71 

72 CALC_MM_RELAX = Gromacs() 

73 CALC_MM_RELAX.set_own_params_runs('force_field', 'oplsaa') 

74 """ 

75 

76 implemented_properties = ['energy', 'forces'] 

77 discard_results_on_any_change = True 

78 

79 default_parameters = dict( 

80 define='-DFLEXIBLE', 

81 integrator='cg', 

82 nsteps='10000', 

83 nstfout='10', 

84 nstlog='10', 

85 nstenergy='10', 

86 nstlist='10', 

87 ns_type='grid', 

88 pbc='xyz', 

89 rlist='1.15', 

90 coulombtype='PME-Switch', 

91 rcoulomb='0.8', 

92 vdwtype='shift', 

93 rvdw='0.8', 

94 rvdw_switch='0.75', 

95 DispCorr='Ener') 

96 

97 def __init__(self, restart=None, 

98 ignore_bad_restart_file=FileIOCalculator._deprecated, 

99 label='gromacs', atoms=None, 

100 do_qmmm=False, clean=True, 

101 water_model='tip3p', force_field='oplsaa', command=None, 

102 **kwargs): 

103 """Construct GROMACS-calculator object. 

104 

105 Parameters 

106 ========== 

107 label: str 

108 Prefix to use for filenames (label.in, label.txt, ...). 

109 Default is 'gromacs'. 

110 

111 do_qmmm : bool 

112 Is gromacs used as mm calculator for a qm/mm calculation 

113 

114 clean : bool 

115 Remove gromacs backup files 

116 and old gormacs.* files 

117 

118 water_model: str 

119 Water model to be used in gromacs runs (see gromacs manual) 

120 

121 force_field: str 

122 Force field to be used in gromacs runs 

123 

124 command : str 

125 Gromacs executable; if None (default), choose available one from 

126 ('gmx', 'gmx_d', 'gmx_mpi', 'gmx_mpi_d') 

127 """ 

128 

129 gmxes = ('gmx', 'gmx_d', 'gmx_mpi', 'gmx_mpi_d') 

130 if command is not None: 

131 self.command = command 

132 else: 

133 for command in gmxes: 

134 if which(command): 

135 self.command = command 

136 break 

137 else: 

138 self.command = None 

139 self.missing_gmx = 'missing gromacs executable {}'.format(gmxes) 

140 

141 self.do_qmmm = do_qmmm 

142 self.water_model = water_model 

143 self.force_field = force_field 

144 self.clean = clean 

145 self.params_doc = {} 

146 # add comments for gromacs input file 

147 self.params_doc['define'] = \ 

148 'flexible/ rigid water' 

149 self.params_doc['integrator'] = \ 

150 'md: molecular dynamics(Leapfrog), \n' + \ 

151 '; md-vv: molecular dynamics(Velocity Verlet), \n' + \ 

152 '; steep: steepest descent minimization, \n' + \ 

153 '; cg: conjugate cradient minimization \n' 

154 

155 self.positions = None 

156 self.atoms = None 

157 

158 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

159 label, atoms, **kwargs) 

160 self.set(**kwargs) 

161 # default values for runtime parameters 

162 # can be changed by self.set_own_params_runs('key', 'value') 

163 self.params_runs = {} 

164 self.params_runs['index_filename'] = 'index.ndx' 

165 self.params_runs['init_structure'] = self.label + '.pdb' 

166 self.params_runs['water'] = self.water_model 

167 self.params_runs['force_field'] = self.force_field 

168 

169 # these below are required by qm/mm 

170 self.topology_filename = self.label + '.top' 

171 

172 # clean up gromacs backups 

173 if self.clean: 

174 do_clean('gromacs.???') 

175 

176 # write input files for gromacs program energy 

177 self.write_energy_files() 

178 

179 if self.do_qmmm: 

180 self.parameters['integrator'] = 'md' 

181 self.parameters['nsteps'] = '0' 

182 

183 def _get_name(self): 

184 return 'Gromacs' 

185 

186 def _execute_gromacs(self, command): 

187 """ execute gmx command 

188 Parameters 

189 ---------- 

190 command : str 

191 """ 

192 if self.command: 

193 subprocess.check_call(self.command + ' ' + command, shell=True) 

194 else: 

195 raise EnvironmentError(self.missing_gmx) 

196 

197 def generate_g96file(self): 

198 """ from current coordinates (self.structure_file) 

199 write a structure file in .g96 format 

200 """ 

201 # generate structure file in g96 format 

202 write_gromos(self.label + '.g96', self.atoms) 

203 

204 def run_editconf(self): 

205 """ run gromacs program editconf, typically to set a simulation box 

206 writing to the input structure""" 

207 subcmd = 'editconf' 

208 command = ' '.join([ 

209 subcmd, 

210 '-f', self.label + '.g96', 

211 '-o', self.label + '.g96', 

212 self.params_runs.get('extra_editconf_parameters', ''), 

213 '> {}.{}.log 2>&1'.format(self.label, subcmd)]) 

214 self._execute_gromacs(command) 

215 

216 def run_genbox(self): 

217 """Run gromacs program genbox, typically to solvate the system 

218 writing to the input structure 

219 as extra parameter you need to define the file containing the solvent 

220 

221 for instance:: 

222 

223 CALC_MM_RELAX = Gromacs() 

224 CALC_MM_RELAX.set_own_params_runs( 

225 'extra_genbox_parameters', '-cs spc216.gro') 

226 """ 

227 subcmd = 'genbox' 

228 command = ' '.join([ 

229 subcmd, 

230 '-cp', self.label + '.g96', 

231 '-o', self.label + '.g96', 

232 '-p', self.label + '.top', 

233 self.params_runs.get('extra_genbox_parameters', ''), 

234 '> {}.{}.log 2>&1'.format(self.label, subcmd)]) 

235 self._execute_gromacs(command) 

236 

237 def run(self): 

238 """ runs a gromacs-mdrun with the 

239 current atom-configuration """ 

240 

241 # clean up gromacs backups 

242 if self.clean: 

243 do_clean('#*') 

244 

245 subcmd = 'mdrun' 

246 command = [subcmd] 

247 if self.do_qmmm: 

248 command += [ 

249 '-s', self.label + '.tpr', 

250 '-o', self.label + '.trr', 

251 '-e', self.label + '.edr', 

252 '-g', self.label + '.log', 

253 '-rerun', self.label + '.g96', 

254 self.params_runs.get('extra_mdrun_parameters', ''), 

255 '> QMMM.log 2>&1'] 

256 command = ' '.join(command) 

257 self._execute_gromacs(command) 

258 else: 

259 command += [ 

260 '-s', self.label + '.tpr', 

261 '-o', self.label + '.trr', 

262 '-e', self.label + '.edr', 

263 '-g', self.label + '.log', 

264 '-c', self.label + '.g96', 

265 self.params_runs.get('extra_mdrun_parameters', ''), 

266 '> MM.log 2>&1'] 

267 command = ' '.join(command) 

268 self._execute_gromacs(command) 

269 

270 atoms = read_gromos(self.label + '.g96') 

271 self.atoms = atoms.copy() 

272 

273 def generate_topology_and_g96file(self): 

274 """ from coordinates (self.label.+'pdb') 

275 and gromacs run input file (self.label + '.mdp) 

276 generate topology (self.label+'top') 

277 and structure file in .g96 format (self.label + '.g96') 

278 """ 

279 # generate structure and topology files 

280 # In case of predefinded topology file this is not done 

281 subcmd = 'pdb2gmx' 

282 command = ' '.join([ 

283 subcmd, 

284 '-f', self.params_runs['init_structure'], 

285 '-o', self.label + '.g96', 

286 '-p', self.label + '.top', 

287 '-ff', self.params_runs['force_field'], 

288 '-water', self.params_runs['water'], 

289 self.params_runs.get('extra_pdb2gmx_parameters', ''), 

290 '> {}.{}.log 2>&1'.format(self.label, subcmd)]) 

291 self._execute_gromacs(command) 

292 

293 atoms = read_gromos(self.label + '.g96') 

294 self.atoms = atoms.copy() 

295 

296 def generate_gromacs_run_file(self): 

297 """ Generates input file for a gromacs mdrun 

298 based on structure file and topology file 

299 resulting file is self.label + '.tpr 

300 """ 

301 

302 # generate gromacs run input file (gromacs.tpr) 

303 try: 

304 os.remove(self.label + '.tpr') 

305 except OSError: 

306 pass 

307 

308 subcmd = 'grompp' 

309 command = ' '.join([ 

310 subcmd, 

311 '-f', self.label + '.mdp', 

312 '-c', self.label + '.g96', 

313 '-p', self.label + '.top', 

314 '-o', self.label + '.tpr', 

315 '-maxwarn', '100', 

316 self.params_runs.get('extra_grompp_parameters', ''), 

317 '> {}.{}.log 2>&1'.format(self.label, subcmd)]) 

318 self._execute_gromacs(command) 

319 

320 def write_energy_files(self): 

321 """write input files for gromacs force and energy calculations 

322 for gromacs program energy""" 

323 filename = 'inputGenergy.txt' 

324 with open(filename, 'w') as output: 

325 output.write('Potential \n') 

326 output.write(' \n') 

327 output.write(' \n') 

328 

329 filename = 'inputGtraj.txt' 

330 with open(filename, 'w') as output: 

331 output.write('System \n') 

332 output.write(' \n') 

333 output.write(' \n') 

334 

335 def set_own_params(self, key, value, docstring=""): 

336 """Set own gromacs parameter with doc strings.""" 

337 self.parameters[key] = value 

338 self.params_doc[key] = docstring 

339 

340 def set_own_params_runs(self, key, value): 

341 """Set own gromacs parameter for program parameters 

342 Add spaces to avoid errors """ 

343 self.params_runs[key] = ' ' + value + ' ' 

344 

345 def write_input(self, atoms=None, properties=None, system_changes=None): 

346 """Write input parameters to input file.""" 

347 

348 FileIOCalculator.write_input(self, atoms, properties, system_changes) 

349 # print self.parameters 

350 with open(self.label + '.mdp', 'w') as myfile: 

351 for key, val in self.parameters.items(): 

352 if val is not None: 

353 docstring = self.params_doc.get(key, '') 

354 myfile.write('%-35s = %s ; %s\n' 

355 % (key, val, ';' + docstring)) 

356 

357 def update(self, atoms): 

358 """ set atoms and do the calculation """ 

359 # performs an update of the atoms 

360 self.atoms = atoms.copy() 

361 # must be g96 format for accuracy, alternatively binary formats 

362 write_gromos(self.label + '.g96', atoms) 

363 # does run to get forces and energies 

364 self.calculate() 

365 

366 def calculate(self, atoms=None, properties=['energy', 'forces'], 

367 system_changes=all_changes): 

368 """ runs a gromacs-mdrun and 

369 gets energy and forces 

370 rest below is to make gromacs calculator 

371 compactible with ase-Calculator class 

372 

373 atoms: Atoms object 

374 Contains positions, unit-cell, ... 

375 properties: list of str 

376 List of what needs to be calculated. Can be any combination 

377 of 'energy', 'forces' 

378 system_changes: list of str 

379 List of what has changed since last calculation. Can be 

380 any combination of these five: 'positions', 'numbers', 'cell', 

381 'pbc', 'initial_charges' and 'initial_magmoms'. 

382 

383 """ 

384 

385 self.run() 

386 if self.clean: 

387 do_clean('#*') 

388 # get energy 

389 try: 

390 os.remove('tmp_ene.del') 

391 except OSError: 

392 pass 

393 

394 subcmd = 'energy' 

395 command = ' '.join([ 

396 subcmd, 

397 '-f', self.label + '.edr', 

398 '-o', self.label + '.Energy.xvg', 

399 '< inputGenergy.txt', 

400 '> {}.{}.log 2>&1'.format(self.label, subcmd)]) 

401 self._execute_gromacs(command) 

402 with open(self.label + '.Energy.xvg') as fd: 

403 lastline = fd.readlines()[-1] 

404 energy = float(lastline.split()[1]) 

405 # We go for ASE units ! 

406 self.results['energy'] = energy * units.kJ / units.mol 

407 # energies are about 100 times bigger in Gromacs units 

408 # when compared to ase units 

409 

410 subcmd = 'traj' 

411 command = ' '.join([ 

412 subcmd, 

413 '-f', self.label + '.trr', 

414 '-s', self.label + '.tpr', 

415 '-of', self.label + '.Force.xvg', 

416 '< inputGtraj.txt', 

417 '> {}.{}.log 2>&1'.format(self.label, subcmd)]) 

418 self._execute_gromacs(command) 

419 with open(self.label + '.Force.xvg', 'r') as fd: 

420 lastline = fd.readlines()[-1] 

421 forces = np.array([float(f) for f in lastline.split()[1:]]) 

422 # We go for ASE units !gromacsForce.xvg 

423 tmp_forces = forces / units.nm * units.kJ / units.mol 

424 tmp_forces = np.reshape(tmp_forces, (-1, 3)) 

425 self.results['forces'] = tmp_forces