Coverage for biobb_amber/leap/leap_add_ions.py: 55%
237 statements
« prev ^ index » next coverage.py v7.6.10, created at 2025-01-28 08:28 +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 Mol/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([
540 # self.stage_io_dict.get("unique_dir", ""),
541 str(self.tmp_folder), "leap.log"
542 ])
543 self.remove_tmp_files()
545 self.check_arguments(output_files_created=True, raise_exception=False)
547 return self.return_code
550def leap_add_ions(
551 input_pdb_path: str,
552 output_pdb_path: str,
553 output_top_path: str,
554 output_crd_path: str,
555 input_lib_path: Optional[str] = None,
556 input_frcmod_path: Optional[str] = None,
557 input_params_path: Optional[str] = None,
558 input_prep_path: Optional[str] = None,
559 input_source_path: Optional[str] = None,
560 properties: Optional[dict] = None,
561 **kwargs,
562) -> int:
563 """Create :class:`LeapAddIons <leap.leap_add_ions.LeapAddIons>`leap.leap_add_ions.LeapAddIons class and
564 execute :meth:`launch() <leap.leap_add_ions.LeapAddIons.launch>` method"""
566 return LeapAddIons(
567 input_pdb_path=input_pdb_path,
568 input_lib_path=input_lib_path,
569 input_frcmod_path=input_frcmod_path,
570 input_params_path=input_params_path,
571 input_prep_path=input_prep_path,
572 input_source_path=input_source_path,
573 output_pdb_path=output_pdb_path,
574 output_top_path=output_top_path,
575 output_crd_path=output_crd_path,
576 properties=properties,
577 ).launch()
579 leap_add_ions.__doc__ = LeapAddIons.__doc__
582def main():
583 parser = argparse.ArgumentParser(
584 description="Adds counterions to a system box for an AMBER MD system using tLeap.",
585 formatter_class=lambda prog: argparse.RawTextHelpFormatter(prog, width=99999),
586 )
587 parser.add_argument("--config", required=False, help="Configuration file")
589 # Specific args
590 required_args = parser.add_argument_group("required arguments")
591 required_args.add_argument(
592 "--input_pdb_path",
593 required=True,
594 help="Input 3D structure PDB file. Accepted formats: pdb.",
595 )
596 required_args.add_argument(
597 "--input_lib_path",
598 required=False,
599 help="Input ligand library parameters file. Accepted formats: lib, zip.",
600 )
601 required_args.add_argument(
602 "--input_frcmod_path",
603 required=False,
604 help="Input ligand frcmod parameters file. Accepted formats: frcmod, zip.",
605 )
606 required_args.add_argument(
607 "--input_params_path",
608 required=False,
609 help="Additional leap parameter files to load with loadAmberParams Leap command. Accepted formats: leapin, in, txt, zip.",
610 )
611 required_args.add_argument(
612 "--input_prep_path",
613 required=False,
614 help="Additional leap parameter files to load with loadAmberPrep Leap command. Accepted formats: leapin, in, txt, zip.",
615 )
616 required_args.add_argument(
617 "--input_source_path",
618 required=False,
619 help="Additional leap command files to load with source Leap command. Accepted formats: leapin, in, txt, zip.",
620 )
621 required_args.add_argument(
622 "--output_pdb_path",
623 required=True,
624 help="Output 3D structure PDB file matching the topology file. Accepted formats: pdb.",
625 )
626 required_args.add_argument(
627 "--output_top_path",
628 required=True,
629 help="Output topology file (AMBER ParmTop). Accepted formats: top.",
630 )
631 required_args.add_argument(
632 "--output_crd_path",
633 required=True,
634 help="Output coordinates file (AMBER crd). Accepted formats: crd.",
635 )
637 args = parser.parse_args()
638 config = args.config if args.config else None
639 properties = settings.ConfReader(config=config).get_prop_dic()
641 # Specific call
642 leap_add_ions(
643 input_pdb_path=args.input_pdb_path,
644 input_lib_path=args.input_lib_path,
645 input_frcmod_path=args.input_frcmod_path,
646 input_params_path=args.input_params_path,
647 input_prep_path=args.input_prep_path,
648 input_source_path=args.input_source_path,
649 output_pdb_path=args.output_pdb_path,
650 output_top_path=args.output_top_path,
651 output_crd_path=args.output_crd_path,
652 properties=properties,
653 )
656if __name__ == "__main__":
657 main()