Coverage for biobb_amber/leap/leap_solvate.py: 71%
160 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-06 14:11 +0000
1#!/usr/bin/env python3
3"""Module containing the LeapSolvate class and the command line interface."""
5import os
6import re
7from pathlib import PurePath
8from typing import List, Optional
10from biobb_common.generic.biobb_object import BiobbObject
11from biobb_common.tools import file_utils as fu
12from biobb_common.tools.file_utils import launchlogger
14from biobb_amber.leap.common import _from_string_to_list
17class LeapSolvate(BiobbObject):
18 """
19 | biobb_amber LeapSolvate
20 | Wrapper of the `AmberTools (AMBER MD Package) leap tool <https://ambermd.org/AmberTools.php>`_ module.
21 | Creates and solvates a system box for an AMBER MD system using tLeap tool from the AmberTools MD package.
23 Args:
24 input_pdb_path (str): Input 3D structure PDB file. File type: input. `Sample file <https://github.com/bioexcel/biobb_amber/raw/master/biobb_amber/test/data/leap/structure.leap.pdb>`_. Accepted formats: pdb (edam:format_1476).
25 input_lib_path (str) (Optional): Input ligand library parameters file. File type: input. `Sample file <https://github.com/bioexcel/biobb_amber/raw/master/biobb_amber/test/data/leap/ligand.lib>`_. Accepted formats: lib (edam:format_3889), zip (edam:format_3987).
26 input_frcmod_path (str) (Optional): Input ligand frcmod parameters file. File type: input. `Sample file <https://github.com/bioexcel/biobb_amber/raw/master/biobb_amber/test/data/leap/ligand.frcmod>`_. Accepted formats: frcmod (edam:format_3888), zip (edam:format_3987).
27 input_params_path (str) (Optional): Additional leap parameter files to load with loadAmberParams Leap command. File type: input. `Sample file <https://github.com/bioexcel/biobb_amber/raw/master/biobb_amber/test/data/leap/frcmod.ionsdang_spce.txt>`_. Accepted formats: in (edam:format_2330), leapin (edam:format_2330), txt (edam:format_2330), zip (edam:format_3987).
28 input_prep_path (str) (Optional): Additional leap parameter files to load with loadAmberPrep Leap command. File type: input. Accepted formats: in (edam:format_2330), leapin (edam:format_2330), txt (edam:format_2330), zip (edam:format_3987).
29 input_source_path (str) (Optional): Additional leap command files to load with source Leap command. File type: input. `Sample file <https://github.com/bioexcel/biobb_amber/raw/master/biobb_amber/test/data/leap/leaprc.water.spce.txt>`_. Accepted formats: in (edam:format_2330), leapin (edam:format_2330), txt (edam:format_2330), zip (edam:format_3987).
30 output_pdb_path (str): Output 3D structure PDB file matching the topology file. File type: output. `Sample file <https://github.com/bioexcel/biobb_amber/raw/master/biobb_amber/test/reference/leap/structure.solv.pdb>`_. Accepted formats: pdb (edam:format_1476).
31 output_top_path (str): Output topology file (AMBER ParmTop). File type: output. `Sample file <https://github.com/bioexcel/biobb_amber/raw/master/biobb_amber/test/reference/leap/structure.solv.top>`_. Accepted formats: top (edam:format_3881), parmtop (edam:format_3881), prmtop (edam:format_3881).
32 output_crd_path (str): Output coordinates file (AMBER crd). File type: output. `Sample file <https://github.com/bioexcel/biobb_amber/raw/master/biobb_amber/test/reference/leap/structure.solv.crd>`_. Accepted formats: crd (edam:format_3878), mdcrd (edam:format_3878), inpcrd (edam:format_3878).
33 properties (dic - Python dictionary object containing the tool parameters, not input/output files):
34 * **forcefield** (*list*) - (["protein.ff14SB","DNA.bsc1","gaff"]) Forcefields to be used for the structure generation. Each item should be either a path to a leaprc file or a string with the leaprc file name if the force field is included with Amber (e.g. "/path/to/leaprc.protein.ff14SB" or "protein.ff14SB"). Default values: ["protein.ff14SB","DNA.bsc1","gaff"].
35 * **water_type** (*str*) - ("TIP3PBOX") Water molecule parameters to be used for the topology. Values: POL3BOX, QSPCFWBOX, SPCBOX, SPCFWBOX, TIP3PBOX, TIP3PFBOX, TIP4PBOX, TIP4PEWBOX, OPCBOX, OPC3BOX, TIP5PBOX.
36 * **box_type** (*str*) - ("truncated_octahedron") Type for the MD system box. Values: cubic, truncated_octahedron.
37 * **ions_type** (*str*) - ("ionsjc_tip3p") Ions type. Values: ionsjc_tip3p, ionsjc_spce, ionsff99_tip3p, ions_charmm22, ionsjc_tip4pew, None.
38 * **neutralise** (*bool*) - ("False") Energetically neutralise the system adding the necessary counterions.
39 * **iso** (*bool*) - ("False") Make the box isometric.
40 * **positive_ions_number** (*int*) - (0) Number of additional positive ions to include in the system box.
41 * **negative_ions_number** (*int*) - (0) Number of additional negative ions to include in the system box.
42 * **positive_ions_type** (*str*) - ("Na+") Type of additional positive ions to include in the system box. Values: Na+,K+.
43 * **negative_ions_type** (*str*) - ("Cl-") Type of additional negative ions to include in the system box. Values: Cl-.
44 * **distance_to_molecule** (*float*) - ("8.0") Size for the MD system box -in Angstroms-, defined such as the minimum distance between any atom originally present in solute and the edge of the periodic box is given by this distance parameter.
45 * **closeness** (*float*) - ("1.0") How close, in Å, solvent ATOMs may come to solute ATOMs.
46 * **binary_path** (*str*) - ("tleap") Path to the tleap executable binary.
47 * **remove_tmp** (*bool*) - (True) [WF property] Remove temporal files.
48 * **restart** (*bool*) - (False) [WF property] Do not execute if output files exist.
49 * **sandbox_path** (*str*) - ("./") [WF property] Parent path to the sandbox directory.
50 * **container_path** (*str*) - (None) Container path definition.
51 * **container_image** (*str*) - ('afandiadib/ambertools:serial') Container image definition.
52 * **container_volume_path** (*str*) - ('/tmp') Container volume path definition.
53 * **container_working_dir** (*str*) - (None) Container working directory definition.
54 * **container_user_id** (*str*) - (None) Container user_id definition.
55 * **container_shell_path** (*str*) - ('/bin/bash') Path to default shell inside the container.
57 Examples:
58 This is a use example of how to use the building block from Python::
59 from biobb_amber.leap.leap_solvate import leap_solvate
60 prop = {
61 'forcefield': ['protein.ff14SB'],
62 'water_type': 'TIP3PBOX',
63 'box_type': 'truncated_octahedron',
64 'neutralise' : True
65 }
66 leap_solvate(input_pdb_path='/path/to/structure.pdb',
67 output_pdb_path='/path/to/newStructure.pdb',
68 output_top_path='/path/to/newTopology.top',
69 output_crd_path='/path/to/newCoordinates.crd',
70 properties=prop)
72 Info:
73 * wrapped_software:
74 * name: AmberTools tLeap
75 * version: >20.9
76 * license: LGPL 2.1
77 * ontology:
78 * name: EDAM
79 * schema: http://edamontology.org/EDAM.owl
81 """
83 def __init__(
84 self,
85 input_pdb_path: str,
86 output_pdb_path: str,
87 output_top_path: str,
88 output_crd_path: str,
89 input_lib_path: Optional[str] = None,
90 input_frcmod_path: Optional[str] = None,
91 input_params_path: Optional[str] = None,
92 input_prep_path: Optional[str] = None,
93 input_source_path: Optional[str] = None,
94 properties: Optional[dict] = None,
95 **kwargs,
96 ):
97 properties = properties or {}
99 # Call parent class constructor
100 super().__init__(properties)
101 self.locals_var_dict = locals().copy()
103 # Input/Output files
104 self.io_dict = {
105 "in": {
106 "input_pdb_path": input_pdb_path,
107 "input_lib_path": input_lib_path,
108 "input_frcmod_path": input_frcmod_path,
109 "input_params_path": input_params_path,
110 "input_prep_path": input_prep_path,
111 "input_source_path": input_source_path,
112 },
113 "out": {
114 "output_pdb_path": output_pdb_path,
115 "output_top_path": output_top_path,
116 "output_crd_path": output_crd_path,
117 },
118 }
120 # # Ligand Parameter lists
121 # self.ligands_lib_list = []
122 # if input_lib_path:
123 # self.ligands_lib_list.append(input_lib_path)
124 #
125 # self.ligands_frcmod_list = []
126 # if input_frcmod_path:
127 # self.ligands_frcmod_list.append(input_frcmod_path)
129 # Properties specific for BB
130 self.properties = properties
132 # Set default forcefields
133 if self.container_path:
134 self.forcefield = _from_string_to_list(
135 properties.get("forcefield", ['leaprc.protein.ff14SB', 'leaprc.DNA.bsc1', 'leaprc.gaff'])
136 )
137 else:
138 amber_home_path = os.getenv("AMBERHOME")
139 protein_ff14SB_path = os.path.join(amber_home_path, 'dat', 'leap', 'cmd', 'leaprc.protein.ff14SB')
140 dna_bsc1_path = os.path.join(amber_home_path, 'dat', 'leap', 'cmd', 'leaprc.DNA.bsc1')
141 gaff_path = os.path.join(amber_home_path, 'dat', 'leap', 'cmd', 'leaprc.gaff')
143 self.forcefield = _from_string_to_list(
144 properties.get("forcefield", [protein_ff14SB_path, dna_bsc1_path, gaff_path])
145 )
147 # Find the paths of the leaprc files if only the force field names are provided
148 self.forcefield = self.find_leaprc_paths(self.forcefield)
150 self.water_type = properties.get("water_type", "TIP3PBOX")
151 self.box_type = properties.get("box_type", "truncated_octahedron")
152 self.ions_type = properties.get("ions_type", "ionsjc_tip3p")
153 self.neutralise = properties.get("neutralise", False)
154 self.iso = properties.get("iso", False)
155 self.positive_ions_number = properties.get("positive_ions_number", 0)
156 self.positive_ions_type = properties.get("positive_ions_type", "Na+")
157 self.negative_ions_number = properties.get("negative_ions_number", 0)
158 self.negative_ions_type = properties.get("negative_ions_type", "Cl-")
159 self.distance_to_molecule = properties.get("distance_to_molecule", 8.0)
160 self.closeness = properties.get("closeness", 1.0)
161 self.binary_path = properties.get("binary_path", "tleap")
163 # Check the properties
164 self.check_properties(properties)
165 self.check_arguments()
167 def find_leaprc_paths(self, forcefields: List[str]) -> List[str]:
168 """
169 Find the leaprc paths for the force fields provided.
171 For each item in the forcefields list, the function checks if the str is a path to an existing file.
172 If not, it tries to find the file in the $AMBERHOME/dat/leap/cmd/ directory or the $AMBERHOME/dat/leap/cmd/oldff/
173 directory with and without the leaprc prefix.
175 Args:
176 forcefields (List[str]): List of force fields to find the leaprc files for.
178 Returns:
179 List[str]: List of leaprc file paths.
180 """
182 leaprc_paths = []
184 for forcefield in forcefields:
186 num_paths = len(leaprc_paths)
188 # Check if the forcefield is a path to an existing file
189 if os.path.exists(forcefield):
190 leaprc_paths.append(forcefield)
191 continue
193 # Check if the forcefield is in the leaprc directory
194 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', f"leaprc.{forcefield}")
195 if os.path.exists(leaprc_path):
196 leaprc_paths.append(leaprc_path)
197 continue
199 # Check if the forcefield is in the oldff directory
200 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', 'oldff', f"leaprc.{forcefield}")
201 if os.path.exists(leaprc_path):
202 leaprc_paths.append(leaprc_path)
203 continue
205 # Check if the forcefield is in the leaprc directory without the leaprc prefix
206 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', f"{forcefield}")
207 if os.path.exists(leaprc_path):
208 leaprc_paths.append(leaprc_path)
209 continue
211 # Check if the forcefield is in the oldff directory without the leaprc prefix
212 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', 'oldff', f"{forcefield}")
213 if os.path.exists(leaprc_path):
214 leaprc_paths.append(leaprc_path)
215 continue
217 new_num_paths = len(leaprc_paths)
219 if new_num_paths == num_paths:
220 raise ValueError(f"Force field {forcefield} not found. Check the $AMBERHOME/dat/leap/cmd/ directory for available force fields or provide the path to an existing leaprc file.")
222 return leaprc_paths
224 # def check_data_params(self, out_log, err_log):
225 # """ Checks input/output paths correctness """
227 # # Check input(s)
228 # self.io_dict["in"]["input_pdb_path"] = check_input_path(self.io_dict["in"]["input_pdb_path"], "input_pdb_path", False, out_log, self.__class__.__name__)
229 # self.io_dict["in"]["input_lib_path"] = check_input_path(self.io_dict["in"]["input_lib_path"], "input_lib_path", True, out_log, self.__class__.__name__)
230 # self.io_dict["in"]["input_frcmod_path"] = check_input_path(self.io_dict["in"]["input_frcmod_path"], "input_frcmod_path", True, out_log, self.__class__.__name__)
231 # # self.io_dict["in"]["input_params_path"] = check_input_path(self.io_dict["in"]["input_params_path"], "input_params_path", True, out_log, self.__class__.__name__)
232 # # self.io_dict["in"]["input_source_path"] = check_input_path(self.io_dict["in"]["input_source_path"], "input_source_path", True, out_log, self.__class__.__name__)
234 # # Check output(s)
235 # self.io_dict["out"]["output_pdb_path"] = check_output_path(self.io_dict["out"]["output_pdb_path"], "output_pdb_path", False, out_log, self.__class__.__name__)
236 # self.io_dict["out"]["output_top_path"] = check_output_path(self.io_dict["out"]["output_top_path"], "output_top_path", False, out_log, self.__class__.__name__)
237 # self.io_dict["out"]["output_crd_path"] = check_output_path(self.io_dict["out"]["output_crd_path"], "output_crd_path", False, out_log, self.__class__.__name__)
239 def _resolve_input_files(self, io_key: str, tmp_folder: Optional[str]) -> List[str]:
240 """Resolve an optional input (zip or single file) to the list of paths the tleap
241 script must load.
243 In container mode the zipball is staged into the host sandbox (unique_dir), which is
244 bind-mounted at ``container_volume_path``; it is unpacked on the host and the entries
245 are referenced by their container path. In local mode the zipball is unpacked into a
246 temporary folder and referenced by its host path.
247 """
248 file_path = self.io_dict["in"][io_key]
249 if file_path is None:
250 return []
251 if not file_path.endswith(".zip"):
252 return [self.stage_io_dict["in"][io_key]]
253 if self.container_path:
254 unique_dir = self.stage_io_dict["unique_dir"]
255 host_zip = os.path.join(unique_dir, os.path.basename(file_path))
256 extracted = fu.unzip_list(host_zip, dest_dir=unique_dir, out_log=self.out_log)
257 return [
258 os.path.join(self.container_volume_path, os.path.relpath(p, unique_dir))
259 for p in extracted
260 ]
261 return fu.unzip_list(
262 self.stage_io_dict["in"][io_key], dest_dir=tmp_folder, out_log=self.out_log
263 )
265 @launchlogger
266 def launch(self):
267 """Launches the execution of the LeapSolvate module."""
269 # check input/output paths and parameters
270 # self.check_data_params(self.out_log, self.err_log)
272 # Setup Biobb
273 if self.check_restart():
274 return 0
275 self.stage_files()
277 box_command = "solvateOct"
278 if self.box_type == "cubic":
279 box_command = "solvateBox"
281 # Forcefield
282 # source_ff_command = "source leaprc." + self.forcefield
284 # Water Type
285 # leaprc.water.tip4pew, tip4pd, tip3p, spceb, spce, opc, fb4, fb3
286 # Values: POL3BOX, QSPCFWBOX, SPCBOX, SPCFWBOX, TIP3PBOX, TIP3PFBOX, TIP4PBOX, TIP4PEWBOX, OPCBOX, OPC3BOX, TIP5PBOX.
287 source_wat_command = "source leaprc.water.tip3p"
288 if self.water_type == "TIP4PEWBOX":
289 source_wat_command = "leaprc.water.tip4pew"
290 if self.water_type == "TIP4PBOX":
291 source_wat_command = "leaprc.water.tip4pd"
292 if re.match(r"SPC", self.water_type):
293 source_wat_command = "source leaprc.water.spce"
294 if re.match(r"OPC", self.water_type):
295 source_wat_command = "source leaprc.water.opc"
297 # Counterions
298 ions_command = ""
299 if self.neutralise:
300 ions_command = (
301 ions_command + "addions mol " + self.negative_ions_type + " 0 \n"
302 )
303 ions_command = (
304 ions_command + "addions mol " + self.positive_ions_type + " 0 \n"
305 )
307 if self.negative_ions_number != 0:
308 ions_command = (
309 ions_command + "addions mol " + self.negative_ions_type + " " + str(self.negative_ions_number) + " \n"
310 )
311 if self.positive_ions_number != 0:
312 ions_command = (
313 ions_command + "addions mol " + self.positive_ions_type + " " + str(self.positive_ions_number) + " \n"
314 )
316 # Creating temporary folder & Leap configuration (instructions) file
317 if self.container_path:
318 instructions_file = str(
319 PurePath(self.stage_io_dict["unique_dir"]).joinpath("leap.in")
320 )
321 instructions_file_path = str(
322 PurePath(self.container_volume_path).joinpath("leap.in")
323 )
324 tmp_folder = None
325 else:
326 tmp_folder = fu.create_unique_dir()
327 instructions_file = str(PurePath(tmp_folder).joinpath("leap.in"))
328 fu.log("Creating %s temporary folder" % tmp_folder, self.out_log)
329 instructions_file_path = instructions_file
331 ligands_lib_list = self._resolve_input_files("input_lib_path", tmp_folder)
332 ligands_frcmod_list = self._resolve_input_files("input_frcmod_path", tmp_folder)
333 amber_params_list = self._resolve_input_files("input_params_path", tmp_folder)
334 amber_prep_list = self._resolve_input_files("input_prep_path", tmp_folder)
335 leap_source_list = self._resolve_input_files("input_source_path", tmp_folder)
337 with open(instructions_file, "w") as leapin:
338 # Forcefields loaded by default:
339 # Protein: ff14SB (PARM99 + frcmod.ff99SB + frcmod.parmbsc0 + OL3 for RNA)
340 # leapin.write("source leaprc.protein.ff14SB \n")
341 # DNA: parmBSC1 (ParmBSC1 (ff99 + bsc0 + bsc1) for DNA. Ivani et al. Nature Methods 13: 55, 2016)
342 # leapin.write("source leaprc.DNA.bsc1 \n")
343 # Ligands: GAFF (General Amber Force field, J. Comput. Chem. 2004 Jul 15;25(9):1157-74)
344 # leapin.write("source leaprc.gaff \n")
346 # Forcefields loaded from input forcefield property
347 for t in self.forcefield:
348 if self.container_path:
349 leapin.write("source leaprc.{}\n".format(t))
350 else:
351 leapin.write("source {}\n".format(t))
353 # Additional Leap commands
354 for leap_commands in leap_source_list:
355 leapin.write("source " + leap_commands + "\n")
357 # Ions Type
358 if self.ions_type != "None":
359 leapin.write("loadamberparams frcmod." + self.ions_type + "\n")
361 # Additional Amber parameters
362 for amber_params in amber_params_list:
363 leapin.write("loadamberparams " + amber_params + "\n")
365 # Additional Amber prep files
366 for amber_prep in amber_prep_list:
367 leapin.write("loadamberprep " + amber_prep + "\n")
369 # Water Model loaded from input water_model property
370 leapin.write(source_wat_command + " \n")
372 # Ligand(s) libraries (if any)
373 for amber_lib in ligands_lib_list:
374 leapin.write("loadOff " + amber_lib + "\n")
375 for amber_frcmod in ligands_frcmod_list:
376 leapin.write("loadamberparams " + amber_frcmod + "\n")
378 # Loading PDB file
379 leapin.write(
380 "mol = loadpdb " + self.stage_io_dict["in"]["input_pdb_path"] + " \n"
381 )
383 # Generating box + adding water molecules
384 leapin.write(
385 box_command + " mol " + self.water_type + " " + str(self.distance_to_molecule)
386 )
387 leapin.write(
388 " iso " + str(self.closeness) + "\n"
389 ) if self.iso else leapin.write(" " + str(self.closeness) + "\n")
391 # Adding counterions
392 leapin.write(ions_command)
394 # Saving output PDB file, coordinates and topology
395 leapin.write(
396 "savepdb mol " + self.stage_io_dict["out"]["output_pdb_path"] + " \n"
397 )
398 leapin.write(
399 "saveAmberParm mol " + self.stage_io_dict["out"]["output_top_path"] + " " + self.stage_io_dict["out"]["output_crd_path"] + "\n"
400 )
401 leapin.write("quit \n")
403 # Command line
404 self.cmd = [self.binary_path, "-f", instructions_file_path]
406 # Run Biobb block
407 self.run_biobb()
409 # Copy files to host
410 self.copy_to_host()
412 # Saving octahedron box with all decimals in PDB file. Needed for the add_ions BB.
414 # Getting octahedron box from generated crd file
415 with open(self.io_dict["out"]["output_crd_path"], "r") as file:
416 for line in file:
417 pass
419 # Adding box as a first line in the generated pdb file with OCTBOX tag
420 octbox = "OCTBOX " + line
421 with open(self.io_dict["out"]["output_pdb_path"], "r+") as f:
422 content = f.read()
423 f.seek(0, 0)
424 f.write(octbox + content)
426 # remove temporary folder(s)
427 self.tmp_files.extend([str(tmp_folder), "leap.log"])
428 self.remove_tmp_files()
430 self.check_arguments(output_files_created=True, raise_exception=False)
432 return self.return_code
435def leap_solvate(
436 input_pdb_path: str,
437 output_pdb_path: str,
438 output_top_path: str,
439 output_crd_path: str,
440 input_lib_path: Optional[str] = None,
441 input_frcmod_path: Optional[str] = None,
442 input_params_path: Optional[str] = None,
443 input_prep_path: Optional[str] = None,
444 input_source_path: Optional[str] = None,
445 properties: Optional[dict] = None,
446 **kwargs,
447) -> int:
448 """Create the :class:`LeapSolvate <leap.leap_solvate.LeapSolvate>` class and
449 execute the :meth:`launch() <leap.leap_solvate.LeapSolvate.launch>` method."""
450 return LeapSolvate(**dict(locals())).launch()
453leap_solvate.__doc__ = LeapSolvate.__doc__
454main = LeapSolvate.get_main(leap_solvate, "Generating and solvating a system box for an AMBER MD system. using tLeap program from AmberTools MD package.")
456if __name__ == "__main__":
457 main()