Penetrating potential evaluation of cyclotriphosphazenes target blood-brain barrier via computational toxicology methods.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches

[1] § Results › BBB penetrating mechanisms of CTPs › Stable binding of HPCTP with two SLC transporters ↔ GMXMMPBSA/amber_outputs.py, lines 1780–1899 · score 0.59 · van der Waals, polar solvation, electrostatic, Binding
[2] § STAR★Methods › Method details › Molecular dynamics simulations ↔ GMXMMPBSA/radii.py, lines 1–47 · score 0.54 · TIP3P, force field, solvation, Amber

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 2,028 lines · 99 KB · GPL-3.0 · 1 match

"""
This module contains all the classes and code to collect data and calculate
statistics from the output files of various calculation types. Each calculation
type needs its own class.
All data is stored in a special class derived from the list.
"""
# ##############################################################################
# GPLv3 LICENSE INFO #
# #
# Copyright (C) 2020 Mario S. Valdés-Tresanco and Mario E. Valdés-Tresanco #
# Copyright (C) 2014 Jason Swails, Bill Miller III, and Dwight McGee #
# #
# Project: https://github.com/Valdes-Tresanco-MS/gmx_MMPBSA #
# #
# This program is free software; you can redistribute it and/or modify it #
# under the terms of the GNU General Public License version 3 as published #
# by the Free Software Foundation. #
# #
# This program is distributed in the hope that it will be useful, but #
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY #
# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License #
# for more details. #
# ##############################################################################
import logging
from copy import deepcopy
from math import sqrt
from GMXMMPBSA.exceptions import (OutputError, LengthError, DecompError, GMXMMPBSA_ERROR)
from GMXMMPBSA.utils import EnergyVector, block_statistics, get_std
from types import SimpleNamespace
import numpy as np
import sys
import re
from csv import writer
idecompString = ['idecomp = 0: No decomposition analysis',
'idecomp = 1: Per-residue decomp adding 1-4 interactions to Internal.',
'idecomp = 2: Per-residue decomp adding 1-4 interactions to EEL and VDW.',
'idecomp = 3: Pairwise decomp adding 1-4 interactions to Internal.',
'idecomp = 4: Pairwise decomp adding 1-4 interactions to EEL and VDW.']
# Keep the summary rule aligned with the current eight-column statistics rows.
sep = '-' * 101
data_key_owner = {'BOND': ['GGAS', 'TOTAL'], 'ANGLE': ['GGAS', 'TOTAL'], 'DIHED': ['GGAS', 'TOTAL'],
'VDWAALS': ['GGAS', 'TOTAL'], 'EEL': ['GGAS', 'TOTAL'], '1-4 VDW': ['GGAS', 'TOTAL'],
'1-4 EEL': ['GGAS', 'TOTAL'],
# charmm
'UB': ['GGAS', 'TOTAL'], 'IMP': ['GGAS', 'TOTAL'], 'CMAP': ['GGAS', 'TOTAL'],
# non lineal PB
'EEL+EPB': ['TOTAL'],
# PB
'EPB': ['GSOLV', 'TOTAL'], 'ENPOLAR': ['GSOLV', 'TOTAL'], 'EDISPER': ['GSOLV', 'TOTAL'],
# GB
'EGB': ['GSOLV', 'TOTAL'], 'ESURF': ['GSOLV', 'TOTAL'],
# QM/GB
'ESCF': ['GGAS', 'TOTAL'],
# RISM
'POLAR SOLV': ['GSOLV', 'TOTAL'], 'APOLAR SOLV': ['GSOLV', 'TOTAL'],
'ERISM': ['GSOLV', 'TOTAL'],
# NMODE and QH
'TRANSLATIONAL': ['TOTAL'], 'ROTATIONAL': ['TOTAL'], 'VIBRATIONAL': ['TOTAL']
}
def _vector_statistics(vector):
"""Return legacy and block statistics for one energy vector."""
block_sd = block_sem = float('nan')
_, _, block_sd, block_sem = block_statistics(vector)
return (float(vector.mean()), float(vector.stdev()), float(vector.std()),
float(vector.semp()), float(vector.sem()), block_sd, block_sem)
class AmberOutput(dict):
"""
Base Amber output class. It takes a basename as a file name and parses
through all of the thread-specific output files (assumed to have the suffix
.# where # spans from 0 to num_files - 1
"""
print_levels = {'BOND': 2, 'ANGLE': 2, 'DIHED': 2, 'VDWAALS': 1, 'EEL': 1, '1-4 VDW': 2, '1-4 EEL': 2, 'EPOL': 1}
def __init__(self, mol: str, INPUT, chamber=False, **kwargs):
super(AmberOutput, self).__init__(**kwargs)
self.numframes = None
self.mol = mol
self.INPUT = INPUT
self.chamber = chamber
self.basename = None
self.num_files = None
self.frame_idx = 0
self.extraframe_idx = 0
self.is_read = False
self.apbs = INPUT['pb']['sander_apbs']
# This variable is used to get if the nmode calculation hasn't at least one frame
self.no_nmode_convergence = False
self.data_keys = ['BOND', 'ANGLE', 'DIHED', 'VDWAALS', 'EEL', '1-4 VDW', '1-4 EEL']
self.chamber_keys = ['CMAP', 'IMP', 'UB']
self.composite_keys = ['GGAS', 'GSOLV', 'TOTAL']
if self.chamber:
for key in self.chamber_keys:
if key not in self.data_keys:
self.data_keys.insert(3, key)
def parse_from_file(self, basename, num_files=1, numframes=1):
self.num_files = num_files
self.basename = basename
self.temperature = self.INPUT['general']['temperature']
self.numframes = numframes
for key in self.data_keys:
self[key] = EnergyVector(numframes)
for key in self.composite_keys:
self[key] = EnergyVector(numframes)
AmberOutput._read(self)
self._fill_composite_terms()
def _print_vectors(self, csvwriter):
""" Prints the energy vectors to a CSV file for easy viewing
in spreadsheets
"""
print_keys = list(self.data_keys)
# Add on the composite keys
print_keys += self.composite_keys
# write the header
csvwriter.writerow(['Frame #'] + print_keys)
# write out each frame
c = self.INPUT['nmode']['nmstartframe'] if self.__class__ == NMODEout else self.INPUT['general']['startframe']
for i in range(self.numframes):
csvwriter.writerow([c] + [round(self[key][i], 2) for key in print_keys])
c += self.INPUT['nmode']['nminterval'] if self.__class__ == NMODEout else self.INPUT['general']['interval']
def set_frame_range(self, start=None, end=None, interval=None):
d = deepcopy(self)
for key in d.data_keys:
d[key] = d[key][start:end:interval]
d._fill_composite_terms()
return d
def summary_output(self):
if not self.is_read:
raise OutputError('Cannot print summary before reading output files')
text = [f'{self.mol.capitalize()}:']
summary = self.summary()
for c, row in enumerate(summary, start=1):
key, avg, stdev, std, semp, sem, block_sd, block_sem = row
if key in ['GGAS', 'TOTAL']:
text.append('')
if isinstance(avg, str):
text.extend(
(
f'{key:16s} {avg:>13s} {stdev:>13s} {std:>10s} {semp:>12s} {sem:>10s} '
f'{block_sd:>10s} {block_sem:>10s}',
sep,
)
)
else:
text.append(f'{key:16s} {avg:13.2f} {stdev:13.2f} {std:10.2f} {semp:12.2f} {sem:10.2f} '
f'{block_sd:10.2f} {block_sem:10.2f}')
return '\n'.join(text) + '\n\n'
def summary(self):
""" Returns a formatted string that can be printed directly to the
output file
"""
if not self.is_read:
raise OutputError('Cannot print summary before reading output files')
comp_name = 'Entropy Component' if self.__class__ in [NMODEout, QHout] else 'Energy Component'
summary_list = [[comp_name, 'Average', 'SD(Prop.)', 'SD', 'SEM(Prop.)', 'SEM', 'Block SD', 'Block SEM']]
for key in self.data_keys:
# Skip the composite terms, since we print those at the end
if key in self.composite_keys:
continue
avg = float(self[key].mean())
stdev = float(self[key].stdev())
semp = float(self[key].semp())
std = float(self[key].std())
sem = float(self[key].sem())
_, _, block_sd, block_sem = block_statistics(self[key])
summary_list.append([key, avg, stdev, std, semp, sem, block_sd, block_sem])
for key in self.composite_keys:
# Now print out the composite terms
avg = float(self[key].mean())
stdev = float(self[key].stdev())
semp = float(self[key].semp())
std = float(self[key].std())
sem = float(self[key].sem())
_, _, block_sd, block_sem = block_statistics(self[key])
summary_list.append([key, avg, stdev, std, semp, sem, block_sd, block_sem])
return summary_list
def _read(self):
""" Internal reading function. This should be called at the end of __init__"""
if self.is_read:
return None # don't read through them twice
# Loop through all filenames
for fileno in range(self.num_files):
with open('%s.%d' % (self.basename, fileno)) as output_file:
self._get_energies(output_file)
self._extra_reading(fileno)
self._fill_nmode_values()
self.is_read = True
def _get_energies(self, output_file):
pass
def _extra_reading(self, fileno):
pass
def _fill_nmode_values(self):
pass
def _fill_composite_terms(self):
"""
Fills in the composite terms WITHOUT adding in terms we're not printing.
This should be called after the final verbosity level has been set (based
on whether or not certain terms need to be added in)
"""
for key in self.composite_keys:
self[key] = EnergyVector(self.numframes)
for key in self.data_keys:
for component in data_key_owner[key]:
self[component] = self[key] + self[component]
class IEout(dict):
"""
Interaction Entropy output
"""
def __init__(self, INPUT, method, **kwargs):
super(IEout, self).__init__(**kwargs)
self.INPUT = INPUT
self.method = method
self.block_analysis = []
def parse_from_dict(self, d: dict):
values = dict(d)
self.block_analysis = values.pop('block_analysis', [])
self.update(values)
def parse_from_file(self, filename, numframes=1):
self['data'] = EnergyVector(numframes)
with open(filename) as of:
c = 0
f = 0
while line := of.readline():
f += 1
if line.startswith('| BLOCK IE '):
fields = line.split()
self.block_analysis.append({
'block_size': int(fields[3]),
'nblocks': int(fields[4]),
'used_frames': int(fields[5]),
'mean': float(fields[6]),
'std': float(fields[7]),
'sem': float(fields[8]),
'p025': float(fields[9]),
'p975': float(fields[10]),
})
continue
if line.startswith('|') or not line.split():
# Legacy 1.6.x wrote "| Interaction Entropy (-TΔS): mean +/- sd"
# on a comment-style line; still parse that primary value.
stripped = line.lstrip('| ').strip()
if stripped.startswith('Interaction Entropy (-TΔS):') and 'ie_value' not in self:
tokens = stripped.split(':', 1)[1].split()
self['ie_value'] = float(tokens[0])
if len(tokens) >= 3 and tokens[1] == '+/-':
self['tail_std'] = float(tokens[2])
continue
if line.startswith('IE-frames:'):
self['ieframes'] = int(line.strip('\n').split()[-1])
elif line.startswith('Internal Energy SD (sigma):'):
self['sigma'] = float(line.strip('\n').split()[-1])
elif line.startswith('Full-ensemble Interaction Entropy (-TΔS):'):
self['ie_value'] = float(line.split()[-1])
elif line.startswith('Interaction Entropy (-TΔS):'):
# Legacy header without the leading comment marker.
tokens = line.split(':', 1)[1].split()
self['ie_value'] = float(tokens[0])
if len(tokens) >= 3 and tokens[1] == '+/-':
self['tail_std'] = float(tokens[2])
elif line.startswith('Tail convergence mean'):
self['tail_mean'] = float(line.split()[-1])
elif line.startswith('Tail convergence SD'):
self['tail_std'] = float(line.split()[-1])
elif line.startswith('Block diagnostic:'):
fields = line.split()
self['block_size'] = int(fields[3])
self['block_nblocks'] = int(fields[5])
self['block_std'] = float(fields[7])
self['block_sem'] = float(fields[9])
elif line.startswith('Frame'):
continue
else:
frame, value = line.strip('\n').split()
self['data'][c] = float(value)
c += 1
f += 1
self['iedata'] = self['data'][-self['ieframes']:]
# Prefer the explicit primary value. For legacy files without
# Full-ensemble / header text, restore the 1.6.x primary (tail mean).
self['ie_value'] = self.get('ie_value', float(self['iedata'].mean()))
self['tail_mean'] = self.get('tail_mean', float(self['iedata'].mean()))
self['tail_std'] = self.get('tail_std', float(self['iedata'].std()))
self['block_size'] = self.get('block_size', 0)
self['block_nblocks'] = self.get('block_nblocks', 0)
self['block_std'] = self.get('block_std', self['tail_std'])
self['block_sem'] = self.get(
'block_sem', self['tail_std'] / sqrt(self['ieframes'])
)
def _print_vectors(self, csvwriter):
""" Prints the energy vectors to a CSV file for easy viewing
in spreadsheets
"""
csvwriter.writerow(['Frame #', 'Interaction Entropy'])
f = self.INPUT['general']['startframe']
for d in self['data']:
csvwriter.writerow([f] + [round(d, 2)])
f += self.INPUT['general']['interval']
csvwriter.writerow([])
def summary_output(self):
summary = self.summary()
text = []
for row in summary:
met, key, sigma, avg, std, sem = row
if isinstance(avg, str):
text.extend((
f'{met:16s} {key:>13s} {sigma:>13s} {avg:>10s} {std:>12s} {sem:>10s}',
sep,
))
else:
text.append(f'{met:16s} {key:>13s} {sigma:13.2f} {avg:10.2f} {std:12.2f} {sem:10.2f}')
if self.get('block_size'):
text.append(
f"Block diagnostic: {self['block_nblocks']} nonoverlapping blocks "
f"of {self['block_size']} frames"
)
text.append(
f"Tail convergence diagnostic: {self['tail_mean']:.2f} +/- {self['tail_std']:.2f} "
f"over the last {self['ieframes']} prefixes"
)
return '\n'.join(text) + '\n\n'
def summary(self):
""" Formatted summary of Interaction Entropy results """
avg = float(self.get('ie_value', self['iedata'].mean()))
stdev = float(self.get('block_std', self['data'][-self['ieframes']:].stdev()))
sem = float(self.get('block_sem', self['data'][-self['ieframes']:].sem()))
return [
[
'Energy Method',
'Entropy',
'σ(Int. Energy)',
'Full IE',
'Block SD',
'Block SEM'
],
[self.method.upper(), 'IE', self['sigma'], avg, stdev, sem]
]
class C2out(dict):
"""
C2 Entropy output
"""
def __init__(self, method, **kwargs):
super(C2out, self).__init__(**kwargs)
self.method = method
self.block_analysis = []
def parse_from_dict(self, d):
values = dict(d)
self.block_analysis = values.pop('block_analysis', [])
self.update(values)
def parse_from_file(self, filename):
with open(filename) as of:
while line := of.readline():
if line.startswith('| BLOCK C2 '):
fields = line.split()
self.block_analysis.append({
'block_size': int(fields[3]),
'nblocks': int(fields[4]),
'used_frames': int(fields[5]),
'mean': float(fields[6]),
'std': float(fields[7]),
'sem': float(fields[8]),
'p025': float(fields[9]),
'p975': float(fields[10]),
})
continue
if line.startswith('|') or not line:
continue
if line.startswith('C2 Entropy (-TΔS):'):
self['c2data'] = float(line.strip('\n').split()[-1])
elif line.startswith(('C2 Entropy SD:', 'C2 Block SD:')):
self['c2_std'] = float(line.strip('\n').split()[-1])
elif line.startswith('C2 Block SEM:'):
self['c2_sem'] = float(line.strip('\n').split()[-1])
elif line.startswith('Internal Energy SD (sigma):'):
self['sigma'] = float(line.strip('\n').split()[-1])
elif line.startswith(('C2 Entropy CI:', 'C2 Block P2.5-P97.5:')):
self['c2_ci'] = [float(line.strip('\n').split()[-2]), float(line.strip('\n').split()[-1])]
elif line.startswith('Block diagnostic:'):
fields = line.split()
self['block_size'] = int(fields[3])
self['block_nblocks'] = int(fields[5])
self['c2_sem'] = self.get('c2_sem', self['c2_std'])
self['block_size'] = self.get('block_size', 0)
self['block_nblocks'] = self.get('block_nblocks', 0)
def summary_output(self):
summary = self.summary()
text = []
for row in summary:
met, key, sigma, avg, std, sem, ci = row
if isinstance(avg, str):
text.extend((f'{met:16s} {key:>13s} {sigma:>13s} {avg:>10s} {std:>8s} {sem:>8s} {ci:>14s}', sep))
else:
text.append(f"{met:16s} {key:>13s} {sigma:13.2f} {avg:10.2f} {std:8.2f} {sem:8.2f} {ci:>14s}")
if self.get('block_size'):
text.append(
f"Block diagnostic: {self['block_nblocks']} nonoverlapping blocks "
f"of {self['block_size']} frames"
)
return '\n'.join(text) + '\n\n'
def summary(self):
""" Formatted summary of C2 Entropy results """
return [
[
'Energy Method',
'Entropy',
'σ(Int. Energy)',
'C2 Value',
'Block SD',
'Block SEM',
'Block P2.5-P97.5'
],
[self.method.upper(), 'C2', float(self['sigma']), float(self['c2data']), float(self['c2_std']),
float(self.get('c2_sem', self['c2_std'])), f"{self['c2_ci'][0]:.2f}-{self['c2_ci'][1]:.2f}",]
]
class QHout(dict):
""" Quasi-harmonic output file class. QH output files are strange so we won't
derive from AmberOutput
"""
def __init__(self, filename=None, temp=298.15, **kwargs):
super(QHout, self).__init__(**kwargs)
self.filename = filename
self.temperature = temp
self.stability = False
self._read()
def summary_output(self):
text = [' TRANSLATIONAL ROTATIONAL VIBRATIONAL TOTAL',
'Complex %13.4f %15.4f %16.4f %15.4f' % (self['complex']['TRANSLATIONAL'],
self['complex']['ROTATIONAL'],
self['complex']['VIBRATIONAL'],
self['complex']['TOTAL'],)]
if not self.stability:
text.extend(['Receptor %13.4f %15.4f %16.4f %15.4f' % (self['receptor']['TRANSLATIONAL'],
self['receptor']['ROTATIONAL'],
self['receptor']['VIBRATIONAL'],
self['receptor']['TOTAL']),
'Ligand %13.4f %15.4f %16.4f %15.4f' % (self['ligand']['TRANSLATIONAL'],
self['ligand']['ROTATIONAL'],
self['ligand']['VIBRATIONAL'],
self['ligand']['TOTAL']),
'',
'Delta %13.4f %15.4f %16.4f %15.4f' % (self['delta']['TRANSLATIONAL'],
self['delta']['ROTATIONAL'],
self['delta']['VIBRATIONAL'],
self['delta']['TOTAL'])])
return '\n'.join(text) + '\n\n'
def summary(self):
""" Formatted summary of quasi-harmonic results """
summry_list = [['', 'TRANSLATIONAL', 'ROTATIONAL', 'VIBRATIONAL', 'TOTAL'],
['Complex', self['complex']['TRANSLATIONAL'], self['complex']['ROTATIONAL'],
self['complex']['VIBRATIONAL'], self['complex']['TOTAL']]]
if not self.stability:
summry_list.extend([['Receptor',
self['receptor']['TRANSLATIONAL'],
self['receptor']['ROTATIONAL'],
self['receptor']['VIBRATIONAL'],
self['receptor']['TOTAL']],
['Ligand',
self['ligand']['TRANSLATIONAL'],
self['ligand']['ROTATIONAL'],
self['ligand']['VIBRATIONAL'],
self['ligand']['TOTAL']],
['-TΔS',
self['delta']['TRANSLATIONAL'],
self['delta']['ROTATIONAL'],
self['delta']['VIBRATIONAL'],
self['delta']['TOTAL']]])
return summry_list
def _read(self):
""" Parses the output files and fills the data arrays """
with open(self.filename, 'r') as output:
rawline = output.readline()
self['complex'] = {}
self['receptor'] = {'TOTAL': 0}
self['ligand'] = {'TOTAL': 0}
self['delta'] = {}
comdone = False # if we've done the complex yet (filled in self.com)
recdone = False # if we've done the receptor yet (filled in self.rec)
# Try to fill in all found entropy values. If we can only find 1 set,
# we're doing stability calculations
while rawline:
if rawline[:6] == " Total":
if not comdone:
self['complex']['TOTAL'] = (float(rawline.split()[3]) * self.temperature / 1000 * -1)
self['complex']['TRANSLATIONAL'] = (
float(output.readline().split()[3]) * self.temperature / 1000 * -1)
self['complex']['ROTATIONAL'] = (
float(output.readline().split()[3]) * self.temperature / 1000 * -1)
self['complex']['VIBRATIONAL'] = (
float(output.readline().split()[3]) * self.temperature / 1000 * -1)
comdone = True
elif not recdone:
self['receptor']['TOTAL'] = (float(rawline.split()[3]) * self.temperature / 1000 * -1)
self['receptor']['TRANSLATIONAL'] = (
float(output.readline().split()[3]) * self.temperature / 1000 * -1)
self['receptor']['ROTATIONAL'] = (
float(output.readline().split()[3]) * self.temperature / 1000 * -1)
self['receptor']['VIBRATIONAL'] = (
float(output.readline().split()[3]) * self.temperature / 1000 * -1)
recdone = True
else:
self['ligand']['TOTAL'] = (float(rawline.split()[3]) * self.temperature / 1000 * -1)
self['ligand']['TRANSLATIONAL'] = (
float(output.readline().split()[3]) * self.temperature / 1000 * -1)
self['ligand']['ROTATIONAL'] = (
float(output.readline().split()[3]) * self.temperature / 1000 * -1)
self['ligand']['VIBRATIONAL'] = (
float(output.readline().split()[3]) * self.temperature / 1000 * -1)
break
rawline = output.readline()
# end while rawline
self.stability = not recdone
# fill the delta if not stability
if not self.stability:
self['delta']['TOTAL'] = self['complex']['TOTAL'] - self['receptor']['TOTAL'] - self['ligand']['TOTAL']
self['delta']['TRANSLATIONAL'] = (self['complex']['TRANSLATIONAL'] - self['receptor']['TRANSLATIONAL'] -
self['ligand']['TRANSLATIONAL'])
self['delta']['ROTATIONAL'] = (self['complex']['ROTATIONAL'] - self['receptor']['ROTATIONAL'] -
self['ligand']['ROTATIONAL'])
self['delta']['VIBRATIONAL'] = (self['complex']['VIBRATIONAL'] - self['receptor']['VIBRATIONAL'] -
self['ligand']['VIBRATIONAL'])
class DeltaDeltaQH(dict):
def __init__(self, mut, norm, **kwargs):
super(DeltaDeltaQH, self).__init__(**kwargs)
self.mut = mut
self.norm = norm
self._delta()
def _delta(self):
for key, mut_values in self.mut.items():
norm_values = self.norm.get(key, {})
self[key] = {
term: value - norm_values[term]
for term, value in mut_values.items()
if term in norm_values
}
def summary_output(self):
text = [' TRANSLATIONAL ROTATIONAL VIBRATIONAL TOTAL',
'Complex %13.4f %15.4f %16.4f %15.4f' % (self['complex']['TRANSLATIONAL'],
self['complex']['ROTATIONAL'],
self['complex']['VIBRATIONAL'],
self['complex']['TOTAL'],)]
if not self.norm.stability:
text.extend(['Receptor %13.4f %15.4f %16.4f %15.4f' % (self['receptor']['TRANSLATIONAL'],
self['receptor']['ROTATIONAL'],
self['receptor']['VIBRATIONAL'],
self['receptor']['TOTAL']),
'Ligand %13.4f %15.4f %16.4f %15.4f' % (self['ligand']['TRANSLATIONAL'],
self['ligand']['ROTATIONAL'],
self['ligand']['VIBRATIONAL'],
self['ligand']['TOTAL']),
'',
'Delta %13.4f %15.4f %16.4f %15.4f' % (self['delta']['TRANSLATIONAL'],
self['delta']['ROTATIONAL'],
self['delta']['VIBRATIONAL'],
self['delta']['TOTAL'])])
return '\n'.join(text) + '\n\n'
def summary(self):
""" Formatted summary of quasi-harmonic results """
summry_list = [['', 'TRANSLATIONAL', 'ROTATIONAL', 'VIBRATIONAL', 'TOTAL'],
['Complex', self['complex']['TRANSLATIONAL'], self['complex']['ROTATIONAL'],
self['complex']['VIBRATIONAL'], self['complex']['TOTAL']]]
if not self.norm.stability:
summry_list.extend([['Receptor',
self['receptor']['TRANSLATIONAL'],
self['receptor']['ROTATIONAL'],
self['receptor']['VIBRATIONAL'],
self['receptor']['TOTAL']],
['Ligand',
self['ligand']['TRANSLATIONAL'],
self['ligand']['ROTATIONAL'],
self['ligand']['VIBRATIONAL'],
self['ligand']['TOTAL']],
['-TΔS',
self['delta']['TRANSLATIONAL'],
self['delta']['ROTATIONAL'],
self['delta']['VIBRATIONAL'],
self['delta']['TOTAL']]])
return summry_list
class NMODEout(AmberOutput):
""" Normal mode entropy approximation output class """
print_levels = {'TRANSLATIONAL': 1, 'ROTATIONAL': 1, 'VIBRATIONAL': 1, 'TOTAL': 1}
def __init__(self, mol: str, INPUT, chamber=False, **kwargs):
super(NMODEout, self).__init__(mol, INPUT, chamber, **kwargs)
# Ordered list of keys in the data dictionary
self.data_keys = ['TRANSLATIONAL', 'ROTATIONAL', 'VIBRATIONAL']
# Other aspects of AmberOutputs, which are just blank arrays
self.composite_keys = ['TOTAL']
def _get_energies(self, outfile):
""" Parses the energy terms from the output file. This will parse 1 line
at a time in order to minimize the memory requirements (we should only
have to store a single line at a time in addition to the arrays of
data)
"""
while rawline := outfile.readline():
if "|---- Entropy not Calculated---|" in rawline:
self['TOTAL'][self.frame_idx] = np.nan
self['TRANSLATIONAL'][self.frame_idx] = np.nan
self['ROTATIONAL'][self.frame_idx] = np.nan
self['VIBRATIONAL'][self.frame_idx] = np.nan
self.frame_idx += 1
if rawline[:6] == 'Total:':
self['TOTAL'][self.frame_idx] = float(rawline.split()[3]) * self.temperature / 1000 * -1
self['TRANSLATIONAL'][self.frame_idx] = (float(outfile.readline().split()[3]) * self.temperature /
1000 * -1)
self['ROTATIONAL'][self.frame_idx] = (float(outfile.readline().split()[3]) * self.temperature / 1000
* -1)
self['VIBRATIONAL'][self.frame_idx] = (float(outfile.readline().split()[3]) * self.temperature / 1000
* -1)
self.frame_idx += 1
def _fill_nmode_values(self):
"""Leave unconverged NMODE frames as NaN; do not impute values.
Averages and uncertainties omit non-finite frames. If every frame
failed minimization, NMODE is disabled for the run.
"""
nonfinite_frames = np.zeros(self.numframes, dtype=bool)
for term in self.data_keys + ['TOTAL']:
nonfinite_frames |= ~np.isfinite(self[term])
unconverged = int(nonfinite_frames.sum())
if unconverged == 0:
return
if np.isnan(self['TOTAL']).all():
logging.warning(
f'{self.mol.capitalize()}: {unconverged} of {self.numframes} NMODE frames did not satisfy the '
'minimized-energy-gradient convergence criterion. No converged values are available. '
'Increase drms or maxcyc in the NMODE settings.')
self.no_nmode_convergence = True
return
logging.warning(
f'{self.mol.capitalize()}: {unconverged} of {self.numframes} NMODE frames did not satisfy the '
'minimized-energy-gradient convergence criterion. Those frames remain NaN and are omitted from '
'NMODE averages and uncertainties. Increase drms or maxcyc if more frames should converge.')
def conv_float(word):
if '*' in word:
GMXMMPBSA_ERROR('Some energy terms are undefined. Please, check the input structure and trajectory. Check this '
'section the docs for more info '
'https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/Q%26A/calculations/#possible-solutions')
elif 'nan' in word.lower():
GMXMMPBSA_ERROR('Some energy terms are undefined. Please, check the input structure and trajectory. Check this '
'section the docs for more info '
'https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/Q%26A/calculations/#possible-solutions')
else:
return float(word)
class GBout(AmberOutput):
""" Amber output class for normal generalized Born simulations """
print_levels = {'BOND': 2, 'ANGLE': 2, 'DIHED': 2, 'VDWAALS': 1, 'EEL': 1, '1-4 VDW': 2, '1-4 EEL': 2, 'EGB': 1,
'ESURF': 1}
# Ordered list of keys in the data dictionary
def __init__(self, mol, INPUT, chamber=False, **kwargs):
AmberOutput.__init__(self, mol, INPUT, chamber, **kwargs)
self.data_keys.extend(['EGB', 'ESURF'])
def _get_energies(self, outfile):
""" Parses the mdout files for the GB potential terms """
while rawline := outfile.readline():
if rawline[:5] == ' BOND':
words = rawline.split()
self['BOND'][self.frame_idx] = conv_float(words[2])
self['ANGLE'][self.frame_idx] = conv_float(words[5])
self['DIHED'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
if self.chamber:
self['UB'][self.frame_idx] = conv_float(words[2])
self['IMP'][self.frame_idx] = conv_float(words[5])
self['CMAP'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
self['VDWAALS'][self.frame_idx] = conv_float(words[2])
self['EEL'][self.frame_idx] = conv_float(words[5])
self['EGB'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
self['1-4 VDW'][self.frame_idx] = conv_float(words[3])
self['1-4 EEL'][self.frame_idx] = conv_float(words[7])
self.frame_idx += 1
def _extra_reading(self, fileno):
# Load the ESURF data from the cpptraj output
fname = '%s.%d' % (self.basename, fileno)
fname = fname.replace('gb.mdout', 'gb_surf.dat')
surf_data = _get_cpptraj_surf(fname)
for sd in surf_data:
self['ESURF'][self.extraframe_idx] = sd * self.INPUT['gb']['surften'] + self.INPUT['gb']['surfoff']
self.extraframe_idx += 1
class GBNSR6out(AmberOutput):
""" Amber output class for normal generalized Born simulations """
print_levels = {'BOND': 2, 'ANGLE': 2, 'DIHED': 2, 'VDWAALS': 1, 'EEL': 1, '1-4 VDW': 2, '1-4 EEL': 2, 'EGB': 1,
'ESURF': 1}
# Ordered list of keys in the data dictionary
def __init__(self, mol, INPUT, chamber=False, **kwargs):
AmberOutput.__init__(self, mol, INPUT, chamber, **kwargs)
# As the MM terms will be updated, in order to maintain order, we need to initialize these keys
self.data_keys.extend(['EGB', 'ESURF'])
def _get_energies(self, outfile):
""" Parses the mdout files for the GB potential terms """
while rawline := outfile.readline():
if rawline[:5] == ' BOND':
words = rawline.split()
self['BOND'][self.frame_idx] = conv_float(words[2])
self['ANGLE'][self.frame_idx] = conv_float(words[5])
self['DIHED'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
if self.chamber:
self['UB'][self.frame_idx] = conv_float(words[2])
self['IMP'][self.frame_idx] = conv_float(words[5])
self['CMAP'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
self['VDWAALS'][self.frame_idx] = conv_float(words[2])
self['EEL'][self.frame_idx] = conv_float(words[5])
self['EGB'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
self['1-4 VDW'][self.frame_idx] = conv_float(words[3])
self['1-4 EEL'][self.frame_idx] = conv_float(words[7])
words = outfile.readline().split()
self['ESURF'][self.frame_idx] = conv_float(words[2])
self.frame_idx += 1
class MMout(AmberOutput):
""" Amber output class for normal MM simulations """
print_levels = {'BOND': 2, 'ANGLE': 2, 'DIHED': 2, 'VDWAALS': 1, 'EEL': 1, '1-4 VDW': 2, '1-4 EEL': 2}
# Ordered list of keys in the data dictionary
def __init__(self, mol, INPUT, chamber=False, **kwargs):
AmberOutput.__init__(self, mol, INPUT, chamber, **kwargs)
def _get_energies(self, outfile):
""" Parses the mdout files for the GB potential terms """
while rawline := outfile.readline():
if rawline[:5] == ' BOND':
words = rawline.split()
self['BOND'][self.frame_idx] = conv_float(words[2])
self['ANGLE'][self.frame_idx] = conv_float(words[5])
self['DIHED'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
if self.chamber:
self['UB'][self.frame_idx] = conv_float(words[2])
self['IMP'][self.frame_idx] = conv_float(words[5])
self['CMAP'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
self['VDWAALS'][self.frame_idx] = conv_float(words[2])
self['EEL'][self.frame_idx] = conv_float(words[5])
words = outfile.readline().split()
self['1-4 VDW'][self.frame_idx] = conv_float(words[3])
self['1-4 EEL'][self.frame_idx] = conv_float(words[7])
class PBout(AmberOutput):
# What the value of verbosity must be to print out this data
print_levels = {'BOND': 2, 'ANGLE': 2, 'DIHED': 2, 'VDWAALS': 1, 'EEL': 1,
'1-4 VDW': 2, '1-4 EEL': 2, 'EPB': 1, 'ENPOLAR': 1, 'EDISPER': 1}
def __init__(self, mol, INPUT, chamber=False, **kwargs):
AmberOutput.__init__(self, mol, INPUT, chamber, **kwargs)
# EPB is still parsed when eneopt=1/P3M (incl. Amber-forced NLPB): Amber reports EPB≈0 and
# folds RF+Coulomb into EEL. GGAS/GSOLV labels are then not a separable partition; TOTAL is.
self.data_keys.extend(['EPB', 'ENPOLAR', 'EDISPER'])
def _get_energies(self, outfile):
""" Parses the energy values from the output files """
while rawline := outfile.readline():
if rawline[:5] == ' BOND':
words = rawline.split()
self['BOND'][self.frame_idx] = conv_float(words[2])
self['ANGLE'][self.frame_idx] = conv_float(words[5])
self['DIHED'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
if self.chamber:
self['UB'][self.frame_idx] = conv_float(words[2])
self['IMP'][self.frame_idx] = conv_float(words[5])
self['CMAP'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
self['VDWAALS'][self.frame_idx] = conv_float(words[2])
self['EEL'][self.frame_idx] = conv_float(words[5])
self['EPB'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
self['1-4 VDW'][self.frame_idx] = conv_float(words[3])
self['1-4 EEL'][self.frame_idx] = conv_float(words[7])
words = outfile.readline().split()
# electrostatic solvation free energy will not report in amber ouput
# when `ipb=0`. in such case, set it as 0 would be reasonable
if self.INPUT['pb']['ipb'] == 0:
self['ENPOLAR'][self.frame_idx] = 0
else:
self['ENPOLAR'][self.frame_idx] = conv_float(words[2])
if self.INPUT['pb']['inp'] == 2 and not self.apbs:
self['EDISPER'][self.frame_idx] = conv_float(words[5])
self.frame_idx += 1
class RISMout(AmberOutput):
# Which of those keys belong to the gas phase energy contributions
print_levels = {'BOND': 2, 'ANGLE': 2, 'DIHED': 2, 'VDWAALS': 1, 'EEL': 1,
'1-4 VDW': 2, '1-4 EEL': 2, 'ERISM': 1}
def __init__(self, mol, INPUT, chamber=False, solvtype=0, **kwargs):
AmberOutput.__init__(self, mol, INPUT, chamber)
self.solvtype = solvtype
self.data_keys.extend(['ERISM'])
def _get_energies(self, outfile):
""" Parses the RISM output file for energy terms """
# Getting the RISM solvation energies requires some decision-making.
# There are 2 possibilities (right now):
#
# 1. Standard free energy (solvtype==0)
# 2. GF free energy (solvtype==1)
# 3. PC+ GF free energy (solvtype==2)
while rawline := outfile.readline():
if re.match(r'(solute_epot|solutePotentialEnergy)', rawline):
words = rawline.split()
self['VDWAALS'][self.frame_idx] = conv_float(words[2])
self['EEL'][self.frame_idx] = conv_float(words[3])
self['BOND'][self.frame_idx] = conv_float(words[4])
self['ANGLE'][self.frame_idx] = conv_float(words[5])
self['DIHED'][self.frame_idx] = conv_float(words[6])
self['1-4 VDW'][self.frame_idx] = conv_float(words[7])
self['1-4 EEL'][self.frame_idx] = conv_float(words[8])
elif self.solvtype == 0 and re.match(r'(rism_exchem|rism_excessChemicalPotential)\s', rawline):
self['ERISM'][self.frame_idx] = conv_float(rawline.split()[1])
self.frame_idx += 1
elif self.solvtype == 1 and re.match(r'(rism_exchGF|rism_excessChemicalPotentialGF)\s', rawline):
self['ERISM'][self.frame_idx] = conv_float(rawline.split()[1])
self.frame_idx += 1
elif self.solvtype == 2 and re.match(r'(rism_exchPCPLUS|rism_excessChemicalPotentialPCPLUS)\s', rawline):
self['ERISM'][self.frame_idx] = conv_float(rawline.split()[1])
self.frame_idx += 1
class RISM_std_Out(RISMout):
""" No polar decomp RISM output file for standard free energy """
def __init__(self, mol, INPUT, chamber=False):
RISMout.__init__(self, mol, INPUT, chamber, 0)
class RISM_gf_Out(RISMout):
""" No polar decomp RISM output file for Gaussian Fluctuation free energy """
def __init__(self, mol, INPUT, chamber=False):
RISMout.__init__(self, mol, INPUT, chamber, 1)
class RISM_pcplus_Out(RISMout):
""" No polar decomp RISM output file for PC+ free energy """
def __init__(self, mol, INPUT, chamber=False):
RISMout.__init__(self, mol, INPUT, chamber, 2)
class PolarRISMout(RISMout):
# Which of those keys belong to the gas phase energy contributions
print_levels = {'BOND': 2, 'ANGLE': 2, 'DIHED': 2, 'VDWAALS': 1, 'EEL': 1,
'1-4 VDW': 2, '1-4 EEL': 2, 'POLAR SOLV': 1, 'APOLAR SOLV': 1}
def __init__(self, mol, INPUT, chamber=False, solvtype=0, **kwargs):
AmberOutput.__init__(self, mol, INPUT, chamber)
self.solvtype = solvtype
self.data_keys.extend(['POLAR SOLV', 'APOLAR SOLV'])
def _get_energies(self, outfile):
""" Parses the RISM output file for energy terms """
# Getting the RISM solvation energies requires some decision-making.
# There are 2 possibilities (right now):
#
# 1. Standard free energy (solvtype==0)
# 2. GF free energy (solvtype==1)
while rawline := outfile.readline():
if re.match(r'(solute_epot|solutePotentialEnergy)', rawline):
words = rawline.split()
self['VDWAALS'][self.frame_idx] = conv_float(words[2])
self['EEL'][self.frame_idx] = conv_float(words[3])
self['BOND'][self.frame_idx] = conv_float(words[4])
self['ANGLE'][self.frame_idx] = conv_float(words[5])
self['DIHED'][self.frame_idx] = conv_float(words[6])
self['1-4 VDW'][self.frame_idx] = conv_float(words[8])
self['1-4 EEL'][self.frame_idx] = conv_float(words[8])
elif self.solvtype == 0 and re.match(
r'(rism_polar|rism_polarExcessChemicalPotential)\s', rawline):
self['POLAR SOLV'][self.frame_idx] = conv_float(rawline.split()[1])
elif self.solvtype == 0 and re.match(
r'(rism_apolar|rism_apolarExcessChemicalPotential)\s', rawline):
self['APOLAR SOLV'][self.frame_idx] = conv_float(rawline.split()[1])
self.frame_idx += 1
elif self.solvtype == 1 and re.match(
r'(rism_polGF|rism_polarExcessChemicalPotentialGF)\s', rawline):
self['POLAR SOLV'][self.frame_idx] = conv_float(rawline.split()[1])
elif self.solvtype == 1 and re.match(
r'(rism_apolGF|rism_apolarExcessChemicalPotentialGF)\s', rawline):
self['APOLAR SOLV'][self.frame_idx] = conv_float(rawline.split()[1])
self.frame_idx += 1
elif self.solvtype == 2 and re.match(
r'(rism_polPCPLUS|rism_polarExcessChemicalPotentialPCPLUS)\s', rawline):
self['POLAR SOLV'][self.frame_idx] = conv_float(rawline.split()[1])
elif self.solvtype == 2 and re.match(
r'(rism_apolPCPLUS|rism_apolarExcessChemicalPotentialPCPLUS)\s', rawline):
self['APOLAR SOLV'][self.frame_idx] = conv_float(rawline.split()[1])
self.frame_idx += 1
class PolarRISM_std_Out(PolarRISMout):
""" Polar decomp RISM output file for standard free energy """
def __init__(self, mol, INPUT, chamber=False):
PolarRISMout.__init__(self, mol, INPUT, chamber, 0)
class PolarRISM_gf_Out(PolarRISMout):
""" Polar decomp RISM output file for Gaussian Fluctuation free energy """
def __init__(self, mol, INPUT, chamber=False):
PolarRISMout.__init__(self, mol, INPUT, chamber, 1)
class PolarRISM_pcplus_Out(PolarRISMout):
""" Polar decomp RISM output file for PC+ free energy """
def __init__(self, mol, INPUT, chamber=False):
PolarRISMout.__init__(self, mol, INPUT, chamber, 2)
class QMMMout(GBout):
""" Class for QM/MM GBSA output files """
# What the value of verbosity must be to print out this data
print_levels = {'BOND': 2, 'ANGLE': 2, 'DIHED': 2, 'VDWAALS': 1, 'EEL': 1,
'1-4 VDW': 2, '1-4 EEL': 2, 'EGB': 1, 'ESURF': 1, 'ESCF': 1}
def __init__(self, mol, INPUT, chamber=False, **kwargs):
GBout.__init__(self, mol, INPUT, chamber, **kwargs)
self.data_keys.extend(['ESCF'])
def _get_energies(self, outfile):
""" Parses the energies from a QM/MM output file. NOTE, however, that a
QMMMout *could* just be a GBout with ESCF==0 if the QM region lies
entirely outside this system
"""
while rawline := outfile.readline():
if rawline[:5] == ' BOND':
words = rawline.split()
self['BOND'][self.frame_idx] = conv_float(words[2])
self['ANGLE'][self.frame_idx] = conv_float(words[5])
self['DIHED'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
if self.chamber:
self['UB'][self.frame_idx] = conv_float(words[2])
self['IMP'][self.frame_idx] = conv_float(words[5])
self['CMAP'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
self['VDWAALS'][self.frame_idx] = conv_float(words[2])
self['EEL'][self.frame_idx] = conv_float(words[5])
self['EGB'][self.frame_idx] = conv_float(words[8])
words = outfile.readline().split()
self['1-4 VDW'][self.frame_idx] = conv_float(words[3])
self['1-4 EEL'][self.frame_idx] = conv_float(words[7])
words = outfile.readline().split()
# This is where ESCF will be. Since ESCF can differ based on which
# qmtheory was chosen, we just check to see if it's != ESURF:
if words[0] == 'minimization':
continue
elif words[0].endswith('='):
self['ESCF'][self.frame_idx] = conv_float(words[1])
else:
self['ESCF'][self.frame_idx] = conv_float(words[2])
self.frame_idx += 1
class BindingStatistics(dict):
""" Base class for compiling the binding statistics """
st_null = ['BOND', 'ANGLE', 'DIHED', '1-4 VDW', '1-4 EEL']
def __init__(self, com, rec, lig, chamber=False, traj_protocol='STP', **kwargs):
super(BindingStatistics, self).__init__(**kwargs)
self.com = com
self.rec = rec
self.lig = lig
self.mol = 'delta'
self.numframes = self.com.numframes
self.INPUT = self.com.INPUT
self.chamber = chamber
self.traj_protocol = traj_protocol
self.inconsistent = False
self.missing_terms = False
self.data_keys = self.com.data_keys
self.composite_keys = []
try:
self._delta()
self.missing_terms = False
except LengthError:
self._delta2()
self.missing_terms = True
def _delta(self):
"""
Calculates the delta statistics. Should check for any consistencies that
would cause verbosity levels to change, and it should change them
accordingly in the child classes
"""
# First thing we do is check to make sure that all of the terms that
# should *not* be printed actually cancel out (i.e. bonded terms)
if not isinstance(self.com, NMODEout) and self.traj_protocol == 'STP':
TINY = 0.005
for key in self.st_null:
diff = self.com[key] - self.rec[key] - self.lig[key]
if diff.abs_gt(TINY):
self.inconsistent = True
logging.warning(f"{key} component is reported as inconsistent. Please, check the output file "
f"for more details")
break
for key in self.com.data_keys:
if self.traj_protocol == 'STP':
temp = self.com[key].corr_sub(self.rec[key])
self[key] = temp.corr_sub(self.lig[key])
else:
self[key] = self.com[key] - self.rec[key] - self.lig[key]
for key in self.com.composite_keys:
self[key] = EnergyVector(self.numframes)
self.composite_keys.append(key)
for key in self.com.data_keys:
if self.traj_protocol == 'STP' and key in self.st_null:
continue
for component in data_key_owner[key]:
self[component] = self[key] + self[component]
def _print_vectors(self, csvwriter):
""" Output all of the energy terms including the differences if we're
doing a single trajectory simulation and there are no missing terms
"""
term_text = 'Entropy' if isinstance(self.com, NMODEout) else 'Energy'
csvwriter.writerow([f'Complex {term_text} Terms'])
self.com._print_vectors(csvwriter)
csvwriter.writerow([])
csvwriter.writerow([f'Receptor {term_text} Terms'])
self.rec._print_vectors(csvwriter)
csvwriter.writerow([])
csvwriter.writerow([f'Ligand {term_text} Terms'])
self.lig._print_vectors(csvwriter)
csvwriter.writerow([])
csvwriter.writerow([f'Delta {term_text} Terms'])
print_keys = list(self.data_keys)
# Add on the composite keys
print_keys += self.composite_keys
# write the header
csvwriter.writerow(['Frame #'] + print_keys)
# write out each frame
c = self.com.INPUT['nmode']['nmstartframe'] if isinstance(self.com, NMODEout) else self.com.INPUT['general'][
'startframe']
for i in range(self.numframes):
csvwriter.writerow([c] + [round(self[key][i], 2) for key in print_keys])
c += self.com.INPUT['nmode']['nminterval'] if isinstance(self.com, NMODEout) else self.com.INPUT['general'][
'interval']
csvwriter.writerow([])
def report_inconsistency(self, output_format: str = 'ascii'):
_output_format = 0 if output_format == 'ascii' else 1
text = []
if _output_format:
text.append(['WARNING: INCONSISTENCIES EXIST WITHIN INTERNAL POTENTIAL TERMS AND\n'
'THE VALIDITY OF THESE RESULTS ARE HIGHLY QUESTIONABLE!\n'
'\n'
'Some absolute differences in the internal potential terms are greater than 0.005.\n'
'This should not happen when using Single Trajectory Protocol!\n'
'\n'
'You can generate a detailed *.csv file with all the terms and differences as follows:\n'
'gmx_MMPBSA --rewrite-output -eo energy_terms_differences.csv'])
else:
text.append('WARNING: INCONSISTENCIES EXIST WITHIN INTERNAL POTENTIAL TERMS AND\n'
'THE VALIDITY OF THESE RESULTS ARE HIGHLY QUESTIONABLE!\n'
'\n'
'Some absolute differences in the internal potential terms are greater than 0.005.\n'
'This should not happen when using Single Trajectory Protocol!\n'
'\n'
'You can generate a detailed *.csv file with all the terms and differences as follows:\n'
'gmx_MMPBSA --rewrite-output -eo energy_terms_differences.csv')
return text if _output_format else '\n'.join(text) + '\n'
def summary_output(self, output_format: str = 'ascii'):
_output_format = 0 if output_format == 'ascii' else 1
summary = self.summary()
text = []
if _output_format:
text.append(['Delta (Complex - Receptor - Ligand):'])
else:
text.append('Delta (Complex - Receptor - Ligand):')
for c, row in enumerate(summary, start=1):
# Skip the composite terms, since we print those at the end
if _output_format:
text.append(row)
else:
key, avg, stdev, std, semp, sem, block_sd, block_sem = row
if key in ['GGAS', 'TOTAL']:
text.append('')
if isinstance(avg, str):
text.extend(
(
f'{key:16s} {avg:>13s} {stdev:>13s} {std:>10s} {semp:>12s} {sem:>10s} '
f'{block_sd:>10s} {block_sem:>10s}',
sep,
)
)
else:
text.append(f'{f"Δ{key}":16s} {avg:13.2f} {stdev:13.2f} {std:10.2f} {semp:12.2f} {sem:10.2f} '
f'{block_sd:10.2f} {block_sem:10.2f}')
return text if _output_format else '\n'.join(text) + '\n'
def summary(self):
""" Returns a string printing the summary of the binding statistics """
if isinstance(self.com, NMODEout):
col_name = '%-16s' % 'Entropy Term'
else:
col_name = '%-16s' % 'Energy Component'
summary_list = [
[col_name] + ['Average', 'SD(Prop.)', 'SD', 'SEM(Prop.)', 'SEM', 'Block SD', 'Block SEM']
]
for key in self.data_keys:
# Skip the composite terms, since we print those at the end
if key in self.composite_keys:
continue
# Now print out the stats
stdev = float(self[key].stdev())
avg = float(self[key].mean())
std = float(self[key].std())
semp = float(self[key].semp())
sem = float(self[key].sem())
_, _, block_sd, block_sem = block_statistics(self[key])
summary_list.append([key, avg, stdev, std, semp, sem, block_sd, block_sem])
for key in self.composite_keys:
# Now print out the composite terms
stdev = float(self[key].stdev())
avg = float(self[key].mean())
std = float(self[key].std())
semp = float(self[key].semp())
sem = float(self[key].sem())
_, _, block_sd, block_sem = block_statistics(self[key])
summary_list.append([key, avg, stdev, std, semp, sem, block_sd, block_sem])
return summary_list
class DeltaDeltaStatistics(dict):
""" Base class for compiling the binding statistics. Include """
st_null = ['BOND', 'ANGLE', 'DIHED', '1-4 VDW', '1-4 EEL']
def __init__(self, mut, norm, **kwargs):
super(DeltaDeltaStatistics, self).__init__(**kwargs)
self.mut = mut
self.norm = norm
self.numframes = self.norm.numframes
self.mol = self.norm.mol
self.data_keys = self.norm.data_keys
self.composite_keys = ['GGAS', 'GSOLV', 'TOTAL']
self.term_text = 'Entropy' if 'ROTATIONAL' in self.data_keys else 'Energy'
self._delta()
def _delta(self):
"""
Calculates the delta statistics. Should check for any consistencies that
would cause verbosity levels to change, and it should change them
accordingly in the child classes
"""
for key in self.norm:
if key in self.composite_keys:
continue
self[key] = self.mut[key].corr_sub(self.norm[key])
for key in self.composite_keys:
self[key] = EnergyVector(self.numframes)
# self.composite_keys.append(key)
for key in self.data_keys:
for component in data_key_owner[key]:
self[component] = self[key] + self[component]
def _print_vectors(self, csvwriter):
""" Output all of the energy terms including the differences if we're
doing a single trajectory simulation and there are no missing terms
"""
csvwriter.writerow([f'Delta Delta {self.term_text} Terms (Mutant - Normal)'])
# write the header
csvwriter.writerow(['Frame #'] + list(self.keys()))
# write out each frame
c = self.norm.INPUT['nmode']['nmstartframe'] if isinstance(self.norm, NMODEout) else \
self.norm.INPUT['general']['startframe']
for i in range(self.numframes):
csvwriter.writerow([c] + [round(self[key][i], 2) for key in self])
c += self.norm.INPUT['nmode']['nminterval'] if isinstance(self.norm, NMODEout) else \
self.norm.INPUT['general']['interval']
csvwriter.writerow([])
def summary_output(self):
summary = self.summary()
text = ['Delta Delta (Mutant - Normal):']
for c, row in enumerate(summary, start=1):
# Skip the composite terms, since we print those at the end
key, avg, stdev, std, semp, sem, block_sd, block_sem = row
if key in ['GGAS', 'TOTAL']:
text.append('')
if isinstance(avg, str):
text.extend(
(
f'{key:16s} {avg:>13s} {stdev:>13s} {std:>10s} {semp:>12s} {sem:>10s} '
f'{block_sd:>10s} {block_sem:>10s}',
sep,
)
)
else:
text.append(f'{f"ΔΔ{key}":16s} {avg:13.2f} {stdev:13.2f} {std:10.2f} {semp:12.2f} {sem:10.2f} '
f'{block_sd:10.2f} {block_sem:10.2f}')
return '\n'.join(text) + '\n'
def summary(self):
""" Returns a string printing the summary of the binding statistics """
summary_list = [
[f'{self.term_text} Component'] + ['Average', 'SD(Prop.)', 'SD', 'SEM(Prop.)', 'SEM', 'Block SD', 'Block SEM']
]
for key in self.norm.data_keys:
# Now print out the stats
stdev = float(self[key].stdev())
avg = float(self[key].mean())
std = float(self[key].std())
sem = float(self[key].sem())
semp = float(self[key].semp())
_, _, block_sd, block_sem = block_statistics(self[key])
summary_list.append([key, avg, stdev, std, semp, sem, block_sd, block_sem])
for key in self.composite_keys:
# Now print out the composite terms
stdev = float(self[key].stdev())
avg = float(self[key].mean())
std = float(self[key].std())
sem = float(self[key].sem())
semp = float(self[key].semp())
_, _, block_sd, block_sem = block_statistics(self[key])
summary_list.append([key, avg, stdev, std, semp, sem, block_sd, block_sem])
return summary_list
class DeltaIEC2Statistic(dict):
def __init__(self, mut, norm, **kwargs):
super(DeltaIEC2Statistic, self).__init__(**kwargs)
self.mut = mut
self.norm = norm
self.meth = 'C2' if 'c2data' in self.norm else 'IE'
self._delta()
def _delta(self):
"""
Calculates the delta statistics. Should check for any consistencies that
would cause verbosity levels to change, and it should change them
accordingly in the child classes
"""
for key in self.norm:
if key == 'sigma':
self[key] = (self.mut[key] + self.norm[key]) / 2
elif key == 'ieframes':
self[key] = self.norm[key]
elif key in ('ie_value', 'tail_mean'):
self[key] = self.mut[key] - self.norm[key]
elif key in ('tail_std', 'block_std', 'block_sem'):
self[key] = get_std(self.mut[key], self.norm[key])
elif key in ('block_size', 'block_nblocks'):
self[key] = min(self.mut[key], self.norm[key])
elif key == 'c2data':
self[key] = self.mut[key] - self.norm[key]
elif key in ('c2_std', 'c2_sem'):
self[key] = get_std(self.mut[key], self.norm[key])
elif key == 'c2_ci':
continue
else:
self[key] = self.mut[key].corr_sub(self.norm[key])
def _print_vectors(self, csvwriter):
""" Output all of the energy terms including the differences if we're
doing a single trajectory simulation and there are no missing terms
"""
csvwriter.writerow(['Delta Delta Entropy Terms (Mutant - Normal)'])
# write the header
csvwriter.writerow(['Frame #'] + list(self.keys()))
# write out each frame
c = self.norm.INPUT['nmode']['nmstartframe'] if isinstance(self.norm, NMODEout) else\
self.norm.INPUT['general']['startframe']
for i in range(self.numframes):
csvwriter.writerow([c] + [round(self[key][i], 2) for key in self])
c += self.norm.INPUT['nmode']['nminterval'] if isinstance(self.norm, NMODEout) else\
self.norm.INPUT['general']['interval']
csvwriter.writerow([])
def summary_output(self):
summary = self.summary()
text = []
for row in summary:
key, sigma, avg, std, sem = row
if isinstance(avg, str):
text.extend((f'{key:15s} {sigma:>14s} {avg:>16s} {std:>14s} {sem:>14s}', sep))
else:
text.append(f"Δ{key:14s} {sigma:14.2f} {avg:16.2f} {std:14.2f} {sem:14.2f}")
return '\n'.join(text) + '\n\n'
def summary(self):
""" Returns a string printing the summary of the binding statistics """
if self.meth == 'C2':
return [
[
'Method',
'σ(Int. Energy)',
'C2 Value',
'Block SD',
'Block SEM'
],
['C2', float(self['sigma']), float(self['c2data']), float(self['c2_std']),
float(self.get('c2_sem', self['c2_std']))]
]
avg = float(self.get('ie_value', self['iedata'].mean()))
stdev = float(self.get('block_std', self['data'][-self['ieframes']:].stdev()))
return [
[
'Method',
'σ(Int. Energy)',
'Full IE',
'Block SD',
'Block SEM'
],
['IE', self['sigma'], avg, stdev,
float(self.get('block_sem', self['data'][-self['ieframes']:].sem()))]
]
class DecompOut(dict):
""" Class for decomposition output file to collect statistics and output them """
indicator = " PRINT DECOMP - TOTAL ENERGIES"
descriptions = {'TDC': 'Total Energy Decomposition:',
'SDC': 'Sidechain Energy Decomposition:',
'BDC': 'Backbone Energy Decomposition:'}
def __init__(self, mol: str, **kwargs):
super(DecompOut, self).__init__(**kwargs)
self.numframes = None
self.mut = None
self.mol = mol
self.decfile = None
self.num_terms = None
self.allowed_tokens = ('TDC',)
self.verbose = None
self.num_files = None
self.resl = None
self.basename = None
self.csvwriter = None
self.surften = None
self.frame_idx = 0
self.current_file = 0 # File counter
def set_frame_range(self, start=0, end=None, interval=1):
frames_updated = False
for term in self.allowed_tokens:
for res in self[term]:
for et in self[term][res]:
self[term][res][et] = self[term][res][et][start:end:interval]
if not frames_updated:
self.numframes = len(self[term][res][et])
frames_updated = True
self._fill_composite_terms()
def parse_from_file(self, basename, resl, INPUT, surften, num_files=1, numframes=1, mut=False):
self.basename = basename # base name of output files
self.resl = resl
self.mut = mut
self.numframes = numframes
self.num_files = num_files # how many MPI files we created
self.INPUT = INPUT
self.verbose = INPUT['decomp']['dec_verbose']
self.surften = surften # explicitly defined since is for GB and PB models
if self.verbose in [1, 3]:
self.allowed_tokens = 'TDC', 'SDC', 'BDC'
try:
self.num_terms = int(self._get_num_terms())
except TypeError:
raise OutputError('DecompOut: Not a decomp output file')
for token in self.allowed_tokens:
self[token] = {}
self._read()
self._fill_composite_terms()
def _get_num_terms(self):
""" Gets the number of terms in the output file """
with open('%s.%d' % (self.basename, 0), 'r') as decfile:
lines = decfile.readlines()
num_terms = 0
flag = False
for line in lines:
if line[:3] == 'TDC':
num_terms += 1
flag = True
elif flag:
break
# We've now gotten to the end of the Total Decomp Contribution,
# so we know how many terms we have
if not flag:
raise TypeError(f"{self.basename}.0 have 0 TDC starts")
return num_terms
def _read(self):
"""
Internal reading function. This should be called at the end of __init__.
It loops through all of the output files to populate the arrays
"""
for fileno in range(self.num_files):
with open('%s.%d' % (self.basename, fileno)) as output_file:
self._get_decomp_energies(output_file)
def _get_decomp_energies(self, outfile):
while line := outfile.readline():
if self.frame_idx == self.numframes:
self.frame_idx = 0
if line[:3] in self.allowed_tokens:
if self.mut and self.resl[int(line[4:10])].is_mutant():
resnum = self.resl[int(line[4:10])].mutant_string
else:
resnum = self.resl[int(line[4:10])].string
internal = float(line[11:20])
vdw = float(line[21:30])
eel = float(line[31:40])
pol = float(line[41:50])
sas = float(line[51:60]) * self.surften
if resnum not in self[line[:3]]:
self[line[:3]][resnum] = {}
for term in ['int', 'vdw', 'eel', 'pol', 'sas']:
self[line[:3]][resnum][term] = EnergyVector(self.numframes)
self[line[:3]][resnum]['int'][self.frame_idx] = internal
self[line[:3]][resnum]['vdw'][self.frame_idx] = vdw
self[line[:3]][resnum]['eel'][self.frame_idx] = eel
self[line[:3]][resnum]['pol'][self.frame_idx] = pol
self[line[:3]][resnum]['sas'][self.frame_idx] = sas
if line[:3] == self.allowed_tokens[-1] and resnum == list(self.resl.values())[-1].string:
self.frame_idx += 1
def _print_vectors(self, csvwriter):
tokens = {'TDC': 'Total Decomposition Contribution (TDC)',
'SDC': 'Sidechain Decomposition Contribution (SDC)',
'BDC': 'Backbone Decomposition Contribution (BDC)'}
for term in self.allowed_tokens:
csvwriter.writerow([tokens[term]])
csvwriter.writerow(['Frame #', 'Residue', 'Internal', 'van der Waals', 'Electrostatic', 'Polar Solvation',
'Non-Polar Solv.', 'TOTAL'])
c = self.INPUT['general']['startframe']
for i in range(self.numframes):
for res in self[term]:
csvwriter.writerow([c, res] + [round(self[term][res][key][i], 2) for key in self[term][res]])
c += self.INPUT['general']['interval']
def _fill_composite_terms(self):
for term in self:
for res in self[term]:
item = self[term][res][list(self[term][res].keys())[0]]
# pair decomp scheme
if isinstance(item, dict):
for res2 in self[term][res]:
tot = EnergyVector(self.numframes)
for e in self[term][res][res2]:
tot = tot + self[term][res][res2][e]
self[term][res][res2]['tot'] = tot
else:
tot = EnergyVector(self.numframes)
for e in self[term][res]:
tot = tot + self[term][res][e]
self[term][res]['tot'] = tot
def summary(self, output_format: str = 'ascii'):
""" Writes the summary in ASCII format to and open output_file """
_output_format = 0 if output_format == 'ascii' else 1
text = []
if _output_format:
text.append([f'{self.mol.capitalize()}:'])
else:
text.append(f'{self.mol.capitalize()}:')
for term in self:
if _output_format:
text.extend([[self.descriptions[term]],
['Residue', 'Internal', '', '', 'van der Waals', '',
'', 'Electrostatic', '', '', 'Polar Solvation', '',
'', 'Non-Polar Solv.', '', '', 'TOTAL', '', ''],
[''] + ['Avg.', 'SD(Prop.)', 'SD', 'SEM(Prop.)', 'SEM', 'Block SD', 'Block SEM'] * 6])
else:
text.extend([self.descriptions[term],
'Residue | Internal Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| van der Waals Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| Electrostatic Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| Polar Solvation Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| Non-Polar Solv. Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| TOTAL Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD',
'-------------------------------------------------------------------------------------'
'-------------------------------------------------------------'])
for res in self[term]:
int_avg, int_stdev, int_std, int_semp, int_sem, int_block_sd, int_block_sem = _vector_statistics(self[term][res]['int'])
vdw_avg, vdw_stdev, vdw_std, vdw_semp, vdw_sem, vdw_block_sd, vdw_block_sem = _vector_statistics(self[term][res]['vdw'])
eel_avg, eel_stdev, eel_std, eel_semp, eel_sem, eel_block_sd, eel_block_sem = _vector_statistics(self[term][res]['eel'])
pol_avg, pol_stdev, pol_std, pol_semp, pol_sem, pol_block_sd, pol_block_sem = _vector_statistics(self[term][res]['pol'])
sas_avg, sas_stdev, sas_std, sas_semp, sas_sem, sas_block_sd, sas_block_sem = _vector_statistics(self[term][res]['sas'])
tot_avg, tot_stdev, tot_std, tot_semp, tot_sem, tot_block_sd, tot_block_sem = _vector_statistics(self[term][res]['tot'])
if _output_format:
text.append([res,
int_avg, int_stdev, int_std, int_semp, int_sem, int_block_sd, int_block_sem,
vdw_avg, vdw_stdev, vdw_std, vdw_semp, vdw_sem, vdw_block_sd, vdw_block_sem,
eel_avg, eel_stdev, eel_std, eel_semp, eel_sem, eel_block_sd, eel_block_sem,
pol_avg, pol_stdev, pol_std, pol_semp, pol_sem, pol_block_sd, pol_block_sem,
sas_avg, sas_stdev, sas_std, sas_semp, sas_sem, sas_block_sd, sas_block_sem,
tot_avg, tot_stdev, tot_std, tot_semp, tot_sem, tot_block_sd, tot_block_sem])
else:
text.append(f"{res:14s} "
f"|{int_avg:9.3f} +/- {int_block_sem:6.3f} [{int_stdev:6.3f}/{int_std:6.3f}/{int_sem:6.3f}] / {int_block_sd:6.3f} "
f"|{vdw_avg:9.3f} +/- {vdw_block_sem:6.3f} [{vdw_stdev:6.3f}/{vdw_std:6.3f}/{vdw_sem:6.3f}] / {vdw_block_sd:6.3f} "
f"|{eel_avg:9.3f} +/- {eel_block_sem:6.3f} [{eel_stdev:6.3f}/{eel_std:6.3f}/{eel_sem:6.3f}] / {eel_block_sd:6.3f} "
f"|{pol_avg:9.3f} +/- {pol_block_sem:6.3f} [{pol_stdev:6.3f}/{pol_std:6.3f}/{pol_sem:6.3f}] / {pol_block_sd:6.3f} "
f"|{sas_avg:9.3f} +/- {sas_block_sem:6.3f} [{sas_stdev:6.3f}/{sas_std:6.3f}/{sas_sem:6.3f}] / {sas_block_sd:6.3f} "
f"|{tot_avg:9.3f} +/- {tot_block_sem:6.3f} [{tot_stdev:6.3f}/{tot_std:6.3f}/{tot_sem:6.3f}] / {tot_block_sd:6.3f}")
if _output_format:
text.append([])
else:
text.append('')
return text if _output_format else '\n'.join(text)
class PairDecompOut(DecompOut):
""" Same as DecompOut, but for Pairwise decomposition """
indicator = " PRINT PAIR DECOMP - TOTAL ENERGIES"
def set_frame_range(self, start=0, end=None, interval=1):
frames_updated = False
for term in self.allowed_tokens:
for res in self[term]:
for res2 in self[term][res]:
for et in self[term][res][res2]:
self[term][res][res2][et] = self[term][res][res2][et][start:end:interval]
if not frames_updated:
self.numframes = len(self[term][res][res2][et])
frames_updated = True
self._fill_composite_terms()
def _get_decomp_energies(self, outfile):
while line := outfile.readline():
if line[:3] in self.allowed_tokens:
if self.mut and self.resl[int(line[4:11])].is_mutant():
resnum = self.resl[int(line[4:11])].mutant_string
else:
resnum = self.resl[int(line[4:11])].string
if self.mut and self.resl[int(line[13:20])].is_mutant():
resnum2 = self.resl[int(line[13:20])].mutant_string
else:
resnum2 = self.resl[int(line[13:20])].string
internal = float(line[21:33])
vdw = float(line[34:46])
eel = float(line[47:59])
pol = float(line[60:72])
sas = float(line[73:85]) * self.surften
if resnum not in self[line[:3]]:
self[line[:3]][resnum] = {}
if resnum2 not in self[line[:3]][resnum]:
self[line[:3]][resnum][resnum2] = {}
for term in ['int', 'vdw', 'eel', 'pol', 'sas']:
self[line[:3]][resnum][resnum2][term] = EnergyVector(self.numframes)
self[line[:3]][resnum][resnum2]['int'][self.frame_idx] = internal
self[line[:3]][resnum][resnum2]['vdw'][self.frame_idx] = vdw
self[line[:3]][resnum][resnum2]['eel'][self.frame_idx] = eel
self[line[:3]][resnum][resnum2]['pol'][self.frame_idx] = pol
self[line[:3]][resnum][resnum2]['sas'][self.frame_idx] = sas
if line[:3] == self.allowed_tokens[-1] and resnum == list(self.resl.values())[-1].string == resnum2:
self.frame_idx += 1
def _print_vectors(self, csvwriter):
tokens = {'TDC': 'Total Decomposition Contribution (TDC)',
'SDC': 'Sidechain Decomposition Contribution (SDC)',
'BDC': 'Backbone Decomposition Contribution (BDC)'}
for term in self.allowed_tokens:
csvwriter.writerow([tokens[term]])
csvwriter.writerow(['Frame #', 'Resid 1', 'Resid 2', 'Internal', 'van der Waals', 'Electrostatic',
'Polar Solvation', 'Non-Polar Solv.', 'TOTAL'])
c = self.INPUT['general']['startframe']
for i in range(self.numframes):
for res in self[term]:
for res2 in self[term][res]:
csvwriter.writerow([c, res, res2] +
[round(self[term][res][res2][key][i], 2) for key in self[term][res][res2]])
c += self.INPUT['general']['interval']
def summary(self, output_format: str = 'ascii'):
""" Writes the summary in ASCII format to and open output_file """
_output_format = 0 if output_format == 'ascii' else 1
text = []
if _output_format:
text.extend([[f'{self.mol.capitalize()}:']])
else:
text.extend([f'{self.mol.capitalize()}:'])
for term in self:
if _output_format:
text.extend([[self.descriptions[term]],
['Resid 1', 'Resid 2', 'Internal', '', '', 'van der Waals', '', '', 'Electrostatic',
'', '', 'Polar Solvation', '', '', 'Non-Polar Solv.', '', '', 'TOTAL', '', ''],
[''] * 2 + ['Avg.', 'SD(Prop.)', 'SD', 'SEM(Prop.)', 'SEM', 'Block SD', 'Block SEM'] * 6])
else:
text.append(self.descriptions[term] + '\n' +
'Resid 1 | Resid 2 | Internal Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD | '
'van der Waals Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD | Electrostatic Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| Polar Solvation Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD | Non-Polar Solv. Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| TOTAL Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD\n' +
'-----------------------------------------------------------------------------'
'--------------------------------------------------------------------------------------')
for res in self[term]:
for res2 in self[term][res]:
int_avg, int_stdev, int_std, int_semp, int_sem, int_block_sd, int_block_sem = _vector_statistics(self[term][res][res2]['int'])
vdw_avg, vdw_stdev, vdw_std, vdw_semp, vdw_sem, vdw_block_sd, vdw_block_sem = _vector_statistics(self[term][res][res2]['vdw'])
eel_avg, eel_stdev, eel_std, eel_semp, eel_sem, eel_block_sd, eel_block_sem = _vector_statistics(self[term][res][res2]['eel'])
pol_avg, pol_stdev, pol_std, pol_semp, pol_sem, pol_block_sd, pol_block_sem = _vector_statistics(self[term][res][res2]['pol'])
sas_avg, sas_stdev, sas_std, sas_semp, sas_sem, sas_block_sd, sas_block_sem = _vector_statistics(self[term][res][res2]['sas'])
tot_avg, tot_stdev, tot_std, tot_semp, tot_sem, tot_block_sd, tot_block_sem = _vector_statistics(self[term][res][res2]['tot'])
if _output_format:
text.append([res, res2,
int_avg, int_stdev, int_std, int_semp, int_sem, int_block_sd, int_block_sem,
vdw_avg, vdw_stdev, vdw_std, vdw_semp, vdw_sem, vdw_block_sd, vdw_block_sem,
eel_avg, eel_stdev, eel_std, eel_semp, eel_sem, eel_block_sd, eel_block_sem,
pol_avg, pol_stdev, pol_std, pol_semp, pol_sem, pol_block_sd, pol_block_sem,
sas_avg, sas_stdev, sas_std, sas_semp, sas_sem, sas_block_sd, sas_block_sem,
tot_avg, tot_stdev, tot_std, tot_semp, tot_sem, tot_block_sd, tot_block_sem])
else:
text.append(f"{res:14s} | {res2:14s} "
f"|{int_avg:9.3f} +/- {int_block_sem:6.3f} [{int_stdev:6.3f}/{int_std:6.3f}/{int_sem:6.3f}] / {int_block_sd:6.3f} "
f"|{vdw_avg:9.3f} +/- {vdw_block_sem:6.3f} [{vdw_stdev:6.3f}/{vdw_std:6.3f}/{vdw_sem:6.3f}] / {vdw_block_sd:6.3f} "
f"|{eel_avg:9.3f} +/- {eel_block_sem:6.3f} [{eel_stdev:6.3f}/{eel_std:6.3f}/{eel_sem:6.3f}] / {eel_block_sd:6.3f} "
f"|{pol_avg:9.3f} +/- {pol_block_sem:6.3f} [{pol_stdev:6.3f}/{pol_std:6.3f}/{pol_sem:6.3f}] / {pol_block_sd:6.3f} "
f"|{sas_avg:9.3f} +/- {sas_block_sem:6.3f} [{sas_stdev:6.3f}/{sas_std:6.3f}/{sas_sem:6.3f}] / {sas_block_sd:6.3f} "
f"|{tot_avg:9.3f} +/- {tot_block_sem:6.3f} [{tot_stdev:6.3f}/{tot_std:6.3f}/{tot_sem:6.3f}] / {tot_block_sd:6.3f}")
if _output_format:
text.append([])
else:
text.append('')
return text if _output_format else '\n'.join(text) + '\n\n'
class DecompBinding(dict):
""" Class for decomposition binding (per-residue) """
def __init__(self, com, rec, lig, INPUT, desc=None, **kwargs):
"""
output should be an open file and csvfile should be a csv.writer class. If
the output format is specified as csv, then output should be a csv.writer
class as well.
"""
super(DecompBinding, self).__init__(**kwargs)
self.com, self.rec, self.lig = com, rec, lig
self.num_terms = self.com.num_terms
self.desc = desc # Description
self.INPUT = INPUT
self.idecomp = INPUT['decomp']['idecomp']
self.verbose = INPUT['decomp']['dec_verbose']
# Set up the data for the DELTAs
if self.verbose in [1, 3]:
self.allowed_tokens = 'TDC', 'SDC', 'BDC'
else:
self.allowed_tokens = ('TDC',)
for token in self.allowed_tokens:
self[token] = {}
# Parse everything
self._parse_all_begin()
def _print_vectors(self, csvwriter):
tokens = {'TDC': 'Total Decomposition Contribution (TDC)',
'SDC': 'Sidechain Decomposition Contribution (SDC)',
'BDC': 'Backbone Decomposition Contribution (BDC)'}
for term in self.allowed_tokens:
csvwriter.writerow([tokens[term]])
csvwriter.writerow(['Frame #', 'Residue', 'Internal', 'van der Waals', 'Electrostatic', 'Polar Solvation',
'Non-Polar Solv.', 'TOTAL'])
c = self.INPUT['general']['startframe']
for i in range(self.com.numframes):
for res in self[term]:
csvwriter.writerow([c, res] + [round(self[term][res][key][i], 2) for key in self[term][res]])
c += self.INPUT['general']['interval']
def _parse_all_begin(self):
""" Parses through all of the terms in all of the frames, but doesn't
do any printing
"""
# For per-residue decomp, we need terms to match up
if self.com.num_terms != (self.rec.num_terms + self.lig.num_terms):
raise DecompError('Mismatch in number of decomp terms!')
for term in self.com:
for res in self.com[term]:
self[term][res] = {}
other_token = self.rec[term][res] if res.startswith('R') else self.lig[term][res]
for e in self.com[term][res]:
self[term][res][e] = self.com[term][res][e] - other_token[e]
def summary(self, output_format: str = 'ascii'):
_output_format = 0 if output_format == 'ascii' else 1
text = []
if _output_format:
text.extend([[idecompString[self.idecomp]], [self.desc], []])
else:
text.extend((idecompString[self.idecomp], self.desc))
if self.verbose > 1:
if _output_format:
text.extend(self.com.summary(output_format))
text.extend(self.rec.summary(output_format))
text.extend(self.lig.summary(output_format))
else:
text.extend((self.com.summary(), self.rec.summary(), self.lig.summary()))
if _output_format:
text.append(['DELTAS:'])
else:
text.append('DELTAS:')
for term in self:
if _output_format:
text.extend([[DecompOut.descriptions[term]],
['Residue', 'Internal', '', '', 'van der Waals', '', '', 'Electrostatic',
'', '', 'Polar Solvation', '', '', 'Non-Polar Solv.', '', '', 'TOTAL', '', ''],
[''] + ['Avg.', 'SD(Prop.)', 'SD', 'SEM(Prop.)', 'SEM', 'Block SD', 'Block SEM'] * 6])
else:
text.extend([DecompOut.descriptions[term],
'Residue | Internal Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| van der Waals Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| Electrostatic Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| Polar Solvation Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| Non-Polar Solv. Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| TOTAL Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD',
'----------------------------------------------------------------------------------'
'----------------------------------------------------------------'])
for res in self[term]:
int_avg, int_stdev, int_std, int_semp, int_sem, int_block_sd, int_block_sem = _vector_statistics(self[term][res]['int'])
vdw_avg, vdw_stdev, vdw_std, vdw_semp, vdw_sem, vdw_block_sd, vdw_block_sem = _vector_statistics(self[term][res]['vdw'])
eel_avg, eel_stdev, eel_std, eel_semp, eel_sem, eel_block_sd, eel_block_sem = _vector_statistics(self[term][res]['eel'])
pol_avg, pol_stdev, pol_std, pol_semp, pol_sem, pol_block_sd, pol_block_sem = _vector_statistics(self[term][res]['pol'])
sas_avg, sas_stdev, sas_std, sas_semp, sas_sem, sas_block_sd, sas_block_sem = _vector_statistics(self[term][res]['sas'])
tot_avg, tot_stdev, tot_std, tot_semp, tot_sem, tot_block_sd, tot_block_sem = _vector_statistics(self[term][res]['tot'])
if _output_format:
text.append([res,
int_avg, int_stdev, int_std, int_semp, int_sem, int_block_sd, int_block_sem,
vdw_avg, vdw_stdev, vdw_std, vdw_semp, vdw_sem, vdw_block_sd, vdw_block_sem,
eel_avg, eel_stdev, eel_std, eel_semp, eel_sem, eel_block_sd, eel_block_sem,
pol_avg, pol_stdev, pol_std, pol_semp, pol_sem, pol_block_sd, pol_block_sem,
sas_avg, sas_stdev, sas_std, sas_semp, sas_sem, sas_block_sd, sas_block_sem,
tot_avg, tot_stdev, tot_std, tot_semp, tot_sem, tot_block_sd, tot_block_sem])
else:
text.append(f"{res:14s} "
f"|{int_avg:9.3f} +/- {int_block_sem:6.3f} [{int_stdev:6.3f}/{int_std:6.3f}/{int_sem:6.3f}] / {int_block_sd:6.3f} "
f"|{vdw_avg:9.3f} +/- {vdw_block_sem:6.3f} [{vdw_stdev:6.3f}/{vdw_std:6.3f}/{vdw_sem:6.3f}] / {vdw_block_sd:6.3f} "
f"|{eel_avg:9.3f} +/- {eel_block_sem:6.3f} [{eel_stdev:6.3f}/{eel_std:6.3f}/{eel_sem:6.3f}] / {eel_block_sd:6.3f} "
f"|{pol_avg:9.3f} +/- {pol_block_sem:6.3f} [{pol_stdev:6.3f}/{pol_std:6.3f}/{pol_sem:6.3f}] / {pol_block_sd:6.3f} "
f"|{sas_avg:9.3f} +/- {sas_block_sem:6.3f} [{sas_stdev:6.3f}/{sas_std:6.3f}/{sas_sem:6.3f}] / {sas_block_sd:6.3f} "
f"|{tot_avg:9.3f} +/- {tot_block_sem:6.3f} [{tot_stdev:6.3f}/{tot_std:6.3f}/{tot_sem:6.3f}] / {tot_block_sd:6.3f}")
if _output_format:
text.append([])
else:
text.append('')
return text if _output_format else '\n'.join(text)
class PairDecompBinding(DecompBinding):
""" Class for decomposition binding (pairwise) """
def _parse_all_begin(self):
""" Parses through all of the terms in all of the frames, but doesn't
do any printing
"""
for term in self.com:
for res in self.com[term]:
self[term][res] = {}
for res2 in self.com[term][res]:
self[term][res][res2] = {}
if res.startswith('R') and res2.startswith('R'):
# Both residues are in the receptor -- pull the next one
other_token = self.rec[term][res][res2]
elif (
res.startswith('R')
and not res2.startswith('R')
or not res.startswith('R')
and not res2.startswith('L')
):
other_token = {}
else:
# Both residues are in the ligand -- pull the next one
other_token = self.lig[term][res][res2]
for e in self.com[term][res][res2]:
if other_token:
self[term][res][res2][e] = self.com[term][res][res2][e] - other_token[e]
else:
self[term][res][res2][e] = self.com[term][res][res2][e]
def _print_vectors(self, csvwriter):
tokens = {'TDC': 'Total Decomposition Contribution (TDC)',
'SDC': 'Sidechain Decomposition Contribution (SDC)',
'BDC': 'Backbone Decomposition Contribution (BDC)'}
for term in self.allowed_tokens:
csvwriter.writerow([tokens[term]])
csvwriter.writerow(['Frame #', 'Resid 1', 'Resid 2', 'Internal', 'van der Waals', 'Electrostatic',
'Polar Solvation', 'Non-Polar Solv.', 'TOTAL'])
c = self.com.INPUT['general']['startframe']
for i in range(self.com.numframes):
for res in self[term]:
for res2 in self[term][res]:
csvwriter.writerow([c, res, res2] +
[round(self[term][res][res2][key][i], 2) for key in self[term][res][res2]])
c += self.com.INPUT['general']['interval']
def summary(self, output_format: str = 'ascii'):
_output_format = 0 if output_format == 'ascii' else 1
text = []
if _output_format:
text.extend([[idecompString[self.idecomp]], [self.desc], []])
else:
text.extend((idecompString[self.idecomp], self.desc, ''))
if self.verbose > 1:
if _output_format:
text.extend(self.com.summary(output_format))
text.extend(self.rec.summary(output_format))
text.extend(self.lig.summary(output_format))
else:
text.extend((self.com.summary(), self.rec.summary(), self.lig.summary()))
if _output_format:
text.append(['DELTAS:'])
else:
text.append('DELTAS:')
for term in self:
if _output_format:
text.extend([[DecompOut.descriptions[term]],
['Resid 1', 'Resid 2', 'Internal', '', '', 'van der Waals', '', '', 'Electrostatic', '',
'', 'Polar Solvation', '', '', 'Non-Polar Solv.', '', '', 'TOTAL', '', ''],
['', ''] + ['Avg.', 'SD(Prop.)', 'SD', 'SEM(Prop.)', 'SEM', 'Block SD', 'Block SEM'] * 6])
else:
text.extend([DecompOut.descriptions[term],
'Resid 1 | Resid 2 | Internal Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD | '
'van der Waals Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD | Electrostatic Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| Polar Solvation Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD | Non-Polar Solv. Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD '
'| TOTAL Avg +/- Block SEM [SD(Prop.)/SD/SEM] / Block SD',
'-----------------------------------------------------------------------------------------'
'--------------------------------------------------------------------------'])
for res in self[term]:
for res2 in self[term][res]:
int_avg, int_stdev, int_std, int_semp, int_sem, int_block_sd, int_block_sem = _vector_statistics(self[term][res][res2]['int'])
vdw_avg, vdw_stdev, vdw_std, vdw_semp, vdw_sem, vdw_block_sd, vdw_block_sem = _vector_statistics(self[term][res][res2]['vdw'])
eel_avg, eel_stdev, eel_std, eel_semp, eel_sem, eel_block_sd, eel_block_sem = _vector_statistics(self[term][res][res2]['eel'])
pol_avg, pol_stdev, pol_std, pol_semp, pol_sem, pol_block_sd, pol_block_sem = _vector_statistics(self[term][res][res2]['pol'])
sas_avg, sas_stdev, sas_std, sas_semp, sas_sem, sas_block_sd, sas_block_sem = _vector_statistics(self[term][res][res2]['sas'])
tot_avg, tot_stdev, tot_std, tot_semp, tot_sem, tot_block_sd, tot_block_sem = _vector_statistics(self[term][res][res2]['tot'])
if _output_format:
text.append([res, res2,
int_avg, int_stdev, int_std, int_semp, int_sem, int_block_sd, int_block_sem,
vdw_avg, vdw_stdev, vdw_std, vdw_semp, vdw_sem, vdw_block_sd, vdw_block_sem,
eel_avg, eel_stdev, eel_std, eel_semp, eel_sem, eel_block_sd, eel_block_sem,
pol_avg, pol_stdev, pol_std, pol_semp, pol_sem, pol_block_sd, pol_block_sem,
sas_avg, sas_stdev, sas_std, sas_semp, sas_sem, sas_block_sd, sas_block_sem,
tot_avg, tot_stdev, tot_std, tot_semp, tot_sem, tot_block_sd, tot_block_sem])
else:
text.append(f"{res:14s} | {res2:14s} "
f"|{int_avg:9.3f} +/- {int_block_sem:6.3f} [{int_stdev:6.3f}/{int_std:6.3f}/{int_sem:6.3f}] / {int_block_sd:6.3f} "
f"|{vdw_avg:9.3f} +/- {vdw_block_sem:6.3f} [{vdw_stdev:6.3f}/{vdw_std:6.3f}/{vdw_sem:6.3f}] / {vdw_block_sd:6.3f} "
f"|{eel_avg:9.3f} +/- {eel_block_sem:6.3f} [{eel_stdev:6.3f}/{eel_std:6.3f}/{eel_sem:6.3f}] / {eel_block_sd:6.3f} "
f"|{pol_avg:9.3f} +/- {pol_block_sem:6.3f} [{pol_stdev:6.3f}/{pol_std:6.3f}/{pol_sem:6.3f}] / {pol_block_sd:6.3f} "
f"|{sas_avg:9.3f} +/- {sas_block_sem:6.3f} [{sas_stdev:6.3f}/{sas_std:6.3f}/{sas_sem:6.3f}] / {sas_block_sd:6.3f} "
f"|{tot_avg:9.3f} +/- {tot_block_sem:6.3f} [{tot_stdev:6.3f}/{tot_std:6.3f}/{tot_sem:6.3f}] / {tot_block_sd:6.3f}")
if _output_format:
text.append([])
else:
text.append('')
return text if _output_format else '\n'.join(text) + '\n\n'
def _get_cpptraj_surf(fname):
"""
This function will parse out the surface areas printed out by cpptraj in a
standard data file and return it as an EnergyVector instance.
"""
f = open(fname, 'r')
vec = EnergyVector()
for line in f:
if line.startswith('#'):
continue
vec = vec.append(float(line.split()[1]))
return vec

amber_outputs.py at commit 6db310d, under GPL-3.0 · at the source

Overview

Authors: Qiyu Yang¹, Jinming Zhang¹, Yuxuan Dao¹, Zhengrui He¹, Rui Lu¹, Rummana Jaman¹, Yaole Wu¹, Shuqi Liu¹, Conglin Zhang², Zhibi Zhang³, Jiaqi Zhou⁴

ORCID iDs: Jiaqi Zhou

Institutions: Dali University (China); Kunming Medical University (China)

Journal: iScience, volume 29, issue 9, article 117147

Dates: received 25 January 2026; accepted 25 July 2026; published online 13 August 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1016/j.isci.2026.117147 · PMID 42633227 · PMCID PMC13499207 · OpenAlex W7202370640

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: methods / tools (subfield)

Methods: Statistics

Keywords: cyclotriphosphazenes, toxicological evaluation, blood-brain barrier, computational toxicology methods

Topic: Cholinesterase and Neurodegenerative Diseases (Pharmacology, Medicine), according to OpenAlex

Funding: Dali University (KY2396120340); Yunnan Provincial Department of Education (KY2414107040, S202510679099)

Citations: not cited yet (Europe PMC); 57 references in the paper

Abstract

Cyclotriphosphazenes (CTPs), widely used as lithium-ion battery electrolyte additives, are increasingly recognized as environmental contaminants with potential neurotoxicity. We developed an integrated computational framework to evaluate the blood-brain barrier (BBB) penetration potential and neurotoxic mechanisms of CTPs. Six representative CTPs were predicted to penetrate the BBB via passive diffusion and impair BBB function. Among them, HPCTP, PFPCTP, HCCTP, and HFCTP primarily disrupted BBB function through inhibition of ATP-binding cassette transporters. HPCTP exhibited the strongest BBB-disrupting potential, primarily through interactions with SLC2A1 and SLC6A3 and through the modulation of BBB-related gene expression. Experimental validation in hCMEC/D3 cells showed that HPCTP impaired glucose transport and utilization and significantly downregulated SLC2A1 expression. These findings provide a strategy for assessing BBB penetration of exogenous chemicals and highlight the potential CNS risks associated with environmental CTP exposure.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

Valdes-Tresanco-MS/gmx_MMPBSA

License: GPL-3.0

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 6db310d5d48074d442ebac6348ccb135567d95a7, 13 September 2026

Languages: Python (114), JavaScript (4), Shell (2), Jupyter (2)

Size: 802 files, 122 scripts

Software Heritage: not archived

Found in: the text, “Key resources table”

Holds: README, license file, environment (setup.cfg, setup.py, docs/env.yml, docs/requirements.txt), tests, continuous integration, documentation, 2 notebooks

Not found: CITATION.cff

Tools: NumPy (15 files), pandas (11 files), Matplotlib (7 files), seaborn (5 files), SciPy (3 files)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

124 files

GMXMMPBSA/API.py, Python, 1,175 lines
GMXMMPBSA/GMXMMPBSA.sh, Shell, 111 lines
GMXMMPBSA/__init__.py, Python, 99 lines
GMXMMPBSA/_version.py, Python, 683 lines
GMXMMPBSA/alamdcrd.py, Python, 375 lines
GMXMMPBSA/amber_outputs.py, Python, 2,028 lines, 1 match
GMXMMPBSA/analyzer/__init__.py, Python, 19 lines
GMXMMPBSA/analyzer/chartsettings.py, Python, 1,254 lines
GMXMMPBSA/analyzer/customitem.py, Python, 588 lines
GMXMMPBSA/analyzer/dialogs.py, Python, 702 lines
GMXMMPBSA/analyzer/gui.py, Python, 1,421 lines
GMXMMPBSA/analyzer/items_delegate.py, Python, 61 lines
GMXMMPBSA/analyzer/parametertree/Parameter.py, Python, 758 lines
GMXMMPBSA/analyzer/parametertree/ParameterItem.py, Python, 133 lines
GMXMMPBSA/analyzer/parametertree/ParameterTree.py, Python, 125 lines
GMXMMPBSA/analyzer/parametertree/__init__.py, Python, 25 lines
GMXMMPBSA/analyzer/parametertree/parameterTypes.py, Python, 785 lines
GMXMMPBSA/analyzer/plots.py, Python, 992 lines
GMXMMPBSA/analyzer/style/__init__.py, Python, 99 lines
GMXMMPBSA/analyzer/style/app_theme.py, Python, 428 lines
GMXMMPBSA/analyzer/utils.py, Python, 416 lines
GMXMMPBSA/app.py, Python, 237 lines
GMXMMPBSA/calculation.py, Python, 1,391 lines
GMXMMPBSA/commandlineparser.py, Python, 584 lines
GMXMMPBSA/createinput.py, Python, 826 lines
GMXMMPBSA/error_bundle.py, Python, 338 lines
GMXMMPBSA/exceptions.py, Python, 185 lines
GMXMMPBSA/fake_mpi.py, Python, 70 lines
GMXMMPBSA/gbnsr6_topology.py, Python, 172 lines
GMXMMPBSA/infofile.py, Python, 286 lines
GMXMMPBSA/input_parser.py, Python, 748 lines
GMXMMPBSA/logging_utils.py, Python, 122 lines
GMXMMPBSA/main.py, Python, 1,774 lines
GMXMMPBSA/make_top.py, Python, 2,297 lines
GMXMMPBSA/make_top_amber.py, Python, 1,692 lines
GMXMMPBSA/make_trajs.py, Python, 857 lines
GMXMMPBSA/membrane.py, Python, 400 lines
GMXMMPBSA/output_file.py, Python, 691 lines
GMXMMPBSA/parm_setup.py, Python, 886 lines
GMXMMPBSA/progress.py, Python, 353 lines
GMXMMPBSA/qmmm_diagnostics.py, Python, 114 lines
GMXMMPBSA/radii.py, Python, 408 lines, 1 match
GMXMMPBSA/test_manifest.py, Python, 324 lines
GMXMMPBSA/tester.py, Python, 246 lines
GMXMMPBSA/timer.py, Python, 127 lines
GMXMMPBSA/topology_preprocess.py, Python, 86 lines
GMXMMPBSA/utils.py, Python, 1,235 lines
docs/assets/js/custom.js, JavaScript, 107 lines
docs/assets/js/news.js, JavaScript, 40 lines
docs/assets/js/termynal.js, JavaScript, 264 lines
docs/overrides/assets/javascripts/bundle.de8129cf.min.js, JavaScript, 18 lines
examples/API/extract_api_data.py, Python, 50 lines
examples/psf_dcd/protein_protein/script.py, Python, 43 lines
notebooks/gmx_MMPBSA_Colab.ipynb, Jupyter, 526 lines
notebooks/gmx_MMPBSA_Local.ipynb, Jupyter, 261 lines
scripts/bootstrap_test_manifest.py, Python, 152 lines
scripts/conda_pip_install.sh, Shell, 31 lines
scripts/example_readme_sync.py, Python, 57 lines
scripts/run_example_docs_symlink_spike.py, Python, 164 lines
scripts/sync_example_docs.py, Python, 79 lines
scripts/validate_example_readme_parity.py, Python, 67 lines
scripts/validate_gmx_MMPBSA_test_docs.py, Python, 233 lines
scripts/validation/__init__.py, Python, 1 line
scripts/validation/common.py, Python, 538 lines
scripts/validation/compare_api_results.py, Python, 211 lines
scripts/validation/compare_combination_matrix.py, Python, 111 lines
scripts/validation/compare_legacy_matrix.py, Python, 169 lines
scripts/validation/export_validation_csv.py, Python, 335 lines
scripts/validation/plan_combination_matrix.py, Python, 114 lines
scripts/validation/run_analyzer_validation.py, Python, 78 lines
scripts/validation/run_api_validation.py, Python, 184 lines
scripts/validation/run_calculation_matrix.py, Python, 65 lines
scripts/validation/run_cleanup_validation.py, Python, 107 lines
scripts/validation/run_combination_matrix.py, Python, 242 lines
scripts/validation/run_concurrency_validation.py, Python, 155 lines
scripts/validation/run_disk_analyzer_validation.py, Python, 54 lines
scripts/validation/run_mpi_validation.py, Python, 113 lines
scripts/validation/run_negative_validation.py, Python, 77 lines
scripts/validation/run_packaging_validation.py, Python, 93 lines
setup.py, Python, 64 lines
tests/test_amber_baseline.py, Python, 312 lines
tests/test_amber_mtp.py, Python, 124 lines
tests/test_analyzer.py, Python, 58 lines
tests/test_api_stability.py, Python, 145 lines
tests/test_calculation.py, Python, 308 lines
tests/test_cleanup.py, Python, 86 lines
tests/test_commandlineparser.py, Python, 220 lines
tests/test_composite_mutations.py, Python, 396 lines
tests/test_createinput.py, Python, 160 lines
tests/test_energy_parse_golden.py, Python, 96 lines
tests/test_error_bundle.py, Python, 105 lines
tests/test_error_handling.py, Python, 173 lines
tests/test_example_docs_symlink_spike.py, Python, 42 lines
tests/test_explicit_waters.py, Python, 571 lines
tests/test_gbnsr6_topology.py, Python, 83 lines
tests/test_input_parser.py, Python, 153 lines
tests/test_logging.py, Python, 141 lines
tests/test_make_top.py, Python, 203 lines
tests/test_make_trajs.py, Python, 174 lines
tests/test_membrane.py, Python, 100 lines
tests/test_output_file.py, Python, 231 lines
tests/test_pb_inp1_nonpolar.py, Python, 129 lines
tests/test_phase10_docs.py, Python, 332 lines
tests/test_phase4_logging_and_validation.py, Python, 325 lines
tests/test_phase5_severity.py, Python, 75 lines
tests/test_phase6_messages.py, Python, 50 lines
tests/test_phase7_topology_messages.py, Python, 37 lines
tests/test_phase9_errors.py, Python, 63 lines
tests/test_progress.py, Python, 231 lines
tests/test_qmmm_convergence.py, Python, 160 lines
tests/test_qmmmgbsa_support.py, Python, 62 lines
tests/test_radius_compatibility.py, Python, 64 lines
tests/test_radius_provenance.py, Python, 228 lines
tests/test_reference_matching.py, Python, 17 lines
tests/test_require_complex_topology.py, Python, 76 lines
tests/test_statistics.py, Python, 89 lines
tests/test_test_manifest.py, Python, 97 lines
tests/test_tester.py, Python, 57 lines
tests/test_utils_create_input_args.py, Python, 39 lines
tests/test_validate_gmx_MMPBSA_test_docs.py, Python, 50 lines
tests/test_validation_scripts.py, Python, 216 lines
versioneer.py, Python, 2,277 lines
LICENSE.txt, License, 674 lines
README.md, Text, 58 lines

At the source: github.com/Valdes-Tresanco-MS/gmx_MMPBSA

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
122 scripts, each with its path and the digest of its content;
2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data and code availability

All data generated or analyzed during this study are included in this article and its supplemental information. Additional raw data supporting the findings of this study are available from the corresponding authors upon reasonable request. No original code was developed for this study; all computational analyses were performed using publicly available software and online databases as described in the STAR Methods section.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 3, 28 September 2026

Authors: added Jiaqi Zhou (0009-0004-5802-5697); removed Jiaqi Zhou

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 11 authors, 4 keywords, 2 funders, 57 references.

Cite

This paper

Yang, Q., Zhang, J., Dao, Y., He, Z., Lu, R., Jaman, R., Wu, Y., Liu, S., Zhang, C., Zhang, Z., & Zhou, J. (2026). Penetrating potential evaluation of cyclotriphosphazenes target blood-brain barrier via computational toxicology methods. iScience, 29(9), 117147. https://doi.org/10.1016/j.isci.2026.117147

BibTeX

@article{yang2026penetrating, author = {Yang, Qiyu and Zhang, Jinming and Dao, Yuxuan and He, Zhengrui and Lu, Rui and Jaman, Rummana and Wu, Yaole and Liu, Shuqi and Zhang, Conglin and Zhang, Zhibi and Zhou, Jiaqi}, title = {{Penetrating potential evaluation of cyclotriphosphazenes target blood-brain barrier via computational toxicology methods}}, journal = {iScience}, year = {2026}, month = aug, volume = {29}, number = {9}, pages = {117147}, publisher = {Elsevier}, issn = {2589-0042}, doi = {10.1016/j.isci.2026.117147}, url = {https://doi.org/10.1016/j.isci.2026.117147}, pmid = {42633227}, pmcid = {PMC13499207} }

RIS

TY - JOUR AU - Yang, Qiyu AU - Zhang, Jinming AU - Dao, Yuxuan AU - He, Zhengrui AU - Lu, Rui AU - Jaman, Rummana AU - Wu, Yaole AU - Liu, Shuqi AU - Zhang, Conglin AU - Zhang, Zhibi AU - Zhou, Jiaqi TI - Penetrating potential evaluation of cyclotriphosphazenes target blood-brain barrier via computational toxicology methods T2 - iScience J2 - iScience PY - 2026 DA - 2026/08/13 VL - 29 IS - 9 SP - 117147 SN - 2589-0042 PB - Elsevier DO - 10.1016/j.isci.2026.117147 UR - https://doi.org/10.1016/j.isci.2026.117147 LA - en ER -

CSL-JSON

{ "id": "10.1016/j.isci.2026.117147", "type": "article-journal", "title": "Penetrating potential evaluation of cyclotriphosphazenes target blood-brain barrier via computational toxicology methods", "container-title": "iScience", "author": [ { "family": "Yang", "given": "Qiyu" }, { "family": "Zhang", "given": "Jinming" }, { "family": "Dao", "given": "Yuxuan" }, { "family": "He", "given": "Zhengrui" }, { "family": "Lu", "given": "Rui" }, { "family": "Jaman", "given": "Rummana" }, { "family": "Wu", "given": "Yaole" }, { "family": "Liu", "given": "Shuqi" }, { "family": "Zhang", "given": "Conglin" }, { "family": "Zhang", "given": "Zhibi" }, { "family": "Zhou", "given": "Jiaqi" } ], "container-title-short": "iScience", "volume": "29", "issue": "9", "page": "117147", "DOI": "10.1016/j.isci.2026.117147", "PMID": "42633227", "PMCID": "PMC13499207", "ISSN": "2589-0042", "publisher": "Elsevier", "URL": "https://doi.org/10.1016/j.isci.2026.117147", "language": "en", "issued": { "date-parts": [ [ 2026, 8, 13 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.