Coverage for biobb_amber / leap / leap_add_ions.py: 59%
219 statements
« prev ^ index » next coverage.py v7.13.0, created at 2025-12-22 10:27 +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.solv.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 # Set default forcefields
130 amber_home_path = os.getenv("AMBERHOME")
131 protein_ff14SB_path = os.path.join(amber_home_path, 'dat', 'leap', 'cmd', 'leaprc.protein.ff14SB')
132 dna_bsc1_path = os.path.join(amber_home_path, 'dat', 'leap', 'cmd', 'leaprc.DNA.bsc1')
133 gaff_path = os.path.join(amber_home_path, 'dat', 'leap', 'cmd', 'leaprc.gaff')
135 # Properties specific for BB
136 self.properties = properties
137 self.forcefield = _from_string_to_list(
138 properties.get("forcefield", [protein_ff14SB_path, dna_bsc1_path, gaff_path])
139 )
140 # Find the paths of the leaprc files if only the force field names are provided
141 self.forcefield = self.find_leaprc_paths(self.forcefield)
142 self.water_type = properties.get("water_type", "TIP3PBOX")
143 self.box_type = properties.get("box_type", "truncated_octahedron")
144 self.ions_type = properties.get("ions_type", "ionsjc_tip3p")
145 self.neutralise = properties.get("neutralise", True)
146 self.ionic_concentration = properties.get("ionic_concentration", 50)
147 self.positive_ions_number = properties.get("positive_ions_number", 0)
148 self.positive_ions_type = properties.get("positive_ions_type", "Na+")
149 self.negative_ions_number = properties.get("negative_ions_number", 0)
150 self.negative_ions_type = properties.get("negative_ions_type", "Cl-")
151 self.binary_path = properties.get("binary_path", "tleap")
153 # Check the properties
154 self.check_properties(properties)
155 self.check_arguments()
157 def find_leaprc_paths(self, forcefields: List[str]) -> List[str]:
158 """
159 Find the leaprc paths for the force fields provided.
161 For each item in the forcefields list, the function checks if the str is a path to an existing file.
162 If not, it tries to find the file in the $AMBERHOME/dat/leap/cmd/ directory or the $AMBERHOME/dat/leap/cmd/oldff/
163 directory with and without the leaprc prefix.
165 Args:
166 forcefields (List[str]): List of force fields to find the leaprc files for.
168 Returns:
169 List[str]: List of leaprc file paths.
170 """
172 leaprc_paths = []
174 for forcefield in forcefields:
176 num_paths = len(leaprc_paths)
178 # Check if the forcefield is a path to an existing file
179 if os.path.exists(forcefield):
180 leaprc_paths.append(forcefield)
181 continue
183 # Check if the forcefield is in the leaprc directory
184 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', f"leaprc.{forcefield}")
185 if os.path.exists(leaprc_path):
186 leaprc_paths.append(leaprc_path)
187 continue
189 # Check if the forcefield is in the oldff directory
190 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', 'oldff', f"leaprc.{forcefield}")
191 if os.path.exists(leaprc_path):
192 leaprc_paths.append(leaprc_path)
193 continue
195 # Check if the forcefield is in the leaprc directory without the leaprc prefix
196 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', f"{forcefield}")
197 if os.path.exists(leaprc_path):
198 leaprc_paths.append(leaprc_path)
199 continue
201 # Check if the forcefield is in the oldff directory without the leaprc prefix
202 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', 'oldff', f"{forcefield}")
203 if os.path.exists(leaprc_path):
204 leaprc_paths.append(leaprc_path)
205 continue
207 new_num_paths = len(leaprc_paths)
209 if new_num_paths == num_paths:
210 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.")
212 return leaprc_paths
214 # def check_data_params(self, out_log, err_log):
215 # """ Checks input/output paths correctness """
217 # # Check input(s)
218 # 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__)
219 # 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__)
220 # 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__)
221 # # 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__)
222 # # 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__)
224 # # Check output(s)
225 # 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__)
226 # 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__)
227 # 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__)
229 def find_out_number_of_ions(self):
230 """Computes the number of positive and negative ions from the input ionic
231 concentration and the number of water molecules in the system box."""
233 # Finding number of water molecules in the box
234 cnt = 0
235 with open(self.io_dict["in"]["input_pdb_path"]) as fp:
236 for line in fp:
237 # Water molecules are identified with different resids
238 # by different forcefields / MD packages
239 pq = re.compile((r"WAT|HOH|TP3|TIP3|SOL"), re.M)
240 if pq.search(line):
241 # Incrementing number of water molecules just in the water
242 # oxygen atom, ignoring hydrogen or dummy atoms
243 pq2 = re.compile(r"H|EPW|EP1|EP2", re.M)
244 if not pq2.search(line):
245 cnt = cnt + 1
247 # To get to X mM ions we need
248 # n(NaCl) = #waters / 55 Mol * X M
249 # where X M = X mM / 1000
250 self.nio = int((cnt / 55) * (self.ionic_concentration / 1000))
252 @launchlogger
253 def launch(self):
254 """Launches the execution of the LeapAddIons module."""
256 # check input/output paths and parameters
257 # self.check_data_params(self.out_log, self.err_log)
259 # Setup Biobb
260 if self.check_restart():
261 return 0
262 self.stage_files()
264 # Water Type
265 # leaprc.water.tip4pew, tip4pd, tip3p, spceb, spce, opc, fb4, fb3
266 # Values: POL3BOX, QSPCFWBOX, SPCBOX, SPCFWBOX, TIP3PBOX, TIP3PFBOX, TIP4PBOX, TIP4PEWBOX, OPCBOX, OPC3BOX, TIP5PBOX.
267 source_wat_command = "source leaprc.water.tip3p"
268 if self.water_type == "TIP4PEWBOX":
269 source_wat_command = "leaprc.water.tip4pew"
270 if self.water_type == "TIP4PBOX":
271 source_wat_command = "leaprc.water.tip4pd"
272 if re.match(r"SPC", self.water_type):
273 source_wat_command = "source leaprc.water.spce"
274 if re.match(r"OPC", self.water_type):
275 source_wat_command = "source leaprc.water.opc"
277 # Counterions
278 ions_command = ""
279 if self.neutralise:
280 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " 0 \n"
281 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " 0 \n"
282 ions_command = (
283 ions_command + "addionsRand mol " + self.negative_ions_type + " 0 \n"
284 )
285 ions_command = (
286 ions_command + "addionsRand mol " + self.positive_ions_type + " 0 \n"
287 )
289 if (
290 self.ionic_concentration and self.negative_ions_number == 0 and self.positive_ions_number == 0
291 ):
292 self.find_out_number_of_ions()
293 nneg = self.nio # Update with function
294 npos = self.nio # Update with function
295 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " " + str(nneg) + " \n"
296 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " " + str(npos) + " \n"
297 ions_command = (
298 ions_command + "addionsRand mol " + self.negative_ions_type + " " + str(nneg) + " \n"
299 )
300 ions_command = (
301 ions_command + "addionsRand mol " + self.positive_ions_type + " " + str(npos) + " \n"
302 )
303 else:
304 if self.negative_ions_number != 0:
305 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " " + str(self.negative_ions_number) + " \n"
306 ions_command = (
307 ions_command + "addionsRand mol " + self.negative_ions_type + " " + str(self.negative_ions_number) + " \n"
308 )
309 if self.positive_ions_number != 0:
310 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " " + str(self.positive_ions_number) + " \n"
311 ions_command = (
312 ions_command + "addionsRand mol " + self.positive_ions_type + " " + str(self.positive_ions_number) + " \n"
313 )
315 # Creating temporary folder & Leap configuration (instructions) file
316 if self.container_path:
317 instructions_file = str(
318 PurePath(self.stage_io_dict["unique_dir"]).joinpath("leap.in")
319 )
320 instructions_file_path = str(
321 PurePath(self.container_volume_path).joinpath("leap.in")
322 )
323 tmp_folder = None
324 else:
325 tmp_folder = fu.create_unique_dir()
326 instructions_file = str(PurePath(tmp_folder).joinpath("leap.in"))
327 fu.log("Creating %s temporary folder" % tmp_folder, self.out_log)
328 instructions_file_path = instructions_file
330 ligands_lib_list = []
331 if self.io_dict["in"]["input_lib_path"] is not None:
332 if self.io_dict["in"]["input_lib_path"].endswith(".zip"):
333 ligands_lib_list = fu.unzip_list(
334 self.stage_io_dict["in"]["input_lib_path"],
335 dest_dir=tmp_folder,
336 out_log=self.out_log,
337 )
338 else:
339 ligands_lib_list.append(self.stage_io_dict["in"]["input_lib_path"])
341 ligands_frcmod_list = []
342 if self.io_dict["in"]["input_frcmod_path"] is not None:
343 if self.io_dict["in"]["input_frcmod_path"].endswith(".zip"):
344 ligands_frcmod_list = fu.unzip_list(
345 self.stage_io_dict["in"]["input_frcmod_path"],
346 dest_dir=tmp_folder,
347 out_log=self.out_log,
348 )
349 else:
350 ligands_frcmod_list.append(
351 self.stage_io_dict["in"]["input_frcmod_path"]
352 )
354 amber_params_list = []
355 if self.io_dict["in"]["input_params_path"] is not None:
356 if self.io_dict["in"]["input_params_path"].endswith(".zip"):
357 amber_params_list = fu.unzip_list(
358 self.stage_io_dict["in"]["input_params_path"],
359 dest_dir=tmp_folder,
360 out_log=self.out_log,
361 )
362 else:
363 amber_params_list.append(self.stage_io_dict["in"]["input_params_path"])
365 amber_prep_list = []
366 if self.io_dict["in"]["input_prep_path"] is not None:
367 if self.io_dict["in"]["input_prep_path"].endswith(".zip"):
368 amber_prep_list = fu.unzip_list(
369 self.stage_io_dict["in"]["input_prep_path"],
370 dest_dir=tmp_folder,
371 out_log=self.out_log,
372 )
373 else:
374 amber_prep_list.append(self.stage_io_dict["in"]["input_prep_path"])
376 leap_source_list = []
377 if self.io_dict["in"]["input_source_path"] is not None:
378 if self.io_dict["in"]["input_source_path"].endswith(".zip"):
379 leap_source_list = fu.unzip_list(
380 self.stage_io_dict["in"]["input_source_path"],
381 dest_dir=tmp_folder,
382 out_log=self.out_log,
383 )
384 else:
385 leap_source_list.append(self.stage_io_dict["in"]["input_source_path"])
387 # instructions_file = str(PurePath(tmp_folder).joinpath("leap.in"))
388 with open(instructions_file, "w") as leapin:
389 # Forcefields loaded by default:
390 # Protein: ff14SB (PARM99 + frcmod.ff99SB + frcmod.parmbsc0 + OL3 for RNA)
391 # leapin.write("source leaprc.protein.ff14SB \n")
392 # DNA: parmBSC1 (ParmBSC1 (ff99 + bsc0 + bsc1) for DNA. Ivani et al. Nature Methods 13: 55, 2016)
393 # leapin.write("source leaprc.DNA.bsc1 \n")
394 # Ligands: GAFF (General Amber Force field, J. Comput. Chem. 2004 Jul 15;25(9):1157-74)
395 # leapin.write("source leaprc.gaff \n")
397 # Forcefields loaded from input forcefield property
398 for t in self.forcefield:
399 leapin.write("source {}\n".format(t))
401 # Additional Leap commands
402 for leap_commands in leap_source_list:
403 leapin.write("source " + leap_commands + "\n")
405 # Water Model loaded from input water_model property
406 leapin.write(source_wat_command + " \n")
408 # Ions Type
409 if self.ions_type != "None":
410 leapin.write("loadamberparams frcmod." + self.ions_type + "\n")
412 # Additional Amber parameters
413 for amber_params in amber_params_list:
414 leapin.write("loadamberparams " + amber_params + "\n")
416 # Additional Amber prep files
417 for amber_prep in amber_prep_list:
418 leapin.write("loadamberprep " + amber_prep + "\n")
420 # Ligand(s) libraries (if any)
421 for amber_lib in ligands_lib_list:
422 leapin.write("loadOff " + amber_lib + "\n")
423 for amber_frcmod in ligands_frcmod_list:
424 leapin.write("loadamberparams " + amber_frcmod + "\n")
426 # Loading PDB file
427 leapin.write(
428 "mol = loadpdb " + self.stage_io_dict["in"]["input_pdb_path"] + " \n"
429 )
431 # Adding ions
432 leapin.write(ions_command)
434 # Generating box
435 leapin.write("setBox mol vdw \n")
437 # Saving output PDB file, coordinates and topology
438 leapin.write(
439 "savepdb mol " + self.stage_io_dict["out"]["output_pdb_path"] + " \n"
440 )
441 leapin.write(
442 "saveAmberParm mol " + self.stage_io_dict["out"]["output_top_path"] + " " + self.stage_io_dict["out"]["output_crd_path"] + "\n"
443 )
444 leapin.write("quit \n")
446 # Command line
447 self.cmd = [self.binary_path, "-f", instructions_file_path]
449 # Run Biobb block
450 self.run_biobb()
452 # Copy files to host
453 self.copy_to_host()
455 if self.box_type != "cubic":
456 fu.log(
457 "Fixing truncated octahedron Box in the topology and coordinates files",
458 self.out_log,
459 self.global_log,
460 )
462 # Taking box info from input PDB file, CRYST1 tag (first line)
463 with open(self.io_dict["in"]["input_pdb_path"]) as file:
464 lines = file.readlines()
465 pdb_line = lines[0]
467 if "OCTBOX" not in pdb_line:
468 fu.log(
469 "WARNING: box info not found in input PDB file (OCTBOX). Needed to correctly assign the octahedron box. Assuming cubic box.",
470 self.out_log,
471 self.global_log,
472 )
473 else:
474 # PDB info: CRYST1 86.316 86.316 86.316 109.47 109.47 109.47 P 1
475 # PDB info: OCTBOX 86.1942924 86.1942924 86.1942924 109.4712190 109.4712190 109.4712190
476 # regex_box = 'CRYST1\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*P 1'
477 regex_box = r"OCTBOX\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)"
478 box = re.findall(regex_box, pdb_line)[0]
479 box_line = ""
480 for coord in box:
481 box_line += "{:12.7f}".format(float(coord))
483 # PRMTOP info: 1.09471219E+02 8.63157502E+01 8.63157502E+01 8.63157502E+01
484 top_box_line = ""
485 top_box_line += " %.8E" % Decimal(float(box[3]))
486 top_box_line += " %.8E" % Decimal(float(box[0]))
487 top_box_line += " %.8E" % Decimal(float(box[1]))
488 top_box_line += " %.8E" % Decimal(float(box[2]))
490 # Removing box generated by tleap from the crd file (last line)
491 with open(self.io_dict["out"]["output_crd_path"]) as file:
492 lines = file.readlines()
493 crd_lines = lines[:-1]
495 # Adding old box coordinates (taken from the input pdb)
496 crd_lines.append(box_line)
498 with open(self.io_dict["out"]["output_crd_path"], "w") as file:
499 for line in crd_lines:
500 file.write(str(line))
501 file.write("\n")
503 # Now fixing IFBOX param in prmtop.
504 box_flag = False
505 ifbox_flag = 0
506 # %FLAG BOX_DIMENSIONS
507 # %FORMAT(5E16.8)
508 # 1.09471219E+02 8.63157502E+01 8.63157502E+01 8.63157502E+01
510 tmp_parmtop = str(
511 PurePath(str(tmp_folder)).joinpath("top_temp.parmtop")
512 )
513 shutil.copyfile(self.io_dict["out"]["output_top_path"], tmp_parmtop)
515 with open(self.io_dict["out"]["output_top_path"], "w") as new_top:
516 with open(tmp_parmtop, "r") as old_top:
517 for line in old_top:
518 if "BOX_DIMENSIONS" in line:
519 box_flag = True
520 new_top.write(line)
521 elif box_flag and "FORMAT" not in line:
522 new_top.write(top_box_line + "\n")
523 box_flag = False
524 elif (
525 "FLAG POINTERS" in line or ifbox_flag == 1 or ifbox_flag == 2 or ifbox_flag == 3
526 ):
527 ifbox_flag += 1
528 new_top.write(line)
529 elif ifbox_flag == 4:
530 # new_top.write(top_box_line + "\n")
531 new_top.write(line[:56] + " 2" + line[64:])
532 ifbox_flag += 1
533 else:
534 new_top.write(line)
536 # remove temporary folder(s)
537 self.tmp_files.extend([str(tmp_folder), "leap.log"])
538 self.remove_tmp_files()
540 self.check_arguments(output_files_created=True, raise_exception=False)
542 return self.return_code
545def leap_add_ions(
546 input_pdb_path: str,
547 output_pdb_path: str,
548 output_top_path: str,
549 output_crd_path: str,
550 input_lib_path: Optional[str] = None,
551 input_frcmod_path: Optional[str] = None,
552 input_params_path: Optional[str] = None,
553 input_prep_path: Optional[str] = None,
554 input_source_path: Optional[str] = None,
555 properties: Optional[dict] = None,
556 **kwargs,
557) -> int:
558 """Create the :class:`LeapAddIons <leap.leap_add_ions.LeapAddIons>` class and
559 execute the :meth:`launch() <leap.leap_add_ions.LeapAddIons.launch>` method."""
560 return LeapAddIons(**dict(locals())).launch()
563leap_add_ions.__doc__ = LeapAddIons.__doc__
564main = LeapAddIons.get_main(leap_add_ions, "Adds counterions to a system box for an AMBER MD system using tLeap.")
566if __name__ == "__main__":
567 main()