adding force halving

Author: skfegan

CodeEntropy/GeometricFunctions.py

@@ -200,7 +200,7 @@ def get_sphCoord_axes(arg_r):
     return spherical_basis
 # END
 
-def get_weighted_forces(data_container, bead, trans_axes):
+def get_weighted_forces(data_container, bead, trans_axes, highest_level, force_partitioning=0.5):
     """
     Function to calculate the mass weighted forces for a given bead.
 
@@ -222,6 +222,12 @@ def get_weighted_forces(data_container, bead, trans_axes):
         forces_local = nmp.matmul(trans_axes, data_container.atoms[atom.index].force)
         forces_trans += forces_local
 
+    if highest_level:
+        # multiply by the force_partitioning parameter to avoid double counting
+        # of the forces on weakly correlated atoms
+        # the default value of force_partitioning is 0.5 (dividing by two)
+        forces_trans = force_partitioning * forces_trans
+
     # divide the sum of forces by the mass of the bead to get the weighted forces
     mass = bead.total_mass()
 
@@ -230,7 +236,7 @@ def get_weighted_forces(data_container, bead, trans_axes):
     return weighted_force
 #END
 
-def get_weighted_torques(data_container, bead, rot_axes):
+def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5):
     """
     Function to calculate the moment of inertia weighted torques for a given bead.
 
@@ -256,6 +262,11 @@ def get_weighted_torques(data_container, bead, rot_axes):
         # update local forces in rotational frame
         forces_rot = nmp.matmul(rot_axes, data_container.atoms[atom.index].force)
 
+        # multiply by the force_partitioning parameter to avoid double counting
+        # of the forces on weakly correlated atoms
+        # the default value of force_partitioning is 0.5 (dividing by two)
+        forces_rot = force_partitioning * forces_rot
+
         # define torques (cross product of coordinates and forces) in rotational axes
         torques_local = nmp.cross(coords_rot, forces_rot)
         torques += torques_local

CodeEntropy/LevelFunctions.py

@@ -45,7 +45,7 @@ def select_levels(data_container, verbose):
     return number_molecules, levels
 #END get_levels
 
-def get_matrices(data_container, level, verbose, start, end, step, number_frames):
+def get_matrices(data_container, level, verbose, start, end, step, number_frames, highest_level):
     """
     Function to create the force matrix needed for the transvibrational entropy calculation
     and the torque matrix for the rovibrational entropy calculation.
@@ -84,7 +84,7 @@ def get_matrices(data_container, level, verbose, start, end, step, number_frames
             trans_axes, rot_axes = GF.get_axes(data_container, level, bead_index)
 
             ## Sort out coordinates, forces, and torques for each atom in the bead
-            weighted_forces[bead_index][timestep.frame] = GF.get_weighted_forces(data_container, list_of_beads[bead_index], trans_axes)
+            weighted_forces[bead_index][timestep.frame] = GF.get_weighted_forces(data_container, list_of_beads[bead_index], trans_axes, highest_level)
             weighted_torques[bead_index][timestep.frame] = GF.get_weighted_torques(data_container, list_of_beads[bead_index], rot_axes)
 
     ## Make covariance matrices - looping over pairs of beads

CodeEntropy/main_mcc.py

@@ -1,10 +1,10 @@
+import math
 import MDAnalysis as mda
 import numpy as nmp
 import pandas as pd
 from CodeEntropy import LevelFunctions as LF
 from CodeEntropy import EntropyFunctions as EF
 from CodeEntropy import MDAUniverseHelper as MDAHelper
-from CodeEntropy import poseidon
 
 def main(arg_dict):
     """
@@ -23,13 +23,13 @@ def main(arg_dict):
 
     # Define trajectory slicing from inputs
     start = arg_dict['start']
-    if start == None:
+    if start is None:
         start = 0
     end = arg_dict['end']
-    if end == None:
+    if end is None:
         end = -1
     step = arg_dict['step']
-    if step == None:
+    if step is None:
         step = 1
     # Count number of frames, easy if not slicing
     if start == 0 and end == -1 and step == 1:
@@ -100,7 +100,7 @@ def main(arg_dict):
 
                     ## Vibrational entropy at every level
                     # Get the force and torque matrices for the beads at the relevant level
-                    force_matrix, torque_matrix = LF.get_matrices(residue_container, level, arg_dict['verbose'], start, end, step, number_frames)
+                    force_matrix, torque_matrix = LF.get_matrices(residue_container, level, arg_dict['verbose'], start, end, step, number_frames, highest_level)
 
                     # Calculate the entropy from the diagonalisation of the matrices
                     S_trans_residue = EF.vibrational_entropy(force_matrix, "force", arg_dict['temper'],highest_level)
@@ -177,7 +177,7 @@ def main(arg_dict):
             if level in ('polymer', 'residue'):
                 ## Vibrational entropy at every level
                 # Get the force and torque matrices for the beads at the relevant level
-                force_matrix, torque_matrix = LF.get_matrices(molecule_container, level, arg_dict['verbose'], start, end, step, number_frames)
+                force_matrix, torque_matrix = LF.get_matrices(molecule_container, level, arg_dict['verbose'], start, end, step, number_frames, highest_level)
 
                 # Calculate the entropy from the diagonalisation of the matrices
                 S_trans = EF.vibrational_entropy(force_matrix, "force", arg_dict['temper'],highest_level)