Coverage for /builds/debichem-team/python-ase/ase/calculators/gromacs.py: 86.11%
144 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
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
23import numpy as np
25from ase import units
26from ase.calculators.calculator import (
27 CalculatorSetupError,
28 FileIOCalculator,
29 all_changes,
30)
31from ase.io.gromos import read_gromos, write_gromos
34def parse_gromacs_version(output):
35 import re
36 match = re.search(r'GROMACS version\:\s*(\S+)', output, re.M)
37 return match.group(1)
40def get_gromacs_version(executable):
41 output = subprocess.check_output([executable, '--version'],
42 encoding='utf-8')
43 return parse_gromacs_version(output)
46def do_clean(name='#*'):
47 """ remove files matching wildcards """
48 myfiles = glob(name)
49 for myfile in myfiles:
50 try:
51 os.remove(myfile)
52 except OSError:
53 pass
56class Gromacs(FileIOCalculator):
57 """Class for doing GROMACS calculations.
58 Before running a gromacs calculation you must prepare the input files
59 separately (pdb2gmx and grompp for instance.)
61 Input parameters for gromacs runs (the .mdp file)
62 are given in self.params and can be set when initializing the calculator
63 or by method set_own.
64 for example::
66 CALC_MM_RELAX = Gromacs()
67 CALC_MM_RELAX.set_own_params('integrator', 'steep',
68 'use steepest descent')
70 Run command line arguments for gromacs related programs:
71 pdb2gmx, grompp, mdrun, energy, traj. These can be given as::
73 CALC_MM_RELAX = Gromacs()
74 CALC_MM_RELAX.set_own_params_runs('force_field', 'oplsaa')
75 """
77 implemented_properties = ['energy', 'forces']
78 discard_results_on_any_change = True
80 default_parameters = dict(
81 define='-DFLEXIBLE',
82 integrator='cg',
83 nsteps='10000',
84 nstfout='10',
85 nstlog='10',
86 nstenergy='10',
87 nstlist='10',
88 ns_type='grid',
89 pbc='xyz',
90 rlist='1.15',
91 coulombtype='PME-Switch',
92 rcoulomb='0.8',
93 vdwtype='shift',
94 rvdw='0.8',
95 rvdw_switch='0.75',
96 DispCorr='Ener')
98 def __init__(self, restart=None,
99 ignore_bad_restart_file=FileIOCalculator._deprecated,
100 label='gromacs', atoms=None,
101 do_qmmm=False, clean=True,
102 water_model='tip3p', force_field='oplsaa', command=None,
103 **kwargs):
104 """Construct GROMACS-calculator object.
106 Parameters
107 ==========
108 label: str
109 Prefix to use for filenames (label.in, label.txt, ...).
110 Default is 'gromacs'.
112 do_qmmm : bool
113 Is gromacs used as mm calculator for a qm/mm calculation
115 clean : bool
116 Remove gromacs backup files
117 and old gormacs.* files
119 water_model: str
120 Water model to be used in gromacs runs (see gromacs manual)
122 force_field: str
123 Force field to be used in gromacs runs
125 command : str
126 Gromacs executable; if None (default), choose available one from
127 ('gmx', 'gmx_d', 'gmx_mpi', 'gmx_mpi_d')
128 """
130 self.do_qmmm = do_qmmm
131 self.water_model = water_model
132 self.force_field = force_field
133 self.clean = clean
134 self.params_doc = {}
135 # add comments for gromacs input file
136 self.params_doc['define'] = \
137 'flexible/ rigid water'
138 self.params_doc['integrator'] = \
139 'md: molecular dynamics(Leapfrog), \n' + \
140 '; md-vv: molecular dynamics(Velocity Verlet), \n' + \
141 '; steep: steepest descent minimization, \n' + \
142 '; cg: conjugate cradient minimization \n'
144 self.positions = None
145 self.atoms = None
147 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
148 label, atoms, command=command,
149 **kwargs)
150 self.set(**kwargs)
151 # default values for runtime parameters
152 # can be changed by self.set_own_params_runs('key', 'value')
153 self.params_runs = {}
154 self.params_runs['index_filename'] = 'index.ndx'
155 self.params_runs['init_structure'] = self.label + '.pdb'
156 self.params_runs['water'] = self.water_model
157 self.params_runs['force_field'] = self.force_field
159 # these below are required by qm/mm
160 self.topology_filename = self.label + '.top'
162 # clean up gromacs backups
163 if self.clean:
164 do_clean('gromacs.???')
166 # write input files for gromacs program energy
167 self.write_energy_files()
169 if self.do_qmmm:
170 self.parameters['integrator'] = 'md'
171 self.parameters['nsteps'] = '0'
173 def _get_name(self):
174 return 'Gromacs'
176 def _execute_gromacs(self, command):
177 """ execute gmx command
178 Parameters
179 ----------
180 command : str
181 """
182 if self.command:
183 subprocess.check_call(self.command + ' ' + command, shell=True)
184 else:
185 raise CalculatorSetupError('Missing gromacs executable')
187 def generate_g96file(self):
188 """ from current coordinates (self.structure_file)
189 write a structure file in .g96 format
190 """
191 # generate structure file in g96 format
192 write_gromos(self.label + '.g96', self.atoms)
194 def run_editconf(self):
195 """ run gromacs program editconf, typically to set a simulation box
196 writing to the input structure"""
197 subcmd = 'editconf'
198 command = ' '.join([
199 subcmd,
200 '-f', self.label + '.g96',
201 '-o', self.label + '.g96',
202 self.params_runs.get('extra_editconf_parameters', ''),
203 f'> {self.label}.{subcmd}.log 2>&1'])
204 self._execute_gromacs(command)
206 def run_genbox(self):
207 """Run gromacs program genbox, typically to solvate the system
208 writing to the input structure
209 as extra parameter you need to define the file containing the solvent
211 for instance::
213 CALC_MM_RELAX = Gromacs()
214 CALC_MM_RELAX.set_own_params_runs(
215 'extra_genbox_parameters', '-cs spc216.gro')
216 """
217 subcmd = 'genbox'
218 command = ' '.join([
219 subcmd,
220 '-cp', self.label + '.g96',
221 '-o', self.label + '.g96',
222 '-p', self.label + '.top',
223 self.params_runs.get('extra_genbox_parameters', ''),
224 f'> {self.label}.{subcmd}.log 2>&1'])
225 self._execute_gromacs(command)
227 def run(self):
228 """ runs a gromacs-mdrun with the
229 current atom-configuration """
231 # clean up gromacs backups
232 if self.clean:
233 do_clean('#*')
235 subcmd = 'mdrun'
236 command = [subcmd]
237 if self.do_qmmm:
238 command += [
239 '-s', self.label + '.tpr',
240 '-o', self.label + '.trr',
241 '-e', self.label + '.edr',
242 '-g', self.label + '.log',
243 '-rerun', self.label + '.g96',
244 self.params_runs.get('extra_mdrun_parameters', ''),
245 '> QMMM.log 2>&1']
246 command = ' '.join(command)
247 self._execute_gromacs(command)
248 else:
249 command += [
250 '-s', self.label + '.tpr',
251 '-o', self.label + '.trr',
252 '-e', self.label + '.edr',
253 '-g', self.label + '.log',
254 '-c', self.label + '.g96',
255 self.params_runs.get('extra_mdrun_parameters', ''),
256 '> MM.log 2>&1']
257 command = ' '.join(command)
258 self._execute_gromacs(command)
260 atoms = read_gromos(self.label + '.g96')
261 self.atoms = atoms.copy()
263 def generate_topology_and_g96file(self):
264 """ from coordinates (self.label.+'pdb')
265 and gromacs run input file (self.label + '.mdp)
266 generate topology (self.label+'top')
267 and structure file in .g96 format (self.label + '.g96')
268 """
269 # generate structure and topology files
270 # In case of predefinded topology file this is not done
271 subcmd = 'pdb2gmx'
272 command = ' '.join([
273 subcmd,
274 '-f', self.params_runs['init_structure'],
275 '-o', self.label + '.g96',
276 '-p', self.label + '.top',
277 '-ff', self.params_runs['force_field'],
278 '-water', self.params_runs['water'],
279 self.params_runs.get('extra_pdb2gmx_parameters', ''),
280 f'> {self.label}.{subcmd}.log 2>&1'])
281 self._execute_gromacs(command)
283 atoms = read_gromos(self.label + '.g96')
284 self.atoms = atoms.copy()
286 def generate_gromacs_run_file(self):
287 """ Generates input file for a gromacs mdrun
288 based on structure file and topology file
289 resulting file is self.label + '.tpr
290 """
292 # generate gromacs run input file (gromacs.tpr)
293 try:
294 os.remove(self.label + '.tpr')
295 except OSError:
296 pass
298 subcmd = 'grompp'
299 command = ' '.join([
300 subcmd,
301 '-f', self.label + '.mdp',
302 '-c', self.label + '.g96',
303 '-p', self.label + '.top',
304 '-o', self.label + '.tpr',
305 '-maxwarn', '100',
306 self.params_runs.get('extra_grompp_parameters', ''),
307 f'> {self.label}.{subcmd}.log 2>&1'])
308 self._execute_gromacs(command)
310 def write_energy_files(self):
311 """write input files for gromacs force and energy calculations
312 for gromacs program energy"""
313 filename = 'inputGenergy.txt'
314 with open(filename, 'w') as output:
315 output.write('Potential \n')
316 output.write(' \n')
317 output.write(' \n')
319 filename = 'inputGtraj.txt'
320 with open(filename, 'w') as output:
321 output.write('System \n')
322 output.write(' \n')
323 output.write(' \n')
325 def set_own_params(self, key, value, docstring=""):
326 """Set own gromacs parameter with doc strings."""
327 self.parameters[key] = value
328 self.params_doc[key] = docstring
330 def set_own_params_runs(self, key, value):
331 """Set own gromacs parameter for program parameters
332 Add spaces to avoid errors """
333 self.params_runs[key] = ' ' + value + ' '
335 def write_input(self, atoms=None, properties=None, system_changes=None):
336 """Write input parameters to input file."""
338 FileIOCalculator.write_input(self, atoms, properties, system_changes)
339 # print self.parameters
340 with open(self.label + '.mdp', 'w') as myfile:
341 for key, val in self.parameters.items():
342 if val is not None:
343 docstring = self.params_doc.get(key, '')
344 myfile.write('%-35s = %s ; %s\n'
345 % (key, val, ';' + docstring))
347 def update(self, atoms):
348 """ set atoms and do the calculation """
349 # performs an update of the atoms
350 self.atoms = atoms.copy()
351 # must be g96 format for accuracy, alternatively binary formats
352 write_gromos(self.label + '.g96', atoms)
353 # does run to get forces and energies
354 self.calculate()
356 def calculate(self, atoms=None, properties=['energy', 'forces'],
357 system_changes=all_changes):
358 """ runs a gromacs-mdrun and
359 gets energy and forces
360 rest below is to make gromacs calculator
361 compactible with ase-Calculator class
363 atoms: Atoms object
364 Contains positions, unit-cell, ...
365 properties: list of str
366 List of what needs to be calculated. Can be any combination
367 of 'energy', 'forces'
368 system_changes: list of str
369 List of what has changed since last calculation. Can be
370 any combination of these five: 'positions', 'numbers', 'cell',
371 'pbc', 'initial_charges' and 'initial_magmoms'.
373 """
375 self.run()
376 if self.clean:
377 do_clean('#*')
378 # get energy
379 try:
380 os.remove('tmp_ene.del')
381 except OSError:
382 pass
384 subcmd = 'energy'
385 command = ' '.join([
386 subcmd,
387 '-f', self.label + '.edr',
388 '-o', self.label + '.Energy.xvg',
389 '< inputGenergy.txt',
390 f'> {self.label}.{subcmd}.log 2>&1'])
391 self._execute_gromacs(command)
392 with open(self.label + '.Energy.xvg') as fd:
393 lastline = fd.readlines()[-1]
394 energy = float(lastline.split()[1])
395 # We go for ASE units !
396 self.results['energy'] = energy * units.kJ / units.mol
397 # energies are about 100 times bigger in Gromacs units
398 # when compared to ase units
400 subcmd = 'traj'
401 command = ' '.join([
402 subcmd,
403 '-f', self.label + '.trr',
404 '-s', self.label + '.tpr',
405 '-of', self.label + '.Force.xvg',
406 '< inputGtraj.txt',
407 f'> {self.label}.{subcmd}.log 2>&1'])
408 self._execute_gromacs(command)
409 with open(self.label + '.Force.xvg') as fd:
410 lastline = fd.readlines()[-1]
411 forces = np.array([float(f) for f in lastline.split()[1:]])
412 # We go for ASE units !gromacsForce.xvg
413 tmp_forces = forces / units.nm * units.kJ / units.mol
414 tmp_forces = np.reshape(tmp_forces, (-1, 3))
415 self.results['forces'] = tmp_forces