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.

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]