Coverage for /builds/debichem-team/python-ase/ase/calculators/gromacs.py : 83.01%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""This module defines an ASE interface to GROMACS.
3http://www.gromacs.org/
4It is VERY SLOW compared to standard Gromacs
5(due to slow formatted io required here).
7Mainly intended to be the MM part in the ase QM/MM
9Markus.Kaukonen@iki.fi
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
17"""
19import os
20import subprocess
21from glob import glob
22from shutil import which
24import numpy as np
26from ase import units
27from ase.calculators.calculator import (EnvironmentError,
28 FileIOCalculator,
29 all_changes)
30from ase.io.gromos import read_gromos, write_gromos
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)
39def get_gromacs_version(executable):
40 output = subprocess.check_output([executable, '--version'],
41 encoding='utf-8')
42 return parse_gromacs_version(output)
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
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.)
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::
65 CALC_MM_RELAX = Gromacs()
66 CALC_MM_RELAX.set_own_params('integrator', 'steep',
67 'use steepest descent')
69 Run command line arguments for gromacs related programs:
70 pdb2gmx, grompp, mdrun, energy, traj. These can be given as::
72 CALC_MM_RELAX = Gromacs()
73 CALC_MM_RELAX.set_own_params_runs('force_field', 'oplsaa')
74 """
76 implemented_properties = ['energy', 'forces']
77 discard_results_on_any_change = True
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')
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.
105 Parameters
106 ==========
107 label: str
108 Prefix to use for filenames (label.in, label.txt, ...).
109 Default is 'gromacs'.
111 do_qmmm : bool
112 Is gromacs used as mm calculator for a qm/mm calculation
114 clean : bool
115 Remove gromacs backup files
116 and old gormacs.* files
118 water_model: str
119 Water model to be used in gromacs runs (see gromacs manual)
121 force_field: str
122 Force field to be used in gromacs runs
124 command : str
125 Gromacs executable; if None (default), choose available one from
126 ('gmx', 'gmx_d', 'gmx_mpi', 'gmx_mpi_d')
127 """
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)
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'
155 self.positions = None
156 self.atoms = None
157 # storage for energy and forces
158 #self.energy = None
159 #self.forces = None
161 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
162 label, atoms, **kwargs)
163 self.set(**kwargs)
164 # default values for runtime parameters
165 # can be changed by self.set_own_params_runs('key', 'value')
166 self.params_runs = {}
167 self.params_runs['index_filename'] = 'index.ndx'
168 self.params_runs['init_structure'] = self.label + '.pdb'
169 self.params_runs['water'] = self.water_model
170 self.params_runs['force_field'] = self.force_field
172 # these below are required by qm/mm
173 self.topology_filename = self.label + '.top'
174 self.name = 'Gromacs'
176 # clean up gromacs backups
177 if self.clean:
178 do_clean('gromacs.???')
180 # write input files for gromacs program energy
181 self.write_energy_files()
183 if self.do_qmmm:
184 self.parameters['integrator'] = 'md'
185 self.parameters['nsteps'] = '0'
187 def _execute_gromacs(self, command):
188 """ execute gmx command
189 Parameters
190 ----------
191 command : str
192 """
193 if self.command:
194 subprocess.check_call(self.command + ' ' + command, shell=True)
195 else:
196 raise EnvironmentError(self.missing_gmx)
198 def generate_g96file(self):
199 """ from current coordinates (self.structure_file)
200 write a structure file in .g96 format
201 """
202 # generate structure file in g96 format
203 write_gromos(self.label + '.g96', self.atoms)
205 def run_editconf(self):
206 """ run gromacs program editconf, typically to set a simulation box
207 writing to the input structure"""
208 subcmd = 'editconf'
209 command = ' '.join([
210 subcmd,
211 '-f', self.label + '.g96',
212 '-o', self.label + '.g96',
213 self.params_runs.get('extra_editconf_parameters', ''),
214 '> {}.{}.log 2>&1'.format(self.label, subcmd)])
215 self._execute_gromacs(command)
217 def run_genbox(self):
218 """Run gromacs program genbox, typically to solvate the system
219 writing to the input structure
220 as extra parameter you need to define the file containing the solvent
222 for instance::
224 CALC_MM_RELAX = Gromacs()
225 CALC_MM_RELAX.set_own_params_runs(
226 'extra_genbox_parameters', '-cs spc216.gro')
227 """
228 subcmd = 'genbox'
229 command = ' '.join([
230 subcmd,
231 '-cp', self.label + '.g96',
232 '-o', self.label + '.g96',
233 '-p', self.label + '.top',
234 self.params_runs.get('extra_genbox_parameters', ''),
235 '> {}.{}.log 2>&1'.format(self.label, subcmd)])
236 self._execute_gromacs(command)
238 def run(self):
239 """ runs a gromacs-mdrun with the
240 current atom-configuration """
242 # clean up gromacs backups
243 if self.clean:
244 do_clean('#*')
246 subcmd = 'mdrun'
247 command = [subcmd]
248 if self.do_qmmm:
249 command += [
250 '-s', self.label + '.tpr',
251 '-o', self.label + '.trr',
252 '-e', self.label + '.edr',
253 '-g', self.label + '.log',
254 '-rerun', self.label + '.g96',
255 self.params_runs.get('extra_mdrun_parameters', ''),
256 '> QMMM.log 2>&1']
257 command = ' '.join(command)
258 self._execute_gromacs(command)
259 else:
260 command += [
261 '-s', self.label + '.tpr',
262 '-o', self.label + '.trr',
263 '-e', self.label + '.edr',
264 '-g', self.label + '.log',
265 '-c', self.label + '.g96',
266 self.params_runs.get('extra_mdrun_parameters', ''),
267 '> MM.log 2>&1']
268 command = ' '.join(command)
269 self._execute_gromacs(command)
271 atoms = read_gromos(self.label + '.g96')
272 self.atoms = atoms.copy()
274 def generate_topology_and_g96file(self):
275 """ from coordinates (self.label.+'pdb')
276 and gromacs run input file (self.label + '.mdp)
277 generate topology (self.label+'top')
278 and structure file in .g96 format (self.label + '.g96')
279 """
280 #generate structure and topology files
281 # In case of predefinded topology file this is not done
282 subcmd = 'pdb2gmx'
283 command = ' '.join([
284 subcmd,
285 '-f', self.params_runs['init_structure'],
286 '-o', self.label + '.g96',
287 '-p', self.label + '.top',
288 '-ff', self.params_runs['force_field'],
289 '-water', self.params_runs['water'],
290 self.params_runs.get('extra_pdb2gmx_parameters', ''),
291 '> {}.{}.log 2>&1'.format(self.label, subcmd)])
292 self._execute_gromacs(command)
294 atoms = read_gromos(self.label + '.g96')
295 self.atoms = atoms.copy()
297 def generate_gromacs_run_file(self):
298 """ Generates input file for a gromacs mdrun
299 based on structure file and topology file
300 resulting file is self.label + '.tpr
301 """
303 #generate gromacs run input file (gromacs.tpr)
304 try:
305 os.remove(self.label + '.tpr')
306 except OSError:
307 pass
309 subcmd = 'grompp'
310 command = ' '.join([
311 subcmd,
312 '-f', self.label + '.mdp',
313 '-c', self.label + '.g96',
314 '-p', self.label + '.top',
315 '-o', self.label + '.tpr',
316 '-maxwarn', '100',
317 self.params_runs.get('extra_grompp_parameters', ''),
318 '> {}.{}.log 2>&1'.format(self.label, subcmd)])
319 self._execute_gromacs(command)
321 def write_energy_files(self):
322 """write input files for gromacs force and energy calculations
323 for gromacs program energy"""
324 filename = 'inputGenergy.txt'
325 with open(filename, 'w') as output:
326 output.write('Potential \n')
327 output.write(' \n')
328 output.write(' \n')
330 filename = 'inputGtraj.txt'
331 with open(filename, 'w') as output:
332 output.write('System \n')
333 output.write(' \n')
334 output.write(' \n')
336 def set_own_params(self, key, value, docstring=""):
337 """Set own gromacs parameter with doc strings."""
338 self.parameters[key] = value
339 self.params_doc[key] = docstring
341 def set_own_params_runs(self, key, value):
342 """Set own gromacs parameter for program parameters
343 Add spaces to avoid errors """
344 self.params_runs[key] = ' ' + value + ' '
346 def write_input(self, atoms=None, properties=None, system_changes=None):
347 """Write input parameters to input file."""
349 FileIOCalculator.write_input(self, atoms, properties, system_changes)
350 #print self.parameters
351 with open(self.label + '.mdp', 'w') as myfile:
352 for key, val in self.parameters.items():
353 if val is not None:
354 docstring = self.params_doc.get(key, '')
355 myfile.write('%-35s = %s ; %s\n'
356 % (key, val, ';' + docstring))
358 def update(self, atoms):
359 """ set atoms and do the calculation """
360 # performs an update of the atoms
361 self.atoms = atoms.copy()
362 #must be g96 format for accuracy, alternatively binary formats
363 write_gromos(self.label + '.g96', atoms)
364 # does run to get forces and energies
365 self.calculate()
367 def calculate(self, atoms=None, properties=['energy', 'forces'],
368 system_changes=all_changes):
369 """ runs a gromacs-mdrun and
370 gets energy and forces
371 rest below is to make gromacs calculator
372 compactible with ase-Calculator class
374 atoms: Atoms object
375 Contains positions, unit-cell, ...
376 properties: list of str
377 List of what needs to be calculated. Can be any combination
378 of 'energy', 'forces'
379 system_changes: list of str
380 List of what has changed since last calculation. Can be
381 any combination of these five: 'positions', 'numbers', 'cell',
382 'pbc', 'initial_charges' and 'initial_magmoms'.
384 """
386 self.run()
387 if self.clean:
388 do_clean('#*')
389 # get energy
390 try:
391 os.remove('tmp_ene.del')
392 except OSError:
393 pass
395 subcmd = 'energy'
396 command = ' '.join([
397 subcmd,
398 '-f', self.label + '.edr',
399 '-o', self.label + '.Energy.xvg',
400 '< inputGenergy.txt',
401 '> {}.{}.log 2>&1'.format(self.label, subcmd)])
402 self._execute_gromacs(command)
403 with open(self.label + '.Energy.xvg') as fd:
404 lastline = fd.readlines()[-1]
405 energy = float(lastline.split()[1])
406 #We go for ASE units !
407 #self.energy = energy * units.kJ / units.mol
408 self.results['energy'] = energy * units.kJ / units.mol
409 # energies are about 100 times bigger in Gromacs units
410 # when compared to ase units
412 subcmd = 'traj'
413 command = ' '.join([
414 subcmd,
415 '-f', self.label + '.trr',
416 '-s', self.label + '.tpr',
417 '-of', self.label + '.Force.xvg',
418 '< inputGtraj.txt',
419 '> {}.{}.log 2>&1'.format(self.label, subcmd)])
420 self._execute_gromacs(command)
421 with open(self.label + '.Force.xvg', 'r') as fd:
422 lastline = fd.readlines()[-1]
423 forces = np.array([float(f) for f in lastline.split()[1:]])
424 #We go for ASE units !gromacsForce.xvg
425 #self.forces = np.array(forces)/ units.nm * units.kJ / units.mol
426 #self.forces = np.reshape(self.forces, (-1, 3))
427 tmp_forces = forces / units.nm * units.kJ / units.mol
428 tmp_forces = np.reshape(tmp_forces, (-1, 3))
429 self.results['forces'] = tmp_forces
430 #self.forces = np.array(forces)