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