Coverage for biobb_amber / leap / leap_add_ions.py: 59%
227 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-12 10:52 +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 # 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 @launchlogger
261 def launch(self):
262 """Launches the execution of the LeapAddIons module."""
264 # check input/output paths and parameters
265 # self.check_data_params(self.out_log, self.err_log)
267 # Setup Biobb
268 if self.check_restart():
269 return 0
270 self.stage_files()
272 # Water Type
273 # leaprc.water.tip4pew, tip4pd, tip3p, spceb, spce, opc, fb4, fb3
274 # Values: POL3BOX, QSPCFWBOX, SPCBOX, SPCFWBOX, TIP3PBOX, TIP3PFBOX, TIP4PBOX, TIP4PEWBOX, OPCBOX, OPC3BOX, TIP5PBOX.
275 source_wat_command = "source leaprc.water.tip3p"
276 if self.water_type == "TIP4PEWBOX":
277 source_wat_command = "leaprc.water.tip4pew"
278 if self.water_type == "TIP4PBOX":
279 source_wat_command = "leaprc.water.tip4pd"
280 if re.match(r"SPC", self.water_type):
281 source_wat_command = "source leaprc.water.spce"
282 if re.match(r"OPC", self.water_type):
283 source_wat_command = "source leaprc.water.opc"
285 # Counterions
286 ions_command = ""
287 if self.neutralise:
288 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " 0 \n"
289 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " 0 \n"
290 ions_command = (
291 ions_command + "addionsRand mol " + self.negative_ions_type + " 0 \n"
292 )
293 ions_command = (
294 ions_command + "addionsRand mol " + self.positive_ions_type + " 0 \n"
295 )
297 if (
298 self.ionic_concentration and self.negative_ions_number == 0 and self.positive_ions_number == 0
299 ):
300 self.find_out_number_of_ions()
301 nneg = self.nio # Update with function
302 npos = self.nio # Update with function
303 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " " + str(nneg) + " \n"
304 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " " + str(npos) + " \n"
305 ions_command = (
306 ions_command + "addionsRand mol " + self.negative_ions_type + " " + str(nneg) + " \n"
307 )
308 ions_command = (
309 ions_command + "addionsRand mol " + self.positive_ions_type + " " + str(npos) + " \n"
310 )
311 else:
312 if self.negative_ions_number != 0:
313 # ions_command = ions_command + "addions mol " + self.negative_ions_type + " " + str(self.negative_ions_number) + " \n"
314 ions_command = (
315 ions_command + "addionsRand mol " + self.negative_ions_type + " " + str(self.negative_ions_number) + " \n"
316 )
317 if self.positive_ions_number != 0:
318 # ions_command = ions_command + "addions mol " + self.positive_ions_type + " " + str(self.positive_ions_number) + " \n"
319 ions_command = (
320 ions_command + "addionsRand mol " + self.positive_ions_type + " " + str(self.positive_ions_number) + " \n"
321 )
323 # Creating temporary folder & Leap configuration (instructions) file
324 if self.container_path:
325 instructions_file = str(
326 PurePath(self.stage_io_dict["unique_dir"]).joinpath("leap.in")
327 )
328 instructions_file_path = str(
329 PurePath(self.container_volume_path).joinpath("leap.in")
330 )
331 tmp_folder = None
332 else:
333 tmp_folder = fu.create_unique_dir()
334 instructions_file = str(PurePath(tmp_folder).joinpath("leap.in"))
335 fu.log("Creating %s temporary folder" % tmp_folder, self.out_log)
336 instructions_file_path = instructions_file
338 ligands_lib_list = []
339 if self.io_dict["in"]["input_lib_path"] is not None:
340 if self.io_dict["in"]["input_lib_path"].endswith(".zip"):
341 ligands_lib_list = fu.unzip_list(
342 self.stage_io_dict["in"]["input_lib_path"],
343 dest_dir=tmp_folder,
344 out_log=self.out_log,
345 )
346 else:
347 ligands_lib_list.append(self.stage_io_dict["in"]["input_lib_path"])
349 ligands_frcmod_list = []
350 if self.io_dict["in"]["input_frcmod_path"] is not None:
351 if self.io_dict["in"]["input_frcmod_path"].endswith(".zip"):
352 ligands_frcmod_list = fu.unzip_list(
353 self.stage_io_dict["in"]["input_frcmod_path"],
354 dest_dir=tmp_folder,
355 out_log=self.out_log,
356 )
357 else:
358 ligands_frcmod_list.append(
359 self.stage_io_dict["in"]["input_frcmod_path"]
360 )
362 amber_params_list = []
363 if self.io_dict["in"]["input_params_path"] is not None:
364 if self.io_dict["in"]["input_params_path"].endswith(".zip"):
365 amber_params_list = fu.unzip_list(
366 self.stage_io_dict["in"]["input_params_path"],
367 dest_dir=tmp_folder,
368 out_log=self.out_log,
369 )
370 else:
371 amber_params_list.append(self.stage_io_dict["in"]["input_params_path"])
373 amber_prep_list = []
374 if self.io_dict["in"]["input_prep_path"] is not None:
375 if self.io_dict["in"]["input_prep_path"].endswith(".zip"):
376 amber_prep_list = fu.unzip_list(
377 self.stage_io_dict["in"]["input_prep_path"],
378 dest_dir=tmp_folder,
379 out_log=self.out_log,
380 )
381 else:
382 amber_prep_list.append(self.stage_io_dict["in"]["input_prep_path"])
384 leap_source_list = []
385 if self.io_dict["in"]["input_source_path"] is not None:
386 if self.io_dict["in"]["input_source_path"].endswith(".zip"):
387 leap_source_list = fu.unzip_list(
388 self.stage_io_dict["in"]["input_source_path"],
389 dest_dir=tmp_folder,
390 out_log=self.out_log,
391 )
392 else:
393 leap_source_list.append(self.stage_io_dict["in"]["input_source_path"])
395 # instructions_file = str(PurePath(tmp_folder).joinpath("leap.in"))
396 with open(instructions_file, "w") as leapin:
397 # Forcefields loaded by default:
398 # Protein: ff14SB (PARM99 + frcmod.ff99SB + frcmod.parmbsc0 + OL3 for RNA)
399 # leapin.write("source leaprc.protein.ff14SB \n")
400 # DNA: parmBSC1 (ParmBSC1 (ff99 + bsc0 + bsc1) for DNA. Ivani et al. Nature Methods 13: 55, 2016)
401 # leapin.write("source leaprc.DNA.bsc1 \n")
402 # Ligands: GAFF (General Amber Force field, J. Comput. Chem. 2004 Jul 15;25(9):1157-74)
403 # leapin.write("source leaprc.gaff \n")
405 # Forcefields loaded from input forcefield property
406 for t in self.forcefield:
407 if self.container_path:
408 leapin.write("source leaprc.{}\n".format(t))
409 else:
410 leapin.write("source {}\n".format(t))
412 # Additional Leap commands
413 for leap_commands in leap_source_list:
414 leapin.write("source " + leap_commands + "\n")
416 # Water Model loaded from input water_model property
417 leapin.write(source_wat_command + " \n")
419 # Ions Type
420 if self.ions_type != "None":
421 leapin.write("loadamberparams frcmod." + self.ions_type + "\n")
423 # Additional Amber parameters
424 for amber_params in amber_params_list:
425 leapin.write("loadamberparams " + amber_params + "\n")
427 # Additional Amber prep files
428 for amber_prep in amber_prep_list:
429 leapin.write("loadamberprep " + amber_prep + "\n")
431 # Ligand(s) libraries (if any)
432 for amber_lib in ligands_lib_list:
433 leapin.write("loadOff " + amber_lib + "\n")
434 for amber_frcmod in ligands_frcmod_list:
435 leapin.write("loadamberparams " + amber_frcmod + "\n")
437 # Loading PDB file
438 leapin.write(
439 "mol = loadpdb " + self.stage_io_dict["in"]["input_pdb_path"] + " \n"
440 )
442 # Adding ions
443 leapin.write(ions_command)
445 # Generating box
446 leapin.write("setBox mol vdw \n")
448 # Saving output PDB file, coordinates and topology
449 leapin.write(
450 "savepdb mol " + self.stage_io_dict["out"]["output_pdb_path"] + " \n"
451 )
452 leapin.write(
453 "saveAmberParm mol " + self.stage_io_dict["out"]["output_top_path"] + " " + self.stage_io_dict["out"]["output_crd_path"] + "\n"
454 )
455 leapin.write("quit \n")
457 # Container path
458 if self.container_path:
459 if not self.container_working_dir:
460 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)
461 self.container_working_dir = self.container_volume_path
463 # Command line
464 self.cmd = [self.binary_path, "-f", instructions_file_path]
466 # Run Biobb block
467 self.run_biobb()
469 # Copy files to host
470 self.copy_to_host()
472 if self.box_type != "cubic":
473 fu.log(
474 "Fixing truncated octahedron Box in the topology and coordinates files",
475 self.out_log,
476 self.global_log,
477 )
479 # Taking box info from input PDB file, CRYST1 tag (first line)
480 with open(self.io_dict["in"]["input_pdb_path"]) as file:
481 lines = file.readlines()
482 pdb_line = lines[0]
484 if "OCTBOX" not in pdb_line:
485 fu.log(
486 "WARNING: box info not found in input PDB file (OCTBOX). Needed to correctly assign the octahedron box. Assuming cubic box.",
487 self.out_log,
488 self.global_log,
489 )
490 else:
491 # PDB info: CRYST1 86.316 86.316 86.316 109.47 109.47 109.47 P 1
492 # PDB info: OCTBOX 86.1942924 86.1942924 86.1942924 109.4712190 109.4712190 109.4712190
493 # regex_box = 'CRYST1\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*P 1'
494 regex_box = r"OCTBOX\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)\s*(\d+\.\d+)"
495 box = re.findall(regex_box, pdb_line)[0]
496 box_line = ""
497 for coord in box:
498 box_line += "{:12.7f}".format(float(coord))
500 # PRMTOP info: 1.09471219E+02 8.63157502E+01 8.63157502E+01 8.63157502E+01
501 top_box_line = ""
502 top_box_line += " %.8E" % Decimal(float(box[3]))
503 top_box_line += " %.8E" % Decimal(float(box[0]))
504 top_box_line += " %.8E" % Decimal(float(box[1]))
505 top_box_line += " %.8E" % Decimal(float(box[2]))
507 # Removing box generated by tleap from the crd file (last line)
508 with open(self.io_dict["out"]["output_crd_path"]) as file:
509 lines = file.readlines()
510 crd_lines = lines[:-1]
512 # Adding old box coordinates (taken from the input pdb)
513 crd_lines.append(box_line)
515 with open(self.io_dict["out"]["output_crd_path"], "w") as file:
516 for line in crd_lines:
517 file.write(str(line))
518 file.write("\n")
520 # Now fixing IFBOX param in prmtop.
521 box_flag = False
522 ifbox_flag = 0
523 # %FLAG BOX_DIMENSIONS
524 # %FORMAT(5E16.8)
525 # 1.09471219E+02 8.63157502E+01 8.63157502E+01 8.63157502E+01
527 tmp_parmtop = str(
528 #PurePath(str(tmp_folder)).joinpath("top_temp.parmtop")
529 PurePath(self.stage_io_dict["unique_dir"]).joinpath("top_temp.parmtop")
530 )
531 shutil.copyfile(self.io_dict["out"]["output_top_path"], tmp_parmtop)
533 with open(self.io_dict["out"]["output_top_path"], "w") as new_top:
534 with open(tmp_parmtop, "r") as old_top:
535 for line in old_top:
536 if "BOX_DIMENSIONS" in line:
537 box_flag = True
538 new_top.write(line)
539 elif box_flag and "FORMAT" not in line:
540 new_top.write(top_box_line + "\n")
541 box_flag = False
542 elif (
543 "FLAG POINTERS" in line or ifbox_flag == 1 or ifbox_flag == 2 or ifbox_flag == 3
544 ):
545 ifbox_flag += 1
546 new_top.write(line)
547 elif ifbox_flag == 4:
548 # new_top.write(top_box_line + "\n")
549 new_top.write(line[:56] + " 2" + line[64:])
550 ifbox_flag += 1
551 else:
552 new_top.write(line)
554 # remove temporary folder(s)
555 self.tmp_files.extend([str(tmp_folder), "leap.log"])
556 self.remove_tmp_files()
558 self.check_arguments(output_files_created=True, raise_exception=False)
560 return self.return_code
563def leap_add_ions(
564 input_pdb_path: str,
565 output_pdb_path: str,
566 output_top_path: str,
567 output_crd_path: str,
568 input_lib_path: Optional[str] = None,
569 input_frcmod_path: Optional[str] = None,
570 input_params_path: Optional[str] = None,
571 input_prep_path: Optional[str] = None,
572 input_source_path: Optional[str] = None,
573 properties: Optional[dict] = None,
574 **kwargs,
575) -> int:
576 """Create the :class:`LeapAddIons <leap.leap_add_ions.LeapAddIons>` class and
577 execute the :meth:`launch() <leap.leap_add_ions.LeapAddIons.launch>` method."""
578 return LeapAddIons(**dict(locals())).launch()
581leap_add_ions.__doc__ = LeapAddIons.__doc__
582main = LeapAddIons.get_main(leap_add_ions, "Adds counterions to a system box for an AMBER MD system using tLeap.")
584if __name__ == "__main__":
585 main()