Coverage for biobb_gromacs/gromacs/common.py: 76%
204 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-25 09:23 +0000
1""" Common functions for package biobb_gromacs.gromacs """
2import re
3import warnings
4from pathlib import Path
5from biobb_common.tools import file_utils as fu
6from biobb_common.command_wrapper import cmd_wrapper
7from typing import Mapping, Optional
10def get_gromacs_version(gmx: str = "gmx", minimum_version: int = 512) -> int:
11 """ Gets the GROMACS installed version and returns it as an int(3) for
12 versions older than 5.1.5 and an int(5) for 20XX versions filling the gaps
13 with '0' digits.
15 Args:
16 gmx (str): ('gmx') Path to the GROMACS binary.
17 minimum_version (int): (512) Minimum GROMACS version required.
19 Returns:
20 int: GROMACS version.
21 """
22 unique_dir = fu.create_unique_dir()
23 out_log, err_log = fu.get_logs(path=unique_dir, can_write_console=False)
24 cmd = [gmx, "-version"]
25 try:
26 cmd_wrapper.CmdWrapper(cmd=cmd, out_log=out_log, err_log=err_log).launch()
27 pattern = re.compile(r"GROMACS version:\s+(.+)")
28 version_str: Optional[re.Match[str]] = None
29 with open(Path(unique_dir).joinpath('log.out')) as log_file:
30 for line in log_file:
31 version_str = pattern.match(line.strip())
32 if version_str:
33 break
34 if not version_str:
35 warnings.warn("GROMACS version not found. Your GROMACS installation Biobb compatibility has not been tested.")
36 return 0
37 version = version_str.group(1).replace(".", "").replace("VERSION", "").strip()
38 version = "".join([c for c in version if c.isdigit()])
39 except Exception:
40 warnings.warn("GROMACS version not found. Your GROMACS installation Biobb compatibility has not been tested.")
41 return 0
42 if version.startswith("2"):
43 while len(version) < 5:
44 version += '0'
45 else:
46 while len(version) < 3:
47 version += '0'
49 fu.rm(unique_dir)
50 if int(version) < minimum_version:
51 warnings.warn(f"GROMACS version should be {minimum_version} or newer {version} detected")
52 return int(version)
55# class GromacsVersionError(Exception):
56# """ Exception Raised when the installed version of GROMACS is not
57# compatible with the current function.
58# """
59# # pass
62def gmx_check(file_a: str, file_b: str, gmx: str = 'gmx') -> bool:
63 print("Comparing GROMACS files:")
64 print("FILE_A: %s" % str(Path(file_a).resolve()))
65 print("FILE_B: %s" % str(Path(file_b).resolve()))
66 check_result = 'check_result.out'
67 cmd = [gmx, 'check']
68 if file_a.endswith(".tpr"):
69 cmd.append('-s1')
70 else:
71 cmd.append('-f')
72 cmd.append(file_a)
73 if file_b.endswith(".tpr"):
74 cmd.append('-s2')
75 else:
76 cmd.append('-f2')
77 cmd.append(file_b)
78 cmd.append('> check_result.out')
79 cmd_wrapper.CmdWrapper(cmd).launch()
80 print("Result file: %s" % str(Path(check_result).resolve()))
81 with open(check_result) as check_file:
82 for line_num, line in enumerate(check_file):
83 if not line.rstrip():
84 continue
85 if line.startswith("Both files read correctly"):
86 continue
87 if not line.startswith('comparing'):
88 print('Discrepance found in line %d: %s' % (line_num, line))
89 return False
90 return True
93def gmx_rms(file_a: str, file_b: str, file_tpr: str, gmx: str = 'gmx', tolerance: float = 0.5):
94 print("Comparing GROMACS files:")
95 print("FILE_A: %s" % str(Path(file_a).resolve()))
96 print("FILE_B: %s" % str(Path(file_b).resolve()))
97 rmsd_result = 'rmsd.xvg'
98 cmd = ['echo', '\"Protein Protein\"', '|',
99 gmx, 'rms', '-s', file_tpr, '-f', file_a, '-f2', file_b, '-xvg', 'none']
100 cmd_wrapper.CmdWrapper(cmd).launch()
101 print("Result file: %s" % str(Path(rmsd_result).resolve()))
102 with open(rmsd_result) as check_file:
103 for line in check_file:
104 time_step, rmsd = tuple(line.strip().split())
105 if float(rmsd) > tolerance:
106 print('RMSD: %s bigger than tolerance %g for time step %s' % (rmsd, tolerance, time_step))
107 return False
108 return True
111def read_mdp(input_mdp_path: str) -> dict[str, str]:
112 # Credit for these two reg exps to:
113 # https://github.com/Becksteinlab/GromacsWrapper/blob/master/gromacs/fileformats/mdp.py
114 parameter_re = re.compile(r"\s*(?P<parameter>[^=]+?)\s*=\s*(?P<value>[^;]*)(?P<comment>\s*;.*)?", re.VERBOSE)
116 mdp_dict = {}
117 with open(input_mdp_path) as mdp_file:
118 for line in mdp_file:
119 re_match = parameter_re.match(line.strip())
120 if re_match:
121 parameter = re_match.group('parameter')
122 value = re_match.group('value')
123 mdp_dict[parameter] = value
125 return mdp_dict
128def mdp_preset(sim_type: str) -> dict[str, str]:
129 mdp_dict: dict[str, str] = {}
130 if not sim_type or sim_type == 'index':
131 return mdp_dict
133 minimization = (sim_type == 'minimization') or (sim_type == 'ions')
134 nvt = (sim_type == 'nvt')
135 npt = (sim_type == 'npt')
136 free = (sim_type == 'free')
137 md = (nvt or npt or free)
139 # Position restrain
140 if not free:
141 mdp_dict['define'] = '-DPOSRES'
143 # Run parameters
144 mdp_dict['nsteps'] = '5000'
145 if minimization:
146 mdp_dict['integrator'] = 'steep'
147 mdp_dict['emtol'] = '1000.0'
148 mdp_dict['emstep'] = '0.01'
149 if md:
150 mdp_dict['integrator'] = 'md'
151 mdp_dict['dt'] = '0.002'
153 # Output control
154 if md:
155 if nvt or npt:
156 mdp_dict['nstxout'] = '500'
157 mdp_dict['nstvout'] = '500'
158 mdp_dict['nstenergy'] = '500'
159 mdp_dict['nstlog'] = '500'
160 mdp_dict['nstcalcenergy'] = '100'
161 mdp_dict['nstcomm'] = '100'
162 mdp_dict['nstxout-compressed'] = '1000'
163 mdp_dict['compressed-x-precision'] = '1000'
164 mdp_dict['compressed-x-grps'] = 'System'
165 if free:
166 mdp_dict['nstcomm'] = '100'
167 mdp_dict['nstxout'] = '5000'
168 mdp_dict['nstvout'] = '5000'
169 mdp_dict['nstenergy'] = '5000'
170 mdp_dict['nstlog'] = '5000'
171 mdp_dict['nstcalcenergy'] = '100'
172 mdp_dict['nstxout-compressed'] = '1000'
173 mdp_dict['compressed-x-grps'] = 'System'
174 mdp_dict['compressed-x-precision'] = '1000'
176 # Bond parameters
177 if md:
178 mdp_dict['constraint-algorithm'] = 'lincs'
179 mdp_dict['constraints'] = 'h-bonds'
180 mdp_dict['lincs-iter'] = '1'
181 mdp_dict['lincs-order'] = '4'
182 if nvt:
183 mdp_dict['continuation'] = 'no'
184 if npt or free:
185 mdp_dict['continuation'] = 'yes'
187 # Neighbour searching
188 mdp_dict['cutoff-scheme'] = 'Verlet'
189 mdp_dict['ns-type'] = 'grid'
190 mdp_dict['rcoulomb'] = '1.0'
191 mdp_dict['vdwtype'] = 'cut-off'
192 mdp_dict['rvdw'] = '1.0'
193 mdp_dict['nstlist'] = '10'
194 mdp_dict['rlist'] = '1'
196 # Electrostatics
197 mdp_dict['coulombtype'] = 'PME'
198 if md:
199 mdp_dict['pme-order'] = '4'
200 mdp_dict['fourierspacing'] = '0.12'
201 mdp_dict['fourier-nx'] = '0'
202 mdp_dict['fourier-ny'] = '0'
203 mdp_dict['fourier-nz'] = '0'
204 mdp_dict['ewald-rtol'] = '1e-5'
206 # Temperature coupling
207 if md:
208 mdp_dict['tcoupl'] = 'V-rescale'
209 mdp_dict['tc-grps'] = 'Protein Non-Protein'
210 mdp_dict['tau-t'] = '0.1 0.1'
211 mdp_dict['ref-t'] = '300 300'
213 # Pressure coupling
214 if md:
215 if nvt:
216 mdp_dict['pcoupl'] = 'no'
217 if npt or free:
218 mdp_dict['pcoupl'] = 'Parrinello-Rahman'
219 mdp_dict['pcoupltype'] = 'isotropic'
220 mdp_dict['tau-p'] = '1.0'
221 mdp_dict['ref-p'] = '1.0'
222 mdp_dict['compressibility'] = '4.5e-5'
223 mdp_dict['refcoord-scaling'] = 'com'
225 # Dispersion correction
226 if md:
227 mdp_dict['DispCorr'] = 'EnerPres'
229 # Velocity generation
230 if md:
231 if nvt:
232 mdp_dict['gen-vel'] = 'yes'
233 mdp_dict['gen-temp'] = '300'
234 mdp_dict['gen-seed'] = '-1'
235 if npt or free:
236 mdp_dict['gen-vel'] = 'no'
238 # Periodic boundary conditions
239 mdp_dict['pbc'] = 'xyz'
241 return mdp_dict
244def write_mdp(output_mdp_path: str, mdp_dict: Mapping[str, str]):
245 mdp_list = []
246 for k, v in mdp_dict.items():
247 config_parameter_key = str(k).strip().replace('_', '-')
248 if config_parameter_key != 'type':
249 mdp_list.append(config_parameter_key + ' = ' + str(v))
251 with open(output_mdp_path, 'w') as mdp_file:
252 for line in mdp_list:
253 mdp_file.write(line + '\n')
255 return output_mdp_path
258def create_mdp(output_mdp_path: str, input_mdp_path: Optional[str] = None,
259 preset_dict: Optional[Mapping[str, str]] = None,
260 mdp_properties_dict: Optional[Mapping[str, str]] = None) -> str:
261 """Creates an MDP file using the following hierarchy mdp_properties_dict > input_mdp_path > preset_dict"""
262 mdp_dict = {}
264 if preset_dict:
265 for k, v in preset_dict.items():
266 mdp_dict[clean_key(k)] = v
267 if input_mdp_path:
268 input_mdp_dict = read_mdp(input_mdp_path)
269 for k, v in input_mdp_dict.items():
270 mdp_dict[clean_key(k)] = v
271 if mdp_properties_dict:
272 for k, v in mdp_properties_dict.items():
273 mdp_dict[clean_key(k)] = v
275 return write_mdp(output_mdp_path, mdp_dict)
278def clean_key(key: str) -> str:
279 """ Cleans a keyword by converting it to lower case and replacing '_' by '-' as gromacs will not be sensitive to these differences and every keyword
280 can only be defined once in the MDP file. """
282 return key.lower().replace('_', '-')