Coverage for biobb_amber/leap/leap_add_ions.py: 56%
236 statements
« prev ^ index » next coverage.py v7.11.3, created at 2025-11-14 09:23 +0000
1#!/usr/bin/env python3
3"""Module containing the LeapAddIons class and the command line interface."""
5import os
6import argparse
7import re
8import shutil
9from decimal import Decimal
10from pathlib import PurePath
11from typing import List, Optional
13from biobb_common.configuration import settings
14from biobb_common.generic.biobb_object import BiobbObject
15from biobb_common.tools import file_utils as fu
16from biobb_common.tools.file_utils import launchlogger
18from biobb_amber.leap.common import _from_string_to_list
21class LeapAddIons(BiobbObject):
22 """
23 | biobb_amber LeapAddIons
24 | Wrapper of the `AmberTools (AMBER MD Package) leap tool <https://ambermd.org/AmberTools.php>`_ module.
25 | Adds counterions to a system box for an AMBER MD system using tLeap tool from the AmberTools MD package.
27 Args:
28 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).
29 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).
30 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).
31 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).
32 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).
33 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).
34 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).
35 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).
36 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).
37 properties (dic - Python dictionary object containing the tool parameters, not input/output files):
38 * **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"].
39 * **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.
40 * **box_type** (*str*) - ("truncated_octahedron") Type for the MD system box. Values: cubic, truncated_octahedron.
41 * **ions_type** (*str*) - ("ionsjc_tip3p") Ions type. Values: ionsjc_tip3p, ionsjc_spce, ionsff99_tip3p, ions_charmm22, ionsjc_tip4pew, None.
42 * **neutralise** (*bool*) - ("True") Energetically neutralise the system adding the necessary counterions.
43 * **ionic_concentration** (*float*) - (50) Additional ionic concentration to include in the system box. Units in mM/L.
44 * **positive_ions_number** (*int*) - (0) Number of additional positive ions to include in the system box.
45 * **negative_ions_number** (*int*) - (0) Number of additional negative ions to include in the system box.
46 * **positive_ions_type** (*str*) - ("Na+") Type of additional positive ions to include in the system box. Values: Na+,K+.
47 * **negative_ions_type** (*str*) - ("Cl-") Type of additional negative ions to include in the system box. Values: Cl-.
48 * **binary_path** (*str*) - ("tleap") Path to the tleap executable binary.
49 * **remove_tmp** (*bool*) - (True) [WF property] Remove temporal files.
50 * **restart** (*bool*) - (False) [WF property] Do not execute if output files exist.
51 * **sandbox_path** (*str*) - ("./") [WF property] Parent path to the sandbox directory.
52 * **container_path** (*str*) - (None) Container path definition.
53 * **container_image** (*str*) - ('afandiadib/ambertools:serial') Container image definition.
54 * **container_volume_path** (*str*) - ('/tmp') Container volume path definition.
55 * **container_working_dir** (*str*) - (None) Container working directory definition.
56 * **container_user_id** (*str*) - (None) Container user_id definition.
57 * **container_shell_path** (*str*) - ('/bin/bash') Path to default shell inside the container.
59 Examples:
60 This is a use example of how to use the building block from Python::
62 from biobb_amber.leap.leap_add_ions import leap_add_ions
63 prop = {
64 'forcefield': ['protein.ff14SB'],
65 'water_type': 'TIP3PBOX',
66 'neutralise' : True
67 }
68 leap_add_ions(input_pdb_path='/path/to/structure.pdb',
69 output_pdb_path='/path/to/newStructure.pdb',
70 output_top_path='/path/to/newTopology.top',
71 output_crd_path='/path/to/newCoordinates.crd',
72 properties=prop)
74 Info:
75 * wrapped_software:
76 * name: AmberTools tLeap
77 * version: >20.9
78 * license: LGPL 2.1
79 * ontology:
80 * name: EDAM
81 * schema: http://edamontology.org/EDAM.owl
83 """
85 def __init__(
86 self,
87 input_pdb_path: str,
88 output_pdb_path: str,
89 output_top_path: str,
90 output_crd_path: str,
91 input_lib_path: Optional[str] = None,
92 input_frcmod_path: Optional[str] = None,
93 input_params_path: Optional[str] = None,
94 input_prep_path: Optional[str] = None,
95 input_source_path: Optional[str] = None,
96 properties: Optional[dict] = None,
97 **kwargs,
98 ):
99 properties = properties or {}
101 # Call parent class constructor
102 super().__init__(properties)
103 self.locals_var_dict = locals().copy()
105 # Input/Output files
106 self.io_dict = {
107 "in": {
108 "input_pdb_path": input_pdb_path,
109 "input_lib_path": input_lib_path,
110 "input_frcmod_path": input_frcmod_path,
111 "input_params_path": input_params_path,
112 "input_prep_path": input_prep_path,
113 "input_source_path": input_source_path,
114 },
115 "out": {
116 "output_pdb_path": output_pdb_path,
117 "output_top_path": output_top_path,
118 "output_crd_path": output_crd_path,
119 },
120 }
122 # Ligand Parameter lists
123 self.ligands_lib_list = []
124 if input_lib_path:
125 self.ligands_lib_list.append(input_lib_path)
127 self.ligands_frcmod_list = []
128 if input_frcmod_path:
129 self.ligands_frcmod_list.append(input_frcmod_path)
131 # Set default forcefields
132 amber_home_path = os.getenv("AMBERHOME")
133 protein_ff14SB_path = os.path.join(amber_home_path, 'dat', 'leap', 'cmd', 'leaprc.protein.ff14SB')
134 dna_bsc1_path = os.path.join(amber_home_path, 'dat', 'leap', 'cmd', 'leaprc.DNA.bsc1')
135 gaff_path = os.path.join(amber_home_path, 'dat', 'leap', 'cmd', 'leaprc.gaff')
137 # Properties specific for BB
138 self.properties = properties
139 self.forcefield = _from_string_to_list(
140 properties.get("forcefield", [protein_ff14SB_path, dna_bsc1_path, gaff_path])
141 )
142 # Find the paths of the leaprc files if only the force field names are provided
143 self.forcefield = self.find_leaprc_paths(self.forcefield)
144 self.water_type = properties.get("water_type", "TIP3PBOX")
145 self.box_type = properties.get("box_type", "truncated_octahedron")
146 self.ions_type = properties.get("ions_type", "ionsjc_tip3p")
147 self.neutralise = properties.get("neutralise", True)
148 self.ionic_concentration = properties.get("ionic_concentration", 50)
149 self.positive_ions_number = properties.get("positive_ions_number", 0)
150 self.positive_ions_type = properties.get("positive_ions_type", "Na+")
151 self.negative_ions_number = properties.get("negative_ions_number", 0)
152 self.negative_ions_type = properties.get("negative_ions_type", "Cl-")
153 self.binary_path = properties.get("binary_path", "tleap")
155 # Check the properties
156 self.check_properties(properties)
157 self.check_arguments()
159 def find_leaprc_paths(self, forcefields: List[str]) -> List[str]:
160 """
161 Find the leaprc paths for the force fields provided.
163 For each item in the forcefields list, the function checks if the str is a path to an existing file.
164 If not, it tries to find the file in the $AMBERHOME/dat/leap/cmd/ directory or the $AMBERHOME/dat/leap/cmd/oldff/
165 directory with and without the leaprc prefix.
167 Args:
168 forcefields (List[str]): List of force fields to find the leaprc files for.
170 Returns:
171 List[str]: List of leaprc file paths.
172 """
174 leaprc_paths = []
176 for forcefield in forcefields:
178 num_paths = len(leaprc_paths)
180 # Check if the forcefield is a path to an existing file
181 if os.path.exists(forcefield):
182 leaprc_paths.append(forcefield)
183 continue
185 # Check if the forcefield is in the leaprc directory
186 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', f"leaprc.{forcefield}")
187 if os.path.exists(leaprc_path):
188 leaprc_paths.append(leaprc_path)
189 continue
191 # Check if the forcefield is in the oldff directory
192 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', 'oldff', 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 leaprc directory without the leaprc prefix
198 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', f"{forcefield}")
199 if os.path.exists(leaprc_path):
200 leaprc_paths.append(leaprc_path)
201 continue
203 # Check if the forcefield is in the oldff directory without the leaprc prefix
204 leaprc_path = os.path.join(os.environ.get('AMBERHOME', ''), 'dat', 'leap', 'cmd', 'oldff', f"{forcefield}")
205 if os.path.exists(leaprc_path):
206 leaprc_paths.append(leaprc_path)
207 continue
209 new_num_paths = len(leaprc_paths)
211 if new_num_paths == num_paths:
212 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.")
214 return leaprc_paths
216 # def check_data_params(self, out_log, err_log):
217 # """ Checks input/output paths correctness """
219 # # Check input(s)
220 # 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__)
221 # 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__)
222 # 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__)
223 # # 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__)
224 # # 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__)
226 # # Check output(s)
227 # 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__)
228 # 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__)
229 # 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__)
231 def find_out_number_of_ions(self):
232 """Computes the number of positive and negative ions from the input ionic
233 concentration and the number of water molecules in the system box."""
235 # Finding number of water molecules in the box
236 cnt = 0
237 with open(self.io_dict["in"]["input_pdb_path"]) as fp:
238 for line in fp:
239 # Water molecules are identified with different resids
240 # by different forcefields / MD packages
241 pq = re.compile((r"WAT|HOH|TP3|TIP3|SOL"), re.M)
242 if pq.search(line):
243 # Incrementing number of water molecules just in the water
244 # oxygen atom, ignoring hydrogen or dummy atoms
245 pq2 = re.compile(r"H|EPW|EP1|EP2", re.M)
246 if not pq2.search(line):
247 cnt = cnt + 1
249 # To get to X mM ions we need
250 # n(NaCl) = #waters / 55 Mol * X M
251 # where X M = X mM / 1000
252 self.nio = int((cnt / 55) * (self.ionic_concentration / 1000))
254 @launchlogger
255 def launch(self):
256 """Launches the execution of the LeapAddIons module."""
258 # check input/output paths and parameters
259 # self.check_data_params(self.out_log, self.err_log)
261 # Setup Biobb
262 if self.check_restart():
263 return 0
264 self.stage_files()
266 # Water Type
267 # leaprc.water.tip4pew, tip4pd, tip3p, spceb, spce, opc, fb4, fb3
268 # Values: POL3BOX, QSPCFWBOX, SPCBOX, SPCFWBOX, TIP3PBOX, TIP3PFBOX, TIP4PBOX, TIP4PEWBOX, OPCBOX, OPC3BOX, TIP5PBOX.
269 source_wat_command = "source leaprc.water.tip3p"
270 if self.water_type == "TIP4PEWBOX":
271 source_wat_command = "leaprc.water.tip4pew"
272 if self.water_type == "TIP4PBOX":
273 source_wat_command = "leaprc.water.tip4pd"
274 if re.match(r"SPC", self.water_type):
275 source_wat_command = "source leaprc.water.spce"
276 if re.match(r"OPC", self.water_type):
277 source_wat_command = "source leaprc.water.opc"
279 # Counterions
280 ions_command = ""
281 if self.neutralise:
282 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " 0 \n"
283 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " 0 \n"
284 ions_command = (
285 ions_command + "addionsRand mol " + self.negative_ions_type + " 0 \n"
286 )
287 ions_command = (
288 ions_command + "addionsRand mol " + self.positive_ions_type + " 0 \n"
289 )
291 if (
292 self.ionic_concentration and self.negative_ions_number == 0 and self.positive_ions_number == 0
293 ):
294 self.find_out_number_of_ions()
295 nneg = self.nio # Update with function
296 npos = self.nio # Update with function
297 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " " + str(nneg) + " \n"
298 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " " + str(npos) + " \n"
299 ions_command = (
300 ions_command + "addionsRand mol " + self.negative_ions_type + " " + str(nneg) + " \n"
301 )
302 ions_command = (
303 ions_command + "addionsRand mol " + self.positive_ions_type + " " + str(npos) + " \n"
304 )
305 else:
306 if self.negative_ions_number != 0:
307 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " " + str(self.negative_ions_number) + " \n"
308 ions_command = (
309 ions_command + "addionsRand mol " + self.negative_ions_type + " " + str(self.negative_ions_number) + " \n"
310 )
311 if self.positive_ions_number != 0:
312 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " " + str(self.positive_ions_number) + " \n"
313 ions_command = (
314 ions_command + "addionsRand mol " + self.positive_ions_type + " " + str(self.positive_ions_number) + " \n"
315 )
317 # Creating temporary folder & Leap configuration (instructions) file
318 if self.container_path:
319 instructions_file = str(
320 PurePath(self.stage_io_dict["unique_dir"]).joinpath("leap.in")
321 )
322 instructions_file_path = str(
323 PurePath(self.container_volume_path).joinpath("leap.in")
324 )
325 self.tmp_folder = None
326 else:
327 self.tmp_folder = fu.create_unique_dir()
328 instructions_file = str(PurePath(self.tmp_folder).joinpath("leap.in"))
329 fu.log("Creating %s temporary folder" % self.tmp_folder, self.out_log)
330 instructions_file_path = instructions_file
332 ligands_lib_list = []
333 if self.io_dict["in"]["input_lib_path"] is not None:
334 if self.io_dict["in"]["input_lib_path"].endswith(".zip"):
335 ligands_lib_list = fu.unzip_list(
336 self.stage_io_dict["in"]["input_lib_path"],
337 dest_dir=self.tmp_folder,
338 out_log=self.out_log,
339 )
340 else:
341 ligands_lib_list.append(self.stage_io_dict["in"]["input_lib_path"])
343 ligands_frcmod_list = []
344 if self.io_dict["in"]["input_frcmod_path"] is not None:
345 if self.io_dict["in"]["input_frcmod_path"].endswith(".zip"):
346 ligands_frcmod_list = fu.unzip_list(
347 self.stage_io_dict["in"]["input_frcmod_path"],
348 dest_dir=self.tmp_folder,
349 out_log=self.out_log,
350 )
351 else:
352 ligands_frcmod_list.append(
353 self.stage_io_dict["in"]["input_frcmod_path"]
354 )
356 amber_params_list = []
357 if self.io_dict["in"]["input_params_path"] is not None:
358 if self.io_dict["in"]["input_params_path"].endswith(".zip"):
359 amber_params_list = fu.unzip_list(
360 self.stage_io_dict["in"]["input_params_path"],
361 dest_dir=self.tmp_folder,
362 out_log=self.out_log,
363 )
364 else:
365 amber_params_list.append(self.stage_io_dict["in"]["input_params_path"])
367 amber_prep_list = []
368 if self.io_dict["in"]["input_prep_path"] is not None:
369 if self.io_dict["in"]["input_prep_path"].endswith(".zip"):
370 amber_prep_list = fu.unzip_list(
371 self.stage_io_dict["in"]["input_prep_path"],
372 dest_dir=self.tmp_folder,
373 out_log=self.out_log,
374 )
375 else:
376 amber_prep_list.append(self.stage_io_dict["in"]["input_prep_path"])
378 leap_source_list = []
379 if self.io_dict["in"]["input_source_path"] is not None:
380 if self.io_dict["in"]["input_source_path"].endswith(".zip"):
381 leap_source_list = fu.unzip_list(
382 self.stage_io_dict["in"]["input_source_path"],
383 dest_dir=self.tmp_folder,
384 out_log=self.out_log,
385 )
386 else:
387 leap_source_list.append(self.stage_io_dict["in"]["input_source_path"])
389 # instructions_file = str(PurePath(self.tmp_folder).joinpath("leap.in"))
390 with open(instructions_file, "w") as leapin:
391 # Forcefields loaded by default:
392 # Protein: ff14SB (PARM99 + frcmod.ff99SB + frcmod.parmbsc0 + OL3 for RNA)
393 # leapin.write("source leaprc.protein.ff14SB \n")
394 # DNA: parmBSC1 (ParmBSC1 (ff99 + bsc0 + bsc1) for DNA. Ivani et al. Nature Methods 13: 55, 2016)
395 # leapin.write("source leaprc.DNA.bsc1 \n")
396 # Ligands: GAFF (General Amber Force field, J. Comput. Chem. 2004 Jul 15;25(9):1157-74)
397 # leapin.write("source leaprc.gaff \n")
399 # Forcefields loaded from input forcefield property
400 for t in self.forcefield:
401 leapin.write("source {}\n".format(t))
403 # Additional Leap commands
404 for leap_commands in leap_source_list:
405 leapin.write("source " + leap_commands + "\n")
407 # Water Model loaded from input water_model property
408 leapin.write(source_wat_command + " \n")
410 # Ions Type
411 if self.ions_type != "None":
412 leapin.write("loadamberparams frcmod." + self.ions_type + "\n")
414 # Additional Amber parameters
415 for amber_params in amber_params_list:
416 leapin.write("loadamberparams " + amber_params + "\n")
418 # Additional Amber prep files
419 for amber_prep in amber_prep_list:
420 leapin.write("loadamberprep " + amber_prep + "\n")
422 # Ligand(s) libraries (if any)
423 for amber_lib in ligands_lib_list:
424 leapin.write("loadOff " + amber_lib + "\n")
425 for amber_frcmod in ligands_frcmod_list:
426 leapin.write("loadamberparams " + amber_frcmod + "\n")
428 # Loading PDB file
429 leapin.write(
430 "mol = loadpdb " + self.stage_io_dict["in"]["input_pdb_path"] + " \n"
431 )
433 # Adding ions
434 leapin.write(ions_command)
436 # Generating box
437 leapin.write("setBox mol vdw \n")
439 # Saving output PDB file, coordinates and topology
440 leapin.write(
441 "savepdb mol " + self.stage_io_dict["out"]["output_pdb_path"] + " \n"
442 )
443 leapin.write(
444 "saveAmberParm mol " + self.stage_io_dict["out"]["output_top_path"] + " " + self.stage_io_dict["out"]["output_crd_path"] + "\n"
445 )
446 leapin.write("quit \n")
448 # Command line
449 self.cmd = [self.binary_path, "-f", instructions_file_path]
451 # Run Biobb block
452 self.run_biobb()
454 # Copy files to host
455 self.copy_to_host()
457 if self.box_type != "cubic":
458 fu.log(
459 "Fixing truncated octahedron Box in the topology and coordinates files",
460 self.out_log,
461 self.global_log,
462 )
464 # Taking box info from input PDB file, CRYST1 tag (first line)
465 with open(self.io_dict["in"]["input_pdb_path"]) as file:
466 lines = file.readlines()
467 pdb_line = lines[0]
469 if "OCTBOX" not in pdb_line:
470 fu.log(
471 "WARNING: box info not found in input PDB file (OCTBOX). Needed to correctly assign the octahedron box. Assuming cubic box.",
472 self.out_log,
473 self.global_log,
474 )
475 else:
476 # PDB info: CRYST1 86.316 86.316 86.316 109.47 109.47 109.47 P 1
477 # PDB info: OCTBOX 86.1942924 86.1942924 86.1942924 109.4712190 109.4712190 109.4712190
478 # regex_box = 'CRYST1\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*P 1'
479 regex_box = r"OCTBOX\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)"
480 box = re.findall(regex_box, pdb_line)[0]
481 box_line = ""
482 for coord in box:
483 box_line += "{:12.7f}".format(float(coord))
485 # PRMTOP info: 1.09471219E+02 8.63157502E+01 8.63157502E+01 8.63157502E+01
486 top_box_line = ""
487 top_box_line += " %.8E" % Decimal(float(box[3]))
488 top_box_line += " %.8E" % Decimal(float(box[0]))
489 top_box_line += " %.8E" % Decimal(float(box[1]))
490 top_box_line += " %.8E" % Decimal(float(box[2]))
492 # Removing box generated by tleap from the crd file (last line)
493 with open(self.io_dict["out"]["output_crd_path"]) as file:
494 lines = file.readlines()
495 crd_lines = lines[:-1]
497 # Adding old box coordinates (taken from the input pdb)
498 crd_lines.append(box_line)
500 with open(self.io_dict["out"]["output_crd_path"], "w") as file:
501 for line in crd_lines:
502 file.write(str(line))
503 file.write("\n")
505 # Now fixing IFBOX param in prmtop.
506 box_flag = False
507 ifbox_flag = 0
508 # %FLAG BOX_DIMENSIONS
509 # %FORMAT(5E16.8)
510 # 1.09471219E+02 8.63157502E+01 8.63157502E+01 8.63157502E+01
512 tmp_parmtop = str(
513 PurePath(str(self.tmp_folder)).joinpath("top_temp.parmtop")
514 )
515 shutil.copyfile(self.io_dict["out"]["output_top_path"], tmp_parmtop)
517 with open(self.io_dict["out"]["output_top_path"], "w") as new_top:
518 with open(tmp_parmtop, "r") as old_top:
519 for line in old_top:
520 if "BOX_DIMENSIONS" in line:
521 box_flag = True
522 new_top.write(line)
523 elif box_flag and "FORMAT" not in line:
524 new_top.write(top_box_line + "\n")
525 box_flag = False
526 elif (
527 "FLAG POINTERS" in line or ifbox_flag == 1 or ifbox_flag == 2 or ifbox_flag == 3
528 ):
529 ifbox_flag += 1
530 new_top.write(line)
531 elif ifbox_flag == 4:
532 # new_top.write(top_box_line + "\n")
533 new_top.write(line[:56] + " 2" + line[64:])
534 ifbox_flag += 1
535 else:
536 new_top.write(line)
538 # remove temporary folder(s)
539 self.tmp_files.extend([str(self.tmp_folder), "leap.log"])
540 self.remove_tmp_files()
542 self.check_arguments(output_files_created=True, raise_exception=False)
544 return self.return_code
547def leap_add_ions(
548 input_pdb_path: str,
549 output_pdb_path: str,
550 output_top_path: str,
551 output_crd_path: str,
552 input_lib_path: Optional[str] = None,
553 input_frcmod_path: Optional[str] = None,
554 input_params_path: Optional[str] = None,
555 input_prep_path: Optional[str] = None,
556 input_source_path: Optional[str] = None,
557 properties: Optional[dict] = None,
558 **kwargs,
559) -> int:
560 """Create :class:`LeapAddIons <leap.leap_add_ions.LeapAddIons>`leap.leap_add_ions.LeapAddIons class and
561 execute :meth:`launch() <leap.leap_add_ions.LeapAddIons.launch>` method"""
563 return LeapAddIons(
564 input_pdb_path=input_pdb_path,
565 input_lib_path=input_lib_path,
566 input_frcmod_path=input_frcmod_path,
567 input_params_path=input_params_path,
568 input_prep_path=input_prep_path,
569 input_source_path=input_source_path,
570 output_pdb_path=output_pdb_path,
571 output_top_path=output_top_path,
572 output_crd_path=output_crd_path,
573 properties=properties,
574 ).launch()
576 leap_add_ions.__doc__ = LeapAddIons.__doc__
579def main():
580 parser = argparse.ArgumentParser(
581 description="Adds counterions to a system box for an AMBER MD system using tLeap.",
582 formatter_class=lambda prog: argparse.RawTextHelpFormatter(prog, width=99999),
583 )
584 parser.add_argument("--config", required=False, help="Configuration file")
586 # Specific args
587 required_args = parser.add_argument_group("required arguments")
588 required_args.add_argument(
589 "--input_pdb_path",
590 required=True,
591 help="Input 3D structure PDB file. Accepted formats: pdb.",
592 )
593 required_args.add_argument(
594 "--input_lib_path",
595 required=False,
596 help="Input ligand library parameters file. Accepted formats: lib, zip.",
597 )
598 required_args.add_argument(
599 "--input_frcmod_path",
600 required=False,
601 help="Input ligand frcmod parameters file. Accepted formats: frcmod, zip.",
602 )
603 required_args.add_argument(
604 "--input_params_path",
605 required=False,
606 help="Additional leap parameter files to load with loadAmberParams Leap command. Accepted formats: leapin, in, txt, zip.",
607 )
608 required_args.add_argument(
609 "--input_prep_path",
610 required=False,
611 help="Additional leap parameter files to load with loadAmberPrep Leap command. Accepted formats: leapin, in, txt, zip.",
612 )
613 required_args.add_argument(
614 "--input_source_path",
615 required=False,
616 help="Additional leap command files to load with source Leap command. Accepted formats: leapin, in, txt, zip.",
617 )
618 required_args.add_argument(
619 "--output_pdb_path",
620 required=True,
621 help="Output 3D structure PDB file matching the topology file. Accepted formats: pdb.",
622 )
623 required_args.add_argument(
624 "--output_top_path",
625 required=True,
626 help="Output topology file (AMBER ParmTop). Accepted formats: top.",
627 )
628 required_args.add_argument(
629 "--output_crd_path",
630 required=True,
631 help="Output coordinates file (AMBER crd). Accepted formats: crd.",
632 )
634 args = parser.parse_args()
635 config = args.config if args.config else None
636 properties = settings.ConfReader(config=config).get_prop_dic()
638 # Specific call
639 leap_add_ions(
640 input_pdb_path=args.input_pdb_path,
641 input_lib_path=args.input_lib_path,
642 input_frcmod_path=args.input_frcmod_path,
643 input_params_path=args.input_params_path,
644 input_prep_path=args.input_prep_path,
645 input_source_path=args.input_source_path,
646 output_pdb_path=args.output_pdb_path,
647 output_top_path=args.output_top_path,
648 output_crd_path=args.output_crd_path,
649 properties=properties,
650 )
653if __name__ == "__main__":
654 main()