13726590978

Committed 07 Mar 2025 06:07PM UTC coverage: 1.726%. First build

Build # 13726590978

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #45 from CCPBioSim/39-implement-ci-pipeline-pre-commit-hooks

Implement CI pipeline pre-commit hooks

Coverage Stats

27 of 1135 new or added lines in 23 files covered. (2.38%)

53 of 3071 relevant lines covered (1.73%)

0.05 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/CodeEntropy/main_mcc.py

import argparse
import math

import MDAnalysis as mda

# import numpy as np
import pandas as pd

from CodeEntropy import EntropyFunctions as EF
from CodeEntropy import LevelFunctions as LF
from CodeEntropy import MDAUniverseHelper as MDAHelper

# from datetime import datetime


def main():
    """
    Main function for calculating the entropy of a system using the multiscale cell
    correlation method.
    """

    try:
        parser = argparse.ArgumentParser(
            description="""
        CodeEntropy-POSEIDON is a tool to compute entropy using the
        multiscale-cell-correlation (MCC) theory and force/torque covariance
        methods with the ablity to compute solvent entropy.
        Version:
            0.3.1;

        Authors:
            Arghya Chakravorty (arghya90),
            Jas Kalayan (jkalayan),
            Donald Chang,
            Sarah Fegan
            Ioana Papa;

        Output:
            *.csv = results from different calculateion,
            *.pkl - Pickled reduced universe for further analysis,
            *.out - detailed output such as matrix and spectra"""
        )

        parser.add_argument(
            "-f",
            "--top_traj_file",
            required=True,
            dest="filePath",
            action="store",
            nargs="+",
            help="Path to Structure/topology file (AMBER PRMTOP, GROMACS TPR which "
            "contains topology and dihedral information) followed by Trajectory "
            "file(s) (AMBER NETCDF or GROMACS TRR) you will need to output the "
            "coordinates and forces to the same file. Required.",
        )
        parser.add_argument(
            "-l",
            "--selectString",
            action="store",
            dest="selection_string",
            type=str,
            default="all",
            help="Selection string for CodeEntropy such as protein or resid, refer to "
            "MDAnalysis.select_atoms for more information.",
        )
        parser.add_argument(
            "-b",
            "--begin",
            action="store",
            dest="start",
            help="Start analysing the trajectory from this frame index. Defaults to 0",
            default=0,
            type=int,
        )
        parser.add_argument(
            "-e",
            "--end",
            action="store",
            dest="end",
            help="Stop analysing the trajectory at this frame index. Defaults to -1 "
            "(end of trajectory file)",
            default=-1,
            type=int,
        )
        parser.add_argument(
            "-d",
            "--step",
            action="store",
            dest="step",
            help="interval between two consecutive frames to be read index. "
            "Defaults to 1",
            default=1,
            type=int,
        )
        parser.add_argument(
            "-n",
            "--bin_width",
            action="store",
            dest="bin_width",
            default=30,
            type=int,
            help="Bin width in degrees for making the histogram of the dihedral angles "
            "for the conformational entropy. Default: 30",
        )
        parser.add_argument(
            "-k",
            "--tempra",
            action="store",
            dest="temp",
            help="Temperature for entropy calculation (K). Default to 298.0 K",
            default=298.0,
            type=float,
        )
        parser.add_argument(
            "-v",
            "--verbose",
            action="store",
            dest="verbose",
            default=False,
            type=bool,
            help="True/False flag for noisy or quiet output. Default: False",
        )
        parser.add_argument(
            "-t",
            "--thread",
            action="store",
            dest="thread",
            help="How many multiprocess to use. Default 1 for single core execution.",
            default=1,
            type=int,
        )
        parser.add_argument(
            "-o",
            "--out",
            action="store",
            dest="outfile",
            default="outfile.out",
            help="Name of the file where the output will be written. "
            "Default: outfile.out",
        )
        parser.add_argument(
            "-r",
            "--resout",
            action="store",
            dest="resfile",
            default="res_outfile.out",
            help="Name of the file where the residue entropy output will be written. "
            "Default: res_outfile.out",
        )
        parser.add_argument(
            "-m",
            "--mout",
            action="store",
            dest="moutfile",
            default=None,
            help="Name of the file where certain matrices will be written "
            "(default: None).",
        )

        parser.add_argument(
            "-c",
            "--cutShell",
            action="store",
            dest="cutShell",
            default=None,
            type=float,
            help="include cutoff shell analysis, add cutoff distance in angstrom "
            "Default None will ust the RAD Algorithm",
        )
        parser.add_argument(
            "-p",
            "--pureAtomNum",
            action="store",
            dest="puteAtomNum",
            default=1,
            type=int,
            help="Reference molecule resid for system of pure liquid. " "Default to 1",
        )
        parser.add_argument(
            "-x",
            "--excludedResnames",
            dest="excludedResnames",
            action="store",
            nargs="+",
            default=None,
            help="exclude a list of molecule names from nearest non-like analysis. "
            "Default: None. Multiples are gathered into list.",
        )
        parser.add_argument(
            "-w",
            "--water",
            dest="waterResnames",
            action="store",
            default="WAT",
            nargs="+",
            help="resname for water molecules. "
            "Default: WAT. Multiples are gathered into list.",
        )
        parser.add_argument(
            "-s",
            "--solvent",
            dest="solventResnames",
            action="store",
            nargs="+",
            default=None,
            help="include resname of solvent molecules (case-sensitive) "
            "Default: None. Multiples are gathered into list.",
        )
        parser.add_argument(
            "--solContact",
            action="store_true",
            dest="doSolContact",
            default=False,
            help="Do solute contact calculation",
        )

        args = parser.parse_args()
    except argparse.ArgumentError:
        print("Command line arguments are ill-defined, please check the arguments")
        raise

    # REPLACE INPUTS
    print("printing all input")
    for arg in vars(args):
        print(" {} {}".format(arg, getattr(args, arg) or ""))

    # startTime = datetime.now()

    # Get topology and trajectory file names and make universe
    tprfile = args.filePath[0]
    trrfile = args.filePath[1:]
    u = mda.Universe(tprfile, trrfile)

    # Define bin_width for histogram from inputs
    bin_width = args.bin_width

    # Define trajectory slicing from inputs
    start = args.start
    if start is None:
        start = 0
    end = args.end
    if end is None:
        end = -1
    step = args.step
    if step is None:
        step = 1
    # Count number of frames, easy if not slicing
    if start == 0 and end == -1 and step == 1:
        number_frames = len(u.trajectory)
    elif end == -1:
        end = len(u.trajectory)
        number_frames = math.floor((end - start) / step) + 1
    else:
        number_frames = math.floor((end - start) / step) + 1
    print(number_frames)

    # Create pandas data frame for results
    results_df = pd.DataFrame(columns=["Molecule ID", "Level", "Type", "Result"])
    residue_results_df = pd.DataFrame(
        columns=["Molecule ID", "Residue", "Type", "Result"]
    )

    # printing headings for output files
    with open(args.outfile, "a") as out:
        print("Molecule\tLevel\tType\tResult (J/mol/K)\n", file=out)

    with open(args.resfile, "a") as res:
        print("Molecule\tResidue\tType\tResult (J/mol/K)\n", file=res)

    # Reduce number of atoms in MDA universe to selection_string arg
    # (default all atoms included)
    if args.selection_string == "all":
        reduced_atom = u
    else:
        reduced_atom = MDAHelper.new_U_select_atom(u, args.selection_string)
        reduced_atom_name = f"{len(reduced_atom.trajectory)}_frame_dump_atom_selection"
        MDAHelper.write_universe(reduced_atom, reduced_atom_name)

    # Scan system for molecules and select levels (united atom, residue, polymer)
    # for each
    number_molecules, levels = LF.select_levels(reduced_atom, args.verbose)

    # Loop over molecules
    for molecule in range(number_molecules):
        # molecule data container of MDAnalysis Universe type for internal degrees
        # of freedom getting indices of first and last atoms in the molecule
        # assuming atoms are numbered consecutively and all atoms in a given
        # molecule are together
        index1 = reduced_atom.atoms.fragments[molecule].indices[0]
        index2 = reduced_atom.atoms.fragments[molecule].indices[-1]
        selection_string = f"index {index1}:{index2}"
        molecule_container = MDAHelper.new_U_select_atom(reduced_atom, selection_string)

        # Calculate entropy for each relevent level
        for level in levels[molecule]:
            if level == levels[molecule][-1]:
                highest_level = True
            else:
                highest_level = False

            if level == "united_atom":
                # loop over residues, report results per residue + total united atom
                # level. This is done per residue to reduce the size of the matrices -
                # amino acid resiudes have tens of united atoms but a whole protein
                # could have thousands. Doing the calculation per residue allows for
                # comparisons of contributions from different residues
                num_residues = len(molecule_container.residues)
                S_trans = 0
                S_rot = 0
                S_conf = 0
                for residue in range(num_residues):
                    # molecule data container of MDAnalysis Universe type for internal
                    # degrees of freedom getting indices of first and last atoms in the
                    # molecule assuming atoms are numbered consecutively and all atoms
                    # in a given molecule are together
                    index1 = molecule_container.residues[residue].atoms.indices[0]
                    index2 = molecule_container.residues[residue].atoms.indices[-1]
                    selection_string = f"index {index1}:{index2}"
                    residue_container = MDAHelper.new_U_select_atom(
                        molecule_container, selection_string
                    )

                    # Vibrational entropy at every level
                    # Get the force and torque matrices for the beads at the relevant
                    # level
                    force_matrix, torque_matrix = LF.get_matrices(
                        residue_container,
                        level,
                        args.verbose,
                        start,
                        end,
                        step,
                        number_frames,
                        highest_level,
                    )

                    # Calculate the entropy from the diagonalisation of the matrices
                    S_trans_residue = EF.vibrational_entropy(
                        force_matrix, "force", args.temp, highest_level
                    )
                    S_trans += S_trans_residue
                    print(f"S_trans_{level}_{residue} = {S_trans_residue}")
                    new_row = pd.DataFrame(
                        {
                            "Molecule ID": [molecule],
                            "Residue": [residue],
                            "Type": ["Transvibrational (J/mol/K)"],
                            "Result": [S_trans_residue],
                        }
                    )
                    residue_results_df = pd.concat(
                        [residue_results_df, new_row], ignore_index=True
                    )
                    with open(args.resfile, "a") as res:
                        print(
                            molecule,
                            "\t",
                            residue,
                            "\tTransvibration\t",
                            S_trans_residue,
                            file=res,
                        )

                    S_rot_residue = EF.vibrational_entropy(
                        torque_matrix, "torque", args.temp, highest_level
                    )
                    S_rot += S_rot_residue
                    print(f"S_rot_{level}_{residue} = {S_rot_residue}")
                    new_row = pd.DataFrame(
                        {
                            "Molecule ID": [molecule],
                            "Residue": [residue],
                            "Type": ["Rovibrational (J/mol/K)"],
                            "Result": [S_rot_residue],
                        }
                    )
                    residue_results_df = pd.concat(
                        [residue_results_df, new_row], ignore_index=True
                    )
                    with open(args.resfile, "a") as res:
                        #  print(new_row, file=res)
                        print(
                            molecule,
                            "\t",
                            residue,
                            "\tRovibrational \t",
                            S_rot_residue,
                            file=res,
                        )

                    # Conformational entropy based on atom dihedral angle distributions
                    # Gives entropy of conformations within each residue

                    # Get dihedral angle distribution
                    dihedrals = LF.get_dihedrals(residue_container, level)

                    # Calculate conformational entropy
                    S_conf_residue = EF.conformational_entropy(
                        residue_container,
                        dihedrals,
                        bin_width,
                        start,
                        end,
                        step,
                        number_frames,
                    )
                    S_conf += S_conf_residue
                    print(f"S_conf_{level}_{residue} = {S_conf_residue}")
                    new_row = pd.DataFrame(
                        {
                            "Molecule ID": [molecule],
                            "Residue": [residue],
                            "Type": ["Conformational (J/mol/K)"],
                            "Result": [S_conf_residue],
                        }
                    )
                    residue_results_df = pd.concat(
                        [residue_results_df, new_row], ignore_index=True
                    )
                    with open(args.resfile, "a") as res:
                        print(
                            molecule,
                            "\t",
                            residue,
                            "\tConformational\t",
                            S_conf_residue,
                            file=res,
                        )

                # Print united atom level results summed over all residues
                print(f"S_trans_{level} = {S_trans}")
                new_row = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Transvibrational (J/mol/K)"],
                        "Result": [S_trans],
                    }
                )
                with open(args.outfile, "a") as out:
                    print(
                        molecule, "\t", level, "\tTransvibration\t", S_trans, file=out
                    )

                results_df = pd.concat([results_df, new_row], ignore_index=True)

                print(f"S_rot_{level} = {S_rot}")
                new_row = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Rovibrational (J/mol/K)"],
                        "Result": [S_rot],
                    }
                )
                results_df = pd.concat([results_df, new_row], ignore_index=True)
                with open(args.outfile, "a") as out:
                    print(molecule, "\t", level, "\tRovibrational \t", S_rot, file=out)

                print(f"S_conf_{level} = {S_conf}")
                new_row = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Conformational (J/mol/K)"],
                        "Result": [S_conf],
                    }
                )
                results_df = pd.concat([results_df, new_row], ignore_index=True)
                with open(args.outfile, "a") as out:
                    print(molecule, "\t", level, "\tConformational\t", S_conf, file=out)

                # End united atom vibrational and conformational calculations #

            if level in ("polymer", "residue"):
                # Vibrational entropy at every level
                # Get the force and torque matrices for the beads at the relevant level
                force_matrix, torque_matrix = LF.get_matrices(
                    molecule_container,
                    level,
                    args.verbose,
                    start,
                    end,
                    step,
                    number_frames,
                    highest_level,
                )

                # Calculate the entropy from the diagonalisation of the matrices
                S_trans = EF.vibrational_entropy(
                    force_matrix, "force", args.temp, highest_level
                )
                print(f"S_trans_{level} = {S_trans}")
                new_row = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Transvibrational (J/mol/K)"],
                        "Result": [S_trans],
                    }
                )
                results_df = pd.concat([results_df, new_row], ignore_index=True)
                with open(args.outfile, "a") as out:
                    print(
                        molecule, "\t", level, "\tTransvibrational\t", S_trans, file=out
                    )

                S_rot = EF.vibrational_entropy(
                    torque_matrix, "torque", args.temp, highest_level
                )
                print(f"S_rot_{level} = {S_rot}")
                new_row = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Rovibrational (J/mol/K)"],
                        "Result": [S_rot],
                    }
                )
                results_df = pd.concat([results_df, new_row], ignore_index=True)
                with open(args.outfile, "a") as out:
                    print(molecule, "\t", level, "\tRovibrational \t", S_rot, file=out)

                # Note: conformational entropy is not calculated at the polymer level,
                # because there is at most one polymer bead per molecule so no dihedral
                # angles.

            if level == "residue":
                # Conformational entropy based on distributions of dihedral angles
                # of residues. Gives conformational entropy of secondary structure

                # Get dihedral angle distribution
                dihedrals = LF.get_dihedrals(molecule_container, level)
                # Calculate conformational entropy
                S_conf = EF.conformational_entropy(
                    molecule_container,
                    dihedrals,
                    bin_width,
                    start,
                    end,
                    step,
                    number_frames,
                )
                print(f"S_conf_{level} = {S_conf}")
                new_row = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Conformational (J/mol/K)"],
                        "Result": [S_conf],
                    }
                )
                results_df = pd.concat([results_df, new_row], ignore_index=True)
                with open(args.outfile, "a") as out:
                    print(molecule, "\t", level, "\tConformational\t", S_conf, file=out)

            # Orientational entropy based on network of neighbouring molecules,
            #  only calculated at the highest level (whole molecule)
            #    if highest_level:
            #        neigbours = LF.get_neighbours(reduced_atom, molecule)
            #        S_orient = EF.orientational_entropy(neighbours)
            #        print(f"S_orient_{level} = {S_orient}")
            #        new_row = pd.DataFrame({
            #            'Molecule ID': [molecule],
            #            'Level': [level],
            #            'Type':['Orientational (J/mol/K)'],
            #            'Result': [S_orient],})
            #        results_df = pd.concat([results_df, new_row], ignore_index=True)
            #        with open(args.outfile, "a") as out:
            #    print(molecule,
            #          "\t",
            #          level,
            #          "\tOrientational\t",
            #          S_orient,
            #          file=out)

        # Report total entropy for the molecule
        S_molecule = results_df[results_df["Molecule ID"] == molecule]["Result"].sum()
        print(f"S_molecule = {S_molecule}")
        new_row = pd.DataFrame(
            {
                "Molecule ID": [molecule],
                "Level": ["Molecule Total"],
                "Type": ["Molecule Total Entropy "],
                "Result": [S_molecule],
            }
        )
        results_df = pd.concat([results_df, new_row], ignore_index=True)
        with open(args.outfile, "a") as out:
            print(molecule, "\t Molecule\tTotal Entropy\t", S_molecule, file=out)


# END main function

if __name__ == "__main__":

    main()

NEW 1	import argparse	×
2	import math	×
3
4	import MDAnalysis as mda	×
5
6	# import numpy as np
7	import pandas as pd	×
8
9	from CodeEntropy import EntropyFunctions as EF	×
NEW 10	from CodeEntropy import LevelFunctions as LF	×
11	from CodeEntropy import MDAUniverseHelper as MDAHelper	×
12
13	# from datetime import datetime
14
15
16	def main():	×
17	"""
18	Main function for calculating the entropy of a system using the multiscale cell
19	correlation method.
20	"""
21
22	try:	×
NEW 23	parser = argparse.ArgumentParser(	×
24	description="""
25	CodeEntropy-POSEIDON is a tool to compute entropy using the
26	multiscale-cell-correlation (MCC) theory and force/torque covariance
27	methods with the ablity to compute solvent entropy.
28	Version:
29	0.3.1;
30
31	Authors:
32	Arghya Chakravorty (arghya90),
33	Jas Kalayan (jkalayan),
34	Donald Chang,
35	Sarah Fegan
36	Ioana Papa;
37
38	Output:
39	*.csv = results from different calculateion,
40	*.pkl - Pickled reduced universe for further analysis,
41	*.out - detailed output such as matrix and spectra"""
42	)
43
NEW 44	parser.add_argument(	×
45	"-f",
46	"--top_traj_file",
47	required=True,
48	dest="filePath",
49	action="store",
50	nargs="+",
51	help="Path to Structure/topology file (AMBER PRMTOP, GROMACS TPR which "
52	"contains topology and dihedral information) followed by Trajectory "
53	"file(s) (AMBER NETCDF or GROMACS TRR) you will need to output the "
54	"coordinates and forces to the same file. Required.",
55	)
NEW 56	parser.add_argument(	×
57	"-l",
58	"--selectString",
59	action="store",
60	dest="selection_string",
61	type=str,
62	default="all",
63	help="Selection string for CodeEntropy such as protein or resid, refer to "
64	"MDAnalysis.select_atoms for more information.",
65	)
NEW 66	parser.add_argument(	×
67	"-b",
68	"--begin",
69	action="store",
70	dest="start",
71	help="Start analysing the trajectory from this frame index. Defaults to 0",
72	default=0,
73	type=int,
74	)
NEW 75	parser.add_argument(	×
76	"-e",
77	"--end",
78	action="store",
79	dest="end",
80	help="Stop analysing the trajectory at this frame index. Defaults to -1 "
81	"(end of trajectory file)",
82	default=-1,
83	type=int,
84	)
NEW 85	parser.add_argument(	×
86	"-d",
87	"--step",
88	action="store",
89	dest="step",
90	help="interval between two consecutive frames to be read index. "
91	"Defaults to 1",
92	default=1,
93	type=int,
94	)
NEW 95	parser.add_argument(	×
96	"-n",
97	"--bin_width",
98	action="store",
99	dest="bin_width",
100	default=30,
101	type=int,
102	help="Bin width in degrees for making the histogram of the dihedral angles "
103	"for the conformational entropy. Default: 30",
104	)
NEW 105	parser.add_argument(	×
106	"-k",
107	"--tempra",
108	action="store",
109	dest="temp",
110	help="Temperature for entropy calculation (K). Default to 298.0 K",
111	default=298.0,
112	type=float,
113	)
NEW 114	parser.add_argument(	×
115	"-v",
116	"--verbose",
117	action="store",
118	dest="verbose",
119	default=False,
120	type=bool,
121	help="True/False flag for noisy or quiet output. Default: False",
122	)
NEW 123	parser.add_argument(	×
124	"-t",
125	"--thread",
126	action="store",
127	dest="thread",
128	help="How many multiprocess to use. Default 1 for single core execution.",
129	default=1,
130	type=int,
131	)
NEW 132	parser.add_argument(	×
133	"-o",
134	"--out",
135	action="store",
136	dest="outfile",
137	default="outfile.out",
138	help="Name of the file where the output will be written. "
139	"Default: outfile.out",
140	)
NEW 141	parser.add_argument(	×
142	"-r",
143	"--resout",
144	action="store",
145	dest="resfile",
146	default="res_outfile.out",
147	help="Name of the file where the residue entropy output will be written. "
148	"Default: res_outfile.out",
149	)
NEW 150	parser.add_argument(	×
151	"-m",
152	"--mout",
153	action="store",
154	dest="moutfile",
155	default=None,
156	help="Name of the file where certain matrices will be written "
157	"(default: None).",
158	)
159
NEW 160	parser.add_argument(	×
161	"-c",
162	"--cutShell",
163	action="store",
164	dest="cutShell",
165	default=None,
166	type=float,
167	help="include cutoff shell analysis, add cutoff distance in angstrom "
168	"Default None will ust the RAD Algorithm",
169	)
NEW 170	parser.add_argument(	×
171	"-p",
172	"--pureAtomNum",
173	action="store",
174	dest="puteAtomNum",
175	default=1,
176	type=int,
177	help="Reference molecule resid for system of pure liquid. " "Default to 1",
178	)
NEW 179	parser.add_argument(	×
180	"-x",
181	"--excludedResnames",
182	dest="excludedResnames",
183	action="store",
184	nargs="+",
185	default=None,
186	help="exclude a list of molecule names from nearest non-like analysis. "
187	"Default: None. Multiples are gathered into list.",
188	)
NEW 189	parser.add_argument(	×
190	"-w",
191	"--water",
192	dest="waterResnames",
193	action="store",
194	default="WAT",
195	nargs="+",
196	help="resname for water molecules. "
197	"Default: WAT. Multiples are gathered into list.",
198	)
NEW 199	parser.add_argument(	×
200	"-s",
201	"--solvent",
202	dest="solventResnames",
203	action="store",
204	nargs="+",
205	default=None,
206	help="include resname of solvent molecules (case-sensitive) "
207	"Default: None. Multiples are gathered into list.",
208	)
NEW 209	parser.add_argument(	×
210	"--solContact",
211	action="store_true",
212	dest="doSolContact",
213	default=False,
214	help="Do solute contact calculation",
215	)
216
217	args = parser.parse_args()	×
218	except argparse.ArgumentError:	×
NEW 219	print("Command line arguments are ill-defined, please check the arguments")	×
220	raise	×
221
222	# REPLACE INPUTS
223	print("printing all input")	×
224	for arg in vars(args):	×
NEW 225	print(" {} {}".format(arg, getattr(args, arg) or ""))	×
226
227	# startTime = datetime.now()
228
229	# Get topology and trajectory file names and make universe
230	tprfile = args.filePath[0]	×
231	trrfile = args.filePath[1:]	×
232	u = mda.Universe(tprfile, trrfile)	×
233
234	# Define bin_width for histogram from inputs
235	bin_width = args.bin_width	×
236
237	# Define trajectory slicing from inputs
238	start = args.start	×
239	if start is None:	×
240	start = 0	×
241	end = args.end	×
242	if end is None:	×
243	end = -1	×
244	step = args.step	×
245	if step is None:	×
246	step = 1	×
247	# Count number of frames, easy if not slicing
248	if start == 0 and end == -1 and step == 1:	×
249	number_frames = len(u.trajectory)	×
250	elif end == -1:	×
251	end = len(u.trajectory)	×
NEW 252	number_frames = math.floor((end - start) / step) + 1	×
253	else:
NEW 254	number_frames = math.floor((end - start) / step) + 1	×
255	print(number_frames)	×
256
257	# Create pandas data frame for results
NEW 258	results_df = pd.DataFrame(columns=["Molecule ID", "Level", "Type", "Result"])	×
NEW 259	residue_results_df = pd.DataFrame(	×
260	columns=["Molecule ID", "Residue", "Type", "Result"]
261	)
262
263	# printing headings for output files
264	with open(args.outfile, "a") as out:	×
265	print("Molecule\tLevel\tType\tResult (J/mol/K)\n", file=out)	×
266
267	with open(args.resfile, "a") as res:	×
268	print("Molecule\tResidue\tType\tResult (J/mol/K)\n", file=res)	×
269
270	# Reduce number of atoms in MDA universe to selection_string arg
271	# (default all atoms included)
NEW 272	if args.selection_string == "all":	×
273	reduced_atom = u	×
274	else:
275	reduced_atom = MDAHelper.new_U_select_atom(u, args.selection_string)	×
276	reduced_atom_name = f"{len(reduced_atom.trajectory)}_frame_dump_atom_selection"	×
277	MDAHelper.write_universe(reduced_atom, reduced_atom_name)	×
278
279	# Scan system for molecules and select levels (united atom, residue, polymer)
280	# for each
281	number_molecules, levels = LF.select_levels(reduced_atom, args.verbose)	×
282
283	# Loop over molecules
284	for molecule in range(number_molecules):	×
285	# molecule data container of MDAnalysis Universe type for internal degrees
286	# of freedom getting indices of first and last atoms in the molecule
287	# assuming atoms are numbered consecutively and all atoms in a given
288	# molecule are together
289	index1 = reduced_atom.atoms.fragments[molecule].indices[0]	×
290	index2 = reduced_atom.atoms.fragments[molecule].indices[-1]	×
291	selection_string = f"index {index1}:{index2}"	×
292	molecule_container = MDAHelper.new_U_select_atom(reduced_atom, selection_string)	×
293
294	# Calculate entropy for each relevent level
295	for level in levels[molecule]:	×
296	if level == levels[molecule][-1]:	×
297	highest_level = True	×
298	else:
299	highest_level = False	×
300
NEW 301	if level == "united_atom":	×
302	# loop over residues, report results per residue + total united atom
303	# level. This is done per residue to reduce the size of the matrices -
304	# amino acid resiudes have tens of united atoms but a whole protein
305	# could have thousands. Doing the calculation per residue allows for
306	# comparisons of contributions from different residues
307	num_residues = len(molecule_container.residues)	×
308	S_trans = 0	×
309	S_rot = 0	×
310	S_conf = 0	×
311	for residue in range(num_residues):	×
312	# molecule data container of MDAnalysis Universe type for internal
313	# degrees of freedom getting indices of first and last atoms in the
314	# molecule assuming atoms are numbered consecutively and all atoms
315	# in a given molecule are together
316	index1 = molecule_container.residues[residue].atoms.indices[0]	×
317	index2 = molecule_container.residues[residue].atoms.indices[-1]	×
318	selection_string = f"index {index1}:{index2}"	×
NEW 319	residue_container = MDAHelper.new_U_select_atom(	×
320	molecule_container, selection_string
321	)
322
323	# Vibrational entropy at every level
324	# Get the force and torque matrices for the beads at the relevant
325	# level
NEW 326	force_matrix, torque_matrix = LF.get_matrices(	×
327	residue_container,
328	level,
329	args.verbose,
330	start,
331	end,
332	step,
333	number_frames,
334	highest_level,
335	)
336
337	# Calculate the entropy from the diagonalisation of the matrices
NEW 338	S_trans_residue = EF.vibrational_entropy(	×
339	force_matrix, "force", args.temp, highest_level
340	)
341	S_trans += S_trans_residue	×
342	print(f"S_trans_{level}_{residue} = {S_trans_residue}")	×
NEW 343	new_row = pd.DataFrame(	×
344	{
345	"Molecule ID": [molecule],
346	"Residue": [residue],
347	"Type": ["Transvibrational (J/mol/K)"],
348	"Result": [S_trans_residue],
349	}
350	)
NEW 351	residue_results_df = pd.concat(	×
352	[residue_results_df, new_row], ignore_index=True
353	)
354	with open(args.resfile, "a") as res:	×
NEW 355	print(	×
356	molecule,
357	"\t",
358	residue,
359	"\tTransvibration\t",
360	S_trans_residue,
361	file=res,
362	)
363
NEW 364	S_rot_residue = EF.vibrational_entropy(	×
365	torque_matrix, "torque", args.temp, highest_level
366	)
367	S_rot += S_rot_residue	×
368	print(f"S_rot_{level}_{residue} = {S_rot_residue}")	×
NEW 369	new_row = pd.DataFrame(	×
370	{
371	"Molecule ID": [molecule],
372	"Residue": [residue],
373	"Type": ["Rovibrational (J/mol/K)"],
374	"Result": [S_rot_residue],
375	}
376	)
NEW 377	residue_results_df = pd.concat(	×
378	[residue_results_df, new_row], ignore_index=True
379	)
380	with open(args.resfile, "a") as res:	×
381	# print(new_row, file=res)
NEW 382	print(	×
383	molecule,
384	"\t",
385	residue,
386	"\tRovibrational \t",
387	S_rot_residue,
388	file=res,
389	)
390
391	# Conformational entropy based on atom dihedral angle distributions
392	# Gives entropy of conformations within each residue
393
394	# Get dihedral angle distribution
395	dihedrals = LF.get_dihedrals(residue_container, level)	×
396
397	# Calculate conformational entropy
NEW 398	S_conf_residue = EF.conformational_entropy(	×
399	residue_container,
400	dihedrals,
401	bin_width,
402	start,
403	end,
404	step,
405	number_frames,
406	)
407	S_conf += S_conf_residue	×
408	print(f"S_conf_{level}_{residue} = {S_conf_residue}")	×
NEW 409	new_row = pd.DataFrame(	×
410	{
411	"Molecule ID": [molecule],
412	"Residue": [residue],
413	"Type": ["Conformational (J/mol/K)"],
414	"Result": [S_conf_residue],
415	}
416	)
NEW 417	residue_results_df = pd.concat(	×
418	[residue_results_df, new_row], ignore_index=True
419	)
420	with open(args.resfile, "a") as res:	×
NEW 421	print(	×
422	molecule,
423	"\t",
424	residue,
425	"\tConformational\t",
426	S_conf_residue,
427	file=res,
428	)
429
430	# Print united atom level results summed over all residues
431	print(f"S_trans_{level} = {S_trans}")	×
NEW 432	new_row = pd.DataFrame(	×
433	{
434	"Molecule ID": [molecule],
435	"Level": [level],
436	"Type": ["Transvibrational (J/mol/K)"],
437	"Result": [S_trans],
438	}
439	)
440	with open(args.outfile, "a") as out:	×
NEW 441	print(	×
442	molecule, "\t", level, "\tTransvibration\t", S_trans, file=out
443	)
444
445	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
446
447	print(f"S_rot_{level} = {S_rot}")	×
NEW 448	new_row = pd.DataFrame(	×
449	{
450	"Molecule ID": [molecule],
451	"Level": [level],
452	"Type": ["Rovibrational (J/mol/K)"],
453	"Result": [S_rot],
454	}
455	)
456	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
457	with open(args.outfile, "a") as out:	×
NEW 458	print(molecule, "\t", level, "\tRovibrational \t", S_rot, file=out)	×
459
460	print(f"S_conf_{level} = {S_conf}")	×
NEW 461	new_row = pd.DataFrame(	×
462	{
463	"Molecule ID": [molecule],
464	"Level": [level],
465	"Type": ["Conformational (J/mol/K)"],
466	"Result": [S_conf],
467	}
468	)
469	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
470	with open(args.outfile, "a") as out:	×
NEW 471	print(molecule, "\t", level, "\tConformational\t", S_conf, file=out)	×
472
473	# End united atom vibrational and conformational calculations #
474
NEW 475	if level in ("polymer", "residue"):	×
476	# Vibrational entropy at every level
477	# Get the force and torque matrices for the beads at the relevant level
NEW 478	force_matrix, torque_matrix = LF.get_matrices(	×
479	molecule_container,
480	level,
481	args.verbose,
482	start,
483	end,
484	step,
485	number_frames,
486	highest_level,
487	)
488
489	# Calculate the entropy from the diagonalisation of the matrices
NEW 490	S_trans = EF.vibrational_entropy(	×
491	force_matrix, "force", args.temp, highest_level
492	)
493	print(f"S_trans_{level} = {S_trans}")	×
NEW 494	new_row = pd.DataFrame(	×
495	{
496	"Molecule ID": [molecule],
497	"Level": [level],
498	"Type": ["Transvibrational (J/mol/K)"],
499	"Result": [S_trans],
500	}
501	)
502	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
503	with open(args.outfile, "a") as out:	×
NEW 504	print(	×
505	molecule, "\t", level, "\tTransvibrational\t", S_trans, file=out
506	)
507
NEW 508	S_rot = EF.vibrational_entropy(	×
509	torque_matrix, "torque", args.temp, highest_level
510	)
511	print(f"S_rot_{level} = {S_rot}")	×
NEW 512	new_row = pd.DataFrame(	×
513	{
514	"Molecule ID": [molecule],
515	"Level": [level],
516	"Type": ["Rovibrational (J/mol/K)"],
517	"Result": [S_rot],
518	}
519	)
520	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
521	with open(args.outfile, "a") as out:	×
NEW 522	print(molecule, "\t", level, "\tRovibrational \t", S_rot, file=out)	×
523
524	# Note: conformational entropy is not calculated at the polymer level,
525	# because there is at most one polymer bead per molecule so no dihedral
526	# angles.
527
NEW 528	if level == "residue":	×
529	# Conformational entropy based on distributions of dihedral angles
530	# of residues. Gives conformational entropy of secondary structure
531
532	# Get dihedral angle distribution
533	dihedrals = LF.get_dihedrals(molecule_container, level)	×
534	# Calculate conformational entropy
NEW 535	S_conf = EF.conformational_entropy(	×
536	molecule_container,
537	dihedrals,
538	bin_width,
539	start,
540	end,
541	step,
542	number_frames,
543	)
544	print(f"S_conf_{level} = {S_conf}")	×
NEW 545	new_row = pd.DataFrame(	×
546	{
547	"Molecule ID": [molecule],
548	"Level": [level],
549	"Type": ["Conformational (J/mol/K)"],
550	"Result": [S_conf],
551	}
552	)
553	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
554	with open(args.outfile, "a") as out:	×
NEW 555	print(molecule, "\t", level, "\tConformational\t", S_conf, file=out)	×
556
557	# Orientational entropy based on network of neighbouring molecules,
558	# only calculated at the highest level (whole molecule)
559	# if highest_level:
560	# neigbours = LF.get_neighbours(reduced_atom, molecule)
561	# S_orient = EF.orientational_entropy(neighbours)
562	# print(f"S_orient_{level} = {S_orient}")
563	# new_row = pd.DataFrame({
564	# 'Molecule ID': [molecule],
565	# 'Level': [level],
566	# 'Type':['Orientational (J/mol/K)'],
567	# 'Result': [S_orient],})
568	# results_df = pd.concat([results_df, new_row], ignore_index=True)
569	# with open(args.outfile, "a") as out:
570	# print(molecule,
571	# "\t",
572	# level,
573	# "\tOrientational\t",
574	# S_orient,
575	# file=out)
576
577	# Report total entropy for the molecule
NEW 578	S_molecule = results_df[results_df["Molecule ID"] == molecule]["Result"].sum()	×
579	print(f"S_molecule = {S_molecule}")	×
NEW 580	new_row = pd.DataFrame(	×
581	{
582	"Molecule ID": [molecule],
583	"Level": ["Molecule Total"],
584	"Type": ["Molecule Total Entropy "],
585	"Result": [S_molecule],
586	}
587	)
588	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
589	with open(args.outfile, "a") as out:	×
NEW 590	print(molecule, "\t Molecule\tTotal Entropy\t", S_molecule, file=out)	×
591
592
593	# END main function
594
595	if __name__ == "__main__":	×
596
597	main()	×

CCPBioSim / CodeEntropy / 13726590978

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous