diff --git a/CodeEntropy/PoseidonHelper.py b/CodeEntropy/PoseidonHelper.py deleted file mode 100644 index 8cd3e08..0000000 --- a/CodeEntropy/PoseidonHelper.py +++ /dev/null @@ -1,164 +0,0 @@ -import logging -import sys -from datetime import datetime - -# from CodeEntropy.poseidon.analysis.EECalculation import processEE -# from CodeEntropy.poseidon.analysis.helper import memoryInfo, weightingPopulation -# from CodeEntropy.poseidon.analysis.populateClasses import classPopulation -from CodeEntropy.poseidon.extractData.dihedrals import calculateDihedrals -from CodeEntropy.poseidon.extractData.forceTorques import calculateFTMatrix - -# # Energy is not needed -# from CodeEntropy.poseidon.extractData.readFiles import populateEnergy, UAEnergyGroup -from CodeEntropy.poseidon.extractData.HBRAD import HBCalc, UALevelRAD, distCutoffNc -from CodeEntropy.poseidon.extractData.mainClass import clearClass -from CodeEntropy.poseidon.extractData.nearestNonlike2 import ( - getShellAssignment, - moleculePositionRankingRAD, -) - -# from CodeEntropy.poseidon.extractData.outputFiles import moleculeObjectPopulation -from CodeEntropy.poseidon.extractData.readFiles import ( # populateTopology, - getCoordsForces, - getDistArray, -) - - -def frame_iteration( - container, - all_data, - dimensions, - startTime, - verbosePrint, - waterTuple, - cutShell, - excludedResnames, - frame, -): - clearClass(all_data) - print(f"frame = {frame}") - - all_data, dimensions = getCoordsForces( - container, all_data, dimensions, frame, startTime, verbosePrint - ) - # # Energy is not needed - # populateEnergy(container, all_data, - # dimensions, frame, startTime, verbosePrint) - # UAEnergyGroup(all_data) - - calculateDihedrals(all_data, dimensions) - verbosePrint("DIH") - verbosePrint(datetime.now() - startTime) - sys.stdout.flush() - - traj = container.trajectory[frame] - neighbour_coords = None - neighbour_coords = traj.positions - - max_cutoff = 10 - for x in range(0, len(all_data)): - atom = all_data[x] - # start nearest array from solutes - if atom.resname not in waterTuple: - getDistArray( - atom, - all_data, - traj, - max_cutoff, - dimensions, - neighbour_coords, - startTime, - verbosePrint, - ) - # find nearest array for solute neighbours that are solvent - for nearDist in atom.nearest_all_atom_array[0:20]: - neighbour = all_data[nearDist[0]] - if ( - neighbour.resname in waterTuple - and neighbour.nearest_all_atom_array is None - ): - getDistArray( - neighbour, - all_data, - traj, - max_cutoff, - dimensions, - neighbour_coords, - startTime, - verbosePrint, - ) - # find solvent neighbours neighbours nearest array - if neighbour.nearest_all_atom_array is not None: - for nearDist2 in neighbour.nearest_all_atom_array[0:20]: - neighbour2 = all_data[nearDist2[0]] - if ( - neighbour2.resname in waterTuple - and neighbour2.nearest_all_atom_array is None - ): - getDistArray( - neighbour2, - all_data, - traj, - max_cutoff, - dimensions, - neighbour_coords, - startTime, - verbosePrint, - ) - else: - continue - else: - continue - else: - continue - verbosePrint("NEAREST ARRAYS") - verbosePrint(datetime.now() - startTime) - sys.stdout.flush() - - if cutShell is not None: - # Used for fixed cut-off coordination shells - try: - cutoff_dist = float(cutShell) - distCutoffNc(all_data, dimensions, cutoff_dist) - # #don't know what this is it is not defined anywhere - # NcPairs(all_data, dimensions) - verbosePrint("distCutoffNc") - verbosePrint(datetime.now() - startTime) - sys.stdout.flush() - except ValueError: - logging.error("Cutoff distance needs to be a float") - - if cutShell is None: - UALevelRAD(all_data, dimensions) - verbosePrint("RAD") - verbosePrint(datetime.now() - startTime) - - # if inputType != 'pdb': - HBCalc(all_data, waterTuple, dimensions) - verbosePrint("HB") - verbosePrint(datetime.now() - startTime) - sys.stdout.flush() - - getShellAssignment(all_data, excludedResnames, dimensions, startTime, verbosePrint) - verbosePrint("PROX") - verbosePrint(datetime.now() - startTime) - sys.stdout.flush() - - # if (force is not None): - calculateFTMatrix(all_data, dimensions) - verbosePrint("FTMATRIX") - verbosePrint(datetime.now() - startTime) - sys.stdout.flush() - - moleculePositionRankingRAD(all_data, waterTuple, dimensions) - verbosePrint("ORIENTS") - verbosePrint(datetime.now() - startTime) - sys.stdout.flush() - - """ - if iterations == 1: - verbosePrint('Generating .pdb file for frame %s' % (frame+1)) - pdbGenerate(all_data, ('frame_%s' % (frame+1)), dimensions) - """ - # print(all_data, frame, dimensions) - return (all_data, frame, dimensions) diff --git a/CodeEntropy/poseidon/__init__.py b/CodeEntropy/poseidon/__init__.py deleted file mode 100644 index 201e256..0000000 --- a/CodeEntropy/poseidon/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -""" -Empty init file -""" diff --git a/CodeEntropy/poseidon/analysis/EECalculation.py b/CodeEntropy/poseidon/analysis/EECalculation.py deleted file mode 100644 index d3ee34d..0000000 --- a/CodeEntropy/poseidon/analysis/EECalculation.py +++ /dev/null @@ -1,1719 +0,0 @@ -#!/usr/bin/env python - -# import logging -import math - -# import sys -from collections import Counter, defaultdict - -import numpy as np -import pandas as pd -from numpy import linalg as LA - -from CodeEntropy.FunctionCollection import Utils - - -def nested_dict(): - return defaultdict(nested_dict) - - -# create nested dict in one go - - -def processEE( - num_frames, - totFrames, - Aclass, - solvent, - waterTuple, - temperature, - level, - name, - forceUnits, - verbosePrint, -): - """ - Populate and process Smix and Sor into global dictionaries - """ - - verbosePrint("\n\n---ANALYSIS_TYPE_%s---\n" % (name)) - verbosePrint("\n\n---ORIENTATIONAL_ENTROPY_CALCULATIONS---\n") - if any(i in solvent for i in waterTuple) is True: - SorCalculation(Aclass, "Sor_test2", level, waterTuple, verbosePrint) - - verbosePrint("\n\n --UA_VIBRATIONAL_ENTROPY--\n") - SvibCalculations(totFrames, Aclass, temperature, forceUnits, verbosePrint) - verbosePrint("\n\n --WM_VIBRATIONAL_ENTROPY--\n") - WMSvibCalculations(Aclass, totFrames, temperature, forceUnits, verbosePrint) - - verbosePrint("\n\n --WM_CONFORMATIONAL_ENTROPY(DIH)--\n") - process_dihedrals(Aclass, verbosePrint) - - verbosePrint("\n\n--PRINT_ALL_VARIABLES--\n") - # solventData, soluteData = saveAllVariables( - # num_frames, Aclass, level, name, solvent, waterTuple, verbosePrint - # ) - solventData, soluteData = saveAllVariables( - num_frames, Aclass, level, name, solvent, waterTuple, verbosePrint - ) - - if level in ["residLevel_resname", "atomLevel"]: - verbosePrint("\n\n--CONTACTS--\n") - # contactMatrix = contactCalculation(Aclass, level, totFrames, verbosePrint) - contactMatrix = contactCalculation(Aclass, level, totFrames, verbosePrint) - - return { - "solventData": solventData, - "soluteData": soluteData, - "contactMatrix": contactMatrix, - } - - return {"solventData": solventData, "soluteData": soluteData} - - -def contactCalculation(Aclass, level, totFrames, verbosePrint): - """ - resid matrix of RAD contacts between resids that are not water and - the nearest non-like to it. - """ - - verbosePrint("\n*******_RESID_CONTACTS_**********\n") - # data = open('resid_contact_matrix_%s.csv' % (level), 'w') - contact_matrix_name = f"resid_contact_matrix_{level}.csv" - if level == "atomLevel": - # data.write('\n'.join(['centre_resid,neighbour_resid,count,'\ - # 'centre_resname,neighbour_resname,centre_atom,'\ - # 'neighbour_atom']) - # + '\n') - Utils.printOut( - contact_matrix_name, - "centre_resid," - "neighbour_resid," - "count," - "centre_resname," - "neighbour_resname," - "centre_atom," - "neighbour_atom", - ) - contactMatrix = pd.DataFrame( - columns=[ - "centre_resid", - "neighbour_resid", - "count", - "centre_resname", - "neighbour_resname", - "centre_atom", - "neighbour_atom", - ] - ) - else: - # data.write('\n'.join(['centre_resid,neighbour_resid,count,'\ - # 'centre_resname,neighbour_resname']) - # + '\n') - Utils.printOut( - contact_matrix_name, - "centre_resid,neighbour_resid,count,centre_resname,neighbour_resname", - ) - contactMatrix = pd.DataFrame( - columns=[ - "centre_resid", - "neighbour_resid", - "count", - "centre_resname", - "neighbour_resname", - ] - ) - for centre_info, neighbour_resid_key in sorted( - list(Aclass.resid_contact_matrix_dict.items()) - ): - for neighbour_info, count in sorted(list(neighbour_resid_key.items())): - # verbosePrint(centre_info, neighbour_info, count) - if count > 0: - if level == "atomLevel": - centre_atom = centre_info[0] - centre_resname = centre_info[1] - centre_resid = centre_info[2] - neighbour_atom = neighbour_info[0] - neighbour_resname = neighbour_info[1] - neighbour_resid = neighbour_info[2] - count = float(count) / float(totFrames) - # data.write('\n'.join(['%s,%s,%s,%s,%s,%s,%s' - # % (centre_resid, - # neighbour_resid, count, centre_resname, - # neighbour_resname, centre_atom, - # neighbour_atom)]) + '\n') - Utils.printOut( - contact_matrix_name, - ( - f"{centre_resid},{neighbour_resid},{count}," - f"{centre_resname},{neighbour_resname}," - f"{centre_atom},{neighbour_atom}" - ), - ) - newRowContact = pd.DataFrame( - { - "centre_resid": centre_resid, - "neighbour_resid": neighbour_resid, - "count": count, - "centre_resname": centre_resname, - "neighbour_resname": neighbour_resname, - "centre_atom": centre_atom, - "neighbour_atom": neighbour_atom, - }, - index=[0], - ) - contactMatrix = pd.concat( - [contactMatrix, newRowContact], ignore_index=True - ) - else: - centre_resname = centre_info[0] - centre_resid = centre_info[1] - neighbour_resname = neighbour_info[0] - neighbour_resid = neighbour_info[1] - count = float(count) / float(totFrames) - # data.write('\n'.join(['%s,%s,%s,%s,%s' % (centre_resid, - # neighbour_resid, count, centre_resname, - # neighbour_resname)]) + '\n') - Utils.printOut( - contact_matrix_name, - ( - f"{centre_resid},{neighbour_resid}," - f"{count},{centre_resname}," - f"{neighbour_resname}" - ), - ) - newRowContact = pd.DataFrame( - { - "centre_resid": centre_resid, - "neighbour_resid": neighbour_resid, - "count": count, - "centre_resname": centre_resname, - "neighbour_resname": neighbour_resname, - }, - index=[0], - ) - contactMatrix = pd.concat( - [contactMatrix, newRowContact], ignore_index=True - ) - else: - continue - - # data.close() - return contactMatrix - - -def SorCalculation(Aclass, methodType, level, waterTuple, verbosePrint): - """ """ - - # por_dict = nested_dict() - shells_ppD_ppA_dict = nested_dict() - verbosePrint("\n\n__ORIENTATIONAL_ENTROPY_%s" % (methodType)) - num_molecules0 = sum( - count - for atomType in Aclass.RADshell_Sor_dict.values() - # for Nc in atomType.values() - for RADshell_num in atomType.values() - for RAD_str in RADshell_num.values() - for ADs in RAD_str.values() - for key, HBcount in ADs.items() - for HB, count in HBcount.items() - if key == "D" - ) - verbosePrint("\nall molecule count:", num_molecules0) - S1 = 0 - for nearest, assigned_key in sorted(list(Aclass.RADshell_Sor_dict.items())): - num_molecules1 = sum( - count - # for Nc in assigned_key.values() - for RADshell_num in assigned_key.values() - for RAD_str in RADshell_num.values() - for ADs in RAD_str.values() - for AD, HBcount in ADs.items() - for HB, count in HBcount.items() - if AD == "D" - ) - - nearest_resname = nearest.split("_")[0] - verbosePrint( - "\n\nnearest:", nearest, "count:", num_molecules1, "/", num_molecules0 - ) - - S2 = 0 - for assigned, shell_num_key in sorted(list(assigned_key.items())): - num_molecules2 = sum( - count - # for RADshell_num in Nc_key.values() - for RAD_str in shell_num_key.values() - for ADs in RAD_str.values() - for AD, HBcount in ADs.items() - for HB, count in HBcount.items() - if AD == "D" - ) - - # assigned_resname = assigned[0].split("_")[0] - verbosePrint( - "\tassigned:", assigned, "count:", num_molecules2, "/", num_molecules1 - ) - N_H = assigned[1] - verbosePrint("\tN_H", N_H) - S4 = 0 - for shell_num, RADlist_key in sorted(list(shell_num_key.items())): - - num_molecules4 = sum( - count - for ADs in RADlist_key.values() - for AD, HBcount in ADs.items() - for HB, count in HBcount.items() - if AD == "D" - ) - - verbosePrint( - "\t\t\tshell_num:", - shell_num[0], - "count:", - num_molecules4, - "/", - num_molecules2, - ) - S5 = 0 - S_shell = 0 - ave_pbias_ave = 0 - ave_Nc_eff = 0 - - for RADlist, AD_key in sorted(list(RADlist_key.items())): - verbosePrint("\t" * 4, "-----RADshell------") - num_molecules5 = sum( - count - for AD, HBcount in AD_key.items() - for HB, count in HBcount.items() - if AD == "D" - ) - - # dict for p(c) log p(c) - if ( - methodType == "Sor_test2" - and level == "moleculeLevel" - and assigned[0] in waterTuple - and int(shell_num[0]) == 1 - ): - newRADlist = [] - for n in RADlist: - if n == "0": - newRADlist.append(nearest_resname) - else: - newRADlist.append(n) - newRADlist = tuple(sorted(newRADlist)) - try: - if newRADlist not in Aclass.p_coord_dict[len(RADlist)]: - Aclass.p_coord_dict[len(RADlist)][newRADlist] = 0 - Aclass.p_coord_dict[len(RADlist)][ - newRADlist - ] += num_molecules5 - except AttributeError: - pass - - RAD_l = "%s" % (",".join(map(str, RADlist))) - verbosePrint( - "\t\t\t\tneighbours:", - RAD_l, - "Nc:", - len(RADlist), - "count:", - num_molecules5, - "/", - num_molecules4, - ) - S6 = 0 - Sor_list = [] - pbias_ave = 0 - Nc_eff = 0 - - if methodType == "Sor_test2": - D_values = None - A_values = None - for AD, values in sorted(list(AD_key.items())): - if AD == "D": - D_values = values - if AD == "A": - A_values = values - verbosePrint("\t" * 4, "* SEPARATE ADs *") - S7, pbias_ave, Nc_eff, ppD_ppA_dict = RAD_or_entropyCalc_test2( - RADlist, - D_values, - A_values, - N_H, - shell_num[1], - num_molecules5, - nearest_resname, - waterTuple, - verbosePrint, - ) - - # ppD_ppA_dict population - for c in range(0, int(num_molecules5)): - for i, pp in ppD_ppA_dict.items(): - ppD = pp[0] - ppA = pp[1] - if ( - i - not in shells_ppD_ppA_dict[nearest][assigned][ - float(shell_num[0]) - ] - ): - shells_ppD_ppA_dict[nearest][assigned][ - float(shell_num[0]) - ][i] = [0, 0, 0] - ppD_stored = shells_ppD_ppA_dict[nearest][assigned][ - float(shell_num[0]) - ][i][0] - ppA_stored = shells_ppD_ppA_dict[nearest][assigned][ - float(shell_num[0]) - ][i][1] - count_stored = shells_ppD_ppA_dict[nearest][assigned][ - float(shell_num[0]) - ][i][2] - - count_add = count_stored + 1 - ppD_add = (ppD + ppD_stored * count_stored) / float( - count_add - ) - ppA_add = (ppA + ppA_stored * count_stored) / float( - count_add - ) - - shells_ppD_ppA_dict[nearest][assigned][ - float(shell_num[0]) - ][i][0] = ppD_add - shells_ppD_ppA_dict[nearest][assigned][ - float(shell_num[0]) - ][i][1] = ppA_add - shells_ppD_ppA_dict[nearest][assigned][ - float(shell_num[0]) - ][i][2] = count_add - - Sor_list.append([S7, "DA"]) - verbosePrint("\t\t\t\tS7_DA = ", round(S7, 3)) - - Sor_list = sorted(Sor_list) - verbosePrint() - verbosePrint("\t" * 4, "*** Sor-list:") - c = 0 - for s, label in Sor_list: - c += 1 - if c == 1: - verbosePrint( - "\t" * 4, round(s, 3), label, "<-- selected Sor" - ) - else: - verbosePrint("\t" * 4, round(s, 3), label) - # use smallest Sor as selected value - S6 += Sor_list[0][0] * float(num_molecules5) / float(num_molecules4) - S_shell += ( - Sor_list[0][0] * float(num_molecules5) / float(num_molecules4) - ) - ave_pbias_ave += ( - pbias_ave * float(num_molecules5) / float(num_molecules4) - ) - - ave_Nc_eff += Nc_eff * float(num_molecules5) / float(num_molecules4) - - if methodType == "Sor_test2": - Aclass.Sor_reference_dict[nearest][assigned[0]][ - str(shell_num[0]) - ][RADlist] = (Sor_list[0][0] * 8.314) - # print(RADlist, Sor_list[0][0] * 8.314) - - verbosePrint("\t\t\t\tS6 = ", round(S6, 3)) - verbosePrint() - S5 += S6 * float(num_molecules4) / float(num_molecules2) - - Aclass.allVariables_dict[nearest][assigned][float(shell_num[0])][ - "%s" % (methodType) - ] = (S_shell * 8.314, num_molecules4) - - """ - if methodType == 'Sor_test2': - Aclass.allVariables_dict[nearest][assigned]\ - [float(shell_num[0])]\ - ['pbias_ave_%s' % (methodType)] = \ - ave_pbias_ave, num_molecules4 - - Aclass.allVariables_dict[nearest][assigned]\ - [float(shell_num[0])]\ - ['Nc_eff_%s' % (methodType)] = \ - ave_Nc_eff, num_molecules4 - """ - - verbosePrint("\t\t\tS5 = ", round(S5, 3)) - S4 += S5 * float(num_molecules4) / float(num_molecules2) - - verbosePrint("\t\tS4 = ", round(S4, 3)) - S2 += S4 * float(num_molecules1) / float(num_molecules0) - - verbosePrint("\tS2 = ", round(S2, 3)) - S1 += S2 - verbosePrint("S1 = ", round(S1, 3)) - - if methodType == "Sor_test2": - verbosePrint("\n\n__ppX_ORIENTATIONAL_ENTROPY_%s" % (methodType)) - verbosePrint("nearest,assigned,shell_num,i,ppD,ppA,count") - for nearest, assigned_key in sorted(list(shells_ppD_ppA_dict.items())): - for assigned, shell_num_key in sorted(list(assigned_key.items())): - for shell_num, i_key in sorted(list(shell_num_key.items())): - for i, pp in sorted(list(i_key.items())): - ppD, ppA, count = pp[0], pp[1], pp[2] - verbosePrint( - "%s,%s,%s,%s,%s,%s,%s" - % ( - nearest, - assigned[0], - int(shell_num), - i, - round(ppD, 5), - round(ppA, 5), - count, - ) - ) - - -def RAD_or_entropyCalc_test2( - shell, - shellD, - shellA, - N_H, - shell_num, - referenceCount, - nearest, - waterTuple, - verbosePrint, -): - """ - ***** S_or - """ - - # get info about shell - shellCount = Counter([str(N) for N in shell]) - # count constituents in shell - Nc = sum(shellCount.values()) - Nc_solvent_only = 0 - for n in shell: - if n == "X": - Nc_solvent_only += 1 - else: - try: - n = int(n) - if n != 0: - Nc_solvent_only += 1 - else: - continue - except ValueError: - pass - - soluteA_pi, soluteD_pi = None, None - pA_degen_dict, pD_degen_dict = None, None - for AD, shellDs in [["A", shellA], ["D", shellD]]: - verbosePrint("\t" * 5, "__%s__" % (AD)) - # get info about shell - shellCount = Counter([str(N) for N in shell]) - # count constituents in shell - Dcount = Counter([tuple(Ds) for Ds in shellDs]) - num_orients = len(shellDs) - - solute_donor = False - solvent_donor = False - - for ds in shellDs: - for d in ds: - try: - d = int(d) - if d == 0: - solute_donor = True - else: - continue - except ValueError: - if d in ["X", "H"]: - solvent_donor = True - else: - solute_donor = True - - # pi only, degeneracy info for each orientation type observed - i_counts_dict = {} - tot_i_count = 0 - for donors, count in Dcount.items(): - for i in donors: - if i not in i_counts_dict: - i_counts_dict[i] = 0 - i_counts_dict[i] += count - tot_i_count += count - - verbosePrint("\t" * 5, "--i--") - for i, count in i_counts_dict.items(): - verbosePrint("\t" * 5, i, count) - verbosePrint("\t" * 5, "sum i counts:", tot_i_count) - - pi_degen_dict = {} - solute_pis = 0 - largest_pi = 0 - ww_Ds = 0 - for i, count in i_counts_dict.items(): - degeneracy = shellCount[i] - if degeneracy == 0: - degeneracy = 1 - pi = (float(count) / float(tot_i_count)) / float(degeneracy) - pi_degen_dict[i] = [pi, degeneracy, count] - verbosePrint( - "\t" * 6, - "pi", - i, - round(pi, 4), - "degen:", - degeneracy, - "pi_d:", - round((pi * degeneracy), 4), - ) - # verbosePrint(i, pi) - if pi * degeneracy > largest_pi: - largest_pi = pi * degeneracy - - solute_pi, solvent_pi = False, False - try: - i = int(i) - if i == 0 and nearest not in waterTuple: - solute_pi = True - else: - ww_Ds += count - continue - except ValueError: - if i in ["X", "H"]: - solvent_pi = True - ww_Ds += count - else: - solute_pi = True - if solute_pi is True: - solute_pis += pi * degeneracy - - # calculate number of observed orientations - # pi_max = max(map(lambda x: x[0], pi_degen_dict.values())) - # verbosePrint('\t'*6, 'pi_max', round(pi_max, 5)) - - # verbosePrint(ww_Ds, referenceCount) - # per_w_Ds = ww_Ds / float(referenceCount) - - if AD == "A": - pA_degen_dict = pi_degen_dict - if AD == "D": - pD_degen_dict = pi_degen_dict - - # pij, degeneracy info for each orientation type observed - verbosePrint("\t" * 5, "--ij--") - pij_degen_dict = {} - for donors, count in Dcount.items(): - verbosePrint("\t" * 5, donors, count) - donor_constituents = Counter([str(dc) for dc in donors]) - degeneracy = 1 - sum_weight = 0 - for dc, count2 in donor_constituents.items(): - if count2 == N_H and N_H != 1: - degeneracy = 0 - for num in range(shellCount[dc] - (N_H - 1), 0, -1): - degeneracy += num - if N_H == 1: - degeneracy = shellCount[dc] - if count2 != N_H and N_H != 1: - degeneracy *= shellCount[dc] - - if degeneracy == 0: - degeneracy = 1 - pij = (float(count) / float(num_orients)) / float(degeneracy) - pij_degen_dict[donors] = [pij, degeneracy, sum_weight] - - verbosePrint("\t" * 5, "sum pair counts:", num_orients) - - solute_pijs = 0 - largest_pij = 0 - for donors, pij in pij_degen_dict.items(): - verbosePrint( - "\t" * 6, - "pij", - donors, - round(pij[0], 4), - "degen:", - pij[1], - "pij_d:", - round((pij[0] * pij[1]), 4), - ) - solute_pij, solvent_pij = False, False - for d in donors: - try: - d = int(d) - if d == 0: - solute_pij = True - else: - continue - except ValueError: - if d in ["X", "H"]: - solvent_pij = True - else: - solute_pij = True - if solute_pij is True: - solute_pijs += pij[0] * pij[1] - - if pij[0] * pij[1] > largest_pij: - largest_pij = pij[0] * pij[1] - - # calculate number of observed orientations - pij_max = max(map(lambda x: x[0], pij_degen_dict.values())) - verbosePrint("\t" * 6, "pij_max", round(pij_max, 5)) - - por_list = [] - omega_obs = 0 - for donors, pij in pij_degen_dict.items(): - om = (float(pij[0]) / float(pij_max)) * pij[1] - por_list.append((donors, om)) - omega_obs += om - - if AD == "D": - soluteD_pi = solute_pis - if AD == "A": - soluteA_pi = solute_pis - - verbosePrint("\t" * 6, "solute_donor", solute_donor) - verbosePrint("\t" * 6, "Nc, Nc_solvent_only", Nc, Nc_solvent_only) - verbosePrint("\t" * 6, "solute_pijs", round(solute_pijs, 3)) - verbosePrint("\t" * 6, "solute_pis", round(solute_pis, 3)) - verbosePrint("\t" * 6, "largest_pij", round(largest_pij, 3)) - verbosePrint("\t" * 6, "largest_pi", round(largest_pi, 3)) - - S = 0 - ppD_ppA_dict = {} - pbias_ave = 0 - Nc_eff = 0 - if Nc != 0: - pbias_list = [] - for i, N_i in shellCount.items(): - pD_i, pA_i = 0, 0 - pD_count, pA_count, pbias = 0, 0, 0 - if i in pA_degen_dict: - pA_i = pA_degen_dict[i][0] * pA_degen_dict[i][1] - pA_count = pA_degen_dict[i][2] - - if i in pD_degen_dict: - pD_i = pD_degen_dict[i][0] * pD_degen_dict[i][1] - pD_count = pD_degen_dict[i][2] - - verbosePrint( - "\t" * 5, - "i", - i, - "N_i", - N_i, - "pA_count", - round(pA_count, 5), - "pD_count", - round(pD_count, 5), - ) - if pD_count != 0 and pA_count != 0: - ppD_i = pD_i / float(pD_i + pA_i) - ppA_i = pA_i / float(pD_i + pA_i) - N_i_eff = ppD_i * ppA_i / float(0.25) * N_i - Nc_eff += N_i_eff - pbias = ppD_i * ppA_i - for x in range(0, N_i): - pbias_list.append(pbias) - verbosePrint( - "\t" * 6, - "i", - i, - "N_i", - N_i, - "ppA_i", - round(ppA_i, 5), - "ppD_i", - round(ppD_i, 5), - "N_i_eff", - round(N_i_eff, 5), - "pbias_i", - round(pbias, 5), - ) - - ppD_ppA_dict[i] = [ppD_i, ppA_i] - else: - for x in range(0, N_i): - pbias_list.append(pbias) - if pD_count == 0 and pA_count != 0: - ppD_ppA_dict[i] = [0, 1] - if pD_count != 0 and pA_count == 0: - ppD_ppA_dict[i] = [1, 0] - if pD_count == 0 and pA_count == 0: - ppD_ppA_dict[i] = [0, 0] - - verbosePrint("\t" * 5, "Nc_eff", round(Nc_eff, 5)) - - pbias_ave = 0 - if len(pbias_list) != 0: - pbias_ave = sum(pbias_list) / float(len(pbias_list)) - - verbosePrint("\t" * 5, "pbias_ave", round(pbias_ave, 5)) - if Nc_eff != 0: - S = np.log((Nc_eff) ** (3 / float(2)) * np.pi**0.5 * pbias_ave / 2) - verbosePrint("\t" * 5, "S3 (test2)", round(S, 5)) - - if S < 0: - S = 0 - - return S, pbias_ave, Nc_eff, ppD_ppA_dict - - -def forceUnitConversion(forceUnits): - """ - Amber = kcal / mol/ Ang / m^2 - Gromacs = kJ / mol / nm / m^2 - """ - - # for amber - force_half = 0.5**2 - cal_J = 4184**2 - pmol_pmolec = 6.02214086e23**2 - A_m = 1e-10**2 - A = 1e-10 - gmol_kg = 6.02214086e26 - - # - if forceUnits != "kJ": - # top = force_half*cal_J*gmol_kg - top = cal_J * gmol_kg - bottom = pmol_pmolec * A_m - constant = float(top) / float(bottom) - # - - # - if forceUnits == "kJ": - # for gromacs - nm_m = 1e-9**2 - nm = 1e-9 - kJ_J = 1000**2 # sqrd as force is sqred - top = kJ_J * gmol_kg - bottom = pmol_pmolec * nm_m - constant = float(top) / float(bottom) - # constant = 1000 - # - - return constant - - -def SvibCalculations(totFrames, Aclass, temperature, forceUnits, verbosePrint): - """ - Get rotated forces^2 and torques^2 matrices from allMoleculeList and - group into types. - Then for each molecule type (Nc and extended RAD shell), calc S_vib. - - Units: - force: Lammps = Kcal/mole-Angstrom, Amber = kcal/mol-Angstrom, - SI = kg.m/s^2 = N - energy: Lammps = Kcal/mole, Amber = Kcal/mole, SI = kg.m^2/s^2 = J(/mol) - covert Kcal/mole to KJ/mole by * 4.184 - - cal to mol = * 4.184 - - 240.663461 = (1000/6.022E+23)/(1.38E-23*0.5) to get T from KE (kJ/mol) - """ - - constant = forceUnitConversion(forceUnits) - - for nearest, key2 in Aclass.RADshell_dict.items(): - verbosePrint("\nnearest", nearest) - for assigned, key3 in key2.items(): - verbosePrint("\tassigned", assigned) - - for RADshell_num, values in key3.items(): - RADshell_num = float(RADshell_num) - verbosePrint("\n\t\tRADshell_num", RADshell_num, "count", values[5]) - forces = values[0] - torques = values[1] - pe = values[2] - ke = values[3] - Temp = values[4] - count = values[5] - weight = values[6] - weightNorm = weight / float(totFrames) - countNorm = count / float(totFrames) - zpe = 0 - - Strans_qm, Srot_qm, pe_mean, ke_mean, zpe, T = 0, 0, 0, 0, 0, None - - # setting temp here - T = temperature - - try: - force_sqrd_SI = np.multiply(forces, constant) - Strans_qm, zpe_trans, eigenvalues = Svib_calc( - force_sqrd_SI, T, cov=False, out=False, outputEigenValues=True - ) - zpe = zpe_trans - verbosePrint( - "\n\t\t\tave forces:\n\t\t\t%s" - "\n\t\t\t%s\n\t\t\t%s" - % (force_sqrd_SI[0], force_sqrd_SI[1], force_sqrd_SI[2]), - "\n\t\t\tStrans", - Strans_qm, - ) - - verbosePrint("\n\t\t\teigenvalues:\n\t\t\t%s" % (eigenvalues)) - - torque_sqrd_SI = np.multiply(torques, constant) - Srot_qm, zpe_rot, eigenvalues = Svib_calc( - torque_sqrd_SI, T, cov=False, out=False, outputEigenValues=True - ) - zpe += zpe_rot - verbosePrint( - "\n\t\t\tave torques:\n\t\t\t%s" - "\n\t\t\t%s\n\t\t\t%s" - % (torque_sqrd_SI[0], torque_sqrd_SI[1], torque_sqrd_SI[2]), - "\n\t\t\tSrot", - Srot_qm, - ) - - verbosePrint("\n\t\t\teigenvalues:\n\t\t\t%s" % (eigenvalues)) - - except (TypeError, ValueError): # force is a zero - pass - - pe = pe * 4.184 # kcal to kJ - verbosePrint("\n\t\t\tave pe: ", pe) - ke = ke * 4.184 - verbosePrint("\n\t\t\tave ke: ", ke) - - if Strans_qm is not None: - Strans_qm = Strans_qm * 8.314 - if Strans_qm is None: - Strans_qm = 0 - if Srot_qm is not None: - Srot_qm = Srot_qm * 8.314 - if Srot_qm is None: - Srot_qm = 0 - - Aclass.allVariables_dict[nearest][assigned][RADshell_num]["Strans"] = ( - Strans_qm, - count, - ) - Aclass.allVariables_dict[nearest][assigned][RADshell_num]["Srot"] = ( - Srot_qm, - count, - ) - if pe != 0 and ke != 0: - Aclass.allVariables_dict[nearest][assigned][RADshell_num]["PE"] = ( - pe, - count, - ) - Aclass.allVariables_dict[nearest][assigned][RADshell_num]["KE"] = ( - ke, - count, - ) - Aclass.allVariables_dict[nearest][assigned][RADshell_num]["count"] = ( - countNorm, - count, - ) - - -def Svib_calc(matrix_SI, T, cov, out, *args, **kwargs): - """ - Calculate Svib from matrix with SI units. - """ - - outputEigenValues = kwargs.get("outputEigenValues", False) - eigenvalues = kwargs.get("eigenvalues", False) - - kB = 1.38064852e-23 - # if T == None: - # T = 298 #over-write inputted T for now - T = float(T) # uses input T - h = 1.0545718e-34 # value of hbar - # h = 6.62607004e-34 - c = 30000000000 - NA = 6.022e23 - # print(T) - - # verbosePrint(matrix_SI) - - if eigenvalues is False: - if cov is True: - eigenvalues, eigenvectors = LA.eig(matrix_SI) - if cov is False: - eigenvalues = matrix_SI.diagonal() # do this for diagonal only - # CHECK - - if out is True: - if cov is True: - print( - "\n\t\t\teigenvectors:\n\t\t\t%s\n\t\t\t%s\n\t\t\t%s" - % (eigenvectors[0], eigenvectors[1], eigenvectors[2]), - "\n\t\t\teigenvalues:\n\t\t\t%s" % (eigenvalues), - ) - if cov is False: - print("\n\t\t\teigenvalues:\n\t\t\t%s" % (eigenvalues)) - - frequencies = [] - fraction = [] - sums = [] - ZPE = [] - for value in eigenvalues: - w = value / (kB * T) - # print('w = ', (w)) - # Need to avoid square root of negative number -> NaN results - # If w is positive no problems - if w > 0: - w = w**0.5 - # If w is very close to zero but on the negative side, we think this is noise - # and using zero should be fine. - elif w > -0.0001: - w = 0 - else: - raise ValueError("Cannot take square root of negative number") - frequencies.append(w) - zpe = 0.5 * w * h * NA - ZPE.append(zpe) - - for freq in frequencies: - # this is for qm - interim = (h * freq) / (kB * T) - fraction.append(interim) - - for value in fraction: - # verbosePrint(value) - if value != 0: - try: - mid = math.exp(value) - 1 - mid2 = math.log(1 - math.exp(-value)) - mid3 = value / mid - mid2 - sums.append(mid3) - except (ValueError, TypeError, OverflowError): - sums.append(0) - else: - sums.append(0) - - S_qm = sum(sums) - ZPE_qm = sum(ZPE) - # Svib_qm = NA * kB * S_qm - Svib_qm = S_qm - - if outputEigenValues is True: - return Svib_qm, ZPE_qm, eigenvalues - else: - return Svib_qm, ZPE_qm - - -def WMSvibCalculations(Aclass, totFrames, temperature, forceUnits, verbosePrint): - """ - Get forces and torques for molecule and UA levels with - shells as x - """ - - for nearest, resname_Ndih_key in sorted(list(Aclass.WM_FT_shell_dict.items())): - if "Any" in nearest: - verbosePrint("\n\n\n***** NEAREST: %s *****" % (nearest)) - for resname_Ndih, shellNum_key in sorted(list(resname_Ndih_key.items())): - resname = resname_Ndih[0] - for RADshellDist, values in sorted(list(shellNum_key.items())): - verbosePrint( - "\n\n\n***** SHELL_NUM or DIST: %s *****" % (RADshellDist) - ) - RADshellDist = float(RADshellDist) - count = values[7] - # countNorm = count / float(totFrames) - Svib_qm, Strans_qm, Srot_qm, Strans_qm_UA, Srot_qm_UA, T_mean = ( - process_FT( - Aclass, - nearest, - resname_Ndih, - values, - temperature, - forceUnits, - verbosePrint, - ) - ) - - Aclass.allVariables_dict[nearest][resname][RADshellDist][ - "WM_Strans" - ] = (Strans_qm * 8.314, count) - Aclass.allVariables_dict[nearest][resname][RADshellDist][ - "WM_Srot" - ] = (Srot_qm * 8.314, count) - Aclass.allVariables_dict[nearest][resname][RADshellDist][ - "WM_UA_Strans" - ] = (Strans_qm_UA * 8.314, count) - Aclass.allVariables_dict[nearest][resname][RADshellDist][ - "WM_UA_Srot" - ] = (Srot_qm_UA * 8.314, count) - - -def process_FT( - Aclass, nearest, resname_Ndih, values, temperature, forceUnits, verbosePrint -): - """ - Get forces and torques for molecule and UA levels - global - """ - - WM_eigenvalues = [] - WM_UA_eigenvalues = [] - - resname = resname_Ndih[0][0] - Ndih = resname_Ndih[1] - forces = values[0] - torques = values[1] - pe = values[2] - ke = values[3] - T = values[4] - UAs_forces = values[5] - UAs_torques = values[6] - count = values[7] - - verbosePrint("\n\n\t***** MOLECULE LEVEL: %s *****" % (resname), "count:", count) - verbosePrint("\t\t--WM", resname) - - Strans_qm, Srot_qm, pe_mean, ke_mean, zpe, T_mean, eigenvalues = FT_S( - forces, torques, pe, ke, forceUnits, temperature, verbosePrint, print_out=False - ) - verbosePrint("\n\t\t\tTmean", T_mean, "\n") - - for ev in eigenvalues: - e = ev + [resname] - WM_eigenvalues.append(e) - - if not isinstance(UAs_forces, int) and not isinstance(UAs_torques, int): - verbosePrint("\t\t--WM_UAs", resname) - Strans_qm, Srot_qm, pe_mean, ke_mean, zpe, T_mean, eigenvalues = FT_S( - UAs_forces, - UAs_torques, - pe, - ke, - forceUnits, - temperature, - verbosePrint, - print_out=False, - ) - verbosePrint("\n\t\t\tTmean", T_mean, "\n") - - for ev in eigenvalues: - e = ev + [resname] - WM_UA_eigenvalues.append(e) - - # if T_mean != None: - # temperature = T_mean - - verbosePrint("\t" * 2, "--WM eigenvalues") - WM_force_eigenvalues = [] - WM_torque_eigenvalues = [] - WM_eigenvalues.sort(key=lambda x: x[0]) - for ev in WM_eigenvalues: - verbosePrint("\t" * 3, ev) - if ev[1] == "force": - WM_force_eigenvalues.append(ev[0]) - elif ev[1] == "torque": - WM_torque_eigenvalues.append(ev[0]) - else: - continue - - WM_eigenvalues1 = sorted(list(map(lambda x: x[0], WM_eigenvalues))) - Svib_qm, zpe_trans = Svib_calc( - None, temperature, cov=True, out=False, eigenvalues=WM_eigenvalues1 - ) - verbosePrint("\n\t\t\tSvib WM", Svib_qm) - - Strans_qm, zpe_trans = Svib_calc( - None, temperature, cov=True, out=False, eigenvalues=WM_force_eigenvalues - ) - verbosePrint("\n\t\t\tStrans WM", Strans_qm) - - Srot_qm, zpe_trans = Svib_calc( - None, temperature, cov=True, out=False, eigenvalues=WM_torque_eigenvalues - ) - verbosePrint("\n\t\t\tSrot WM", Srot_qm) - - verbosePrint() - - verbosePrint("\t" * 2, "--WM_UAs eigenvalues") - WM_UA_eigenvalues.sort(key=lambda x: x[0]) - for ev in WM_UA_eigenvalues: - verbosePrint("\t" * 3, ev) - - verbosePrint() - verbosePrint("\t***** WM_UA LEVEL: %s *****" % (resname)) - - force_eigenvalues = [] - torque_eigenvalues = [] - for ev in WM_UA_eigenvalues: - if ev[1] == "force": - force_eigenvalues.append(ev[0]) - elif ev[1] == "torque": - torque_eigenvalues.append(ev[0]) - else: - continue - - torque_eigenvalues = sorted(torque_eigenvalues) - force_eigenvalues = sorted(force_eigenvalues) - - verbosePrint("\n\t\t--All WM_UA torque eigenvalues:") - for ev in torque_eigenvalues: - verbosePrint("\t" * 3, ev) - - verbosePrint("\n\t\t--All WM_UA force eigenvalues:") - for ev in force_eigenvalues: - verbosePrint("\t" * 3, ev) - - force_eigenvalues = force_eigenvalues[6:] - # remove 6 smallest eigenvals - # these correspond to WM level evs - if Ndih is not None: # halve Ndih largest remaining eigenvals - verbosePrint("\n\t\tNdih: %s" % (Ndih)) - N_ev = np.array(force_eigenvalues[0:Ndih]) - # lowest Ndih evs are 1/2 - N_ev_halved = [] - for ev_type in N_ev: - ev_half = float(ev_type) / float(4) - N_ev_halved.append(ev_half) - - force_eigenvalues = N_ev_halved + force_eigenvalues[Ndih:] - - verbosePrint("\n\t\t--(-6) WM_UA torque eigenvalues:") - for ev in torque_eigenvalues: - verbosePrint("\t" * 3, ev) - - verbosePrint("\n\t\t--(-6) WM_UA forces eigenvalues") - for ev in force_eigenvalues: - verbosePrint("\t" * 3, ev) - - # WM_UA_eigenvalues = sorted(list(map(lambda x: x[0], - # WM_UA_eigenvalues))) - Strans_qm_UA, zpe_trans_UA = Svib_calc( - None, temperature, cov=True, out=False, eigenvalues=force_eigenvalues - ) - - verbosePrint("\n\t\t\tStrans WM_UA", Strans_qm_UA) - - Srot_qm_UA, zpe_trans_UA = Svib_calc( - None, temperature, cov=True, out=False, eigenvalues=torque_eigenvalues - ) - - verbosePrint("\n\t\t\tSrot WM_UA", Srot_qm_UA) - - # - - return Svib_qm, Strans_qm, Srot_qm, Strans_qm_UA, Srot_qm_UA, T_mean - - -def FT_S(forces, torques, pe, ke, forceUnits, temperature, verbosePrint, print_out): - """ """ - - constant = forceUnitConversion(forceUnits) - forces = np.array(forces) - torques = np.array(torques) - pe = np.array(pe) - ke = np.array(ke) - - Strans_qm, Srot_qm, pe_mean, ke_mean, zpe, T_mean = None, None, None, None, 0, None - - # setting temp here - T_mean = temperature - - DOF = 0 - if not isinstance(forces, int): - for f in forces[0]: - if f != 0: - DOF += 1 - else: - continue - - if not isinstance(torques, int): - for t in torques[0]: - if t != 0: - DOF += 1 - else: - continue - - eigenvalues_list = [] - - if not isinstance(forces, int): - force_sqrd_SI = np.multiply(forces, constant) - Strans_qm, zpe_trans, eigenvalues = Svib_calc( - force_sqrd_SI, T_mean, cov=True, out=False, outputEigenValues=True - ) - zpe = zpe_trans - if print_out is True: - verbosePrint("\n\t\t\tave forces:") - for FF in force_sqrd_SI: - verbosePrint("\n\t\t\t%s" % (FF)) - - verbosePrint("\n\t\t\tStrans", Strans_qm) - - for ev in eigenvalues: - eigenvalues_list.append([ev, "force"]) - - if not isinstance(torques, int): - torque_sqrd_SI = np.multiply(torques, constant) - Srot_qm, zpe_rot, eigenvalues = Svib_calc( - torque_sqrd_SI, T_mean, cov=True, out=False, outputEigenValues=True - ) - zpe += zpe_rot - if print_out is True: - verbosePrint("\n\t\t\tave torques:") - for TT in torque_sqrd_SI: - verbosePrint("\n\t\t\t%s" % (TT)) - verbosePrint("\n\t\t\tSrot", Srot_qm) - - for ev in eigenvalues: - eigenvalues_list.append([ev, "torque"]) - - pe_mean = np.mean(pe) * 4.184 - verbosePrint("\n\t\t\tave pe: ", pe_mean) - - ke_mean = np.mean(ke) * 4.184 - verbosePrint("\n\t\t\tave ke: ", ke_mean) - if DOF != 0: - T_mean = float(ke_mean) / float(DOF) * 240.663461 - else: - T_mean = float(ke_mean) * 240.663461 - - return Strans_qm, Srot_qm, pe_mean, ke_mean, zpe, T_mean, eigenvalues_list - - -def process_dihedrals(Aclass, verbosePrint): - """ """ - - if len(Aclass.adaptive_dih_dict) != 0: - - # adaptive dihedrals - verbosePrint("\nAdaptive Dihedrals p ln p\n") - for nearest, assigned_key in sorted(list(Aclass.adaptive_dih_dict.items())): - for assigned, shell_num_key in sorted(list(assigned_key.items())): - # assigned_count = 0 - for shell_num, dih_atoms_key in sorted(list(shell_num_key.items())): - verbosePrint(nearest, assigned, shell_num) - sum_S = 0 - shell_count = 0 - dih_count = 0 - for dih_atoms, phi_key in sorted(list(dih_atoms_key.items())): - verbosePrint(dih_atoms) - dih_count += 1 - phi_count_list = [] - max_list = [] - for phi, count in sorted(list(phi_key.items())): - phi_count_list.append([phi, count]) - shell_count += count / float(len(dih_atoms_key.keys())) - N = len(phi_count_list) - 1 - # verbosePrint(phi_count_list) - same_dict = nested_dict() - for x in range(0, N + 1): - if phi_count_list[x][1] != 0: - if x == 0: - if ( - phi_count_list[x + 1][1] <= phi_count_list[x][1] - and phi_count_list[N][1] <= phi_count_list[x][1] - ): - max_list.append(phi_count_list[x]) - else: - continue - elif x == N: - if ( - phi_count_list[x - 1][1] <= phi_count_list[x][1] - and phi_count_list[0][1] <= phi_count_list[x][1] - ): - max_list.append(phi_count_list[x]) - else: - continue - else: - if ( - phi_count_list[x - 1][1] <= phi_count_list[x][1] - and phi_count_list[x + 1][1] - <= phi_count_list[x][1] - ): - max_list.append(phi_count_list[x]) - else: - continue - else: - continue - - refined_max_list = [] - excluded = [] - - for max_count in max_list: - for max_count2 in max_list: - if ( - max_count != max_count2 - and len(max_list) > 1 - and max_count not in refined_max_list - ): - diff = max_count[0] - max_count2[0] - - normDeg = diff % 360 - diff2 = min(360 - normDeg, normDeg) - # if less than 30 degrees difference - # then only pick the largest degree to be - # binned the other gets excluded - if diff2 < 60: - if max_count[1] == max_count2[1]: - if ( - max_count[0] > max_count2[0] - and max_count not in refined_max_list - ): - refined_max_list.append(max_count) - excluded.append(max_count2) - else: - continue - elif max_count[1] > max_count2[1]: - if max_count not in refined_max_list: - refined_max_list.append(max_count) - else: - continue - else: - continue - elif ( - max_count == max_count2 - and len(max_list) == 1 - and max_count not in refined_max_list - ): - refined_max_list.append(max_count) - else: - continue - - for max_count in max_list: - for max_count2 in max_list: - if ( - max_count != max_count2 - and len(max_list) > 1 - and max_count not in refined_max_list - ): - diff = max_count[0] - max_count2[0] - normDeg = diff % 360 - diff2 = min(360 - normDeg, normDeg) - if ( - diff2 >= 60 - and max_count not in refined_max_list - and max_count not in excluded - ): - refined_max_list.append(max_count) - - refined_max_dict = {} - for max_count in refined_max_list: - refined_max_dict[max_count[0]] = max_count[1] - max_phis = list(refined_max_dict.keys()) - for phi, count in phi_count_list: - if phi not in refined_max_dict and count != 0: - smallest_diff = 999 - closest_phi = None - for maxPhi in max_phis: - diff = phi - maxPhi - normDeg = diff % 360 - diff2 = min(360 - normDeg, normDeg) - if diff2 < smallest_diff: - smallest_diff = diff2 - closest_phi = maxPhi - refined_max_dict[closest_phi] += count - - # p ln p calcs finally! - sum_count = sum( - count for phi, count in refined_max_dict.items() - ) - topo_S = 0 - for phi, count in refined_max_dict.items(): - p = count / float(sum_count) - S = -p * np.log(p) * 8.314 - topo_S += S - verbosePrint("\t", phi, count, p, S) - verbosePrint("\t\t", topo_S) - verbosePrint() - sum_S += topo_S - - Aclass.allVariables_dict[nearest][ - (assigned[0], assigned[3], assigned[4]) - ][float(shell_num)]["conf_AE"] = sum_S, int(round(shell_count, 0)) - verbosePrint("\nS_adaptive", sum_S, "N_dih:", dih_count, assigned) - verbosePrint() - - -def saveAllVariables( - num_frames, Aclass, level, name, solvent, waterTuple, verbosePrint -): - """ """ - - frames = "" - if num_frames is not None: - frames = num_frames - solventDataName = f"solventVariables{frames}{name}_{level}.csv" - # data = open('solventVariables%s%s_%s.csv' % (frames, name, level), 'w') - # data.write( - # '\n'.join(['nearest,assigned,shell_num,variable,value,count']) - # + '\n') - Utils.printOut(solventDataName, "nearest,assigned,shell_num,variable,value,count") - solventData = pd.DataFrame( - columns=["nearest", "assigned", "shell_num", "variable", "value", "count"] - ) - - soluteDataName = f"soluteVariables{frames}{name}_{level}.csv" - # data2 = open('soluteVariables%s%s_%s.csv' % - # (frames, name, level), 'w') - # data2.write('\n'.join(['resName,variable,value,count']) - # + '\n') - Utils.printOut(soluteDataName, "resName,variable,value,count") - - soluteData = pd.DataFrame(columns=["resName", "variable", "value", "count"]) - - if level == "res_atomLevel": - reduced_dict = nested_dict() - for nearest, assigned_key in sorted(list(Aclass.allVariables_dict.items())): - nearest = nearest.split("_")[0] - for assigned, shellNum_key in sorted(list(assigned_key.items())): - assigned = assigned[0].split("_")[0] - if assigned in waterTuple: - for shellNum, variable_key in sorted(list(shellNum_key.items())): - tot_count = 0 - for nearest2, assigned2_key in Aclass.allVariables_dict.items(): - for assigned2, shellNum2_key in assigned2_key.items(): - for shellNum2, variable2_key in shellNum2_key.items(): - for variable2, count in variable2_key.items(): - nearest3 = nearest2.split("_")[0] - assigned3 = assigned2[0].split("_")[0] - if ( - nearest3 == nearest - and assigned3 == assigned - and shellNum2 == shellNum - and variable2 == "count" - ): - tot_count += count[1] - else: - continue - for variable, value in variable_key.items(): - if ( - variable - not in reduced_dict[nearest][assigned][shellNum] - and value[0] is not None - ): - reduced_dict[nearest][assigned][shellNum][variable] = [ - 0, - 0, - ] - reduced_dict[nearest][assigned][shellNum][variable][0] += ( - value[0] * value[1] / float(tot_count) - ) - reduced_dict[nearest][assigned][shellNum][variable][ - 1 - ] += int(round(value[1], 0)) - else: - continue - - for nearest, assigned_key in sorted(list(reduced_dict.items())): - for assigned, shellNum_key in sorted(list(assigned_key.items())): - for shellNum, variable_key in sorted(list(shellNum_key.items())): - for variable, value in variable_key.items(): - if ( - value[0] is not None - and assigned in solvent - and shellNum == 1 - ): - # data.write('\n'.join(['%s,%s,%s,%s,%s,%s' % - # (nearest, - # assigned, shellNum, variable, value[0], - # value[1])]) + '\n') - Utils.printOut( - solventDataName, - ( - f"{nearest},{assigned}," - f"{shellNum},{variable}," - f"{value[0]},{value[1]}" - ), - ) - newRowSolvent = pd.DataFrame( - { - "nearest": nearest, - "assigned": assigned, - "shell_num": shellNum, - "variable": variable, - "value": value[0], - "count": value[1], - }, - index=[0], - ) - solventData = pd.concat( - [solventData, newRowSolvent], ignore_index=True - ) - else: - continue - - # for solute only - PE_KE_dict = nested_dict() - for nearest, assigned_key in sorted(list(Aclass.allVariables_dict.items())): - for assigned, shellNum_key in sorted(list(assigned_key.items())): - assigned = list(assigned) - if "_" in assigned[0]: - assigned[0] = assigned[0].split("_")[0] - for shellNum, variable_key in sorted(list(shellNum_key.items())): - for variable, value in variable_key.items(): - if assigned[0] not in solvent and shellNum != 0: - if variable in ["PE", "KE"] and value[0] is not None: - new_var = "WM_%s" % (variable) - if new_var not in PE_KE_dict[nearest][assigned[0]]: - PE_KE_dict[nearest][assigned[0]][new_var] = [0, 0] - PE_KE_dict[nearest][assigned[0]][new_var][0] += value[0] - PE_KE_dict[nearest][assigned[0]][new_var][1] = value[1] - - if "WM_" in variable or variable == "conf_AE": - # data2.write('\n'.join(['%s,%s,%s,%s' % - # (assigned[0], - # variable, value[0], - # int(round(value[1], 0)))]) + '\n') - Utils.printOut( - soluteDataName, - ( - f"{assigned[0]},{variable},{value[0]}," - f"{int(round(value[1], 0))}" - ), - ) - newRowSolute = pd.DataFrame( - { - "resName": assigned[0], - "variable": variable, - "value": value[0], - "count": int(round(value[1], 0)), - }, - index=[0], - ) - soluteData = pd.concat( - [soluteData, newRowSolute], ignore_index=True - ) - else: - continue - - for nearest, assigned_key in sorted(list(PE_KE_dict.items())): - for assigned, variable_key in sorted(list(assigned_key.items())): - for variable, value in sorted(list(variable_key.items())): - # data2.write('\n'.join(['%s,%s,%s,%s' % - # (assigned, - # variable, value[0], - # value[1])]) + '\n') - Utils.printOut( - soluteDataName, - f"{assigned[0]},{variable},{value[0]},{value[1]}", - ) - newRowSolute = pd.DataFrame( - { - "resName": assigned, - "variable": variable, - "value": value[0], - "count": value[1], - }, - index=[0], - ) - soluteData = pd.concat( - [soluteData, newRowSolute], ignore_index=True - ) - # - - else: - for nearest, assigned_key in sorted(list(Aclass.allVariables_dict.items())): - for assigned, shellNum_key in sorted(list(assigned_key.items())): - for shellNum, variable_key in sorted(list(shellNum_key.items())): - for variable, value in variable_key.items(): - if ( - assigned[0] in solvent - and shellNum == 1 - and value[0] is not None - ): - # data.write('\n'.join(['%s,%s,%s,%s,%s,%s' % - # (nearest, - # assigned[0], shellNum, - # variable, value[0], - # int(round(value[1], 0)))]) + '\n') - Utils.printOut( - solventDataName, - ( - f"{nearest},{assigned[0]}," - f"{shellNum},{variable},", - f"{value[0]},{int(round(value[1], 0))}", - ), - ) - newRowSolvent = pd.DataFrame( - { - "nearest": nearest, - "assigned": assigned[0], - "shell_num": shellNum, - "variable": variable, - "value": value[0], - "count": int(round(value[1], 0)), - }, - index=[0], - ) - solventData = pd.concat( - [solventData, newRowSolvent], ignore_index=True - ) - else: - continue - - # - # for solute only - PE_KE_dict = nested_dict() - for nearest, assigned_key in sorted(list(Aclass.allVariables_dict.items())): - for assigned, shellNum_key in sorted(list(assigned_key.items())): - for shellNum, variable_key in sorted(list(shellNum_key.items())): - for variable, value in variable_key.items(): - if assigned[0] not in solvent and shellNum != 0: - if variable in ["PE", "KE"] and value[0] is not None: - new_var = "WM_%s" % (variable) - if new_var not in PE_KE_dict[nearest][assigned[0]]: - PE_KE_dict[nearest][assigned[0]][new_var] = [0, 0] - PE_KE_dict[nearest][assigned[0]][new_var][0] += value[0] - PE_KE_dict[nearest][assigned[0]][new_var][1] = value[1] - - if "WM_" in variable or variable == "conf_AE": - # data2.write('\n'.join(['%s,%s,%s,%s' % - # (assigned[0], - # variable, value[0], - # int(round(value[1], 0)))]) + '\n') - Utils.printOut( - soluteDataName, - ( - f"{assigned[0]},{variable}," - f"{value[0]},{int(round(value[1], 0))}" - ), - ) - newRowSolute = pd.DataFrame( - { - "resName": assigned[0], - "variable": variable, - "value": value[0], - "count": int(round(value[1], 0)), - }, - index=[0], - ) - soluteData = pd.concat( - [soluteData, newRowSolute], ignore_index=True - ) - - else: - continue - - for nearest, assigned_key in sorted(list(PE_KE_dict.items())): - for assigned, variable_key in sorted(list(assigned_key.items())): - for variable, value in sorted(list(variable_key.items())): - # data2.write('\n'.join(['%s,%s,%s,%s' % - # (assigned, - # variable, value[0], - # value[1])]) + '\n') - Utils.printOut( - soluteDataName, f"{assigned},{variable},{value[0]},{value[1]}" - ) - newRowSolute = pd.DataFrame( - { - "resName": assigned[0], - "variable": variable, - "value": value[0], - "count": int(round(value[1], 0)), - }, - index=[0], - ) - soluteData = pd.concat( - [soluteData, newRowSolute], ignore_index=True - ) - # - - # data.close() - # data2.close() - return (solventData, soluteData) diff --git a/CodeEntropy/poseidon/analysis/EntropyFunctions.py b/CodeEntropy/poseidon/analysis/EntropyFunctions.py deleted file mode 100644 index 96dc5fb..0000000 --- a/CodeEntropy/poseidon/analysis/EntropyFunctions.py +++ /dev/null @@ -1,55 +0,0 @@ -#!/usr/bin/env python3 - -# import os -import sys - -import numpy as np -from numpy import linalg as la - -# from ast import arg - - -np.set_printoptions(threshold=sys.maxsize) -# from CodeEntropy.ClassCollection import BeadClasses as BC -# from CodeEntropy.ClassCollection import ConformationEntity as CONF -# from CodeEntropy.ClassCollection import ModeClasses -# from CodeEntropy.ClassCollection import CustomDataTypes -# from CodeEntropy.FunctionCollection import CustomFunctions as CF -# from CodeEntropy.FunctionCollection import GeometricFunctions as GF -# from CodeEntropy.FunctionCollection import UnitsAndConversions as UAC -# from CodeEntropy.FunctionCollection import Utils -# from CodeEntropy.IO import Writer -# import multiprocessing as mp -# from functools import partial - -# import MDAnalysis as mda -# import matplotlib.pyplot as plt -# import pandas as pd - -# def frequencies(lambdas): - to add!!! - - -def vibrational_entropies(FTmatrix): - lambdas = la.eig(FTmatrix) - exponent = 1 - exp_positive = np.power(np.e, exponent) - exp_negative = np.power(np.e, -exponent) - print(lambdas) - - -test_matrix = [ - [7, 0, 1, 1, 2, 2], - [5, 5, 0, 0, 4, 5], - [0, 1, 0, 0, 5, 4], - [2, 0, 0, 0, 0, 2], - [ - 1, - 1, - 1, - 1, - 1, - 1, - ], - [0, 0, 0, 0, 0, 1], -] -vibrational_entropies(test_matrix) diff --git a/CodeEntropy/poseidon/analysis/__init__.py b/CodeEntropy/poseidon/analysis/__init__.py deleted file mode 100644 index 201e256..0000000 --- a/CodeEntropy/poseidon/analysis/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -""" -Empty init file -""" diff --git a/CodeEntropy/poseidon/analysis/classes.py b/CodeEntropy/poseidon/analysis/classes.py deleted file mode 100644 index e4bc81d..0000000 --- a/CodeEntropy/poseidon/analysis/classes.py +++ /dev/null @@ -1,45 +0,0 @@ -#!/usr/bin/env python - - -from collections import defaultdict - - -# create nested dict in one go -def nested_dict(): - return defaultdict(nested_dict) - - -class variables_dicts(object): - """ - Main class containing dicts for all - thermodynamic (Energy and Entropy) calculations per proximity - shell in the system. - """ - - def __init__(self, initiate): - - self.initiate = None - - # dict with variables named rather than ordered - self.allVariables_dict = nested_dict() # * - - # get info about weight of bias in biased simulations - self.weighting_list = None - self.weighting_chunks = None - - # For Sor population - self.RADshell_Sor_dict = nested_dict() # for S_or, pijs - self.Sor_reference_dict = nested_dict() - - # For Svib / E population - self.RADshell_dict = nested_dict() - # FT decomposed to prox or dist - self.WM_FT_shell_dict = nested_dict() - - # dicts for dihedrals - self.initiate = None - # bin dihs into 30 degrees bin, uses p ln p (not matrix) - self.adaptive_dih_dict = nested_dict() - - # contacts between resids with central molecule resid using RAD - self.resid_contact_matrix_dict = nested_dict() diff --git a/CodeEntropy/poseidon/analysis/helper.py b/CodeEntropy/poseidon/analysis/helper.py deleted file mode 100644 index 2da7d38..0000000 --- a/CodeEntropy/poseidon/analysis/helper.py +++ /dev/null @@ -1,40 +0,0 @@ -import os - -import psutil - - -def memoryInfo(verbosePrint): - """ - Prints the current process's memory usage in gigabytes. - - Args: - verbosePrint (bool): This parameter is unused in the function. - - Returns: - None - """ - process = psutil.Process(os.getpid()) - bytes_info = float(process.memory_info().rss) - gb = bytes_info * 1e-9 - print("memory use:", round(gb, 3), "GB") # in gbytes - - -def weightingPopulation(weighting): - """ - Save a list of weightings per frame from biased simulations in the - EEclass - """ - dataFile = weighting[0] - try: - frame_chunks = weighting[1] - except IndexError: - frame_chunks = 1 - weighting_list = [] - with open(dataFile) as data: - for line in data: - weighting_list.append(line) - data.close() - - weighting_info = [weighting_list, frame_chunks] - - return weighting_info diff --git a/CodeEntropy/poseidon/analysis/populateClasses.py b/CodeEntropy/poseidon/analysis/populateClasses.py deleted file mode 100644 index 3081355..0000000 --- a/CodeEntropy/poseidon/analysis/populateClasses.py +++ /dev/null @@ -1,527 +0,0 @@ -#!/usr/bin/env python - -import math - -# import os -import sys - -from CodeEntropy.poseidon.analysis.classes import variables_dicts -from CodeEntropy.poseidon.analysis.populateEE import ( - SconfPopulation, - SorPopulation, - SvibPopulation, - contactPopulation, - contactPopulationUA, -) - -# from collections import Counter - -# import numpy as np -# import psutil - - -def classPopulation( - atomList, - entropyEnergy, - level_list, - count, - atom_count, - num_atoms, - waterTuple, - temperature, - solvent, - EEclass, - EEclass_residLevel_resname, - EEclass_soluteContacts, - EEclass_atomLevel, - weighting_info, - verbosePrint, -): - """ - For each flag, populate dictionaries accordingly - The objects read in are used for this. - """ - - # Smix_methods - step_list = None - - # - if count == 1: - - # initiate classes for each flag - - if entropyEnergy: - EEclass = variables_dicts(True) - EEclass_residLevel_resname = variables_dicts(True) - EEclass_soluteContacts = variables_dicts(True) - EEclass_atomLevel = variables_dicts(True) - - if weighting_info is not None: - if EEclass is None: - EEclass = variables_dicts(True) - EEclass.weighting_list = weighting_info[0] - EEclass.weighting_chunks = weighting_info[1] - - # - - for i in atomList: # iterate through each atom - # only need to change this if list order changes (and below) - atom_count += 1 - globalFrame = math.ceil(atom_count / float(num_atoms)) # round UP - frame = i[0] # not used - cumulative_steps = ["All"] - if step_list is not None: - for n in range(1, len(step_list)): - step = step_list[n] - if int(globalFrame) <= step: - cumulative_step = str(step) - # print(globalFrame, cumulative_step) - cumulative_steps.append(cumulative_step) - - atom_name = i[1] - atom_num = i[2] - atom_resname = i[3] - atom_resid = i[4] - N_H = i[5][1] - # molecule_atoms = i[6] #heavy atoms only - molecule_size = len(i[6]) # heavy atoms only - try: - nearest_atomName = i[7][0] - nearest_resname = i[7][1] - nearest_resid = i[7][2] - # dist = float(i[7][3]) - # dist = round(float(i[7][3]), 1) - except IndexError: - nearest_atomName = "unassigned" - nearest_resname = "unassigned" - nearest_resid = 0 - # dist = 0 - RADshell_num = i[8] - if RADshell_num is not None: - RADshell_num = "0" # unassigned to a nearest - RAD_list = sorted(i[9][0]) - As = i[9][1] - Ds = i[9][2] - # Svib_Ecalcs - MweightedForces = i[10] # array - MweightedTorques = i[11] # array - PE = None - KE = None - if i[12] is not None: - PE = i[12][0] - KE = i[12][1] - - UAweightedForces = i[13] - UAweightedTorques = i[14] - # FT POPULATION - try: - dihedral_phi_list = i[15] - Ndih = len(i[15]) - except (IndexError, TypeError): - dihedral_phi_list - Ndih = None - try: - molecule_UA_F = i[16] - molecule_UA_T = i[17] - except IndexError: - Ndih = None - molecule_UA_F = [] - molecule_UA_T = [] - RAD_nAtoms = i[18] # closest atom in each RADshell and distance - - RAD_str = tuple(sorted(RAD_list)) - - # toggle between shells or distances as x variable - RADshell_num_dist = str(RADshell_num) - - # group all solute types together for easier analysis - if atom_resname not in solvent and RADshell_num != "0": - nearest_resname = "Any" - RADshell_num_dist = "1" - - # find out weighting of each frame if simulation is biased - weight = 1 - if weighting_info is not None: - if EEclass.weighting_list is not None: - weight = getWeight(EEclass, globalFrame) - - for level in level_list: - - nearest_resname_list = [] - atom_resnameLevel, nearest_resname_list, waterTuple = levelSelection( - level, - atom_resname, - atom_resid, - nearest_resname, - nearest_resname_list, - atom_name, - nearest_atomName, - nearest_resid, - RAD_nAtoms, - waterTuple, - ) - - for nearest_resnameLevel in nearest_resname_list: - - if "Any" in nearest_resnameLevel: - nearest_resnameLevel = "Any" - - if level in [None, "moleculeLevel"]: - - if entropyEnergy: - - SorPopulation( - EEclass, - atom_name, - atom_resnameLevel, - nearest_resnameLevel, - RAD_str, - RADshell_num, - RADshell_num_dist, - waterTuple, - N_H, - As, - Ds, - weight, - ) - - SvibPopulation( - EEclass, - level, - atom_name, - atom_resnameLevel, - N_H, - atom_num, - nearest_resnameLevel, - RADshell_num_dist, - waterTuple, - UAweightedForces, - UAweightedTorques, - PE, - KE, - MweightedForces, - MweightedTorques, - Ndih, - molecule_UA_F, - molecule_UA_T, - weight, - ) - - SconfPopulation( - EEclass, - nearest_resnameLevel, - atom_name, - atom_resnameLevel, - N_H, - dihedral_phi_list, - molecule_size, - RADshell_num_dist, - weight, - ) - - if level in ["residLevel_resname"]: - if entropyEnergy: - - contactPopulation( - EEclass_residLevel_resname, - RAD_nAtoms, - atom_resnameLevel, - atom_resid, - waterTuple, - weight, - ) - - SorPopulation( - EEclass_residLevel_resname, - atom_name, - atom_resnameLevel, - nearest_resnameLevel, - RAD_str, - RADshell_num, - RADshell_num_dist, - waterTuple, - N_H, - As, - Ds, - weight, - ) - - SvibPopulation( - EEclass_residLevel_resname, - level, - atom_name, - atom_resnameLevel, - N_H, - atom_num, - nearest_resnameLevel, - RADshell_num_dist, - waterTuple, - UAweightedForces, - UAweightedTorques, - PE, - KE, - MweightedForces, - MweightedTorques, - Ndih, - molecule_UA_F, - molecule_UA_T, - weight, - ) - - SconfPopulation( - EEclass_residLevel_resname, - nearest_resnameLevel, - atom_name, - atom_resnameLevel, - N_H, - dihedral_phi_list, - molecule_size, - RADshell_num_dist, - weight, - ) - - if level in ["soluteContacts"]: - - if entropyEnergy: - SorPopulation( - EEclass_soluteContacts, - atom_name, - atom_resnameLevel, - nearest_resnameLevel, - RAD_str, - RADshell_num, - RADshell_num_dist, - waterTuple, - N_H, - As, - Ds, - weight, - ) - - SvibPopulation( - EEclass_soluteContacts, - level, - atom_name, - atom_resnameLevel, - N_H, - atom_num, - nearest_resnameLevel, - RADshell_num_dist, - waterTuple, - UAweightedForces, - UAweightedTorques, - PE, - KE, - MweightedForces, - MweightedTorques, - Ndih, - molecule_UA_F, - molecule_UA_T, - weight, - ) - - SconfPopulation( - EEclass_soluteContacts, - nearest_resnameLevel, - atom_name, - atom_resnameLevel, - N_H, - dihedral_phi_list, - molecule_size, - RADshell_num_dist, - weight, - ) - - if level in ["atomLevel"]: - - if entropyEnergy: - - contactPopulationUA( - EEclass_atomLevel, - RAD_nAtoms, - atom_name, - atom_resnameLevel, - atom_resid, - waterTuple, - weight, - ) - - SorPopulation( - EEclass_atomLevel, - atom_name, - atom_resnameLevel, - nearest_resnameLevel, - RAD_str, - RADshell_num, - RADshell_num_dist, - waterTuple, - N_H, - As, - Ds, - weight, - ) - - SvibPopulation( - EEclass_atomLevel, - level, - atom_name, - atom_resnameLevel, - N_H, - atom_num, - nearest_resnameLevel, - RADshell_num_dist, - waterTuple, - UAweightedForces, - UAweightedTorques, - PE, - KE, - MweightedForces, - MweightedTorques, - Ndih, - molecule_UA_F, - molecule_UA_T, - weight, - ) - - SconfPopulation( - EEclass_atomLevel, - nearest_resnameLevel, - atom_name, - atom_resnameLevel, - N_H, - dihedral_phi_list, - molecule_size, - RADshell_num_dist, - weight, - ) - - sys.stdout.flush() - - return ( - EEclass, - EEclass_residLevel_resname, - EEclass_soluteContacts, - EEclass_atomLevel, - atom_count, - ) - - -def getWeight(EEclass, globalFrame): - """ - Find out what the weighting is for a given simualtion frame - """ - - # chunk = int(EEclass.weighting_chunks) - weighting_list = EEclass.weighting_list - - index = globalFrame - 1 - try: - weight = float(weighting_list[index]) - # print('weight: %s' % (weight)) - except (IndexError, ValueError): - weight = 1 - - # print(globalFrame, index, weight) - - return weight - - -def stripNumbers(atom_name, atom_resname): - """ """ - - atom_name_letters = "".join([i for i in atom_name if not i.isdigit()]) - if atom_name != atom_resname: - atom_name_letters = atom_name_letters[0] - atom_resname = "%s_%s" % (atom_resname, atom_name_letters) - - return atom_resname - - -def levelSelection( - level, - atom_resname, - atom_resid, - nearest_resname, - nearest_resname_list, - atom_name, - nearest_atomName, - nearest_resid, - RAD_nAtoms, - waterTuple, -): - """ - Depending on what level of analysis we want to do, dictionaries are - populated according to how the nearest and assigned are defined. - - For example, we may want to analyse molecule types, the atoms - in each molecule, molecules by residue number or multiple solutes - in a RAD shell as a list of nearest molecules. - """ - - # keep as molecule - if level == "moleculeLevel" or level is None: - atom_resname = atom_resname - nearest_resname_list = [nearest_resname] - - # use resname_atomName for both nearest and assigned - if level == "atomLevel": # res_atomLevel - # - water = False - if atom_resname in waterTuple: - water = True - atom_resname = stripNumbers(atom_name, atom_resname) - if atom_resname not in waterTuple and water: - waterTuple.append(atom_resname) - nearest_resname_list = [stripNumbers(nearest_atomName, nearest_resname)] - - # - - """ - ###Inclusion of atom names, no edits - if level == 'atomLevel': - atom_resname = ('%s_%s' % (atom_resname, atom_name)) - nearest_resname_list = [('%s_%s' % (nearest_resname, - nearest_atomName))] - """ - - # Inclusion of molecule resid - if level == "residLevel_resname": - atom_resname = atom_resname - if atom_resname not in waterTuple: - atom_resname = "%s_%s" % (atom_resname, atom_resid) - nearest_resname_list = ["%s_%s" % (nearest_resname, nearest_resid)] - # print(atom_resname, nearest_resname_list) - - # Inclusion of molecule resid - if level == "residLevel_atom": - atom_resname = "%s_%s" % (atom_resname, atom_name) - nearest_resname_list = ["%s_%s" % (nearest_resname, nearest_resid)] - - if level == "soluteContacts": - atom_resname = atom_resname - if atom_resname in waterTuple: - if RAD_nAtoms is not None: - name_id_list = [] - for nInfo in RAD_nAtoms: - nResid = nInfo[1] - nResname = nInfo[2] - - if nResname not in waterTuple: - name_id = "%s_%s" % (nResname, nResid) - name_id_list.append(name_id) - - name_id_list = sorted(name_id_list) - - if len(name_id_list) > 1: - for solute in name_id_list: - list2 = name_id_list[:] - list2.remove(solute) - for solute2 in list2: - if ( - "%s_%s" % (solute, solute2) - ) not in nearest_resname_list: - nearest_resname_list.append( - ("%s_%s" % (solute, solute2)) - ) - - return atom_resname, nearest_resname_list, waterTuple diff --git a/CodeEntropy/poseidon/analysis/populateEE.py b/CodeEntropy/poseidon/analysis/populateEE.py deleted file mode 100644 index 00b381f..0000000 --- a/CodeEntropy/poseidon/analysis/populateEE.py +++ /dev/null @@ -1,517 +0,0 @@ -#!/usr/bin/env python - -# from itertools import combinations - -import numpy as np - - -def contactPopulation(Aclass, RAD_nAtoms, atom_resname, atom_resid, waterTuple, weight): - """ """ - - # RESID CONTACT MATRIX POP - # contacts defined using QcQn/r^2 between any atoms - # (H or heavy atoms) here! - if RAD_nAtoms is not None: - # RAD_closestRanked = [] - # nRAD_closestRanked_separated = {} - for nInfo in RAD_nAtoms: - nResid = nInfo[1] - nResname = nInfo[2] - if ( - nResname not in waterTuple - and atom_resname not in waterTuple - and nResid != atom_resid - ): - if (nResname, nResid) not in Aclass.resid_contact_matrix_dict[ - (atom_resname, atom_resid) - ]: - Aclass.resid_contact_matrix_dict[(atom_resname, atom_resid)][ - (nResname, nResid) - ] = 0 - Aclass.resid_contact_matrix_dict[(nResname, nResid)][ - (atom_resname, atom_resid) - ] = 0 - Aclass.resid_contact_matrix_dict[(atom_resname, atom_resid)][ - (nResname, nResid) - ] += weight # 1 - - else: - continue - - -def contactPopulationUA( - Aclass, RAD_nAtoms, atom_name, atom_resname, atom_resid, waterTuple, weight -): - """ """ - - # RESID CONTACT MATRIX POP - # contacts defined using QcQn/r^2 between any atoms - # (H or heavy atoms) here! - atom_resname = atom_resname.split("_")[0] - if RAD_nAtoms is not None: - # RAD_closestRanked = [] - # nRAD_closestRanked_separated = {} - for nInfo in RAD_nAtoms: - nAtom = nInfo[0] - nResid = nInfo[1] - nResname = nInfo[2] - if ( - nResname not in waterTuple - and atom_resname not in waterTuple - and nResid != atom_resid - ): - if (nAtom, nResname, nResid) not in Aclass.resid_contact_matrix_dict[ - (atom_name, atom_resname, atom_resid) - ]: - Aclass.resid_contact_matrix_dict[ - (atom_name, atom_resname, atom_resid) - ][(nAtom, nResname, nResid)] = 0 - Aclass.resid_contact_matrix_dict[(nAtom, nResname, nResid)][ - (atom_name, atom_resname, atom_resid) - ] = 0 - Aclass.resid_contact_matrix_dict[(atom_name, atom_resname, atom_resid)][ - (nAtom, nResname, nResid) - ] += weight # 1 - - else: - continue - - -def SorPopulation( - Aclass, - atom_name, - atom_resname, - nearest_resname, - RAD_str, - RADshell_num, - RADshell_num_dist, - waterTuple, - N_H, - As, - Ds, - weight, -): - - Ds = tuple(Ds) - As = tuple(As) - DA_shell = Ds + As - DA_shell = tuple(sorted(DA_shell)) - # print (DA_shell) - - for RR in [[RADshell_num_dist, str(RADshell_num)], ["0", "0"]]: - - RADshell_num_dist2 = RR[0] - RADshell_num2 = RR[1] - - # As and Ds as they are for pijs - for ADs in [["A", As], ["D", Ds]]: - if ( - ADs[1] - not in Aclass.RADshell_Sor_dict[nearest_resname][ - (atom_resname, N_H, atom_name) - ][(RADshell_num_dist2, str(RADshell_num2))][RAD_str][ADs[0]] - ): - - Aclass.RADshell_Sor_dict[nearest_resname][ - (atom_resname, N_H, atom_name) - ][(RADshell_num_dist2, str(RADshell_num2))][RAD_str][ADs[0]][ADs[1]] = 0 - - Aclass.RADshell_Sor_dict[nearest_resname][(atom_resname, N_H, atom_name)][ - (RADshell_num_dist2, str(RADshell_num2)) - ][RAD_str][ADs[0]][ - ADs[1] - ] += weight # 1 - - -def runningWeightedAverageFT( - weight, - weight_stored, - count_stored, - forces_stored, - torques_stored, - T_stored, - weighted_add, - count_add, - forces_add, - torques_add, - T_add, - force, - torque, - KE, -): - """ - Get the running average of force and torque matrices - """ - - try: - if ( - any(elem is None for elem in force) is False - and any(elem is None for elem in torque) is False - ): - try: - if torques_stored == 0 and forces_stored == 0: - weight_add = weight_stored + weight - count_add = count_stored + 1 - forces_add = np.array(force) - torques_add = np.array(torque) - except ValueError: - weight_add = weight_stored + weight - count_add = count_stored + 1 - # rolling_weighted_ave = (V_stored * w_stored + V * w) / - # (w_stored + w) - forces_add = (forces_stored * weight_stored + force * weight) / float( - weight_stored + weight - ) - torques_add = ( - torques_stored * weight_stored + torque * weight - ) / float(weight_stored + weight) - - if KE is not None: - try: - DOF, T = 0, 0 - for f in force[0]: - if f != 0: - DOF += 1 - else: - continue - - for t in torque[0]: - if t != 0: - DOF += 1 - else: - continue - - if DOF != 0: - T = KE / float(DOF) * 240.663461 - if DOF == 0: - T = KE * 240.663461 - T_add = (T + T_stored * count_stored) / float(count_add) - except TypeError: # force is a zero - pass - - except TypeError: - pass - - return weight_add, count_add, forces_add, torques_add, T_add - - -def SvibPopulation( - Aclass, - level, - atom_name, - atom_resname, - N_H, - atom_num, - nearest_resname, - RADshell_num_dist, - waterTuple, - UAweightedForces, - UAweightedTorques, - PE, - KE, - MweightedForces, - MweightedTorques, - Ndih, - molecule_UA_F, - molecule_UA_T, - weight, -): - - # DICTS FOR SVIB_ECALC - - # switch level for solutes only to get WM Svib - if level == "res_atomLevel" and atom_resname not in waterTuple: - level is None - atom_resname = atom_resname.split("_")[0] - if level == "residLevel_resname" and atom_resname not in waterTuple: - level is None - nearest_resname = nearest_resname.split("_")[0] - atom_info = (atom_resname, N_H, atom_name) - - # print(nearest_resname, atom_resname) - # running weighted average of X - - if RADshell_num_dist not in Aclass.RADshell_dict[nearest_resname][atom_info]: - Aclass.RADshell_dict[nearest_resname][atom_info][RADshell_num_dist] = [ - 0, - 0, - 0, - 0, - 0, - 0, - 0, - ] - - if "0" not in Aclass.RADshell_dict[nearest_resname][atom_info]: - Aclass.RADshell_dict[nearest_resname][atom_info]["0"] = [0, 0, 0, 0, 0, 0, 0] - - for RADshell_num_dist2 in [RADshell_num_dist, "0"]: - forces_stored = Aclass.RADshell_dict[nearest_resname][atom_info][ - RADshell_num_dist2 - ][0] - torques_stored = Aclass.RADshell_dict[nearest_resname][atom_info][ - RADshell_num_dist2 - ][1] - PE_stored = Aclass.RADshell_dict[nearest_resname][atom_info][ - RADshell_num_dist2 - ][2] - KE_stored = Aclass.RADshell_dict[nearest_resname][atom_info][ - RADshell_num_dist2 - ][3] - T_stored = Aclass.RADshell_dict[nearest_resname][atom_info][RADshell_num_dist2][ - 4 - ] - count_stored = Aclass.RADshell_dict[nearest_resname][atom_info][ - RADshell_num_dist2 - ][5] - weight_stored = Aclass.RADshell_dict[nearest_resname][atom_info][ - RADshell_num_dist2 - ][6] - - forces_add = forces_stored - torques_add = torques_stored - T_add = T_stored - count_add = count_stored - weight_add = weight_stored - - weight_add, count_add, forces_add, torques_add, T_add = ( - runningWeightedAverageFT( - weight, - weight_stored, - count_stored, - forces_stored, - torques_stored, - T_stored, - weight_add, - count_add, - forces_add, - torques_add, - T_add, - UAweightedForces, - UAweightedTorques, - KE, - ) - ) - - PE_add = PE_stored - KE_add = KE_stored - - if count_add != count_stored and PE is not None: - # PE_add = (PE + PE_stored * count_stored) / float(count_add) - # KE_add = (KE + KE_stored * count_stored) / float(count_add) - - PE_add = (PE_stored * weight_stored + PE * weight) / float( - weight_stored + weight - ) - KE_add = (KE_stored * weight_stored + KE * weight) / float( - weight_stored + weight - ) - - Aclass.RADshell_dict[nearest_resname][atom_info][RADshell_num_dist2][ - 0 - ] = forces_add - Aclass.RADshell_dict[nearest_resname][atom_info][RADshell_num_dist2][ - 1 - ] = torques_add - Aclass.RADshell_dict[nearest_resname][atom_info][RADshell_num_dist2][2] = PE_add - Aclass.RADshell_dict[nearest_resname][atom_info][RADshell_num_dist2][3] = KE_add - Aclass.RADshell_dict[nearest_resname][atom_info][RADshell_num_dist2][4] = T_add - Aclass.RADshell_dict[nearest_resname][atom_info][RADshell_num_dist2][ - 5 - ] = count_add - Aclass.RADshell_dict[nearest_resname][atom_info][RADshell_num_dist2][ - 6 - ] = weight_add - - atom_info = (atom_info, Ndih) - - # PROX SHELL FT ASSIGNMENT - - if type(MweightedForces) is not type(None): # and \ - # atom_resname in waterTuple: #water only for now - - if RADshell_num_dist not in Aclass.WM_FT_shell_dict[nearest_resname][atom_info]: - Aclass.WM_FT_shell_dict[nearest_resname][atom_info][RADshell_num_dist] = [ - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - ] - - forces_stored = Aclass.WM_FT_shell_dict[nearest_resname][atom_info][ - RADshell_num_dist - ][0] - torques_stored = Aclass.WM_FT_shell_dict[nearest_resname][atom_info][ - RADshell_num_dist - ][1] - PE_stored = Aclass.WM_FT_shell_dict[nearest_resname][atom_info][ - RADshell_num_dist - ][2] - KE_stored = Aclass.WM_FT_shell_dict[nearest_resname][atom_info][ - RADshell_num_dist - ][3] - T_stored = Aclass.WM_FT_shell_dict[nearest_resname][atom_info][ - RADshell_num_dist - ][4] - UAforces_stored = Aclass.WM_FT_shell_dict[nearest_resname][atom_info][ - RADshell_num_dist - ][5] - UAtorques_stored = Aclass.WM_FT_shell_dict[nearest_resname][atom_info][ - RADshell_num_dist - ][6] - count_stored = Aclass.WM_FT_shell_dict[nearest_resname][atom_info][ - RADshell_num_dist - ][7] - weight_stored = Aclass.WM_FT_shell_dict[nearest_resname][atom_info][ - RADshell_num_dist - ][8] - - forces_add = forces_stored - torques_add = torques_stored - T_add = T_stored - UAforces_add = UAforces_stored - UAtorques_add = UAtorques_stored - count_add = count_stored - weight_add = weight_stored - - weight_add, count_add, forces_add, torques_add, T_add = ( - runningWeightedAverageFT( - weight, - weight_stored, - count_stored, - forces_stored, - torques_stored, - T_stored, - weight_add, - count_add, - forces_add, - torques_add, - T_add, - MweightedForces, - MweightedTorques, - False, - ) - ) - - if len(molecule_UA_F) != 0 and len(molecule_UA_T) != 0: - weight_add, count_ignore, UAforces_add, UAtorques_add, T_add = ( - runningWeightedAverageFT( - weight, - weight_stored, - count_stored, - UAforces_stored, - UAtorques_stored, - T_stored, - weight_add, - count_add, - UAforces_add, - UAtorques_add, - T_add, - molecule_UA_F, - molecule_UA_T, - KE, - ) - ) - - PE_add = PE_stored - KE_add = KE_stored - - if count_add != count_stored and PE is not None: - # PE_add = (PE + PE_stored * count_stored) / float(count_add) - # KE_add = (KE + KE_stored * count_stored) / float(count_add) - - PE_add = (PE_stored * weight_stored + PE * weight) / float( - weight_stored + weight - ) - KE_add = (KE_stored * weight_stored + KE * weight) / float( - weight_stored + weight - ) - - Aclass.WM_FT_shell_dict[nearest_resname][atom_info][RADshell_num_dist][ - 0 - ] = forces_add - Aclass.WM_FT_shell_dict[nearest_resname][atom_info][RADshell_num_dist][ - 1 - ] = torques_add - Aclass.WM_FT_shell_dict[nearest_resname][atom_info][RADshell_num_dist][ - 2 - ] = PE_add - Aclass.WM_FT_shell_dict[nearest_resname][atom_info][RADshell_num_dist][ - 3 - ] = KE_add - Aclass.WM_FT_shell_dict[nearest_resname][atom_info][RADshell_num_dist][ - 4 - ] = T_add - Aclass.WM_FT_shell_dict[nearest_resname][atom_info][RADshell_num_dist][ - 5 - ] = UAforces_add - Aclass.WM_FT_shell_dict[nearest_resname][atom_info][RADshell_num_dist][ - 6 - ] = UAtorques_add - Aclass.WM_FT_shell_dict[nearest_resname][atom_info][RADshell_num_dist][ - 7 - ] = count_add - Aclass.WM_FT_shell_dict[nearest_resname][atom_info][RADshell_num_dist][ - 8 - ] = weight_add - - -def SconfPopulation( - Aclass, - nearest_resname, - atom_name, - atom_resname, - N_H, - dihedral_phi_list, - molecule_size, - RADshell_num_dist, - weight, -): - """ """ - - def roundPartial(value, resolution): - """ - used for bin widths of 30 degrees - """ - return round(value / resolution) * resolution - - if nearest_resname is not None and dihedral_phi_list is not None: - - assigned = (atom_resname, len(dihedral_phi_list), molecule_size, N_H, atom_name) - - counter = 0 - for dihedral_phi_type in dihedral_phi_list: - # print (dihedral_phi_type) - counter += 1 - # print ('c1:', counter) - dihedral_atoms = (tuple(dihedral_phi_type[0]), counter) - # atomnum_list, dihedral index - phi = dihedral_phi_type[1] # angle, degrees - # dih_index = dihedral_phi_type[2][0] #0, 1, 2 - # dih_type = dihedral_phi_type[2][1] #trans, g-, g+ - - # adaptive method - rounded_phi = roundPartial(phi, 30) - # print (phi, rounded_phi) - all_bins = list(range(-180, 181, 30)) - # print (all_bins) - for b in all_bins: - if ( - b - not in Aclass.adaptive_dih_dict[nearest_resname][assigned][ - RADshell_num_dist - ][dihedral_atoms] - ): - Aclass.adaptive_dih_dict[nearest_resname][assigned][ - RADshell_num_dist - ][dihedral_atoms][b] = 0 - Aclass.adaptive_dih_dict[nearest_resname][assigned][RADshell_num_dist][ - dihedral_atoms - ][ - rounded_phi - ] += weight # 1 diff --git a/CodeEntropy/poseidon/extractData/HBRAD.py b/CodeEntropy/poseidon/extractData/HBRAD.py deleted file mode 100644 index 92bc1f3..0000000 --- a/CodeEntropy/poseidon/extractData/HBRAD.py +++ /dev/null @@ -1,393 +0,0 @@ -#!/usr/bin/env python - -# import logging -import math - -import numpy as np - -from CodeEntropy.poseidon.extractData.generalFunctions import angle - -# import sys -# from datetime import datetime - - -def UALevelRAD(all_data, dimensions): - """ - UA level RAD algorithm implementation. - Iterate through the 20 closest UAs and find which ones - are not blocked by a closer UA. - """ - - for x in range(0, len(all_data)): - i = all_data[x] - if i.mass > 1.1: - RAD = [] - RAD_dist = [] - if i.nearest_sorted_array is not None: - range_limit = len(i.nearest_sorted_array) - if len(i.nearest_sorted_array) > 30: - range_limit = 30 - for y in range(0, range_limit): - j = all_data[i.nearest_sorted_array[y][0]] - rij = float(i.nearest_sorted_array[y][1]) - blocked = False - for z in range(0, y): # only closer UAs can block - k = all_data[i.nearest_sorted_array[z][0]] - rik = float(i.nearest_sorted_array[z][1]) - costheta_jik = float( - angle(j.coords, i.coords, k.coords, dimensions) - ) - if math.isnan(costheta_jik): - break - - LHS = (float(1) / rij) ** 2 - RHS = ((float(1) / rik) ** 2) * costheta_jik - if LHS < RHS: - blocked = True - break - - if blocked is False: - RAD.append(j) - RAD_dist.append((j, rij)) - - i.RAD = RAD - i.RAD_dist = RAD_dist - else: - continue - - """ - Make sure each UA is in the other UAs RAD shell if not, remove - check RAD shell contains resids other than the reference id - (used for systems containing surfaces that are all the same resid) - """ - - for x in range(0, len(all_data)): - atom = all_data[x] - if atom.mass > 1.1: - if len(atom.RAD) != 0: - RAD = [] # symmetric RAD shell here - RAD_resids = [] # checking if neighbours are diff molecules - # make sure atom is in neighbours RAD shell - for neighbour in atom.RAD: - if atom in neighbour.RAD: - RAD.append(neighbour) - RAD_resids.append(neighbour.resid) - else: - continue - - if len(RAD) != 0: - # when RAD shell contains one resid only - if all(n == RAD_resids[0] for n in RAD_resids) is True: - # make sure atom is not surrounded just by itself - if RAD_resids[0] == atom.resid: - atom.RAD = [] - else: - atom.RAD = RAD - else: - atom.RAD = RAD - if len(RAD) == 0: - atom.RAD = RAD - else: - continue - else: - continue - - -def distCutoffNc(all_data, dimensions, cutoff): - """ - Instead of RAD, use a fixed distance cutoff - to define the first coordination shell (Nc) - """ - - for x in range(0, len(all_data)): - i = all_data[x] - if i.mass > 1.1 and i.nearest_sorted_array is not None: - range_limit = len(i.nearest_sorted_array) - if len(i.nearest_sorted_array) > 30: - range_limit = 30 - RAD = [] - RAD_dist = [] - for y in range(0, range_limit): - rij = float(i.nearest_sorted_array[y][1]) - if rij <= cutoff: - near = all_data[i.nearest_sorted_array[y][0]] - RAD.append(near) - RAD_dist.append((near, rij)) - else: - continue - - i.RAD = RAD - i.RAD_dist = RAD_dist - - else: - continue - - -def HBCalc(all_data, waterTuple, dimensions): - """ - HB TYPE CALC METHOD 2 (including all hydrogens (not just - for water) and finding most closest and most negative - atom to donate to - (QA*QD)/r^2 approach) - - ***addition = hydrogen's can only HB to UAs that are in the RAD - shell of their bonded heavy atom. - """ - - for x in range(0, len(all_data)): - atom = all_data[x] - if atom.mass < 1.1 and atom.mass >= 1 and atom.charge > 0.1: - acceptor_charge = 99 - hydrogen = atom - oxygen = None - for bonded in hydrogen.bonded_to_atom_num: - b = all_data[bonded] - if b.mass > 1.1: - oxygen = b - else: - continue - - # Needed for RAD shell HB ranking, HBing inside RAD - # RAD_atom_nums = [neighbour.atom_num for neighbour in oxygen.RAD] - RAD_atom_nums = [] - if oxygen is not None: - for neighbour in oxygen.RAD: - RAD_atom_nums.append(neighbour.atom_num) - for bonded in neighbour.bonded_to_atom_num: - b = all_data[bonded] - if b.mass < 1.1: - RAD_atom_nums.append(b.atom_num) - else: - continue - - if hydrogen.nearest_all_atom_array is not None: - range_limit = len(hydrogen.nearest_all_atom_array) - if len(atom.nearest_all_atom_array) > 50: - range_limit = 50 - for atom_dist in hydrogen.nearest_all_atom_array[:range_limit]: - near = all_data[atom_dist[0]] - HX_dist = atom_dist[1] - - if ( - near.atom_num not in hydrogen.bonded_to_atom_num - and near.atom_num != hydrogen.atom_num - and near.charge != 0 - and HX_dist != 0 - and oxygen is not None - and near.atom_num in RAD_atom_nums - ): - - X = near - QD = hydrogen.charge - QA = X.charge - - r2 = HX_dist**2 - relative_charge = (float(QD) * float(QA)) / float(r2) - cosine_angle = angle( - oxygen.coords, hydrogen.coords, X.coords, dimensions - ) - angle1 = np.arccos(cosine_angle) - OHX_angle = np.degrees(angle1) - - if relative_charge < acceptor_charge and float(OHX_angle) > 90: - acceptor_charge = relative_charge - hydrogen.nearest_atom = X - hydrogen.dist = HX_dist - - else: - continue - else: - continue - else: - continue - - else: - continue - - broken_HBs = False - # find hydrogens with neighbouring eneg atoms and append those - # hydrogens to the nearest_Hs of those eneg atoms - for x in range(0, len(all_data)): - hydrogen = all_data[x] - if hydrogen.nearest_atom is not None: - nearest_eneg = hydrogen.nearest_atom.atom_num - all_data[nearest_eneg].nearest_Hs.append(hydrogen) - - if broken_HBs: - # Dealing with broken HBs - ND, bifurcated and cyclic - # Hard-coded for water resnames only for now - atom = all_data[x] - HHX = False - for b in atom.bonded_to_atom_num: - bonded = all_data[b] - if atom.mass < 1.1 and bonded.mass > 1.1: - if bonded.bondedUA_H is not None: - if bonded.bondedUA_H[0] == 0 and bonded.bondedUA_H[1] == 2: - HHX = True - - if atom.mass < 1.1 and atom.nearest_all_atom_array is not None and HHX: - for atom_dist in atom.nearest_all_atom_array[:30]: - atom2 = all_data[atom_dist[0]] - bonded_overlap = bool( - set(atom.bonded_to_atom_num) & set(atom2.bonded_to_atom_num) - ) - # check if both bonded to same atom - if ( - atom2.atom_num != atom.atom_num - and atom2.resid == atom.resid - and atom.nearest_atom is not None - and atom2.nearest_atom is not None - and bonded_overlap is True - ): - if ( - atom2.nearest_atom.atom_num == atom.nearest_atom.atom_num - and atom2.dist < atom.dist - ): - - atom.broken = [ - "Bifurcated", - atom.atom_num, - atom.atom_name, - atom.nearest_atom.atom_num, - atom.nearest_atom.atom_name, - ] # - - """ - print ('Bifurcated', atom.atom_num, - atom.atom_name, - atom.nearest_atom.atom_num, - atom.nearest_atom.atom_name) - """ - - for b in atom.bonded_to_atom_num: - bonded = all_data[b] - bonded.broken = [ - "Bifurcated", - atom.atom_num, - atom.atom_name, - atom.nearest_atom.atom_num, - atom.nearest_atom.atom_name, - ] # - - # remove this broken donor from acceptor - new_Hlist = [] - for hydrogen in atom.nearest_atom.nearest_Hs: - if hydrogen != atom: - new_Hlist.append(hydrogen) - else: - continue - - atom.nearest_atom.nearest_Hs = new_Hlist - - for b2 in atom.nearest_atom.bonded_to_atom_num: - bonded2 = all_data[b2] - if bonded2.mass < 1.1: - if ( - bonded2.nearest_atom is not None - and bonded2.nearest_atom.resid == atom2.resid - and bonded2.dist < atom.dist - ): - - atom.broken = [ - "ND", - atom.atom_num, - atom.atom_name, - atom.nearest_atom.atom_num, - atom.nearest_atom.atom_name, - ] - # - - """ - print ('ND', atom.atom_num, - atom.atom_name, - atom.nearest_atom.atom_num, - atom.nearest_atom.atom_name) - """ - - for b3 in atom.bonded_to_atom_num: - bonded3 = all_data[b3] - bonded3.broken = [ - "ND", - atom.atom_num, - atom.atom_name, - atom.nearest_atom.atom_num, - atom.nearest_atom.atom_name, - ] - # - - # remove this broken donor from acceptor - new_Hlist = [] - for hydrogen in atom.nearest_atom.nearest_Hs: - if hydrogen != atom: - new_Hlist.append(hydrogen) - else: - continue - - atom.nearest_atom.nearest_Hs = new_Hlist - atom.nearest_atom = None - - else: - continue - else: - continue - - atom.nearest_atom = None - - if ( - atom2.atom_num != atom.atom_num - and atom2.resid == atom.resid - and bonded_overlap is True - and len(atom2.nearest_Hs) != 0 - ): - for hydrogen in atom2.nearest_Hs: - if ( - atom.nearest_atom is not None - and hydrogen.resid == atom.nearest_atom.resid - and hydrogen.dist < atom.dist - ): - - bonded_overlap2 = bool( - set(atom.nearest_atom.bonded_to_atom_num) - & set(hydrogen.bonded_to_atom_num) - ) - # check if both bonded to same atom - - if bonded_overlap2 is True: - - atom.broken = [ - "Cyclic", - atom.atom_num, - atom.atom_name, - atom.nearest_atom.atom_num, - atom.nearest_atom.atom_name, - ] # - - """ - print ('Cyclic', atom.atom_num, - atom.atom_name, - atom.nearest_atom.atom_num, - atom.nearest_atom.atom_name) - """ - - for b4 in atom.bonded_to_atom_num: - bonded4 = all_data[b4] - bonded4.broken = [ - "Cyclic", - atom.atom_num, - atom.atom_name, - atom.nearest_atom.atom_num, - atom.nearest_atom.atom_name, - ] # - - # remove this broken donor from acceptor - new_Hlist = [] - for hydrogen in atom.nearest_atom.nearest_Hs: - if hydrogen != atom: - new_Hlist.append(hydrogen) - else: - continue - - atom.nearest_atom.nearest_Hs = new_Hlist - atom.nearest_atom = None - - else: - continue diff --git a/CodeEntropy/poseidon/extractData/__init__.py b/CodeEntropy/poseidon/extractData/__init__.py deleted file mode 100644 index 201e256..0000000 --- a/CodeEntropy/poseidon/extractData/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -""" -Empty init file -""" diff --git a/CodeEntropy/poseidon/extractData/dihedrals.py b/CodeEntropy/poseidon/extractData/dihedrals.py deleted file mode 100644 index 9cbb450..0000000 --- a/CodeEntropy/poseidon/extractData/dihedrals.py +++ /dev/null @@ -1,151 +0,0 @@ -#!/usr/bin/env python - -# import math -# import sys - -import numpy as np - -from CodeEntropy.poseidon.extractData.generalFunctions import distance, vector - - -def calculateDihedrals(all_data, dimensions): - """ - In-place modification of all_data to calculate dihedrals - For each dihedral in system, find what conformation it is in - Dihedrals: - - - X4 - / - / - X2---X3 - / - / - X1 - - - b1 = vector(X1, X2) - b2 = vector(X2, X3) - b3 = vector(X3, X4) - - source = https://math.stackexchange.com/questions/47059/ - how-do-i-calculate-a-dihedral-angle-given-cartesian-coordinates - - Args: - all_data (list of poseidon.atom_info): Container for all atom info - dimensions (numpy array of size (1,3)): dimension of box of current frame - """ - - rad2deg = float(180) / float(np.pi) - - dihedrals = [] - all_atom_nums_list = [] - - max_atom_num = max(all_data[x].atom_num for x in range(0, len(all_data))) - - for x in range(0, len(all_data)): - atom = all_data[x] - if ( - atom.mass > 1.1 - and len(atom.dihedral_list) != 0 - and atom.atom_num not in all_atom_nums_list - ): - - Mdihedrals = [] # unique dihedrals for each atom in molecule - lowest_atom_num = max_atom_num + 1 - for molecule_atom in atom.molecule_atomNums: - Matom = all_data[molecule_atom] - all_atom_nums_list.append(Matom.atom_num) - # get lowest atom num for renumbering - for dih in Matom.dihedral_list: - for d in dih: - if all_data[d].resid == atom.resid: - if d < lowest_atom_num: - lowest_atom_num = d - else: - continue - - for molecule_atom in atom.molecule_atomNums: - Matom = all_data[molecule_atom] - for dih in Matom.dihedral_list: - - inside_resid = True - for d in dih: - if all_data[d].resid != atom.resid: - inside_resid = False - - if dih not in dihedrals and inside_resid is True: - # only analyse dihedral once - dihedrals.append(dih) - - d_atoms = [] - atom_nums = [] # 1,2,3,4 - for d in dih: # calc dih angles for each dih - atom2 = all_data[d] - d_atoms.append(atom2) - atom_nums.append(atom2.atom_num) - - # re-number atom_nums - renumbered_atom_nums = [] - for num in atom_nums: - renumbered = (num - lowest_atom_num) + 1 - renumbered_atom_nums.append(renumbered) - - # get vectors - b1_vector = vector( - d_atoms[0].coords, d_atoms[1].coords, dimensions - ) - b1_dist = distance( - d_atoms[0].coords, d_atoms[1].coords, dimensions - ) - b1_norm = np.divide(b1_vector, b1_dist) - - b2_vector = vector( - d_atoms[1].coords, d_atoms[2].coords, dimensions - ) - b2_dist = distance( - d_atoms[1].coords, d_atoms[2].coords, dimensions - ) - b2_norm = np.divide(b2_vector, b2_dist) - - b3_vector = vector( - d_atoms[2].coords, d_atoms[3].coords, dimensions - ) - b3_dist = distance( - d_atoms[2].coords, d_atoms[3].coords, dimensions - ) - b3_norm = np.divide(b3_vector, b3_dist) - - # get normals - n1 = np.cross(b1_norm, b2_norm) - n2 = np.cross(b2_norm, b3_norm) - m = np.cross(n1, b2_norm) - - # get dot products - x = np.dot(n1, n2) - y = np.dot(m, n2) - - phi = rad2deg * np.arctan2(y, x) - dih_type = None - - # only needed for fixed-dihedrals - if phi >= 120 or phi < -120: - dih_type = [0, "trans"] - - if phi >= 0 and phi < 120: - dih_type = [1, "g-"] - - if phi >= -120 and phi < 0: - dih_type = [2, "g+"] - - Mdihedrals.append( - [renumbered_atom_nums, int(round(phi, 0)), dih_type] - ) - - else: - continue - - atom.dihedral_phi_type = Mdihedrals - - else: - continue diff --git a/CodeEntropy/poseidon/extractData/forceTorques.py b/CodeEntropy/poseidon/extractData/forceTorques.py deleted file mode 100644 index 78510fb..0000000 --- a/CodeEntropy/poseidon/extractData/forceTorques.py +++ /dev/null @@ -1,709 +0,0 @@ -#!/usr/bin/env python - -# import logging -# import math -# import sys -from collections import defaultdict - -import numpy as np -from numpy import linalg as LA - -from CodeEntropy.poseidon.extractData.generalFunctions import MOI, com, distance, vector - - -# create nested dict in one go -def nested_dict(): - return defaultdict(nested_dict) - - -def calculateFTMatrix(all_data, dimensions): - """ - Consider force/torque calculations as discussed in Ali 2019 paper. - UA level and molecule level. - """ - - resid_list = [] - for x in range(0, len(all_data)): - atom = all_data[x] - if atom.mass > 1.1: - if atom.resid not in resid_list: # iterate though resids once - resid_list.append(atom.resid) - all_molecule_atoms = [] # list of heavy atoms and Hs - UAs_list = [] # list of lists [HA + Hs] - - for num in atom.molecule_atomNums: # atom nums in molecule - UA = all_data[num] - - all_molecule_atoms.append(UA) - heavy_bonded = [] - H_bonded = [] - - for bonded_atom in UA.bonded_to_atom_num: - bonded = all_data[bonded_atom] - if bonded.mass > 1.1 and bonded.atom_num != UA.atom_num: - heavy_bonded.append(bonded) - - if bonded.mass < 1.1 and bonded.atom_num != UA.atom_num: - H_bonded.append(bonded) - all_molecule_atoms.append(bonded) - - UA_atoms_list = [UA] + H_bonded - UAs_list.append(UA_atoms_list) - bonded_atoms_list = [UA] + heavy_bonded + H_bonded - UA.bondedUA_H = [len(heavy_bonded), len(H_bonded)] - - molecule_coords = [] - molecule_masses = [] - for a in all_molecule_atoms: - molecule_coords.append(a.coords) - molecule_masses.append(a.mass) - molecule_COM = com(molecule_coords, molecule_masses) - - # ****** MOLECULE LEVEL ****** - # if molecule is only one UA with bonded Hs, - # this works for one UA molecules - if ( - len(atom.molecule_atomNums) == 1 - and len(atom.bonded_to_atom_num) > 0 - ): - WM_principal_axes, WM_MI_axis = principalAxesMOI( - all_data, all_molecule_atoms, molecule_COM, Hs=True - ) - # use atom list containing Hs for P axes inc. Hs - WM_force, WM_torque = rotateFT( - all_data, - all_molecule_atoms, - WM_principal_axes, - WM_principal_axes, - WM_MI_axis, - molecule_COM, - molecule_COM, - Hs=True, - ) - - # - WM_force = np.outer(WM_force, WM_force) - WM_torque = np.outer(WM_torque, WM_torque) - atom.MweightedForces = np.round(np.divide(WM_force, 4), 3) - atom.MweightedTorques = np.round(np.divide(WM_torque, 4), 3) - atom.WMprincipalAxis = WM_principal_axes, WM_MI_axis, molecule_COM - - # check this with Richard, do you not halve here? - # For Jon's code, you do halve here - # even though molecule == UA level - atom.UAweightedForces = np.round(np.divide(WM_force, 4), 3) - atom.UAweightedTorques = np.round(np.divide(WM_torque, 4), 3) - atom.molecule_UA_Fs = np.round(np.divide(WM_force, 4), 3) - atom.molecule_UA_Ts = np.round(np.divide(WM_torque, 4), 3) - # - - """ - atom.MweightedForces = WM_force - atom.MweightedTorques = WM_torque - atom.UAweightedForces = WM_force - atom.UAweightedTorques = WM_torque - """ - - # if molecule contains more than one UA - # forces work, torques work if Hs = True - if len(atom.molecule_atomNums) > 1: - - WM_principal_axes, WM_MI_axis = principalAxesMOI( - all_data, all_molecule_atoms, molecule_COM, Hs=False - ) - # only use heavy atom list for P axes calc - # we don't consider Hs as they result - # in larger torques as Hs are further away - - WM_force, WM_torque = rotateFT( - all_data, - all_molecule_atoms, - WM_principal_axes, - WM_principal_axes, - WM_MI_axis, - molecule_COM, - molecule_COM, - Hs=True, - ) - - # - # WM_FT = np.concatenate((WM_force, WM_torque), axis=None) - # WM_FT = np.outer(WM_FT, WM_FT) - WM_force = np.outer(WM_force, WM_force) - WM_torque = np.outer(WM_torque, WM_torque) - atom.MweightedForces = np.round(np.divide(WM_force, 4), 3) - atom.MweightedTorques = np.round(np.divide(WM_torque, 4), 3) - atom.WMprincipalAxis = WM_principal_axes, WM_MI_axis, molecule_COM - - atom.molecule_UA_Fs = np.round(np.divide(WM_force, 4), 3) - atom.molecule_UA_Ts = np.round(np.divide(WM_torque, 4), 3) - # - - # atom.MweightedForces = WM_force - # atom.MweightedTorques = WM_torque - - # ****** UNITED-ATOM LEVEL ****** - # works! - UA_F_list = [] - UA_T_list = [] - for UA in UAs_list: - XX_principal_axes, UA_MI_axis = None, None - - if UA[0].mass > 1.1: - bonded_HAs = [UA[0]] # inc itself - bonded_Hs = [] - for b in UA[0].bonded_to_atom_num: - bonded = all_data[b] - if bonded.mass < 1.1: - bonded_Hs.append(bonded) - elif bonded.mass > 1.1: - bonded_HAs.append(bonded) - else: - continue - - if len(bonded_HAs) == 2 and len(bonded_Hs) > 0: - # num HAs including itself - UA_COM = bonded_HAs[0].coords - - R = bonded_HAs[0] - X1 = bonded_HAs[1] - H = bonded_Hs[0] - - RX1_vector = vector(R.coords, X1.coords, dimensions) - RX1_dist = distance(R.coords, X1.coords, dimensions) - Paxis1 = np.divide(RX1_vector, RX1_dist) - - RH_vector = vector(R.coords, H.coords, dimensions) - RH_dist = distance(R.coords, H.coords, dimensions) - RH_vector_norm = np.divide(RH_vector, RH_dist) - - Paxis2 = np.cross(Paxis1, RH_vector_norm) - Paxis3 = np.cross(Paxis1, Paxis2) - - XX_principal_axes = [Paxis1, Paxis2, Paxis3] - - UA_MI_axis, XX_principal_axes = UA_MOI( - all_data, - UA, - UA_COM, - XX_principal_axes, - dimensions, - Hs=True, - ) - - if len(bonded_Hs) == 1: - MIx, MIy, MIz = ( - UA_MI_axis[0], - UA_MI_axis[1], - UA_MI_axis[2], - ) - if MIx < MIy and MIx < MIz: - MIx = 0 - if MIy < MIx and MIy < MIz: - MIy = 0 - if MIz < MIx and MIz < MIy: - MIz = 0 - - UA_MI_axis = [MIx, MIy, MIz] - - elif len(bonded_HAs) > 2: - # num HAs including itself - UA_COM = bonded_HAs[0].coords - - R = bonded_HAs[0] - X1 = bonded_HAs[1] - X2 = bonded_HAs[2] - - Paxis1 = np.zeros(3) - # average of all HA covalent bond vectors - for HA in bonded_HAs[1:]: - Paxis1 += vector(R.coords, HA.coords, dimensions) - - X1X2_vector = vector(X1.coords, X2.coords, dimensions) - X1X2_dist = distance(X1.coords, X2.coords, dimensions) - X1X2_norm = np.divide(X1X2_vector, X1X2_dist) - Paxis2 = np.cross(X1X2_norm, Paxis1) - - Paxis3 = np.cross(Paxis1, Paxis2) - - XX_principal_axes = [Paxis1, Paxis2, Paxis3] - - UA_MI_axis, XX_principal_axes = UA_MOI( - all_data, - UA, - UA_COM, - XX_principal_axes, - dimensions, - Hs=True, - ) - - elif len(bonded_HAs) == 2 and len(bonded_Hs) == 0: - # num use arbitrary coords for 2nd vector - UA_COM = bonded_HAs[0].coords - - R = bonded_HAs[0] - X1 = bonded_HAs[1] - - RX1_vector = vector(R.coords, X1.coords, dimensions) - RX1_dist = distance(R.coords, X1.coords, dimensions) - Paxis1 = np.divide(RX1_vector, RX1_dist) - - RZ_vector = vector(R.coords, [0, 0, 0], dimensions) - RZ_dist = distance(R.coords, [0, 0, 0], dimensions) - RZ_vector_norm = np.divide(RZ_vector, RZ_dist) - - Paxis2 = np.cross(Paxis1, RZ_vector_norm) - Paxis3 = np.cross(Paxis1, Paxis2) - - XX_principal_axes = [Paxis1, Paxis2, Paxis3] - - UA_MI_axis, XX_principal_axes = UA_MOI( - all_data, - UA, - UA_COM, - XX_principal_axes, - dimensions, - Hs=True, - ) - - if len(bonded_Hs) == 0: - MIx, MIy, MIz = ( - UA_MI_axis[0], - UA_MI_axis[1], - UA_MI_axis[2], - ) - if MIx < MIy and MIx < MIz: - MIx = 0 - if MIy < MIx and MIy < MIz: - MIy = 0 - if MIz < MIx and MIz < MIy: - MIz = 0 - - UA_MI_axis = [MIx, MIy, MIz] - - else: - continue - - if UA[0].mass < 1.1: - print("Error: no HA in UA") - - UA_force, UA_torque = rotateFT( - all_data, - UA, - WM_principal_axes, - XX_principal_axes, - UA_MI_axis, - molecule_COM, - UA_COM, - Hs=True, - ) - # torques use HA com, - # forces use molecule COM - - UA_F_list += [UA_force] - UA_T_list += [UA_torque] - - # - UA_force = np.outer(UA_force, UA_force) - UA_torque = np.outer(UA_torque, UA_torque) - UA[0].UAweightedForces = np.round(np.divide(UA_force, 4), 3) - # UA[0].UAweightedForces = UA_force - # 2019 paper, dont halve forces - UA[0].UAweightedTorques = np.round(np.divide(UA_torque, 4), 3) - # - - # UA[0].UAweightedForces = UA_force - # UA[0].UAweightedTorques = UA_torque - - # dont halve forces in molecule UAs (2019 paper), - # divide by 4 to compare with jons output - molecule_UA_Fs = np.concatenate(UA_F_list, axis=None) - molecule_UA_Fs = np.outer(molecule_UA_Fs, molecule_UA_Fs) - # atom.molecule_UA_Fs = np.divide(molecule_UA_Fs, 4) - atom.molecule_UA_Fs = np.round(molecule_UA_Fs) - - # halve torques in molecule UAs - molecule_UA_Ts = np.concatenate(UA_T_list, axis=None) - molecule_UA_Ts = np.outer(molecule_UA_Ts, molecule_UA_Ts) - atom.molecule_UA_Ts = np.round(np.divide(molecule_UA_Ts, 4), 3) - - # if molecule is monoatomic (one UA and no bonded Hs) - if ( - len(atom.molecule_atomNums) == 1 - and len(atom.bonded_to_atom_num) == 0 - ): - mass_sqrt = float(atom.mass**0.5) - forces_sqrt = [ - float(atom.forces[0]) / mass_sqrt, - float(atom.forces[1]) / mass_sqrt, - float(atom.forces[2]) / mass_sqrt, - ] - forces_sqrt = np.sort(forces_sqrt) - - # - UA_force = np.outer(forces_sqrt, forces_sqrt) - UA_torque = np.outer([0, 0, 0], [0, 0, 0]) - atom.UAweightedForces = np.round(np.divide(UA_force, 4), 3) - atom.UAweightedTorques = np.round(UA_torque, 3) # zero - atom.MweightedForces = np.round(np.divide(UA_force, 4), 3) - atom.MweightedTorques = np.round(UA_torque, 3) # zero - atom.WMprincipalAxis = ( - [[0, 0, 0], [0, 0, 0], [0, 0, 0]], - [0, 0, 0], - [0, 0, 0], - ) - - atom.molecule_UA_Fs = np.round(np.divide(UA_force, 4), 3) - atom.molecule_UA_Ts = np.round(np.divide(UA_torque, 4), 3) - # - - """ - atom.UAweightedForces = forces_sqrt - atom.UAweightedTorques = [0, 0, 0] #zero - atom.MweightedForces = forces_sqrt - atom.MweightedTorques = [0, 0, 0] #zero - """ - - -def UA_MOI(all_data, bonded_atoms_list, COM, PIaxes, dimensions, Hs): - """ - Get MOI axis for UA level, where eigenvalues and vectors are - not used - """ - - coord_list = [] - mass_list = [] - - for atom in bonded_atoms_list: - if Hs is True: - # print (atom.atom_name) - coord_list.append(atom.coords) - mass_list.append(atom.mass) - elif Hs is False: - if atom.mass > 1.1: - UA_mass = 0 # add mass of HA and bonded Hs - # print (atom.atom_name) - UA_mass += atom.mass - coord_list.append(atom.coords) # only include coords of UA - for h in atom.bonded_to_atom_num: - H = all_data[h] - if H.mass < 1.1: - UA_mass += H.mass - else: - continue - mass_list.append(UA_mass) - # include mass of Hs on heavy atom mass - else: - continue - - # sorting out PIaxes for MoI for UA fragment - - modPIx = PIaxes[0][0] ** 2 + PIaxes[0][1] ** 2 + PIaxes[0][2] ** 2 - modPIy = PIaxes[1][0] ** 2 + PIaxes[1][1] ** 2 + PIaxes[1][2] ** 2 - modPIz = PIaxes[2][0] ** 2 + PIaxes[2][1] ** 2 + PIaxes[2][2] ** 2 - - PIaxes[0][0] = PIaxes[0][0] / (modPIx**0.5) - PIaxes[0][1] = PIaxes[0][1] / (modPIx**0.5) - PIaxes[0][2] = PIaxes[0][2] / (modPIx**0.5) - - PIaxes[1][0] = PIaxes[1][0] / (modPIy**0.5) - PIaxes[1][1] = PIaxes[1][1] / (modPIy**0.5) - PIaxes[1][2] = PIaxes[1][2] / (modPIy**0.5) - - PIaxes[2][0] = PIaxes[2][0] / (modPIz**0.5) - PIaxes[2][1] = PIaxes[2][1] / (modPIz**0.5) - PIaxes[2][2] = PIaxes[2][2] / (modPIz**0.5) - - # get dot product of Paxis1 and CoM->atom1 vect - # will just be [0,0,0] - RRaxis = vector(coord_list[0], COM, dimensions) - # flip each Paxis if its pointing out of UA - dotProd1 = np.dot(np.array(PIaxes[0]), RRaxis) - if dotProd1 < 0: - PIaxes[0][0] = -PIaxes[0][0] - PIaxes[0][1] = -PIaxes[0][1] - PIaxes[0][2] = -PIaxes[0][2] - - dotProd2 = np.dot(np.array(PIaxes[1]), RRaxis) - if dotProd2 < 0: - PIaxes[1][0] = -PIaxes[1][0] - PIaxes[1][1] = -PIaxes[1][1] - PIaxes[1][2] = -PIaxes[1][2] - - dotProd3 = np.dot(np.array(PIaxes[2]), RRaxis) - if dotProd3 < 0: - PIaxes[2][0] = -PIaxes[2][0] - PIaxes[2][1] = -PIaxes[2][1] - PIaxes[2][2] = -PIaxes[2][2] - - # MI = moment of inertia - axis1MI = 0 # x - axis2MI = 0 # y - axis3MI = 0 # z - # - - for coord, mass in zip(coord_list, mass_list): - dx = coord[0] - COM[0] - dy = coord[1] - COM[1] - dz = coord[2] - COM[2] - - PIaxes_xx = PIaxes[0][1] * dz - PIaxes[0][2] * dy - PIaxes_xy = PIaxes[0][2] * dx - PIaxes[0][0] * dz - PIaxes_xz = PIaxes[0][0] * dy - PIaxes[0][1] * dx - - PIaxes_yx = PIaxes[1][1] * dz - PIaxes[1][2] * dy - PIaxes_yy = PIaxes[1][2] * dx - PIaxes[1][0] * dz - PIaxes_yz = PIaxes[1][0] * dy - PIaxes[1][1] * dx - - PIaxes_zx = PIaxes[2][1] * dz - PIaxes[2][2] * dy - PIaxes_zy = PIaxes[2][2] * dx - PIaxes[2][0] * dz - PIaxes_zz = PIaxes[2][0] * dy - PIaxes[2][1] * dx - - daxisx = (PIaxes_xx**2 + PIaxes_xy**2 + PIaxes_xz**2) ** 0.5 - daxisy = (PIaxes_yx**2 + PIaxes_yy**2 + PIaxes_yz**2) ** 0.5 - daxisz = (PIaxes_zx**2 + PIaxes_zy**2 + PIaxes_zz**2) ** 0.5 - - axis1MI += (daxisx**2) * mass - axis2MI += (daxisy**2) * mass - axis3MI += (daxisz**2) * mass - - UA_moI = [axis1MI, axis2MI, axis3MI] - - return UA_moI, PIaxes - - -def principalAxesMOI(all_data, bonded_atoms_list, COM, Hs): - """ - Calculate the principal axes of the MOI matrix, then calc forces. - If Hs is True, then consider H coords in MOI calcs. - If Hs is False, then don't include H coords, instead add their mass - onto UA heavy atom (HA). - """ - - coord_list = [] - mass_list = [] - - for atom in bonded_atoms_list: - if Hs is True: - # print (atom.atom_name) - coord_list.append(atom.coords) - mass_list.append(atom.mass) - elif Hs is False: - if atom.mass > 1.1: - UA_mass = 0 # add mass of HA and bonded Hs - # print (atom.atom_name) - UA_mass += atom.mass - coord_list.append(atom.coords) # only include coords of UA - for h in atom.bonded_to_atom_num: - H = all_data[h] - if H.mass < 1.1: - UA_mass += H.mass - else: - continue - mass_list.append(UA_mass) - # include mass of Hs on heavy atom mass - else: - continue - - moI = MOI(COM, coord_list, mass_list) - # if Hs == False: - # print ('masses', mass_list) - # print ('moi', moI) - - eigenvalues, eigenvectors = LA.eig(moI) - - # different values generated to Jon's code - # print (eigenvalues) - # print (eigenvectors) - transposed = np.transpose(eigenvectors) # turn columns to rows - - # bonded_atoms_list[0].WMprincipalAxis = transposed, eigenvalues, COM - # print (moI, eigenvalues, eigenvectors) - - min_eigenvalue = abs(eigenvalues[0]) - if eigenvalues[1] < min_eigenvalue: - min_eigenvalue = eigenvalues[1] - if eigenvalues[2] < min_eigenvalue: - min_eigenvalue = eigenvalues[2] - - max_eigenvalue = abs(eigenvalues[0]) - if eigenvalues[1] > max_eigenvalue: - max_eigenvalue = eigenvalues[1] - if eigenvalues[2] > max_eigenvalue: - max_eigenvalue = eigenvalues[2] - - # print (min_eigenvalue * const, max_eigenvalue * const) - # same as Jons - - # - # PA = principal axes - axis1PA = [0, 0, 0] # [x,y,z] - axis2PA = [0, 0, 0] # [x,y,z] - axis3PA = [0, 0, 0] # [x,y,z] - - # MI = moment of inertia - axis1MI = 0 # x - axis2MI = 0 # y - axis3MI = 0 # z - # - - for i in range(0, 3): - if eigenvalues[i] == max_eigenvalue: - axis1PA = transposed[i] - axis1MI = eigenvalues[i] - elif eigenvalues[i] == min_eigenvalue: - axis3PA = transposed[i] - axis3MI = eigenvalues[i] - else: - axis2PA = transposed[i] - axis2MI = eigenvalues[i] - - principal_axes = [axis1PA, axis2PA, axis3PA] - MI_axis = [axis1MI, axis2MI, axis3MI] - - # print ('mI', MI_axis, bonded_atoms_list[0].atom_num) - - return principal_axes, MI_axis - - -def rotateFT( - all_data, - bonded_atoms_list, - WM_principal_axes, - UA_principal_axes, - MI_axis, - molecule_COM, - UA_COM, - Hs, -): - """ """ - - # lists used for torque calcs - coord_list = [] - mass_list = [] - force_list = [] - forces_summed = np.zeros(3) - - for atom in bonded_atoms_list: - if Hs is True: - coord_list.append(atom.coords) - mass_list.append(atom.mass) - force_list.append(atom.forces) - forces_summed += atom.forces - elif Hs is False: - if atom.mass > 1.1: - # print (atom.atom_name) - UA_mass = 0 # add mass of HA and bonded Hs - UA_forces = np.zeros(3) - UA_mass += atom.mass - UA_forces += np.array(atom.forces) - coord_list.append(atom.coords) # only include coords of UA - forces_summed += atom.forces - for h in atom.bonded_to_atom_num: - H = all_data[h] - if H.mass < 1.1: - UA_mass += H.mass - UA_forces += np.array(H.forces) - forces_summed += H.forces - else: - continue - mass_list.append(UA_mass) # include mass of Hs on heavy atom mass - force_list.append(UA_forces) - else: - continue - - torque = [0, 0, 0] # default is no torque - # calc torques here with COM and PA of molecule or UA - if UA_principal_axes is not None and MI_axis is not None: - if Hs is True: - torque = calcTorque( - UA_COM, coord_list, UA_principal_axes, force_list, MI_axis - ) - if Hs is False: - torque = calcTorque( - molecule_COM, coord_list, WM_principal_axes, force_list, MI_axis - ) - - # Rotate forces here - F1 = np.dot(forces_summed, WM_principal_axes[0]) - F2 = np.dot(forces_summed, WM_principal_axes[1]) - F3 = np.dot(forces_summed, WM_principal_axes[2]) - mass_sqrt = sum(mass_list) ** 0.5 - force = np.array( - [ - float(F1) / float(mass_sqrt), - float(F2) / float(mass_sqrt), - float(F3) / float(mass_sqrt), - ] - ) - - return force, torque - - -def calcTorque(cm, coord_list, principal_axes, force_list, MI_axis): - """ """ - x_cm, y_cm, z_cm = cm[0], cm[1], cm[2] - # MI_axis = np.sqrt(MI_axis) #sqrt moi to weight torques - - MI_axis_sqrt = [] - # print (MI_axis) - for coord in MI_axis: - coord_round = round(coord, 10) - if coord_round != 0: - c_sqrt = coord**0.5 - MI_axis_sqrt.append(c_sqrt) - elif coord_round == 0: - MI_axis_sqrt.append(coord) - else: - continue - - MI_axis_sqrt = np.array(MI_axis_sqrt) - # print (MI_axis, MI_axis_sqrt) - - W_torque = np.zeros(3) - for coord, force in zip(coord_list, force_list): - atom_torque = np.array([0, 0, 0]) - # count_atom += 1 - new_coords = [] - new_forces = [] - for axis in principal_axes: - new_c = ( - axis[0] * (coord[0] - x_cm) - + axis[1] * (coord[1] - y_cm) - + axis[2] * (coord[2] - z_cm) - ) - new_f = axis[0] * force[0] + axis[1] * force[1] + axis[2] * force[2] - new_coords.append(new_c) - new_forces.append(new_f) - # print (new_coords) - # print (force) - # CROSS PRODUCT - torquex = float( - new_coords[1] * new_forces[2] - new_coords[2] * new_forces[1] - ) # / float(1e10) - torquey = float( - new_coords[2] * new_forces[0] - new_coords[0] * new_forces[2] - ) # / float(1e10) - torquez = float( - new_coords[0] * new_forces[1] - new_coords[1] * new_forces[0] - ) # / float(1e10) - # print (torquex, torquey, torquez) - - atom_torque = (torquex, torquey, torquez) - # atom_torque = np.divide(atom_torque, MI_axis_sqrt) - - atom_torque2 = [] - for t, MI in zip(atom_torque, MI_axis_sqrt): - if MI != 0: - t_divide = float(t) / float(MI) - atom_torque2.append(t_divide) - elif MI == 0: - atom_torque2.append(0) - else: - continue - - atom_torque2 = np.array(atom_torque2) - - W_torque += atom_torque2 - - return W_torque diff --git a/CodeEntropy/poseidon/extractData/generalFunctions.py b/CodeEntropy/poseidon/extractData/generalFunctions.py deleted file mode 100644 index 2fec5b0..0000000 --- a/CodeEntropy/poseidon/extractData/generalFunctions.py +++ /dev/null @@ -1,608 +0,0 @@ -#!/usr/bin/env python - -import logging -import math - -import numpy as np -from numpy import linalg as LA - -# import sys - - -def distance(x0, x1, dimensions): - """ - calculates distance and accounts for PBCs. - """ - - x0 = np.array(x0) - x1 = np.array(x1) - delta = np.abs(x1 - x0) - delta = np.where(delta > 0.5 * dimensions, delta - dimensions, delta) - - if np.any(delta < 0): - negative_indices = np.where(delta < 0) - negative_values = delta[negative_indices] - raise ValueError( - f"Negative values encountered in 'delta' at indices {negative_indices} " - f"with values {negative_values}. " - f"Cannot take square root of negative values." - ) - - dist = np.sqrt((delta**2).sum(axis=-1)) - return dist - - -def vector(x0, x1, dimensions): - """ - get vector of two coordinates over PBCs. - """ - - x0 = np.array(x0) - x1 = np.array(x1) - dimensions = np.array(dimensions) - delta = x1 - x0 - delta2 = [] - for delt, box in zip(delta, dimensions): - if delt < 0 and delt < 0.5 * box: - delt = delt + box - if delt > 0 and delt > 0.5 * box: - delt = delt - box - - delta2.append(delt) - - delta2 = np.array(delta2) - - return delta2 - - -def getMcoords(x0, x1, dimensions): - """ - Get the middle (M) coords of two coords. - """ - - x0 = np.array(x0) - x1 = np.array(x1) - dimensions = np.array(dimensions) - - delta = x1 - x0 - - delta2 = [] - for delt, box in zip(delta, dimensions): - if delt < 0 and delt < 0.5 * box: - delt = delt + box - if delt > 0 and delt > 0.5 * box: - delt = delt - box - - delta2.append(delt) - - delta2 = np.array(delta2) - - M_coords2 = 0.5 * delta2 + x0 - - return M_coords2 - - -def projectVectorToPlane(vector, plane_normal): - """ - Project a vector onto a plane, need normal of plane as input. - """ - - plane_normal_magnitude = np.linalg.norm(plane_normal) - - project_vector = np.cross(vector, plane_normal) - project_vector = np.divide(project_vector, plane_normal_magnitude) - project_vector = np.cross(plane_normal, project_vector) - projected_vector = np.divide(project_vector, plane_normal_magnitude) - - return projected_vector - - -def normaliseVector(x0, x1, dimensions): - """ - Get vector of two coordinates over PBCs and then normalise by the - magnitude of the vector - """ - - x0 = np.array(x0) - x1 = np.array(x1) - dimensions = np.array(dimensions) - delta = x1 - x0 - - delta2 = [] - for delt, box in zip(delta, dimensions): - if delt < 0 and delt < 0.5 * box: - delt = delt + box - if delt > 0 and delt > 0.5 * box: - delt = delt - box - - delta2.append(delt) - - delta = np.array(delta2) - - delta_magnitude = (delta[0] ** 2 + delta[1] ** 2 + delta[2] ** 2) ** 0.5 - delta_norm = np.divide(delta, delta_magnitude) - - return delta_norm - - -def angleBetweenVectors(x0, x1): - """ - Get the angle between two vectors to determine how orthogonal they are. - """ - - top = np.dot(x0, x1) - - axis1_magnitude = np.linalg.norm(x0) # magnitude/ distance - axis2_magnitude = np.linalg.norm(x1) # magnitude/ distance - bottom = np.multiply(axis1_magnitude, axis2_magnitude) - top_bottom = np.divide(top, bottom) - if top_bottom > 1: - top_bottom = 1 - if top_bottom < -1: - top_bottom = -1 - theta = np.arccos(top_bottom) - theta = np.degrees(theta) - - return theta - - -def angle2(a, b, c, dimensions): - """ - Just used as a check. Doesn't work as periodic boundary conditions - are not considered properly. - """ - - a = np.array(a) # j - b = np.array(b) # i - c = np.array(c) # k - dimensions = np.array(dimensions) - ba = a - b - bc = c - b - ac = c - a - ba = np.where(ba > 0.5 * dimensions, ba - dimensions, ba) - bc = np.where(bc > 0.5 * dimensions, bc - dimensions, bc) - ac = np.where(ac > 0.5 * dimensions, ac - dimensions, ac) - - cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc)) - - angle = np.arccos(cosine_angle) - angle = np.degrees(angle) - - return angle - - -def angle(a, b, c, dimensions): - - a = np.array(a) # j - b = np.array(b) # i - c = np.array(c) # k - dimensions = np.array(dimensions) - ba = np.abs(a - b) - bc = np.abs(c - b) - ac = np.abs(c - a) - ba = np.where(ba > 0.5 * dimensions, ba - dimensions, ba) - bc = np.where(bc > 0.5 * dimensions, bc - dimensions, bc) - ac = np.where(ac > 0.5 * dimensions, ac - dimensions, ac) - - if np.any(ba < 0): - negative_indices = np.where(ba < 0) - negative_values = ba[negative_indices] - raise ValueError( - f"Negative values encountered in 'ba' at indices {negative_indices} " - f"with values {negative_values}. " - f"Cannot take square root of negative values." - ) - - dist_ba = np.sqrt((ba**2).sum(axis=-1)) - - if np.any(bc < 0): - negative_indices = np.where(bc < 0) - negative_values = bc[negative_indices] - raise ValueError( - f"Negative values encountered in 'bc' at indices {negative_indices} " - f"with values {negative_values}. " - f"Cannot take square root of negative values." - ) - - dist_bc = np.sqrt((bc**2).sum(axis=-1)) - - if np.any(ac < 0): - negative_indices = np.where(ac < 0) - negative_values = ac[negative_indices] - raise ValueError( - f"Negative values encountered in 'ac' at indices {negative_indices} " - f"with values {negative_values}. " - f"Cannot take square root of negative values." - ) - - dist_ac = np.sqrt((ac**2).sum(axis=-1)) - - cosine_angle = (dist_ac**2 - dist_bc**2 - dist_ba**2) / (-2 * dist_bc * dist_ba) - - return cosine_angle - - -def com_vectors(coord_list, mass_list, dimensions): - """ - doesnt work - """ - cm = np.zeros(3) - tot_mass = 0 - for coord, mass in zip(coord_list, mass_list): - vec = vector(coord, [0, 0, 0], dimensions) - cm += np.array(vec) * mass - tot_mass += mass - - cm_weighted = np.array(cm) / float(tot_mass) - - return cm_weighted - - -def com(coord_list, mass_list): - """ """ - x_cm = 0 - y_cm = 0 - z_cm = 0 - tot_mass = 0 - for coord, mass in zip(coord_list, mass_list): - x_cm += mass * coord[0] - y_cm += mass * coord[1] - z_cm += mass * coord[2] - tot_mass += mass - x_cm = float(x_cm) / float(tot_mass) - y_cm = float(y_cm) / float(tot_mass) - z_cm = float(z_cm) / float(tot_mass) - - cm = (x_cm, y_cm, z_cm) - - return cm - - -def MOI(cm, coord_list, mass_list): - """ """ - x_cm, y_cm, z_cm = cm[0], cm[1], cm[2] - I_xx = 0 - I_xy = 0 - I_yx = 0 - - I_yy = 0 - I_xz = 0 - I_zx = 0 - - I_zz = 0 - I_yz = 0 - I_zy = 0 - - for coord, mass in zip(coord_list, mass_list): - I_xx += (abs(coord[1] - y_cm) ** 2 + abs(coord[2] - z_cm) ** 2) * mass - I_xy += (coord[0] - x_cm) * (coord[1] - y_cm) * mass - I_yx += (coord[0] - x_cm) * (coord[1] - y_cm) * mass - - I_yy += (abs(coord[0] - x_cm) ** 2 + abs(coord[2] - z_cm) ** 2) * mass - I_xz += (coord[0] - x_cm) * (coord[2] - z_cm) * mass - I_zx += (coord[0] - x_cm) * (coord[2] - z_cm) * mass - - I_zz += (abs(coord[0] - x_cm) ** 2 + abs(coord[1] - y_cm) ** 2) * mass - I_yz += (coord[1] - y_cm) * (coord[2] - z_cm) * mass - I_zy += (coord[1] - y_cm) * (coord[2] - z_cm) * mass - - inertia_tensor = np.array( - [[I_xx, -I_xy, -I_xz], [-I_yx, I_yy, -I_yz], [-I_zx, -I_zy, I_zz]] - ) - return inertia_tensor - - -def UAforceTorque(bonded_atoms_list, F_axes, Fmoi_axes, MI_axes): - """ - Get UA force and torque from previously calculated Faxes and MIaxes - """ - - """ - Faxis_list = (Faxis1, Faxis2, Faxis3) - f_list = (FO, FH1, FH2) - MI_xyz = (axis1MI, axis2MI, axis3MI) - water_torques = torque(c_list, m_list, Faxis_list, f_list, MI_xyz) - """ - - c_list = [] - f_list = [] - forces_summed = np.zeros(3) - m_list = [] - - for atoms in bonded_atoms_list: - c_list.append(atoms.coords) - f_list.append(atoms.forces) - forces_summed += atoms.forces - m_list.append(atoms.mass) - - cm = com(c_list, m_list) - UA_torque = torque_old(cm, c_list, Fmoi_axes, f_list, MI_axes) - - F1 = np.dot(forces_summed, F_axes[0]) - F2 = np.dot(forces_summed, F_axes[1]) - F3 = np.dot(forces_summed, F_axes[2]) - mass_sqrt = sum(m_list) ** 0.5 - UA_force = ( - float(F1) / float(mass_sqrt), - float(F2) / float(mass_sqrt), - float(F3) / float(mass_sqrt), - ) - - return UA_force, UA_torque - - -def torque_old(cm, coord_list, axis_list, force_list, MI_coords): - """ """ - x_cm, y_cm, z_cm = cm[0], cm[1], cm[2] - MI_coords = np.sqrt(MI_coords) # sqrt moi to weight torques - - O_torque = np.array([0, 0, 0]) - H1_torque = np.array([0, 0, 0]) - H2_torque = np.array([0, 0, 0]) - count_atom = 0 - for coord, force in zip(coord_list, force_list): - count_atom += 1 - new_coords = [] - new_forces = [] - for axis in axis_list: - new_c = ( - axis[0] * (coord[0] - x_cm) - + axis[1] * (coord[1] - y_cm) - + axis[2] * (coord[2] - z_cm) - ) - new_f = axis[0] * force[0] + axis[1] * force[1] + axis[2] * force[2] - new_coords.append(new_c) - new_forces.append(new_f) - - torque_list = [] - torquex = float( - new_coords[1] * new_forces[2] - new_coords[2] * new_forces[1] - ) # / float(1e10) - torquey = float( - new_coords[2] * new_forces[0] - new_coords[0] * new_forces[2] - ) # / float(1e10) - torquez = float( - new_coords[0] * new_forces[1] - new_coords[1] * new_forces[0] - ) # / float(1e10) - - if count_atom == 1: - O_torque = (torquex, torquey, torquez) - O_torque = np.divide(O_torque, MI_coords) - elif count_atom == 2: - H1_torque = (torquex, torquey, torquez) - H1_torque = np.divide(H1_torque, MI_coords) - elif count_atom == 3: - H2_torque = (torquex, torquey, torquez) - H2_torque = np.divide(H2_torque, MI_coords) - else: - ("torques outs of range") - break - - W_torque = O_torque + H1_torque + H2_torque - - return W_torque - - -def principalAxesMOI(bonded_atoms_list): - """ - Calculate the principal axes of the MOI matrix, then calc forces - """ - - coord_list = [] - mass_list = [] - forces_summed = np.zeros(3) - - for atom in bonded_atoms_list: - coord_list.append(atom.coords) - mass_list.append(atom.mass) - forces_summed += atom.forces - - cm = com(coord_list, mass_list) - moI = MOI(cm, coord_list, mass_list) - - eigenvalues, eigenvectors = LA.eig(moI) - # different values generated to Jon's code - - transposed = np.transpose(eigenvectors) # turn columns to rows - - min_eigenvalue = abs(eigenvalues[0]) - if eigenvalues[1] < min_eigenvalue: - min_eigenvalue = eigenvalues[1] - if eigenvalues[2] < min_eigenvalue: - min_eigenvalue = eigenvalues[2] - - max_eigenvalue = abs(eigenvalues[0]) - if eigenvalues[1] > max_eigenvalue: - max_eigenvalue = eigenvalues[1] - if eigenvalues[2] > max_eigenvalue: - max_eigenvalue = eigenvalues[2] - - # print (min_eigenvalue * const, max_eigenvalue * const) - # same as Jons - - FMIaxis1 = None - FMIaxis2 = None - FMIaxis3 = None - - axis1MI = None - axis2MI = None - axis3MI = None - - for i in range(0, 3): - if eigenvalues[i] == max_eigenvalue: - FMIaxis1 = transposed[i] - axis1MI = eigenvalues[i] - elif eigenvalues[i] == min_eigenvalue: - FMIaxis3 = transposed[i] - axis3MI = eigenvalues[i] - else: - FMIaxis2 = transposed[i] - axis2MI = eigenvalues[i] - - Fm_axes = [FMIaxis1, FMIaxis2, FMIaxis3] - MI_axes = [axis1MI, axis2MI, axis3MI] - - return Fm_axes, MI_axes - - -def calcAngleWithNearestNonlike(neighbour, nearest, dimensions): - """ - Calc angle in plane and orthog to plane between any molecule - and any other molecule. - """ - - PAxes = np.array(neighbour.WMprincipalAxis[0]) - PAxis = np.array(neighbour.WMprincipalAxis[1]) - COM = np.array(neighbour.WMprincipalAxis[2]) - - # plane_normal = vector(COM, PAxes[0], dimensions) - # plane_normal2 = vector(COM, PAxes[1], dimensions) - # OM_vector = vector(COM, PAxes[2], dimensions) - - plane_normal = PAxes[0] - plane_normal2 = PAxes[1] - OM_vector = PAxes[2] - - OSolute_vector = vector(nearest.coords, COM, dimensions) - - projected_OSolute = projectVectorToPlane(OSolute_vector, plane_normal) - projected_orthog_OSolute = projectVectorToPlane(OSolute_vector, plane_normal2) - - # print ('OM', OM_vector, oxygen.atom_num) - - soluteOMangle = angleBetweenVectors(OM_vector, projected_OSolute) - soluteOorthogAngle = angleBetweenVectors(OM_vector, projected_orthog_OSolute) - - if ( - np.dot(OM_vector, projected_OSolute) > 0 - and np.dot(plane_normal2, projected_OSolute) > 0 - ): - soluteOMangle = (90 - soluteOMangle) + 270 - - if ( - np.dot(-OM_vector, projected_OSolute) > 0 - and np.dot(plane_normal2, projected_OSolute) > 0 - ): - soluteOMangle = (180 - soluteOMangle) + 180 - - if ( - np.dot(-plane_normal, projected_orthog_OSolute) > 0 - and np.dot(-OM_vector, projected_orthog_OSolute) > 0 - ): - soluteOorthogAngle = (180 - soluteOorthogAngle) + 180 - - if ( - np.dot(-plane_normal, projected_orthog_OSolute) > 0 - and np.dot(OM_vector, projected_orthog_OSolute) > 0 - ): - soluteOorthogAngle = (90 - soluteOorthogAngle) + 270 - - # cosAngle = angle(closest_H.coords, oxygen.coords, solute.coords, dimensions) - # arccosAngle = np.arccos(cosAngle) - # angleDegrees = np.degrees(arccosAngle) - # if math.isnan(angleDegrees) is False: - # angleDegrees = int(round(angleDegrees, 0)) - - if math.isnan(soluteOMangle) is False: - soluteOMangle = int(round(soluteOMangle, 0)) - # else: - # print ('Error1: NaN for solute resid: %s, '\ - # 'water resid %s and angle %s.'\ - # % (nearest.resid, neighbour.resid, soluteOMangle)) - # continue - - if math.isnan(soluteOorthogAngle) is False: - soluteOorthogAngle = int(round(soluteOorthogAngle, 0)) - # else: - # print ('Error1: NaN for solute resid: %s, '\ - # 'water resid %s and angle %s.'\ - # % (nearest.resid, neighbour.resid, soluteOorthogAngle)) - # continue - - return soluteOMangle, soluteOorthogAngle - - -def calcWaterSoluteAngle(solute, oxygen, H1, H2, dimensions): - """ - Calc angle in plane and orthog to plane between any water - and any solute atom. - """ - - dist_solute_H1, dist_solute_H2 = None, None - - for atomDist in solute.nearest_all_atom_array: - if atomDist[0] == H1.atom_num: - dist_solute_H1 = atomDist[1] - elif atomDist[0] == H2.atom_num: - dist_solute_H2 = atomDist[1] - else: - continue - - if dist_solute_H1 is None: - dist_solute_H1 = distance(solute.coords, H1.coords, dimensions) - if dist_solute_H2 is None: - dist_solute_H2 = distance(solute.coords, H2.coords, dimensions) - closest_H = None - if dist_solute_H1 < dist_solute_H2: - closest_H = H1 - elif dist_solute_H1 > dist_solute_H2: - closest_H = H2 - else: - closest_H = H1 - - H1O_vector = vector(oxygen.coords, H1.coords, dimensions) - H2O_vector = vector(oxygen.coords, H2.coords, dimensions) - plane_normal = np.cross(H1O_vector, H2O_vector) - M_coords = getMcoords(H1.coords, H2.coords, dimensions) - OM_vector = vector(M_coords, oxygen.coords, dimensions) - plane_normal2 = np.cross(OM_vector, plane_normal) - OSolute_vector = vector(solute.coords, oxygen.coords, dimensions) - projected_OSolute = projectVectorToPlane(OSolute_vector, plane_normal) - projected_orthog_OSolute = projectVectorToPlane(OSolute_vector, plane_normal2) - - # print ('OM', OM_vector, oxygen.atom_num) - - soluteOMangle = angleBetweenVectors(OM_vector, projected_OSolute) - soluteOorthogAngle = angleBetweenVectors(OM_vector, projected_orthog_OSolute) - - if ( - np.dot(OM_vector, projected_OSolute) > 0 - and np.dot(plane_normal2, projected_OSolute) > 0 - ): - soluteOMangle = (90 - soluteOMangle) + 270 - # only go up to 180 degrees rather than 360 as above - # if soluteOMangle > 180: - # soluteOMangle = - - if ( - np.dot(-OM_vector, projected_OSolute) > 0 - and np.dot(plane_normal2, projected_OSolute) > 0 - ): - soluteOMangle = (180 - soluteOMangle) + 180 - - if ( - np.dot(-plane_normal, projected_orthog_OSolute) > 0 - and np.dot(-OM_vector, projected_orthog_OSolute) > 0 - ): - soluteOorthogAngle = (180 - soluteOorthogAngle) + 180 - - if ( - np.dot(-plane_normal, projected_orthog_OSolute) > 0 - and np.dot(OM_vector, projected_orthog_OSolute) > 0 - ): - soluteOorthogAngle = (90 - soluteOorthogAngle) + 270 - - cosAngle = angle(closest_H.coords, oxygen.coords, solute.coords, dimensions) - arccosAngle = np.arccos(cosAngle) - angleDegrees = np.degrees(arccosAngle) - if math.isnan(angleDegrees) is False: - angleDegrees = int(round(angleDegrees, 0)) - soluteOMangle = int(round(soluteOMangle, 0)) - soluteOorthogAngle = int(round(soluteOorthogAngle, 0)) - else: - logging.error( - "NaN for solute resid: %s, water resid %s and angle %s." - % (solute.resid, oxygen.resid, angleDegrees) - ) - - return soluteOMangle, soluteOorthogAngle diff --git a/CodeEntropy/poseidon/extractData/mainClass.py b/CodeEntropy/poseidon/extractData/mainClass.py deleted file mode 100644 index d963887..0000000 --- a/CodeEntropy/poseidon/extractData/mainClass.py +++ /dev/null @@ -1,114 +0,0 @@ -#!/usr/bin/env python - - -class atom_info(object): - """ - Class for all information attributed to atoms and molecules in the system. - """ - - def __init__( - self, - atom_num, - atom_name, - mass, - charge, - resid, - resname, - bonded_to_atom_num, - dihedral_list, - ): # order important - - self.atom_num = atom_num - self.atom_name = atom_name - self.resid = resid - self.resname = resname - self.mass = mass - self.charge = charge - self.bonded_to_atom_num = bonded_to_atom_num - self.dihedral_list = dihedral_list - # list of lists of atom nums in dihedrals, not including Hs - - self.molecule_atomNums = None # list of bonded UA atom nums - self.bonded_to_resids = None # list of resids residue is bonded to - - self.coords = None - self.forces = None - self.nearest_sorted_array = None # list of (atom_nums, distances) - # from nearest to furthest within a set cutoff distance, - # not including Hs - self.nearest_all_atom_array = None # list containing sorted by - # dist (atom_num, dist) for all atoms, including Hs - self.PKenergy = None - self.UA_PKenergy = None - - self.bondedUA_H = None # [num_bondedUAs, num_bonded Hs] - self.nearest_atom = None # HB acceptor all atom info - self.nearest_Hs = [] # HB - self.dist = None # HB distance between donor and acceptor - self.broken = None # HB - self.RAD = None # list of atoms.all_info - self.RAD_dist = None # list of (atoms.all_info, dist from ref) - - self.UAweightedForces = None - self.UAweightedTorques = None - self.MweightedForces = None - self.MweightedTorques = None - self.molecule_UA_Fs = None - self.molecule_UA_Ts = None - self.WMprincipalAxis = None - - self.dihedral_phi_type = None # unique dihedral angles, in degrees - self.nearestAnyNonlike = None # tuple (atom.all_info, dist) - self.RAD_shell_ranked = None # list of 3 tuple - # RAD shell ranked, acceptors ranked, donors ranked by RAD - self.RAD_nAtoms = None - self.hydrationShellRAD_solute = None # nearest solute molecule centric - - def __repr__(self): # return the output below if whole object is printed - return "%s %s %s %s %s %s %s" % ( - self.atom_num, - self.atom_name, - self.mass, - self.charge, - self.resid, - self.resname, - self.bonded_to_atom_num, - ) - - -def clearClass(all_data): - """ - For variables in the class, clear after each frame is analysed. - """ - - def reset(self): - """ - reset attributes in a class. Except those from topology. - This prevents memory overload from each frame. - """ - dic = vars(self) - noEdit = [ - "atom_num", - "atom_name", - "resid", - "resname", - "mass", - "charge", - "bonded_to_atom_num", - "nearest_Hs", - "molecule_atomNums", - "bonded_to_resids", - "dihedral_list", - "bondedUA_H", - ] - for i in dic.keys(): - if i not in noEdit: - dic[i] = None - else: - continue - dic["nearest_Hs"] = [] - - # clear all objects from previous frame - for x in range(0, len(all_data)): - atom = all_data[x] - reset(atom) diff --git a/CodeEntropy/poseidon/extractData/nearestNonlike2.py b/CodeEntropy/poseidon/extractData/nearestNonlike2.py deleted file mode 100644 index efc9b0c..0000000 --- a/CodeEntropy/poseidon/extractData/nearestNonlike2.py +++ /dev/null @@ -1,360 +0,0 @@ -#!/usr/bin/env python - -import logging - -# import math -import sys -from collections import defaultdict # Counter, -from datetime import datetime - -# import MDAnalysis -# import numpy as np -# from MDAnalysis import * - -# from MDAnalysis.analysis import distances - -# from CodeEntropy.poseidon.extractData.generalFunctions import ( -# calcWaterSoluteAngle, -# distance, -# ) - - -def nested_dict(): - return defaultdict(nested_dict) - - -def getShellAssignment(all_data, excludedResnames, dimensions, startTime, verbosePrint): - """ - Assign solvent to solvation shells - """ - - nearestNonlike_dict = {} # need a dictionary of all UAs that are - # assigned as nearest non-likes with values as a list of UAs - - if excludedResnames is None: - excludedResnames = [] - - for x in range(0, len(all_data)): - atom = all_data[x] - if atom.mass > 1.1 and atom.nearest_sorted_array is not None: - # and len(atom.RAD) != 0: - nearestNonlike = None - nearestNonlikeDist = None - for atomDist in atom.nearest_sorted_array: - if nearestNonlikeDist is None: - nonlike = all_data[atomDist[0]] - if ( - nonlike.resid != atom.resid - and nonlike.atom_num not in atom.bonded_to_atom_num - and nonlike.atom_num != atom.atom_num - and nonlike.resname != atom.resname - and nonlike.resname not in excludedResnames - and len(nonlike.RAD) != 0 - ): - nearestNonlike = nonlike - nearestNonlikeDist = atomDist[1] - else: - break - - # When nearest nonlike solute is found - if nearestNonlikeDist is not None: - atom.nearestAnyNonlike = (nearestNonlike, nearestNonlikeDist) - - if nearestNonlike.resid not in nearestNonlike_dict: - nearestNonlike_dict[nearestNonlike.resid] = [[], []] - if nearestNonlike not in nearestNonlike_dict[nearestNonlike.resid][0]: - nearestNonlike_dict[nearestNonlike.resid][0].append(nearestNonlike) - if atom not in nearestNonlike_dict[nearestNonlike.resid][1]: - nearestNonlike_dict[nearestNonlike.resid][1].append(atom) - else: - continue - else: - continue - - verbosePrint("NEAREST NON-LIKE ASSIGNMENT") - verbosePrint(datetime.now() - startTime) - sys.stdout.flush() - - # Part 2: For each nonlike UA with assigned molecules, get prox shell - # This is the slow part of the code - prox shell assignment - - # Find closest nonlike for each molecule, needed below - # so only one atom per molecule is used in prox assignment, rest - # of atoms in molecules are given the same prox as nearest atom - molecule_resid_nearest_dict = {} - for x in range(0, len(all_data)): - atom = all_data[x] - if atom.nearestAnyNonlike is not None: - nearestNonlike = atom.nearestAnyNonlike[0] - dist = atom.nearestAnyNonlike[1] - if atom.resid not in molecule_resid_nearest_dict: - molecule_resid_nearest_dict[atom.resid] = [nearestNonlike.resid, dist] - elif atom.resid in molecule_resid_nearest_dict: - if dist < molecule_resid_nearest_dict[atom.resid][1]: - molecule_resid_nearest_dict[atom.resid] = [ - nearestNonlike.resid, - dist, - ] - else: - continue - else: - continue - - multiple_shellList = [] - multiple_atomNumList = [] - localMols = 0 - - for resid, nonlike_atom_lists in nearestNonlike_dict.items(): - nonlike_list = nonlike_atom_lists[0] - atom_list = nonlike_atom_lists[1] - - for nonlike in nonlike_list: - - if len(nonlike.RAD) != 0: - for neighbour in nonlike.RAD: - X = neighbour - if ( - X in atom_list - and X.atom_num not in multiple_atomNumList - and molecule_resid_nearest_dict[X.resid][0] == resid - ): - for molecule_atom in X.molecule_atomNums: - ma = all_data[molecule_atom] - multiple_shellList.append([ma, 1]) - multiple_atomNumList.append(ma.atom_num) - ma.hydrationShellRAD_solute = 1 - localMols += 1 - else: - continue - - if len(nonlike.RAD) == 0: - logging.warning("no RAD nonlike:", nonlike) - - # Anything beyond the 1st shell is assigned to 2nd - for x in range(0, len(all_data)): - atom = all_data[x] - if atom.hydrationShellRAD_solute is None and atom.mass > 1.1: - atom.hydrationShellRAD_solute = 2 - else: - continue - - -def moleculePositionRankingRAD(all_data, waterTuple, dimensions): - """ - For each UA, rank its RAD shell by proximity to that - references UAs nearest Non-like molecule UA. - num = RAD shell from same molecule type, when nearest nonlike - resid is the same as the reference. - 'X' = when same molecule type has different nearest nonlike resid. - 'RESNAME' = when molecule of different type is in RAD shell. - '0' = closest different type molecule in RAD shell. - (the one its assigned to, its nearest non-like!) - - - Extra: - for each 'X' found in a RAD shell, use this information - to find neighbouring zones. Think voronoi but with zones. - For each molecule that defines a zone, find it's neighbouring zones - from X's in shell type. - - - The following should not occur with new HB defn (HB in RAD shell only): - 'AX' = acceptors not in RAD shell but same molecule type as ref. - 'AS' = acceptors not in RAD shell and diff mol. type to ref. - 'DX' = donor equivalent to 'AX' - 'DS' = donor equivalent to 'AS' - """ - - for x in range(0, len(all_data)): - atom = all_data[x] - if atom.mass > 1.1: - - # part 1: get ranking of RAD Nc according to extended shells. - shell_list_ranked = [] - shell_list_atom_nums = [] - - for neighbour in atom.RAD: - shell_list_atom_nums.append(neighbour.atom_num) - # include Hs in shell_list_atom_nums - for bonded in neighbour.bonded_to_atom_num: - b = all_data[bonded] - if b.mass < 1.1: - shell_list_atom_nums.append(b.atom_num) - else: - continue - - if len(atom.RAD) != 0 and atom.nearestAnyNonlike is not None: - - for neighbour in atom.RAD: - if neighbour.nearestAnyNonlike is not None: - if ( - neighbour.nearestAnyNonlike[0].resid - == atom.nearestAnyNonlike[0].resid - and neighbour.resname == atom.resname - ): - shell_list_ranked.append( - ( - str(neighbour.hydrationShellRAD_solute), - neighbour.atom_num, - ) - ) - - if ( - neighbour.nearestAnyNonlike[0].resid - != atom.nearestAnyNonlike[0].resid - and neighbour.resname == atom.resname - ): - shell_list_ranked.append(("X", neighbour.atom_num)) - - if ( - neighbour.resname != atom.resname - and neighbour.resid == atom.nearestAnyNonlike[0].resid - ): - if neighbour.atom_num == atom.nearestAnyNonlike[0].atom_num: - shell_list_ranked.append(("0", neighbour.atom_num)) - else: - shell_list_ranked.append( - ( - "%s_%s" - % (neighbour.resname, neighbour.atom_name), - neighbour.atom_num, - ) - ) - - if ( - neighbour.resid != atom.nearestAnyNonlike[0].resid - and neighbour.resname != atom.resname - ): - shell_list_ranked.append( - ( - "%s_%s" % (neighbour.resname, neighbour.atom_name), - neighbour.atom_num, - ) - ) - - else: - continue - - # part 2: get HBing info wrt RAD Nc ranking, only for water - - acceptor_list = [] - donor_list = [] - - if atom.hydrationShellRAD_solute == 1 and atom.resname in waterTuple: - - # what Hs is atom accepting from: - if len(atom.nearest_Hs) != 0: - for H in atom.nearest_Hs: - for bonded_atom_num in H.bonded_to_atom_num: - bonded = all_data[bonded_atom_num] - if bonded.mass > 1.1: - X = bonded - if X.atom_num in shell_list_atom_nums: - for rank_num in shell_list_ranked: - rank = rank_num[0] - num = rank_num[1] - if num == X.atom_num: - acceptor_list.append((rank, num)) - else: - continue - elif ( - X.atom_num not in shell_list_atom_nums - and X.resname == atom.resname - ): - acceptor_list.append(("AX", X.atom_num)) - logging.warning("AX") - - elif ( - X.atom_num not in shell_list_atom_nums - and X.resname != atom.resname - ): - acceptor_list.append(("AS", X.atom_num)) - logging.warning("AS") - else: - continue - - # what atoms are bonded Hs donating to: - for bonded_atom_num in atom.bonded_to_atom_num: - H = all_data[bonded_atom_num] - if H.mass < 1.1 and H.nearest_atom is not None: - X = H.nearest_atom - if X.mass > 1.1: - if X.atom_num in shell_list_atom_nums: - for rank_num in shell_list_ranked: - rank = rank_num[0] - num = rank_num[1] - if num == X.atom_num: - # treat a broken HB as a SD for now - if (rank, num) not in donor_list: - donor_list.append((rank, num)) - break - else: - continue - # check donating inside RAD shell - elif ( - X.atom_num not in shell_list_atom_nums - and X.resname == atom.resname - ): - donor_list.append(("DX", X.atom_num)) - logging.warning("DX", atom.atom_name, X.atom_name) - elif ( - X.atom_num not in shell_list_atom_nums - and X.resname != atom.resname - ): - donor_list.append(("DS", X.atom_num)) - logging.warning("DS", atom.atom_name, X.atom_name) - else: - continue - - # dealing with H-H interactions. - elif X.mass < 1.1: - UA_X = None - for bonded_atom_num2 in X.bonded_to_atom_num: - UA_X = all_data[bonded_atom_num2] - if ( - UA_X.mass > 1.1 - and UA_X.atom_num in shell_list_atom_nums - ): - for rank_num in shell_list_ranked: - rank = rank_num[0] - num = rank_num[1] - if num == UA_X.atom_num: - # treat a H-H interaction as SD and - # as unique to H-X interaction. - if ("H", X.atom_num) not in donor_list: - donor_list.append(("H", X.atom_num)) - else: - continue - else: - continue - elif ( - UA_X.mass > 1.1 - and UA_X.atom_num not in shell_list_atom_nums - ): - donor_list.append(("DH", X.atom_num)) - logging.warning( - "DH", atom.atom_name, UA_X.atom_name - ) - else: - continue - - shell_list_ranked.sort(key=lambda x: x[0]) - acceptor_list.sort(key=lambda x: x[0]) - donor_list.sort(key=lambda x: x[0]) - - shell = [i[0] for i in shell_list_ranked] - acceptors = [i[0] for i in acceptor_list] - donors = [i[0] for i in donor_list] - - atom.RAD_shell_ranked = (shell, acceptors, donors) - - # create a list for RAD neighbours and relative - # distances (QQ/r^2) - - RAD_nAtoms = [] - for neighbour in atom.RAD: - RAD_nAtoms.append( - [neighbour.atom_name, neighbour.resid, neighbour.resname] - ) - - atom.RAD_nAtoms = RAD_nAtoms diff --git a/CodeEntropy/poseidon/extractData/outputFiles.py b/CodeEntropy/poseidon/extractData/outputFiles.py deleted file mode 100644 index 3789fa2..0000000 --- a/CodeEntropy/poseidon/extractData/outputFiles.py +++ /dev/null @@ -1,193 +0,0 @@ -#!/usr/bin/env python - -# import math -# import sys - -# import numpy as np - -# from CodeEntropy.poseidon.extractData.generalFunctions import ( -# calcAngleWithNearestNonlike, -# ) - - -def moleculeObjectPopulation(all_data, allMoleculeList, frame, dimensions): - """ - For all molecules in the system, populate a list with required - info for further analysis across any trajectory frame. - """ - - for x in range(0, len(all_data)): - atom = all_data[x] - if atom.mass > 1.1: - nearestInfo = [] - if atom.nearestAnyNonlike is not None: - nearestInfo = [ - atom.nearestAnyNonlike[0].atom_name, - atom.nearestAnyNonlike[0].resname, - atom.nearestAnyNonlike[0].resid, - round(atom.nearestAnyNonlike[1], 3), # dist - ] - - allMoleculeList.append( - [ - frame, # 0 - atom.atom_name, # 1 - atom.atom_num, # 2 - atom.resname, # 3 - atom.resid, # 4 - atom.bondedUA_H, # 5 - atom.molecule_atomNums, # 6 - nearestInfo, # 7 - atom.hydrationShellRAD_solute, # 8 - atom.RAD_shell_ranked, # 9 - atom.MweightedForces, # 10 - atom.MweightedTorques, # 11 - atom.UA_PKenergy, # 12 - atom.UAweightedForces, # 13 - atom.UAweightedTorques, # 14 - atom.dihedral_phi_type, # 15 - atom.molecule_UA_Fs, # 16 - atom.molecule_UA_Ts, # 17 - atom.RAD_nAtoms, # 18 - ] - ) - - else: - continue - - return allMoleculeList - - -def pdbGenerate(all_data, fileName, dimensions): - """ - Output a pdb file to visualise different waters in system. - Might have to also make an equivalent .csv file so that - we have an easier time analysing the files. - But append to one global file rather than each frame, - so need a column of frame number for analysis later. - pdb file pure for producing images with vmd, - but could create these later anyway? - Then we can choose our beta column? - """ - - data = open(fileName + "_solProx_dist.pdb", "w") - # solute resid, int(dist) from solute resid - - data3 = open(fileName + "_RADHBlengths.pdb", "w") - # Nc, N_acceptors/N_donors - - data6 = open(fileName + "_solProx_RAD.pdb", "w") - # solute resid, int(dist) from solute resid - - file_list = [data, data3, data6] - - for x in range(0, len(all_data)): - atom = all_data[x] - - try: - a = int(atom.nearestAnyNonlike[0].resid) - except TypeError: - a = 0 - try: - b = int(atom.nearestAnyNonlike[1]) - except TypeError: - b = 0 - try: - e = len(atom.RAD) - except TypeError: - e = 0 - try: - As = len(atom.RAD_shell_ranked[1]) - Ds = len(atom.RAD_shell_ranked[2]) - DAhh_type = "%s%s" % (Ds, As) - f = str(DAhh_type) - except TypeError: - f = "00" - try: - i = int(atom.hydrationShellRAD_solute) - except (TypeError, ValueError): - i = 0 - - data.write( - "\n".join( - [ - "{:6s}{:5d} {:^4s}{:1s}{:3s}" - " {:1s}{:4d}{:1s} {:8.3f}{:8.3f}{:8.3f}" - "{:6.0f}{:6.0f} {:>2s}{:2s}".format( - "ATOM", - int(str(atom.atom_num)[0:5]), - atom.atom_name, - " ", - atom.resname[0:3], - "A", - int(str(atom.resid)[0:4]), - " ", - float(atom.coords[0]), - float(atom.coords[1]), - float(atom.coords[2]), - a, - b, - " ", - " ", - ) - ] - ) - + "\n" - ) - - data3.write( - "\n".join( - [ - "{:6s}{:5d} {:^4s}{:1s}{:3s}" - " {:1s}{:4d}{:1s} {:8.3f}{:8.3f}{:8.3f}" - "{:6.0f} {:2s} {:>2s}{:2s}".format( - "ATOM", - int(str(atom.atom_num)[0:5]), - atom.atom_name, - " ", - atom.resname[0:3], - "A", - int(str(atom.resid)[0:4]), - " ", - float(atom.coords[0]), - float(atom.coords[1]), - float(atom.coords[2]), - e, - f, - " ", - " ", - ) - ] - ) - + "\n" - ) - - data6.write( - "\n".join( - [ - "{:6s}{:5d} {:^4s}{:1s}{:3s}" - " {:1s}{:4d}{:1s} {:8.3f}{:8.3f}{:8.3f}" - "{:6.0f}{:6.0f} {:>2s}{:2s}".format( - "ATOM", - int(str(atom.atom_num)[0:5]), - atom.atom_name, - " ", - atom.resname[0:3], - "A", - int(str(atom.resid)[0:4]), - " ", - float(atom.coords[0]), - float(atom.coords[1]), - float(atom.coords[2]), - a, - i, - " ", - " ", - ) - ] - ) - + "\n" - ) - - for d in file_list: - d.close() diff --git a/CodeEntropy/poseidon/extractData/readFiles.py b/CodeEntropy/poseidon/extractData/readFiles.py deleted file mode 100644 index ce7c82c..0000000 --- a/CodeEntropy/poseidon/extractData/readFiles.py +++ /dev/null @@ -1,259 +0,0 @@ -#!/usr/bin/env python - -import logging - -# import math -# import operator -import sys -from collections import defaultdict # Counter, -from datetime import datetime - -import MDAnalysis -import numpy as np - -# from CodeEntropy.poseidon.extractData.generalFunctions import * -from CodeEntropy.poseidon.extractData.mainClass import atom_info - -# from numpy import linalg as LA - - -# create nested dict in one go -def nested_dict(): - return defaultdict(nested_dict) - - -def populateTopology(container, all_data, waterTuple): - """ - After reading in topologies using MDAnalysis, - relavant information from topology files are saved into a global - class object. Properties are read in differently depeding on what - file types were read in. - """ - - all_resids = [] - mol = None - - for tp in container.atoms: - bonded_atom_nums = [] - for b in tp.bonds: - for x in b: - if x.index != tp.index: - bonded_atom_nums.append(int(x.index)) - else: - continue - - try: - dihedral_list = [] - if tp.mass > 1.1: - for dih in tp.dihedrals: - diha_list = [] - for di in dih: - if di.mass > 1.1: - diha_list.append(int(di.index)) - else: - break - - if len(diha_list) == 4: - dihedral_list.append(diha_list) - else: - continue - except AttributeError: - dihedral_list = [] - - inf = atom_info( - int(tp.index), - tp.name, - tp.mass, - tp.charge, - int(tp.resid), - tp.resname, - bonded_atom_nums, - dihedral_list, - ) - - all_data.append(inf) - - # Get and populate UA and molecule level information - molecule_dict = nested_dict() - molecule_resids_dict = nested_dict() - for x in range(0, len(all_data)): - atom = all_data[x] - # print(f"x = {x}") - if atom.mass > 1.1: - heavy_bonded = [] - H_bonded = [] - for bonded_atom in atom.bonded_to_atom_num: - # print(f"bonded_atom = {bonded_atom}") - bonded = all_data[bonded_atom] - if bonded.mass > 1.1: - heavy_bonded.append(bonded) - elif bonded.mass < 1.1: - H_bonded.append(bonded) - else: - continue - - bonded_atoms_list = [atom] + heavy_bonded + H_bonded - atom.bondedUA_H = [len(heavy_bonded), len(H_bonded)] - - if atom.resid not in molecule_dict: - molecule_dict[atom.resid] = [] - molecule_resids_dict[atom.resid] = [] - if atom.atom_num not in molecule_dict[atom.resid]: - molecule_dict[atom.resid].append(atom.atom_num) - for bonded2 in heavy_bonded: - if bonded2.resid not in molecule_resids_dict[atom.resid]: - molecule_resids_dict[atom.resid].append(bonded2.resid) - else: - continue - - for x in range(0, len(all_data)): - atom = all_data[x] - if atom.mass > 1.1: - atom.molecule_atomNums = sorted(molecule_dict[atom.resid]) - atom.bonded_to_resids = sorted(molecule_resids_dict[atom.resid]) - else: - continue - - -def getCoordsForces(container, all_data, dimensions, frame, startTime, verbosePrint): - """ - Read in coordinate and force trajectories and populate mainClass. - """ - - t = container.trajectory[frame] - dimensions = np.array(t.dimensions[0:3]) - print(t) - - for x in range(0, len(t)): - crds = np.array(t[x]) - all_data[x].coords = crds - frcs = np.array(t.forces[x]) - all_data[x].forces = frcs - - verbosePrint("COORDS-FORCES") - verbosePrint(datetime.now() - startTime) - sys.stdout.flush() - - return all_data, dimensions - - -# #Energy is not needed -# def populateEnergy(container, all_data, dimensions, frame, startTime, -# verbosePrint): -# ''' -# read in energies from lammps input file. -# ''' - -# for e in container.trajectory: -# if e.frame == frame: -# dimensions = np.array(e.dimensions[0:3]) -# for x in range(0, len(e)): -# energy = np.array(e.velocities[x]) -# all_data[x].PKenergy = [energy[0], energy[1]] -# break -# else: -# continue - - -# verbosePrint('ENERGIES') -# verbosePrint(datetime.now() - startTime) -# sys.stdout.flush() - - -# def UAEnergyGroup(all_data): -# ''' -# For energy on each atom, group together for each UA -# and sum -# ''' - -# for x in range(0, len(all_data)): -# atom = all_data[x] -# if atom.mass > 1.1: -# heavy_bonded = [] -# H_bonded = [] -# for bonded_atom in atom.bonded_to_atom_num: -# bonded = all_data[bonded_atom] -# if bonded.mass > 1.1: -# heavy_bonded.append(bonded) -# elif bonded.mass < 1.1: -# H_bonded.append(bonded) -# else: -# continue - -# bonded_atoms_list = [atom] + heavy_bonded + H_bonded -# atom.bondedUA_H = [len(heavy_bonded), len(H_bonded)] -# UA_atoms_list = [atom] + H_bonded - -# UA_PE_list = [] -# UA_KE_list = [] -# for A in UA_atoms_list: -# if A.PKenergy != None: -# UA_PE_list.append(A.PKenergy[0]) -# UA_KE_list.append(A.PKenergy[1]) -# else: -# continue - -# if len(UA_PE_list) != 0: -# UA_PE = round(sum(UA_PE_list), 3) -# UA_KE = round(sum(UA_KE_list), 3) -# atom.UA_PKenergy = [UA_PE, UA_KE] - - -def getDistArray( - atom, - all_data, - traj, - max_cutoff, - dimensions, - neighbour_coords, - startTime, - verbosePrint, -): - """ - Find the NN list of an atom - Important to use coords directly from MDAnalysis to run NN calc - """ - - atom_coords = traj[atom.atom_num] - - # added a small min cutoff to stop zero distance - array1, array2 = MDAnalysis.lib.distances.capped_distance( - atom_coords, - neighbour_coords, - max_cutoff=max_cutoff, - min_cutoff=None, - box=traj.dimensions, - method=None, - return_distances=True, - ) - - try: - array1, array2 = zip(*sorted(zip(array2, array1), key=lambda x: x[0])) - - except ValueError: - logging.error( - "Bad Arrays for Coordinate/ Atom Number " "Nearest Neighbour Assignments" - ) - - atomNumList = [] - allAtomList = [] - for atoms, dist in zip(array2, array1): - near = atoms[1] - atom_num = all_data[near].atom_num - atom_resid = all_data[near].resid - atom_mass = all_data[near].mass - allAtomList.append((near, dist)) - - # atom_resid != all_data[x].resid removed for quartz - # surface that is all the same resid. - if ( - atom_num != atom.atom_num - and atom_num not in atom.bonded_to_atom_num - and float(atom_mass) > 1.1 - ): - atomNumList.append((near, dist)) - else: - continue - - atom.nearest_sorted_array = atomNumList - atom.nearest_all_atom_array = allAtomList diff --git a/Example/create_new_universe.py b/Example/create_new_universe.py index 12d56ee..b89f1d5 100644 --- a/Example/create_new_universe.py +++ b/Example/create_new_universe.py @@ -69,7 +69,7 @@ def main(): # you can also slice the trajectories select.write("data.trr", frames=u2.trajectory[::2]) # reading data - u_new = mda.Universe("molecules.prmtop", "data.trr") + # u_new = mda.Universe("molecules.prmtop", "data.trr") if __name__ == "__main__": diff --git a/pyproject.toml b/pyproject.toml index d9deb7a..c339c89 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -73,5 +73,4 @@ max-line-length = 88 extend-select = "B950" extend-ignore = [ "E203", # whitespace before `:` - "F841" # local variable is assigned but never used ] \ No newline at end of file