force_halving issue #22

CodeEntropy/GeometricFunctions.py

@@ -33,11 +33,6 @@ def get_beads(data_container, level):
             atom_group = "index " + str(atom.index) + " or (name H* and bonded index " + str(atom.index) +")"
             list_of_beads.append(data_container.select_atoms(atom_group))
 
-    #TODO temporary print
-    print("list of beads")
-    print(level)
-    print(list_of_beads)
-
     return list_of_beads
 #END
 
@@ -206,14 +201,16 @@ 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):
     """
     Function to calculate the mass weighted forces for a given bead.
 
     Input
     -----
+    data_container : information about the system
     bead : the part of the system to be considered
     trans_axes : the axes relative to which the forces are located
+    highest_level : true/false, true if whole molecule level
 
     Output
     ------
@@ -233,6 +230,11 @@ def get_weighted_forces(data_container, bead, trans_axes):
 
     weighted_force = forces_trans / nmp.sqrt(mass)
 
+    # if whole molecule level, divide by 2 to divide the independend forces
+    # between the atoms involved
+    if highest_level:
+        weighted_force = 0.5 * weighted_force
+
     return weighted_force
 #END
 
@@ -261,6 +263,9 @@ def get_weighted_torques(data_container, bead, rot_axes):
         coords_rot = nmp.matmul(rot_axes, coords_rot)
         # update local forces in rotational frame
         forces_rot = nmp.matmul(rot_axes, data_container.atoms[atom.index].force)
+
+        # Divide forces by 2 to avoid double counting of independent motions
+        forces_rot = 0.5 * forces_rot
 
         # define torques (cross product of coordinates and forces) in rotational axes
         torques_local = nmp.cross(coords_rot, forces_rot)

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.
@@ -58,6 +58,8 @@ def get_matrices(data_container, level, verbose, start, end, step, number_frames
         start : int, starting frame, default 0 (first frame)
         end : int, ending frame, default -1 (last frame)
         step : int, step for going through trajectories, default 1
+        number_frames : int, number of frames in trajectory
+        highest_level : true/false, is this the whole molecule level
 
     Returns
     -------
@@ -84,7 +86,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

@@ -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)