biobb_amber/leap/leap_add_ions.py source
1 #!/usr/bin/env python3
2
3 """Module containing the LeapAddIons class and the command line interface."""
4
5 import os
6 import re
7 import shutil
8 from decimal import Decimal
9 from pathlib import PurePath
10 from typing import List, Optional
11
12 from biobb_common.generic.biobb_object import BiobbObject
13 from biobb_common.tools import file_utils as fu
14 from biobb_common.tools.file_utils import launchlogger
15
16 from biobb_amber.leap.common import _from_string_to_list
17
18
19 class 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.
24
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.
56
57 Examples:
58 This is a use example of how to use the building block from Python::
59
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)
71
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
80
81 """
82
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 {}
98
99 # Call parent class constructor
100 super().__init__(properties)
101 self.locals_var_dict = locals().copy()
102
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 }
119
120 # Ligand Parameter lists
121 self.ligands_lib_list = []
122 if input_lib_path:
123 self.ligands_lib_list.append(input_lib_path)
124
125 self.ligands_frcmod_list = []
126 if input_frcmod_path:
127 self.ligands_frcmod_list.append(input_frcmod_path)
128
129 # Properties specific for BB
130 self.properties = properties
131
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')
142
143 self.forcefield = _from_string_to_list(
144 properties.get("forcefield", [protein_ff14SB_path, dna_bsc1_path, gaff_path])
145 )
146
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)
149
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")
160
161 # Check the properties
162 self.check_properties(properties)
163 self.check_arguments()
164
165 def find_leaprc_paths(self, forcefields: List[str]) -> List[str]:
166 """
167 Find the leaprc paths for the force fields provided.
168
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.
172
173 Args:
174 forcefields (List[str]): List of force fields to find the leaprc files for.
175
176 Returns:
177 List[str]: List of leaprc file paths.
178 """
179
180 leaprc_paths = []
181
182 for forcefield in forcefields:
183
184 num_paths = len(leaprc_paths)
185
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
190
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
196
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
202
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
208
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
214
215 new_num_paths = len(leaprc_paths)
216
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.")
219
220 return leaprc_paths
221
222 # def check_data_params(self, out_log, err_log):
223 # """ Checks input/output paths correctness """
224
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__)
231
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__)
236
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."""
240
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
254
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))
259
260 @launchlogger
261 def launch(self):
262 """Launches the execution of the LeapAddIons module."""
263
264 # check input/output paths and parameters
265 # self.check_data_params(self.out_log, self.err_log)
266
267 # Setup Biobb
268 if self.check_restart():
269 return 0
270 self.stage_files()
271
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"
284
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 )
296
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 )
322
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
337
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"])
348
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 )
361
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"])
372
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"])
383
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"])
394
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")
404
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))
411
412 # Additional Leap commands
413 for leap_commands in leap_source_list:
414 leapin.write("source " + leap_commands + "\n")
415
416 # Water Model loaded from input water_model property
417 leapin.write(source_wat_command + " \n")
418
419 # Ions Type
420 if self.ions_type != "None":
421 leapin.write("loadamberparams frcmod." + self.ions_type + "\n")
422
423 # Additional Amber parameters
424 for amber_params in amber_params_list:
425 leapin.write("loadamberparams " + amber_params + "\n")
426
427 # Additional Amber prep files
428 for amber_prep in amber_prep_list:
429 leapin.write("loadamberprep " + amber_prep + "\n")
430
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")
436
437 # Loading PDB file
438 leapin.write(
439 "mol = loadpdb " + self.stage_io_dict["in"]["input_pdb_path"] + " \n"
440 )
441
442 # Adding ions
443 leapin.write(ions_command)
444
445 # Generating box
446 leapin.write("setBox mol vdw \n")
447
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")
456
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
462
463 # Command line
464 self.cmd = [self.binary_path, "-f", instructions_file_path]
465
466 # Run Biobb block
467 self.run_biobb()
468
469 # Copy files to host
470 self.copy_to_host()
471
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 )
478
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]
483
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))
499
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]))
506
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]
511
512 # Adding old box coordinates (taken from the input pdb)
513 crd_lines.append(box_line)
514
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")
519
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
526
527 tmp_parmtop = str(
- E265 Block comment should start with '# '
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)
532
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)
553
554 # remove temporary folder(s)
555 self.tmp_files.extend([str(tmp_folder), "leap.log"])
556 self.remove_tmp_files()
557
558 self.check_arguments(output_files_created=True, raise_exception=False)
559
560 return self.return_code
561
562
563 def 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()
579
580
581 leap_add_ions.__doc__ = LeapAddIons.__doc__
582 main = LeapAddIons.get_main(leap_add_ions, "Adds counterions to a system box for an AMBER MD system using tLeap.")
583
584 if __name__ == "__main__":
585 main()