The dynamics of a protein and the water surrounding it are coupled via nonbonded energy interactions. This coupling can exhibit a complex, nonlinear, and nonstationary nature. The THz frequency spectrum for this interaction energy characterizes both the vibration spectrum of the water hydrogen bond network, and the frequency range of large amplitude modes of proteins. We use a Recurrence Plot based Wiener-Khinchin method RPWK to calculate this spectrum, and the results are compared to those determined using the classical auto-covariance-based Wiener-Khinchin method WK. The frequency spectra for the total nonbonded interaction energy extracted from molecular dynamics simulations between the β-Lactamase Inhibitory Protein BLIP, and water molecules within a 10 Å distance from the protein surface, are calculated at 150, 200, 250, and 310 K, respectively. Similar calculations are also performed for the nonbonded interaction energy between the residues 49ASP, 53TYR, and 142PHE in BLIP, with water molecules within 10 Å from each residue respectively at 150, 200, 250, and 310 K. A comparison of the results shows that RPWK performs better than WK, and is able to detect some frequency data points that WK fails to detect. This points to the importance of using methods capable of taking the complex nature of the protein-solvent energy landscape into consideration, and not to rely on standard linear methods. In general, RPWK can be a valuable addition to the analysis tools for protein molecular dynamics simulations. Proteins 2016; 84:1549-1557. © 2016 Wiley Periodicals, Inc.