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