Coverage for biobb_amber/leap/leap_add_ions.py: 60%
219 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 LeapAddIons class and the command line interface."""
5import os
6import re
7import shutil
8from decimal import Decimal
9from pathlib import PurePath
10from typing import List, Optional
12from biobb_common.generic.biobb_object import BiobbObject
13from biobb_common.tools import file_utils as fu
14from biobb_common.tools.file_utils import launchlogger
16from biobb_amber.leap.common import _from_string_to_list
19class LeapAddIons(BiobbObject):
20 """
21 | biobb_amber LeapAddIons
22 | Wrapper of the `AmberTools (AMBER MD Package) leap tool <https://ambermd.org/AmberTools.php>`_ module.
23 | Adds counterions to a system box for an AMBER MD system using tLeap tool from the AmberTools MD package.
25 Args:
26 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.ions.input.pdb>`_. Accepted formats: pdb (edam:format_1476).
27 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).
28 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).
29 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).
30 input_prep_path (str) (Optional): Additional leap parameter files to load with loadAmberPrep 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).
31 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).
32 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.ions.pdb>`_. Accepted formats: pdb (edam:format_1476).
33 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.ions.top>`_. Accepted formats: top (edam:format_3881), parmtop (edam:format_3881), prmtop (edam:format_3881).
34 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.ions.crd>`_. Accepted formats: crd (edam:format_3878), mdcrd (edam:format_3878), inpcrd (edam:format_3878).
35 properties (dic - Python dictionary object containing the tool parameters, not input/output files):
36 * **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"].
37 * **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.
38 * **box_type** (*str*) - ("truncated_octahedron") Type for the MD system box. Values: cubic, truncated_octahedron.
39 * **ions_type** (*str*) - ("ionsjc_tip3p") Ions type. Values: ionsjc_tip3p, ionsjc_spce, ionsff99_tip3p, ions_charmm22, ionsjc_tip4pew, None.
40 * **neutralise** (*bool*) - ("True") Energetically neutralise the system adding the necessary counterions.
41 * **ionic_concentration** (*float*) - (50) Additional ionic concentration to include in the system box. Units in mM/L.
42 * **positive_ions_number** (*int*) - (0) Number of additional positive ions to include in the system box.
43 * **negative_ions_number** (*int*) - (0) Number of additional negative ions to include in the system box.
44 * **positive_ions_type** (*str*) - ("Na+") Type of additional positive ions to include in the system box. Values: Na+,K+.
45 * **negative_ions_type** (*str*) - ("Cl-") Type of additional negative ions to include in the system box. Values: Cl-.
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::
60 from biobb_amber.leap.leap_add_ions import leap_add_ions
61 prop = {
62 'forcefield': ['protein.ff14SB'],
63 'water_type': 'TIP3PBOX',
64 'neutralise' : True
65 }
66 leap_add_ions(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)
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", True)
154 self.ionic_concentration = properties.get("ionic_concentration", 50)
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.binary_path = properties.get("binary_path", "tleap")
161 # Check the properties
162 self.check_properties(properties)
163 self.check_arguments()
165 def find_leaprc_paths(self, forcefields: List[str]) -> List[str]:
166 """
167 Find the leaprc paths for the force fields provided.
169 For each item in the forcefields list, the function checks if the str is a path to an existing file.
170 If not, it tries to find the file in the $AMBERHOME/dat/leap/cmd/ directory or the $AMBERHOME/dat/leap/cmd/oldff/
171 directory with and without the leaprc prefix.
173 Args:
174 forcefields (List[str]): List of force fields to find the leaprc files for.
176 Returns:
177 List[str]: List of leaprc file paths.
178 """
180 leaprc_paths = []
182 for forcefield in forcefields:
184 num_paths = len(leaprc_paths)
186 # Check if the forcefield is a path to an existing file
187 if os.path.exists(forcefield):
188 leaprc_paths.append(forcefield)
189 continue
191 # Check if the forcefield is in the leaprc directory
192 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', f"leaprc.{forcefield}")
193 if os.path.exists(leaprc_path):
194 leaprc_paths.append(leaprc_path)
195 continue
197 # Check if the forcefield is in the oldff directory
198 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', 'oldff', f"leaprc.{forcefield}")
199 if os.path.exists(leaprc_path):
200 leaprc_paths.append(leaprc_path)
201 continue
203 # Check if the forcefield is in the leaprc directory without the leaprc prefix
204 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', f"{forcefield}")
205 if os.path.exists(leaprc_path):
206 leaprc_paths.append(leaprc_path)
207 continue
209 # Check if the forcefield is in the oldff directory without the leaprc prefix
210 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', 'oldff', f"{forcefield}")
211 if os.path.exists(leaprc_path):
212 leaprc_paths.append(leaprc_path)
213 continue
215 new_num_paths = len(leaprc_paths)
217 if new_num_paths == num_paths:
218 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.")
220 return leaprc_paths
222 # def check_data_params(self, out_log, err_log):
223 # """ Checks input/output paths correctness """
225 # # Check input(s)
226 # 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__)
227 # 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__)
228 # 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__)
229 # # 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__)
230 # # 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__)
232 # # Check output(s)
233 # 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__)
234 # 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__)
235 # 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__)
237 def find_out_number_of_ions(self):
238 """Computes the number of positive and negative ions from the input ionic
239 concentration and the number of water molecules in the system box."""
241 # Finding number of water molecules in the box
242 cnt = 0
243 with open(self.io_dict["in"]["input_pdb_path"]) as fp:
244 for line in fp:
245 # Water molecules are identified with different resids
246 # by different forcefields / MD packages
247 pq = re.compile((r"WAT|HOH|TP3|TIP3|SOL"), re.M)
248 if pq.search(line):
249 # Incrementing number of water molecules just in the water
250 # oxygen atom, ignoring hydrogen or dummy atoms
251 pq2 = re.compile(r"H|EPW|EP1|EP2", re.M)
252 if not pq2.search(line):
253 cnt = cnt + 1
255 # To get to X mM ions we need
256 # n(NaCl) = #waters / 55 Mol * X M
257 # where X M = X mM / 1000
258 self.nio = int((cnt / 55) * (self.ionic_concentration / 1000))
260 def _resolve_input_files(self, io_key: str, tmp_folder: Optional[str]) -> List[str]:
261 """Resolve an optional input (zip or single file) to the list of paths the tleap
262 script must load.
264 In container mode the zipball is staged into the host sandbox (unique_dir), which is
265 bind-mounted at ``container_volume_path``; it is unpacked on the host and the entries
266 are referenced by their container path. In local mode the zipball is unpacked into a
267 temporary folder and referenced by its host path.
268 """
269 file_path = self.io_dict["in"][io_key]
270 if file_path is None:
271 return []
272 if not file_path.endswith(".zip"):
273 return [self.stage_io_dict["in"][io_key]]
274 if self.container_path:
275 unique_dir = self.stage_io_dict["unique_dir"]
276 host_zip = os.path.join(unique_dir, os.path.basename(file_path))
277 extracted = fu.unzip_list(host_zip, dest_dir=unique_dir, out_log=self.out_log)
278 return [
279 os.path.join(self.container_volume_path, os.path.relpath(p, unique_dir))
280 for p in extracted
281 ]
282 return fu.unzip_list(
283 self.stage_io_dict["in"][io_key], dest_dir=tmp_folder, out_log=self.out_log
284 )
286 @launchlogger
287 def launch(self):
288 """Launches the execution of the LeapAddIons module."""
290 # check input/output paths and parameters
291 # self.check_data_params(self.out_log, self.err_log)
293 # Setup Biobb
294 if self.check_restart():
295 return 0
296 self.stage_files()
298 # Water Type
299 # leaprc.water.tip4pew, tip4pd, tip3p, spceb, spce, opc, fb4, fb3
300 # Values: POL3BOX, QSPCFWBOX, SPCBOX, SPCFWBOX, TIP3PBOX, TIP3PFBOX, TIP4PBOX, TIP4PEWBOX, OPCBOX, OPC3BOX, TIP5PBOX.
301 source_wat_command = "source leaprc.water.tip3p"
302 if self.water_type == "TIP4PEWBOX":
303 source_wat_command = "leaprc.water.tip4pew"
304 if self.water_type == "TIP4PBOX":
305 source_wat_command = "leaprc.water.tip4pd"
306 if re.match(r"SPC", self.water_type):
307 source_wat_command = "source leaprc.water.spce"
308 if re.match(r"OPC", self.water_type):
309 source_wat_command = "source leaprc.water.opc"
311 # Counterions
312 ions_command = ""
313 if self.neutralise:
314 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " 0 \n"
315 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " 0 \n"
316 ions_command = (
317 ions_command + "addionsRand mol " + self.negative_ions_type + " 0 \n"
318 )
319 ions_command = (
320 ions_command + "addionsRand mol " + self.positive_ions_type + " 0 \n"
321 )
323 if (
324 self.ionic_concentration and self.negative_ions_number == 0 and self.positive_ions_number == 0
325 ):
326 self.find_out_number_of_ions()
327 nneg = self.nio # Update with function
328 npos = self.nio # Update with function
329 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " " + str(nneg) + " \n"
330 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " " + str(npos) + " \n"
331 ions_command = (
332 ions_command + "addionsRand mol " + self.negative_ions_type + " " + str(nneg) + " \n"
333 )
334 ions_command = (
335 ions_command + "addionsRand mol " + self.positive_ions_type + " " + str(npos) + " \n"
336 )
337 else:
338 if self.negative_ions_number != 0:
339 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " " + str(self.negative_ions_number) + " \n"
340 ions_command = (
341 ions_command + "addionsRand mol " + self.negative_ions_type + " " + str(self.negative_ions_number) + " \n"
342 )
343 if self.positive_ions_number != 0:
344 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " " + str(self.positive_ions_number) + " \n"
345 ions_command = (
346 ions_command + "addionsRand mol " + self.positive_ions_type + " " + str(self.positive_ions_number) + " \n"
347 )
349 # Creating temporary folder & Leap configuration (instructions) file
350 if self.container_path:
351 instructions_file = str(
352 PurePath(self.stage_io_dict["unique_dir"]).joinpath("leap.in")
353 )
354 instructions_file_path = str(
355 PurePath(self.container_volume_path).joinpath("leap.in")
356 )
357 tmp_folder = None
358 else:
359 tmp_folder = fu.create_unique_dir()
360 instructions_file = str(PurePath(tmp_folder).joinpath("leap.in"))
361 fu.log("Creating %s temporary folder" % tmp_folder, self.out_log)
362 instructions_file_path = instructions_file
364 ligands_lib_list = self._resolve_input_files("input_lib_path", tmp_folder)
365 ligands_frcmod_list = self._resolve_input_files("input_frcmod_path", tmp_folder)
366 amber_params_list = self._resolve_input_files("input_params_path", tmp_folder)
367 amber_prep_list = self._resolve_input_files("input_prep_path", tmp_folder)
368 leap_source_list = self._resolve_input_files("input_source_path", tmp_folder)
370 # instructions_file = str(PurePath(tmp_folder).joinpath("leap.in"))
371 with open(instructions_file, "w") as leapin:
372 # Forcefields loaded by default:
373 # Protein: ff14SB (PARM99 + frcmod.ff99SB + frcmod.parmbsc0 + OL3 for RNA)
374 # leapin.write("source leaprc.protein.ff14SB \n")
375 # DNA: parmBSC1 (ParmBSC1 (ff99 + bsc0 + bsc1) for DNA. Ivani et al. Nature Methods 13: 55, 2016)
376 # leapin.write("source leaprc.DNA.bsc1 \n")
377 # Ligands: GAFF (General Amber Force field, J. Comput. Chem. 2004 Jul 15;25(9):1157-74)
378 # leapin.write("source leaprc.gaff \n")
380 # Forcefields loaded from input forcefield property
381 for t in self.forcefield:
382 if self.container_path:
383 leapin.write("source leaprc.{}\n".format(t))
384 else:
385 leapin.write("source {}\n".format(t))
387 # Additional Leap commands
388 for leap_commands in leap_source_list:
389 leapin.write("source " + leap_commands + "\n")
391 # Water Model loaded from input water_model property
392 leapin.write(source_wat_command + " \n")
394 # Ions Type
395 if self.ions_type != "None":
396 leapin.write("loadamberparams frcmod." + self.ions_type + "\n")
398 # Additional Amber parameters
399 for amber_params in amber_params_list:
400 leapin.write("loadamberparams " + amber_params + "\n")
402 # Additional Amber prep files
403 for amber_prep in amber_prep_list:
404 leapin.write("loadamberprep " + amber_prep + "\n")
406 # Ligand(s) libraries (if any)
407 for amber_lib in ligands_lib_list:
408 leapin.write("loadOff " + amber_lib + "\n")
409 for amber_frcmod in ligands_frcmod_list:
410 leapin.write("loadamberparams " + amber_frcmod + "\n")
412 # Loading PDB file
413 leapin.write(
414 "mol = loadpdb " + self.stage_io_dict["in"]["input_pdb_path"] + " \n"
415 )
417 # Adding ions
418 leapin.write(ions_command)
420 # Generating box
421 leapin.write("setBox mol vdw \n")
423 # Saving output PDB file, coordinates and topology
424 leapin.write(
425 "savepdb mol " + self.stage_io_dict["out"]["output_pdb_path"] + " \n"
426 )
427 leapin.write(
428 "saveAmberParm mol " + self.stage_io_dict["out"]["output_top_path"] + " " + self.stage_io_dict["out"]["output_crd_path"] + "\n"
429 )
430 leapin.write("quit \n")
432 # Container path
433 if self.container_path:
434 if not self.container_working_dir:
435 fu.log('WARNING: container_working_dir property was not set. Defining it with the same value as container_volume_path', self.out_log, self.global_log)
436 self.container_working_dir = self.container_volume_path
438 # Command line
439 self.cmd = [self.binary_path, "-f", instructions_file_path]
441 # Run Biobb block
442 self.run_biobb()
444 # Copy files to host
445 self.copy_to_host()
447 if self.box_type != "cubic":
448 fu.log(
449 "Fixing truncated octahedron Box in the topology and coordinates files",
450 self.out_log,
451 self.global_log,
452 )
454 # Taking box info from input PDB file, CRYST1 tag (first line)
455 with open(self.io_dict["in"]["input_pdb_path"]) as file:
456 lines = file.readlines()
457 pdb_line = lines[0]
459 if "OCTBOX" not in pdb_line:
460 fu.log(
461 "WARNING: box info not found in input PDB file (OCTBOX). Needed to correctly assign the octahedron box. Assuming cubic box.",
462 self.out_log,
463 self.global_log,
464 )
465 else:
466 # PDB info: CRYST1 86.316 86.316 86.316 109.47 109.47 109.47 P 1
467 # PDB info: OCTBOX 86.1942924 86.1942924 86.1942924 109.4712190 109.4712190 109.4712190
468 # regex_box = 'CRYST1\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*P 1'
469 regex_box = r"OCTBOX\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)"
470 box = re.findall(regex_box, pdb_line)[0]
471 box_line = ""
472 for coord in box:
473 box_line += "{:12.7f}".format(float(coord))
475 # PRMTOP info: 1.09471219E+02 8.63157502E+01 8.63157502E+01 8.63157502E+01
476 top_box_line = ""
477 top_box_line += " %.8E" % Decimal(float(box[3]))
478 top_box_line += " %.8E" % Decimal(float(box[0]))
479 top_box_line += " %.8E" % Decimal(float(box[1]))
480 top_box_line += " %.8E" % Decimal(float(box[2]))
482 # Removing box generated by tleap from the crd file (last line)
483 with open(self.io_dict["out"]["output_crd_path"]) as file:
484 lines = file.readlines()
485 crd_lines = lines[:-1]
487 # Adding old box coordinates (taken from the input pdb)
488 crd_lines.append(box_line)
490 with open(self.io_dict["out"]["output_crd_path"], "w") as file:
491 for line in crd_lines:
492 file.write(str(line))
493 file.write("\n")
495 # Now fixing IFBOX param in prmtop.
496 box_flag = False
497 ifbox_flag = 0
498 # %FLAG BOX_DIMENSIONS
499 # %FORMAT(5E16.8)
500 # 1.09471219E+02 8.63157502E+01 8.63157502E+01 8.63157502E+01
502 tmp_parmtop = str(
503 # PurePath(str(tmp_folder)).joinpath("top_temp.parmtop")
504 PurePath(self.stage_io_dict["unique_dir"]).joinpath("top_temp.parmtop")
505 )
506 shutil.copyfile(self.io_dict["out"]["output_top_path"], tmp_parmtop)
508 with open(self.io_dict["out"]["output_top_path"], "w") as new_top:
509 with open(tmp_parmtop, "r") as old_top:
510 for line in old_top:
511 if "BOX_DIMENSIONS" in line:
512 box_flag = True
513 new_top.write(line)
514 elif box_flag and "FORMAT" not in line:
515 new_top.write(top_box_line + "\n")
516 box_flag = False
517 elif (
518 "FLAG POINTERS" in line or ifbox_flag == 1 or ifbox_flag == 2 or ifbox_flag == 3
519 ):
520 ifbox_flag += 1
521 new_top.write(line)
522 elif ifbox_flag == 4:
523 # new_top.write(top_box_line + "\n")
524 new_top.write(line[:56] + " 2" + line[64:])
525 ifbox_flag += 1
526 else:
527 new_top.write(line)
529 # remove temporary folder(s)
530 self.tmp_files.extend([str(tmp_folder), "leap.log"])
531 self.remove_tmp_files()
533 self.check_arguments(output_files_created=True, raise_exception=False)
535 return self.return_code
538def leap_add_ions(
539 input_pdb_path: str,
540 output_pdb_path: str,
541 output_top_path: str,
542 output_crd_path: str,
543 input_lib_path: Optional[str] = None,
544 input_frcmod_path: Optional[str] = None,
545 input_params_path: Optional[str] = None,
546 input_prep_path: Optional[str] = None,
547 input_source_path: Optional[str] = None,
548 properties: Optional[dict] = None,
549 **kwargs,
550) -> int:
551 """Create the :class:`LeapAddIons <leap.leap_add_ions.LeapAddIons>` class and
552 execute the :meth:`launch() <leap.leap_add_ions.LeapAddIons.launch>` method."""
553 return LeapAddIons(**dict(locals())).launch()
556leap_add_ions.__doc__ = LeapAddIons.__doc__
557main = LeapAddIons.get_main(leap_add_ions, "Adds counterions to a system box for an AMBER MD system using tLeap.")
559if __name__ == "__main__":
560 main()