Complete remake. Better configuration parsing, results handling, clea…

…ner scripts.

build.sh

@@ -0,0 +1 @@
+$PYTHON setup.py install

lib/config_parser.py

@@ -0,0 +1,147 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+#Marcelo D. Poleto
+#MAR 2023
+
+
+import sys
+from configparser import ConfigParser
+from configparser import NoOptionError
+import numpy as np
+
+
+##############################################################
+# Parsing a configuration file and returning a Config object
+#
+# - configuration parameters are returned as object attributes
+# - exception errors are thrown for wrong values or wrong
+#   section headers are provided
+#
+##############################################################
+
+class Configuration:
+    def __init__(self, configinput):
+        self.config = ConfigParser(allow_no_value=True, inline_comment_prefixes="#")
+        self.configinput = configinput
+        self.parse_config()
+
+
+    def parse_config(self):
+        if self.configinput is None:
+            raise ValueError(">>> Configuration file is a necessary input. Use -template to generate one.\n>>> Exiting...\n")
+        else:
+            self.config.read(self.configinput)
+
+        allowed_sections = ["Environment Selection", "Probe Selection", "Solvent", "Time", "Box Info"]
+        for section in self.config.sections():
+            if section not in allowed_sections:
+                raise ValueError(""">>> ERROR: Unknown section ('""" + section + """') in configuration file!\n""")
+
+        self.parse_environment_selection()
+        self.parse_probe_selection()
+        self.parse_solvent()
+        self.parse_time()
+        self.parse_box_info()
+
+
+    def parse_environment_selection(self):
+        try:
+            self.sele_elecfield = self.config['Environment Selection']["sele_environment"].strip('"')
+        except:
+            raise ValueError("\n>>> ERROR: sele_environment must be defined with a valid MDanalysis selection!\n")
+
+    def parse_probe_selection(self):
+        try:
+            self.mode = self.config['Probe Selection']["mode"].strip('"').lower()
+            if self.mode == "atom":
+                try:
+                    self.selatom = self.config['Probe Selection']["selatom"].strip('"')
+                except:
+                    raise ValueError("""\n>>> ERROR: in "ATOM" mode, selatom must be defined!\n""")
+            elif self.mode == "bond":
+                try:
+                    self.selbond1 = self.config['Probe Selection']["selbond1"].strip('"')
+                    self.selbond2 = self.config['Probe Selection']["selbond2"].strip('"')
+                except:
+                    raise ValueError(""">>> ERROR: in "BOND" mode, both selbond1 and selbond2 must be defined!\n""")
+            elif self.mode == "coordinate":
+                try:
+                    tmp_probecoordinate = self.config['Probe Selection']["probecoordinate"].strip('[]').split(",")
+                    self.probecoordinate = np.array([float(item) for item in tmp_probecoordinate])
+                except:
+                    raise ValueError(""">>> ERROR: in "COORDINATE" mode, a list of coordinates [X,Y,Z] must be provided!\n""")
+            elif self.mode == "list":
+                try:
+                    self.file_of_coordinates = self.config['Probe Selection']["file_of_coordinates"].strip('"')
+                except:
+                    raise ValueError(""">>> ERROR: in "LIST" mode, "file_of_coordinates" must be defined!\n""")
+            else:
+                raise ValueError(""">>> ERROR: "mode" must be defined as "ATOM", "BOND", "COORDINATE" or "LIST"!\n""")
+        except:
+            raise ValueError(""">>> ERROR: "mode" must be defined as "ATOM", "BOND", "COORDINATE" or "LIST"!\n""")
+
+        if self.mode == "coordinate" or self.mode == "list":
+            try:
+                self.remove_self = self.config['Probe Selection'].getboolean("remove_self")
+                if self.remove_self:
+                    try:
+                        self.remove_cutoff = float(self.config['Probe Selection']["remove_cutoff"])
+                    except KeyError:
+                        print("""\n>>> WARNING: "remove_cutoff" was not provided! Using the default value (1 Angstrom)!...\n""")
+                        self.remove_cutoff = 1
+                else:
+                    self.remove_cutoff = 0
+            except:
+                raise ValueError(""">>> ERROR: "remove_self" must be a boolean!\n)""")
+                self.remove_self = False
+                self.remove_cutoff = 0
+        else:
+            self.remove_self = False
+            self.remove_cutoff = 0
+
+
+    def parse_solvent(self):
+        try:
+            self.include_solvent = str(self.config['Solvent']["include_solvent"].strip('"')).lower()
+            if self.include_solvent == "true":
+                self.include_solvent = True
+                try:
+                    self.solvent_selection = str(self.config['Solvent']["solvent_selection"].strip('"'))
+                except:
+                    raise ValueError(""">>> ERROR: solvent_selection must be a valid MDanalysis selection!\n""")
+                try:
+                    self.solvent_cutoff = float(self.config['Solvent']["solvent_cutoff"])
+                except:
+                    raise ValueError(""">>> ERROR: solvent_cutoff must be a number!\n""")
+            else:
+                self.include_solvent = False
+        except:
+            self.include_solvent = False
+
+    def parse_time(self):
+        try:
+            self.dt = int(self.config['Time']['dt'])
+        except:
+            self.dt = 1
+
+
+    def parse_box_info(self):
+        try:
+            if 'Box Info' in self.config.sections() and 'redefine_box' in self.config['Box Info']:
+                redefine_box = str(self.config['Box Info']['redefine_box'].strip('"')).lower()
+                self.redefine_box = True if redefine_box == 'true' else False
+
+                if self.redefine_box:
+                    boxdimensions_str = str(self.config['Box Info']['boxdimensions'].strip('][')).lower()
+                    self.boxdimensions = [float(item) for item in boxdimensions_str.split(',')]
+
+                    if len(self.boxdimensions) != 6:
+                        raise ValueError('boxdimensions should contain [a,b,c, alpha, beta, gamma]!')
+                else:
+                    self.redefine_box = False
+                    self.boxdimensions = None
+            else:
+                self.redefine_box = False
+                self.boxdimensions = None
+        except:
+            raise ValueError('To redefine the box, boxdimensions must be provided!')

lib/format.py

@@ -0,0 +1,244 @@
+#!/usr/bin/env python3
+# -*- coding:utf-8 -*-
+#Marcelo D. Poleto
+#DEC 2021
+
+import os
+from lib.utils import *
+
+
+##############################################################
+# Opening output files
+#
+# - run_info.log
+# - ElecField.dat
+# - ElecField_per_residue.dat
+# - Spatial_Deviation.dat
+# - probe_info.dat
+#
+# if bond mode:
+# - ElecField_proj_onto_bond.dat
+# - ElecField_alignment.dat
+# - bond_axis_info.dat
+#
+##############################################################
+
+def print_run_info(parsed_configuration, top_file, traj_file, outdir):
+
+	##############################
+	# Defining variables
+	sele_elecfield    = parsed_configuration.sele_elecfield
+
+	mode              = parsed_configuration.mode
+	if mode == "atom":
+	    selatom = parsed_configuration.selatom
+	    remove_cutoff = 0
+	elif mode == "bond":
+	    selbond1 = parsed_configuration.selbond1
+	    selbond2 = parsed_configuration.selbond2
+	    remove_cutoff = 0
+	elif mode == "coordinate":
+	    probecoordinate = parsed_configuration.probecoordinate
+	    remove_self = parsed_configuration.remove_self
+	    remove_cutoff = parsed_configuration.remove_cutoff
+	elif mode == "list":
+	    file_of_coordinates = parsed_configuration.file_of_coordinates
+	    remove_self = parsed_configuration.remove_self
+	    remove_cutoff = parsed_configuration.remove_cutoff
+
+	include_solvent   = parsed_configuration.include_solvent
+	if include_solvent:
+	    solvent_selection = parsed_configuration.solvent_selection
+	    solvent_cutoff = parsed_configuration.solvent_cutoff
+
+	dt = parsed_configuration.dt
+
+	redefine_box = parsed_configuration.redefine_box
+	if redefine_box:
+	    boxdimensions = parsed_configuration.boxdimensions
+
+
+	##############################
+	# Defining probe string for each calculation mode
+	probestr = ""
+	if mode == "atom":
+		str1 = 'selatom             = {}'.format(selatom)
+		probestr += str1
+	elif mode == "bond":
+		str1 = 'selbond1            = {}'.format(selbond1)
+		str2 = 'selbond2            = {}'.format(selbond2)
+		probestr += str1 +'\n'+str2
+	elif mode == "coordinate":
+		str1 = 'probecoordinate     = {}'.format(probecoordinate)
+		probestr += str1
+	elif mode == "list":
+		str1 = 'file_of_coordinates = {}'.format(file_of_coordinates)
+		probestr += str1+'\n'
+		if remove_self == True:
+			str2 = 'remove_self         = {}'.format(remove_self)
+			str3 = 'remove_cutoff       = {}'.format(remove_cutoff)
+			probestr += str2 +'\n' + str3
+	solvstr = ""
+	if include_solvent == True:
+		str1 = 'solvent_selection   = {}'.format(solvent_selection)
+		str2 = 'solvent_cutoff      = {}'.format(solvent_cutoff)
+		solvstr += str1 +"\n"+str2
+	boxstr = ""
+	if redefine_box == True:
+		str1 = 'boxdimensions       = {}'.format(boxdimensions)
+		boxstr += str1
+
+
+	##############################################################
+	# Write the parameters to a file
+
+	if not os.path.exists(outdir):
+		os.system("mkdir " + outdir)
+
+	# Writing run info into a file
+	out = open(os.path.join(outdir,"run_info.log"),'w')
+
+	runinfo = """########################################################
+>>> Parameters used to run TUPÃ:
+
+Topology file       = {top}
+Trajectory file     = {traj}
+
+[Environment Selection]
+sele_environment    = {environment}
+
+[Probe Selection]
+mode                = {mode}
+{probestr}
+
+[Solvent]
+include_solvent     = {include_solv}
+{solvstr}
+
+[Time]
+dt                  = {dt}
+
+[Box Info]
+redefine_box        = {box}
+{boxstr}
+########################################################
+	""".format(top = top_file, traj = traj_file,
+	           environment = sele_elecfield,
+	           mode=mode,
+	           probestr=probestr,
+	           include_solv=include_solvent,
+	           solvstr=solvstr,
+	           dt=dt,
+	           box=redefine_box,
+	           boxstr=boxstr)
+
+	out.write(runinfo)
+	out.close()
+
+	print(runinfo)
+
+
+def open_outputs(outdir):
+	out = open(os.path.join(outdir,"ElecField.dat"),'w')
+	out.write("""@    title "Electric Field"\n@    xaxis  label "Time (ps)"\n@    yaxis  label "MV/cm"\n""")
+	out.write("#time             Magnitude      Efield_X       Efield_Y       Efield_Z\n")
+	out.write("@type xy\n")
+
+	outres = open(os.path.join(outdir,"ElecField_per_residue.dat"),'w')
+	outres.write("""@    title "Magnitude of Electric Field"\n@    xaxis  label "Residues"\n@    yaxis  label "MV/cm"\n""")
+	outres.write("#time              Ef_per_res     Std.Dev        Alignment%     Alignment_Std\n")
+	outres.write("@type xydy\n")
+
+	outangle = open(os.path.join(outdir,"Spatial_Deviation.dat"), "w")
+	outangle.write("#time   Angle(Efield(t), avg_field)   Projection(Efield(t), avg_field)   Alignment(Efield(t), avg_field)\n")
+
+	outprobe = open(os.path.join(outdir,"probe_info.dat"), "w")
+	outprobe.write("#time     probe_coordinates\n")
+
+	return out, outres, outangle, outprobe
+
+
+def open_outputs_bond(outdir):
+	outproj = open(os.path.join(outdir,"ElecField_proj_onto_bond.dat"),'w')
+	outproj.write("""@    title "Efield_Projection"\n@    xaxis  label "Time (ps)"\n@    yaxis  label "MV/cm"\n""")
+	outproj.write("#time             Magnitude       Efield_X       Efield_Y      Efield_Z        Direction\n")
+	outproj.write("@type xy\n")
+
+	outalig = open(os.path.join(outdir,"ElecField_alignment.dat"),'w')
+	outalig.write("""@    title "Efield_Projection"\n@    xaxis  label "Time (ps)"\n@    yaxis  label "Alignment_percentage"\n""")
+	outalig.write("#time             Angle           Alignment%\n")
+	outalig.write("@type xy\n")
+
+	outresbond = open(os.path.join(outdir,"Residue_alignment_to_bond.dat"),'w')
+	outresbond.write("""@    title "Magnitude of Electric Field"\n@    xaxis  label "Residues"\n@    yaxis  label "MV/cm"\n""")
+	outresbond.write("#time              Efproj_res     Std.Dev        AVG_align%     Std_align%\n")
+	outresbond.write("@type xydy\n")
+
+	outrbond = open(os.path.join(outdir,"bond_axis_info.dat"), "w")
+	outrbond.write("#time       rbondvec_in_meters                                   rbond_magnitude             rbond_unit_vector\n")
+
+	return outproj, outalig, outrbond, outresbond
+
+
+##############################################################
+# Formatting output lines
+#
+# - ElecField.dat
+# - ElecField_per_residue.dat
+# - Spatial_Deviation.dat
+#
+# if bond mode:
+# - ElecField_proj_onto_bond.dat
+# - ElecField_alignment.dat
+# - bond_axis_info.dat
+#
+##############################################################
+def format_efield_out_line(field_array, t=None, lastline = False):
+	if not lastline:
+		string = "{a:<12}{b: 15.6f}{c: 15.6f}{d: 15.6f}{e: 15.6f}\n".format(a=t, b=field_array[0], c=field_array[1], d=field_array[2],e=field_array[3])
+	else:
+		string = "{a: 15.6f}{b: 15.6f}{c: 15.6f}{d: 15.6f}\n".format(a=field_array[0], b=field_array[1], c=field_array[2], d=field_array[3])
+	return string
+
+
+def format_res_out_line(field_array):
+	string = "{a:<12}{b: 15.6f}{c: 15.6f}{d: 15.6f}{e: 15.6f}\n".format(a=field_array[0], b=field_array[1], c=field_array[2], d=field_array[3],e=field_array[4])
+	return string
+
+
+def format_spatialdev_out_line(field_array,t=None, lastline=False):
+	if lastline == False:
+		string = "{a:<12}{b: 15.6f}{c: 30.6f}{d: 33.6f}\n".format(a=t, b=field_array[0], c=field_array[1], d=field_array[2])
+	else:
+		string = "{a: 15.6f}{b: 30.6f}{c: 33.6f}\n".format(a=field_array[0], b=field_array[1], c=field_array[2])
+	return string
+
+
+def format_proj_out_line(field_array, t=None, lastline = False, optarr=None):
+	if not lastline:
+		string = "{a:<12}{b: 15.6f}{c: 15.6f}{d: 15.6f}{e: 15.6f}{f: 15.6f}\n".format(a=t, b=field_array[0], c=field_array[1], d=field_array[2], e=field_array[3], f=optarr[0])
+	else:
+		string = "{a: 15.6f}{b: 15.6f}{c: 15.6f}{d: 15.6f}\n".format(a=field_array[0], b=field_array[1], c=field_array[2], d=field_array[3])
+	return string
+
+
+def format_align_out_line(field_array, lastline = False):
+	if not lastline:
+		string = "{a:<12}{b: 15.6f}{c: 15.6f}\n".format(a=field_array[0], b=field_array[1], c=field_array[2])
+	else:
+		string = "{a: 15.6f}{b: 15.6f}\n".format(a=field_array[0], b=field_array[1])
+	return string
+
+
+def format_rbond_info(field_array):
+	rvx, rvy, rvz = field_array[1]
+	rhx, rhy, rhz = field_array[2]
+	list1 = "[" + str(rvx) + "," + str(rvy) + "," + str(rvz) + "]"
+	list2 = "[" + str(rhx) + "," + str(rhy) + "," + str(rhz) + "]"
+	string = "{a:<12}{b:50s}{c: 11.6f}                    {d:<50s}\n".format(a=field_array[0], b=list1, c=field_array[3], d=list2)
+	return string
+
+def format_probe_position(time, array):
+	string = "[" + str(array[0]) + ", " + str(array[1]) + ", " + str(array[2]) + "]"
+	string = str(time).ljust(10,' ') +  string.ljust(60,' ') + "\n"
+	return string