docs: def determine_hbonds_for_drawing(self, analysis_cutoff):
          """
          Since plotting all hydrogen bonds could lead to a messy plot, a cutoff has to be imple-
          mented. In this function the frequency of each hydrogen bond is summated and the total
          compared against analysis cutoff - a fraction multiplied by trajectory count. Those
          hydrogen bonds that are present for longer than analysis cutoff will be plotted in the
          final plot.
  
          Takes:
              * analysis_cutoff * - (user-defined) fraction of time a hydrogen bond has to be
              present for to be plotted (default - 0.3). It is multiplied by number of trajectories
          Output:
              * frequency * - dictionary of hydrogen bond donor-acceptor indices and frequencies
              These hydrogen bonds will be plotted in the final image.
  
          """
          self.frequency = defaultdict(int)
          for traj in self.hbonds_by_type:
              for bond in self.hbonds_by_type[traj]:
              # frequency[(residue_atom_idx,ligand_atom_name,residue_atom_name)]=frequency
              # residue atom name will be used to determine if hydrogen bond is interacting with a sidechain or bakcbone
                  if bond["donor_resnm"]!="LIG":
                      self.frequency[(bond["donor_idx"],bond["acceptor_atom"],bond["donor_atom"],bond["acceptor_idx"])] += bond["frequency"]
                  #check whether ligand is donor or acceptor
                  else:
                      self.frequency[(bond["acceptor_idx"],bond["donor_atom"],bond["acceptor_atom"],bond["donor_idx"])] += bond["frequency"]
  
          #Add the frequency counts
          self.frequency = {i:self.frequency[i] for i in self.frequency if self.frequency[i]>(int(len(self.trajectory))*analysis_cutoff)}
  
          #change the ligand atomname to a heavy atom - required for plot since only heavy atoms shown in final image
          self.hbonds_for_drawing = {}
          for bond in self.frequency:
              atomname = bond[1]
              if atomname.startswith("O",0) or atomname.startswith("N",0):
                  lig_atom=atomname
              else:
                  atomindex = [index for index,atom in enumerate(self.topology_data.universe.ligand.atoms) if atom.name==atomname][0]
                  rdkit_atom = self.topology_data.mol.GetAtomWithIdx(atomindex)
                  for neigh in rdkit_atom.GetNeighbors():
                      neigh_atom_id = neigh.GetIdx()
                  lig_atom = [atom.name for index,atom in enumerate(self.topology_data.universe.ligand.atoms) if index==neigh_atom_id][0]
              self.hbonds_for_drawing[(bond[0],lig_atom,bond[2],bond[3])]=self.frequency[bond]