13726590978

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

Build # 13726590978

Build Type

push

github

Committed by

web-flow

Commit Message

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

Implement CI pipeline pre-commit hooks

Coverage Stats

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

53 of 3071 relevant lines covered (1.73%)

0.05 hits per line

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

0.0

/CodeEntropy/poseidon/extractData/forceTorques.py

#!/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

1	#!/usr/bin/env python
2
3	# import logging
4	# import math
5	# import sys
NEW 6	from collections import defaultdict	×
7
8	import numpy as np	×
9	from numpy import linalg as LA	×
10
NEW 11	from CodeEntropy.poseidon.extractData.generalFunctions import MOI, com, distance, vector	×
12
13
14	# create nested dict in one go
NEW 15	def nested_dict():	×
NEW 16	return defaultdict(nested_dict)	×
17
18
19	def calculateFTMatrix(all_data, dimensions):	×
20	"""
21	Consider force/torque calculations as discussed in Ali 2019 paper.
22	UA level and molecule level.
23	"""
24
25	resid_list = []	×
26	for x in range(0, len(all_data)):	×
27	atom = all_data[x]	×
28	if atom.mass > 1.1:	×
NEW 29	if atom.resid not in resid_list: # iterate though resids once	×
30	resid_list.append(atom.resid)	×
NEW 31	all_molecule_atoms = [] # list of heavy atoms and Hs	×
NEW 32	UAs_list = [] # list of lists [HA + Hs]	×
33
NEW 34	for num in atom.molecule_atomNums: # atom nums in molecule	×
35	UA = all_data[num]	×
36
37	all_molecule_atoms.append(UA)	×
38	heavy_bonded = []	×
39	H_bonded = []	×
40
41	for bonded_atom in UA.bonded_to_atom_num:	×
42	bonded = all_data[bonded_atom]	×
NEW 43	if bonded.mass > 1.1 and bonded.atom_num != UA.atom_num:	×
44	heavy_bonded.append(bonded)	×
45
NEW 46	if bonded.mass < 1.1 and bonded.atom_num != UA.atom_num:	×
47	H_bonded.append(bonded)	×
48	all_molecule_atoms.append(bonded)	×
49
50	UA_atoms_list = [UA] + H_bonded	×
51	UAs_list.append(UA_atoms_list)	×
52	bonded_atoms_list = [UA] + heavy_bonded + H_bonded	×
53	UA.bondedUA_H = [len(heavy_bonded), len(H_bonded)]	×
54
55	molecule_coords = []	×
56	molecule_masses = []	×
57	for a in all_molecule_atoms:	×
58	molecule_coords.append(a.coords)	×
59	molecule_masses.append(a.mass)	×
60	molecule_COM = com(molecule_coords, molecule_masses)	×
61
62	# **** MOLECULE LEVEL ****
63	# if molecule is only one UA with bonded Hs,
64	# this works for one UA molecules
NEW 65	if (	×
66	len(atom.molecule_atomNums) == 1
67	and len(atom.bonded_to_atom_num) > 0
68	):
69	WM_principal_axes, WM_MI_axis = principalAxesMOI(	×
70	all_data, all_molecule_atoms, molecule_COM, Hs=True
71	)
72	# use atom list containing Hs for P axes inc. Hs
NEW 73	WM_force, WM_torque = rotateFT(	×
74	all_data,
75	all_molecule_atoms,
76	WM_principal_axes,
77	WM_principal_axes,
78	WM_MI_axis,
79	molecule_COM,
80	molecule_COM,
81	Hs=True,
82	)
83
84	#
85	WM_force = np.outer(WM_force, WM_force)	×
86	WM_torque = np.outer(WM_torque, WM_torque)	×
87	atom.MweightedForces = np.round(np.divide(WM_force, 4), 3)	×
88	atom.MweightedTorques = np.round(np.divide(WM_torque, 4), 3)	×
NEW 89	atom.WMprincipalAxis = WM_principal_axes, WM_MI_axis, molecule_COM	×
90
91	# check this with Richard, do you not halve here?
92	# For Jon's code, you do halve here
93	# even though molecule == UA level
94	atom.UAweightedForces = np.round(np.divide(WM_force, 4), 3)	×
95	atom.UAweightedTorques = np.round(np.divide(WM_torque, 4), 3)	×
96	atom.molecule_UA_Fs = np.round(np.divide(WM_force, 4), 3)	×
97	atom.molecule_UA_Ts = np.round(np.divide(WM_torque, 4), 3)	×
98	#
99
NEW 100	"""	×
101	atom.MweightedForces = WM_force
102	atom.MweightedTorques = WM_torque
103	atom.UAweightedForces = WM_force
104	atom.UAweightedTorques = WM_torque
105	"""
106
107	# if molecule contains more than one UA
108	# forces work, torques work if Hs = True
109	if len(atom.molecule_atomNums) > 1:	×
110
111	WM_principal_axes, WM_MI_axis = principalAxesMOI(	×
112	all_data, all_molecule_atoms, molecule_COM, Hs=False
113	)
114	# only use heavy atom list for P axes calc
115	# we don't consider Hs as they result
116	# in larger torques as Hs are further away
117
NEW 118	WM_force, WM_torque = rotateFT(	×
119	all_data,
120	all_molecule_atoms,
121	WM_principal_axes,
122	WM_principal_axes,
123	WM_MI_axis,
124	molecule_COM,
125	molecule_COM,
126	Hs=True,
127	)
128
129	#
130	# WM_FT = np.concatenate((WM_force, WM_torque), axis=None)
131	# WM_FT = np.outer(WM_FT, WM_FT)
132	WM_force = np.outer(WM_force, WM_force)	×
133	WM_torque = np.outer(WM_torque, WM_torque)	×
134	atom.MweightedForces = np.round(np.divide(WM_force, 4), 3)	×
135	atom.MweightedTorques = np.round(np.divide(WM_torque, 4), 3)	×
NEW 136	atom.WMprincipalAxis = WM_principal_axes, WM_MI_axis, molecule_COM	×
137
138	atom.molecule_UA_Fs = np.round(np.divide(WM_force, 4), 3)	×
139	atom.molecule_UA_Ts = np.round(np.divide(WM_torque, 4), 3)	×
140	#
141
142	# atom.MweightedForces = WM_force
143	# atom.MweightedTorques = WM_torque
144
145	# **** UNITED-ATOM LEVEL ****
146	# works!
147	UA_F_list = []	×
148	UA_T_list = []	×
149	for UA in UAs_list:	×
NEW 150	XX_principal_axes, UA_MI_axis = None, None	×
151
152	if UA[0].mass > 1.1:	×
NEW 153	bonded_HAs = [UA[0]] # inc itself	×
154	bonded_Hs = []	×
155	for b in UA[0].bonded_to_atom_num:	×
156	bonded = all_data[b]	×
157	if bonded.mass < 1.1:	×
158	bonded_Hs.append(bonded)	×
159	elif bonded.mass > 1.1:	×
160	bonded_HAs.append(bonded)	×
161	else:
162	continue	×
163
164	if len(bonded_HAs) == 2 and len(bonded_Hs) > 0:	×
165	# num HAs including itself
166	UA_COM = bonded_HAs[0].coords	×
167
168	R = bonded_HAs[0]	×
169	X1 = bonded_HAs[1]	×
170	H = bonded_Hs[0]	×
171
NEW 172	RX1_vector = vector(R.coords, X1.coords, dimensions)	×
NEW 173	RX1_dist = distance(R.coords, X1.coords, dimensions)	×
174	Paxis1 = np.divide(RX1_vector, RX1_dist)	×
175
NEW 176	RH_vector = vector(R.coords, H.coords, dimensions)	×
NEW 177	RH_dist = distance(R.coords, H.coords, dimensions)	×
178	RH_vector_norm = np.divide(RH_vector, RH_dist)	×
179
180	Paxis2 = np.cross(Paxis1, RH_vector_norm)	×
181	Paxis3 = np.cross(Paxis1, Paxis2)	×
182
183	XX_principal_axes = [Paxis1, Paxis2, Paxis3]	×
184
NEW 185	UA_MI_axis, XX_principal_axes = UA_MOI(	×
186	all_data,
187	UA,
188	UA_COM,
189	XX_principal_axes,
190	dimensions,
191	Hs=True,
192	)
193
194	if len(bonded_Hs) == 1:	×
NEW 195	MIx, MIy, MIz = (	×
196	UA_MI_axis[0],
197	UA_MI_axis[1],
198	UA_MI_axis[2],
199	)
200	if MIx < MIy and MIx < MIz:	×
201	MIx = 0	×
202	if MIy < MIx and MIy < MIz:	×
203	MIy = 0	×
204	if MIz < MIx and MIz < MIy:	×
205	MIz = 0	×
206
207	UA_MI_axis = [MIx, MIy, MIz]	×
208
209	elif len(bonded_HAs) > 2:	×
210	# num HAs including itself
211	UA_COM = bonded_HAs[0].coords	×
212
213	R = bonded_HAs[0]	×
214	X1 = bonded_HAs[1]	×
215	X2 = bonded_HAs[2]	×
216
NEW 217	Paxis1 = np.zeros(3)	×
218	# average of all HA covalent bond vectors
219	for HA in bonded_HAs[1:]:	×
NEW 220	Paxis1 += vector(R.coords, HA.coords, dimensions)	×
221
NEW 222	X1X2_vector = vector(X1.coords, X2.coords, dimensions)	×
NEW 223	X1X2_dist = distance(X1.coords, X2.coords, dimensions)	×
224	X1X2_norm = np.divide(X1X2_vector, X1X2_dist)	×
225	Paxis2 = np.cross(X1X2_norm, Paxis1)	×
226
227	Paxis3 = np.cross(Paxis1, Paxis2)	×
228
229	XX_principal_axes = [Paxis1, Paxis2, Paxis3]	×
230
NEW 231	UA_MI_axis, XX_principal_axes = UA_MOI(	×
232	all_data,
233	UA,
234	UA_COM,
235	XX_principal_axes,
236	dimensions,
237	Hs=True,
238	)
239
240	elif len(bonded_HAs) == 2 and len(bonded_Hs) == 0:	×
241	# num use arbitrary coords for 2nd vector
242	UA_COM = bonded_HAs[0].coords	×
243
244	R = bonded_HAs[0]	×
245	X1 = bonded_HAs[1]	×
246
NEW 247	RX1_vector = vector(R.coords, X1.coords, dimensions)	×
NEW 248	RX1_dist = distance(R.coords, X1.coords, dimensions)	×
249	Paxis1 = np.divide(RX1_vector, RX1_dist)	×
250
NEW 251	RZ_vector = vector(R.coords, [0, 0, 0], dimensions)	×
NEW 252	RZ_dist = distance(R.coords, [0, 0, 0], dimensions)	×
253	RZ_vector_norm = np.divide(RZ_vector, RZ_dist)	×
254
255	Paxis2 = np.cross(Paxis1, RZ_vector_norm)	×
256	Paxis3 = np.cross(Paxis1, Paxis2)	×
257
258	XX_principal_axes = [Paxis1, Paxis2, Paxis3]	×
259
NEW 260	UA_MI_axis, XX_principal_axes = UA_MOI(	×
261	all_data,
262	UA,
263	UA_COM,
264	XX_principal_axes,
265	dimensions,
266	Hs=True,
267	)
268
269	if len(bonded_Hs) == 0:	×
NEW 270	MIx, MIy, MIz = (	×
271	UA_MI_axis[0],
272	UA_MI_axis[1],
273	UA_MI_axis[2],
274	)
275	if MIx < MIy and MIx < MIz:	×
276	MIx = 0	×
277	if MIy < MIx and MIy < MIz:	×
278	MIy = 0	×
279	if MIz < MIx and MIz < MIy:	×
280	MIz = 0	×
281
282	UA_MI_axis = [MIx, MIy, MIz]	×
283
284	else:
285	continue	×
286
287	if UA[0].mass < 1.1:	×
NEW 288	print("Error: no HA in UA")	×
289
NEW 290	UA_force, UA_torque = rotateFT(	×
291	all_data,
292	UA,
293	WM_principal_axes,
294	XX_principal_axes,
295	UA_MI_axis,
296	molecule_COM,
297	UA_COM,
298	Hs=True,
299	)
300	# torques use HA com,
301	# forces use molecule COM
302
303	UA_F_list += [UA_force]	×
304	UA_T_list += [UA_torque]	×
305
306	#
307	UA_force = np.outer(UA_force, UA_force)	×
308	UA_torque = np.outer(UA_torque, UA_torque)	×
NEW 309	UA[0].UAweightedForces = np.round(np.divide(UA_force, 4), 3)	×
310	# UA[0].UAweightedForces = UA_force
311	# 2019 paper, dont halve forces
NEW 312	UA[0].UAweightedTorques = np.round(np.divide(UA_torque, 4), 3)	×
313	#
314
315	# UA[0].UAweightedForces = UA_force
316	# UA[0].UAweightedTorques = UA_torque
317
318	# dont halve forces in molecule UAs (2019 paper),
319	# divide by 4 to compare with jons output
320	molecule_UA_Fs = np.concatenate(UA_F_list, axis=None)	×
321	molecule_UA_Fs = np.outer(molecule_UA_Fs, molecule_UA_Fs)	×
322	# atom.molecule_UA_Fs = np.divide(molecule_UA_Fs, 4)
323	atom.molecule_UA_Fs = np.round(molecule_UA_Fs)	×
324
325	# halve torques in molecule UAs
326	molecule_UA_Ts = np.concatenate(UA_T_list, axis=None)	×
327	molecule_UA_Ts = np.outer(molecule_UA_Ts, molecule_UA_Ts)	×
328	atom.molecule_UA_Ts = np.round(np.divide(molecule_UA_Ts, 4), 3)	×
329
330	# if molecule is monoatomic (one UA and no bonded Hs)
NEW 331	if (	×
332	len(atom.molecule_atomNums) == 1
333	and len(atom.bonded_to_atom_num) == 0
334	):
NEW 335	mass_sqrt = float(atom.mass**0.5)	×
NEW 336	forces_sqrt = [	×
337	float(atom.forces[0]) / mass_sqrt,
338	float(atom.forces[1]) / mass_sqrt,
339	float(atom.forces[2]) / mass_sqrt,
340	]
341	forces_sqrt = np.sort(forces_sqrt)	×
342
343	#
344	UA_force = np.outer(forces_sqrt, forces_sqrt)	×
345	UA_torque = np.outer([0, 0, 0], [0, 0, 0])	×
346	atom.UAweightedForces = np.round(np.divide(UA_force, 4), 3)	×
NEW 347	atom.UAweightedTorques = np.round(UA_torque, 3) # zero	×
348	atom.MweightedForces = np.round(np.divide(UA_force, 4), 3)	×
NEW 349	atom.MweightedTorques = np.round(UA_torque, 3) # zero	×
NEW 350	atom.WMprincipalAxis = (	×
351	[[0, 0, 0], [0, 0, 0], [0, 0, 0]],
352	[0, 0, 0],
353	[0, 0, 0],
354	)
355
356	atom.molecule_UA_Fs = np.round(np.divide(UA_force, 4), 3)	×
357	atom.molecule_UA_Ts = np.round(np.divide(UA_torque, 4), 3)	×
358	#
359
NEW 360	"""	×
361	atom.UAweightedForces = forces_sqrt
362	atom.UAweightedTorques = [0, 0, 0] #zero
363	atom.MweightedForces = forces_sqrt
364	atom.MweightedTorques = [0, 0, 0] #zero
365	"""
366
367
368	def UA_MOI(all_data, bonded_atoms_list, COM, PIaxes, dimensions, Hs):	×
369	"""
370	Get MOI axis for UA level, where eigenvalues and vectors are
371	not used
372	"""
373
374	coord_list = []	×
375	mass_list = []	×
376
377	for atom in bonded_atoms_list:	×
NEW 378	if Hs is True:	×
379	# print (atom.atom_name)
380	coord_list.append(atom.coords)	×
381	mass_list.append(atom.mass)	×
NEW 382	elif Hs is False:	×
383	if atom.mass > 1.1:	×
NEW 384	UA_mass = 0 # add mass of HA and bonded Hs	×
385	# print (atom.atom_name)
386	UA_mass += atom.mass	×
NEW 387	coord_list.append(atom.coords) # only include coords of UA	×
388	for h in atom.bonded_to_atom_num:	×
389	H = all_data[h]	×
390	if H.mass < 1.1:	×
391	UA_mass += H.mass	×
392	else:
393	continue	×
NEW 394	mass_list.append(UA_mass)	×
395	# include mass of Hs on heavy atom mass
396	else:
397	continue	×
398
399	# sorting out PIaxes for MoI for UA fragment
400
401	modPIx = PIaxes[0][0] 2 + PIaxes[0][1] 2 + PIaxes[0][2] ** 2	×
402	modPIy = PIaxes[1][0] 2 + PIaxes[1][1] 2 + PIaxes[1][2] ** 2	×
403	modPIz = PIaxes[2][0] 2 + PIaxes[2][1] 2 + PIaxes[2][2] ** 2	×
404
NEW 405	PIaxes[0][0] = PIaxes[0][0] / (modPIx**0.5)	×
NEW 406	PIaxes[0][1] = PIaxes[0][1] / (modPIx**0.5)	×
NEW 407	PIaxes[0][2] = PIaxes[0][2] / (modPIx**0.5)	×
408
NEW 409	PIaxes[1][0] = PIaxes[1][0] / (modPIy**0.5)	×
NEW 410	PIaxes[1][1] = PIaxes[1][1] / (modPIy**0.5)	×
NEW 411	PIaxes[1][2] = PIaxes[1][2] / (modPIy**0.5)	×
412
NEW 413	PIaxes[2][0] = PIaxes[2][0] / (modPIz**0.5)	×
NEW 414	PIaxes[2][1] = PIaxes[2][1] / (modPIz**0.5)	×
NEW 415	PIaxes[2][2] = PIaxes[2][2] / (modPIz**0.5)	×
416
417	# get dot product of Paxis1 and CoM->atom1 vect
418	# will just be [0,0,0]
419	RRaxis = vector(coord_list[0], COM, dimensions)	×
420	# flip each Paxis if its pointing out of UA
421	dotProd1 = np.dot(np.array(PIaxes[0]), RRaxis)	×
422	if dotProd1 < 0:	×
423	PIaxes[0][0] = -PIaxes[0][0]	×
424	PIaxes[0][1] = -PIaxes[0][1]	×
425	PIaxes[0][2] = -PIaxes[0][2]	×
426
427	dotProd2 = np.dot(np.array(PIaxes[1]), RRaxis)	×
428	if dotProd2 < 0:	×
429	PIaxes[1][0] = -PIaxes[1][0]	×
430	PIaxes[1][1] = -PIaxes[1][1]	×
431	PIaxes[1][2] = -PIaxes[1][2]	×
432
433	dotProd3 = np.dot(np.array(PIaxes[2]), RRaxis)	×
434	if dotProd3 < 0:	×
435	PIaxes[2][0] = -PIaxes[2][0]	×
436	PIaxes[2][1] = -PIaxes[2][1]	×
437	PIaxes[2][2] = -PIaxes[2][2]	×
438
439	# MI = moment of inertia
NEW 440	axis1MI = 0 # x	×
NEW 441	axis2MI = 0 # y	×
NEW 442	axis3MI = 0 # z	×
443	#
444
445	for coord, mass in zip(coord_list, mass_list):	×
446	dx = coord[0] - COM[0]	×
447	dy = coord[1] - COM[1]	×
448	dz = coord[2] - COM[2]	×
449
450	PIaxes_xx = PIaxes[0][1] * dz - PIaxes[0][2] * dy	×
451	PIaxes_xy = PIaxes[0][2] * dx - PIaxes[0][0] * dz	×
452	PIaxes_xz = PIaxes[0][0] * dy - PIaxes[0][1] * dx	×
453
454	PIaxes_yx = PIaxes[1][1] * dz - PIaxes[1][2] * dy	×
455	PIaxes_yy = PIaxes[1][2] * dx - PIaxes[1][0] * dz	×
456	PIaxes_yz = PIaxes[1][0] * dy - PIaxes[1][1] * dx	×
457
458	PIaxes_zx = PIaxes[2][1] * dz - PIaxes[2][2] * dy	×
459	PIaxes_zy = PIaxes[2][2] * dx - PIaxes[2][0] * dz	×
460	PIaxes_zz = PIaxes[2][0] * dy - PIaxes[2][1] * dx	×
461
NEW 462	daxisx = (PIaxes_xx2 + PIaxes_xy2 + PIaxes_xz2) 0.5	×
NEW 463	daxisy = (PIaxes_yx2 + PIaxes_yy2 + PIaxes_yz2) 0.5	×
NEW 464	daxisz = (PIaxes_zx2 + PIaxes_zy2 + PIaxes_zz2) 0.5	×
465
NEW 466	axis1MI += (daxisx*2) mass	×
NEW 467	axis2MI += (daxisy*2) mass	×
NEW 468	axis3MI += (daxisz*2) mass	×
469
470	UA_moI = [axis1MI, axis2MI, axis3MI]	×
471
472	return UA_moI, PIaxes	×
473
474
475	def principalAxesMOI(all_data, bonded_atoms_list, COM, Hs):	×
476	"""
477	Calculate the principal axes of the MOI matrix, then calc forces.
478	If Hs is True, then consider H coords in MOI calcs.
479	If Hs is False, then don't include H coords, instead add their mass
480	onto UA heavy atom (HA).
481	"""
482
483	coord_list = []	×
484	mass_list = []	×
485
486	for atom in bonded_atoms_list:	×
NEW 487	if Hs is True:	×
488	# print (atom.atom_name)
489	coord_list.append(atom.coords)	×
490	mass_list.append(atom.mass)	×
NEW 491	elif Hs is False:	×
492	if atom.mass > 1.1:	×
NEW 493	UA_mass = 0 # add mass of HA and bonded Hs	×
494	# print (atom.atom_name)
495	UA_mass += atom.mass	×
NEW 496	coord_list.append(atom.coords) # only include coords of UA	×
497	for h in atom.bonded_to_atom_num:	×
498	H = all_data[h]	×
499	if H.mass < 1.1:	×
500	UA_mass += H.mass	×
501	else:
502	continue	×
NEW 503	mass_list.append(UA_mass)	×
504	# include mass of Hs on heavy atom mass
505	else:
506	continue	×
507
508	moI = MOI(COM, coord_list, mass_list)	×
509	# if Hs == False:
510	# print ('masses', mass_list)
511	# print ('moi', moI)
512
NEW 513	eigenvalues, eigenvectors = LA.eig(moI)	×
514
515	# different values generated to Jon's code
516	# print (eigenvalues)
517	# print (eigenvectors)
NEW 518	transposed = np.transpose(eigenvectors) # turn columns to rows	×
519
520	# bonded_atoms_list[0].WMprincipalAxis = transposed, eigenvalues, COM
521	# print (moI, eigenvalues, eigenvectors)
522
523	min_eigenvalue = abs(eigenvalues[0])	×
524	if eigenvalues[1] < min_eigenvalue:	×
525	min_eigenvalue = eigenvalues[1]	×
526	if eigenvalues[2] < min_eigenvalue:	×
527	min_eigenvalue = eigenvalues[2]	×
528
529	max_eigenvalue = abs(eigenvalues[0])	×
530	if eigenvalues[1] > max_eigenvalue:	×
531	max_eigenvalue = eigenvalues[1]	×
532	if eigenvalues[2] > max_eigenvalue:	×
533	max_eigenvalue = eigenvalues[2]	×
534
535	# print (min_eigenvalue * const, max_eigenvalue * const)
536	# same as Jons
537
538	#
539	# PA = principal axes
NEW 540	axis1PA = [0, 0, 0] # [x,y,z]	×
NEW 541	axis2PA = [0, 0, 0] # [x,y,z]	×
NEW 542	axis3PA = [0, 0, 0] # [x,y,z]	×
543
544	# MI = moment of inertia
NEW 545	axis1MI = 0 # x	×
NEW 546	axis2MI = 0 # y	×
NEW 547	axis3MI = 0 # z	×
548	#
549
NEW 550	for i in range(0, 3):	×
551	if eigenvalues[i] == max_eigenvalue:	×
552	axis1PA = transposed[i]	×
553	axis1MI = eigenvalues[i]	×
554	elif eigenvalues[i] == min_eigenvalue:	×
555	axis3PA = transposed[i]	×
556	axis3MI = eigenvalues[i]	×
557	else:
558	axis2PA = transposed[i]	×
559	axis2MI = eigenvalues[i]	×
560
561	principal_axes = [axis1PA, axis2PA, axis3PA]	×
562	MI_axis = [axis1MI, axis2MI, axis3MI]	×
563
564	# print ('mI', MI_axis, bonded_atoms_list[0].atom_num)
565
566	return principal_axes, MI_axis	×
567
568
NEW 569	def rotateFT(	×
570	all_data,
571	bonded_atoms_list,
572	WM_principal_axes,
573	UA_principal_axes,
574	MI_axis,
575	molecule_COM,
576	UA_COM,
577	Hs,
578	):
579	""" """
580
581	# lists used for torque calcs
582	coord_list = []	×
583	mass_list = []	×
584	force_list = []	×
585	forces_summed = np.zeros(3)	×
586
587	for atom in bonded_atoms_list:	×
NEW 588	if Hs is True:	×
589	coord_list.append(atom.coords)	×
590	mass_list.append(atom.mass)	×
591	force_list.append(atom.forces)	×
592	forces_summed += atom.forces	×
NEW 593	elif Hs is False:	×
594	if atom.mass > 1.1:	×
595	# print (atom.atom_name)
NEW 596	UA_mass = 0 # add mass of HA and bonded Hs	×
597	UA_forces = np.zeros(3)	×
598	UA_mass += atom.mass	×
599	UA_forces += np.array(atom.forces)	×
NEW 600	coord_list.append(atom.coords) # only include coords of UA	×
601	forces_summed += atom.forces	×
602	for h in atom.bonded_to_atom_num:	×
603	H = all_data[h]	×
604	if H.mass < 1.1:	×
605	UA_mass += H.mass	×
606	UA_forces += np.array(H.forces)	×
607	forces_summed += H.forces	×
608	else:
609	continue	×
NEW 610	mass_list.append(UA_mass) # include mass of Hs on heavy atom mass	×
611	force_list.append(UA_forces)	×
612	else:
613	continue	×
614
NEW 615	torque = [0, 0, 0] # default is no torque	×
616	# calc torques here with COM and PA of molecule or UA
NEW 617	if UA_principal_axes is not None and MI_axis is not None:	×
NEW 618	if Hs is True:	×
NEW 619	torque = calcTorque(	×
620	UA_COM, coord_list, UA_principal_axes, force_list, MI_axis
621	)
NEW 622	if Hs is False:	×
NEW 623	torque = calcTorque(	×
624	molecule_COM, coord_list, WM_principal_axes, force_list, MI_axis
625	)
626
627	# Rotate forces here
628	F1 = np.dot(forces_summed, WM_principal_axes[0])	×
629	F2 = np.dot(forces_summed, WM_principal_axes[1])	×
630	F3 = np.dot(forces_summed, WM_principal_axes[2])	×
631	mass_sqrt = sum(mass_list) ** 0.5	×
NEW 632	force = np.array(	×
633	[
634	float(F1) / float(mass_sqrt),
635	float(F2) / float(mass_sqrt),
636	float(F3) / float(mass_sqrt),
637	]
638	)
639
640	return force, torque	×
641
642
643	def calcTorque(cm, coord_list, principal_axes, force_list, MI_axis):	×
644	""" """
645	x_cm, y_cm, z_cm = cm[0], cm[1], cm[2]	×
646	# MI_axis = np.sqrt(MI_axis) #sqrt moi to weight torques
647
648	MI_axis_sqrt = []	×
649	# print (MI_axis)
650	for coord in MI_axis:	×
651	coord_round = round(coord, 10)	×
652	if coord_round != 0:	×
NEW 653	c_sqrt = coord**0.5	×
654	MI_axis_sqrt.append(c_sqrt)	×
655	elif coord_round == 0:	×
656	MI_axis_sqrt.append(coord)	×
657	else:
658	continue	×
659
660	MI_axis_sqrt = np.array(MI_axis_sqrt)	×
661	# print (MI_axis, MI_axis_sqrt)
662
663	W_torque = np.zeros(3)	×
664	for coord, force in zip(coord_list, force_list):	×
665	atom_torque = np.array([0, 0, 0])	×
666	# count_atom += 1
667	new_coords = []	×
668	new_forces = []	×
669	for axis in principal_axes:	×
NEW 670	new_c = (	×
671	axis[0] * (coord[0] - x_cm)
672	+ axis[1] * (coord[1] - y_cm)
673	+ axis[2] * (coord[2] - z_cm)
674	)
NEW 675	new_f = axis[0] * force[0] + axis[1] * force[1] + axis[2] * force[2]	×
676	new_coords.append(new_c)	×
677	new_forces.append(new_f)	×
678	# print (new_coords)
679	# print (force)
680	# CROSS PRODUCT
NEW 681	torquex = float(	×
682	new_coords[1] * new_forces[2] - new_coords[2] * new_forces[1]
683	) # / float(1e10)
NEW 684	torquey = float(	×
685	new_coords[2] * new_forces[0] - new_coords[0] * new_forces[2]
686	) # / float(1e10)
NEW 687	torquez = float(	×
688	new_coords[0] * new_forces[1] - new_coords[1] * new_forces[0]
689	) # / float(1e10)
690	# print (torquex, torquey, torquez)
691
692	atom_torque = (torquex, torquey, torquez)	×
693	# atom_torque = np.divide(atom_torque, MI_axis_sqrt)
694
695	atom_torque2 = []	×
696	for t, MI in zip(atom_torque, MI_axis_sqrt):	×
697	if MI != 0:	×
698	t_divide = float(t) / float(MI)	×
699	atom_torque2.append(t_divide)	×
700	elif MI == 0:	×
701	atom_torque2.append(0)	×
702	else:
703	continue	×
704
705	atom_torque2 = np.array(atom_torque2)	×
706
707	W_torque += atom_torque2	×
708
NEW 709	return W_torque	×

CCPBioSim / CodeEntropy / 13726590978

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

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