Accurate methods to predict the free energies of protein-ligand interactions have great potential to assist rational drug design. In this work, we used molecular dynamics simulations with alchemical perturbation to predict the binding of carbohydrate-based ligands to influenza virus neuraminidase (N2). This specific drug target is a challenging test system for alchemical free energy methods because it has flexible binding site motifs. We use a molecular dynamics protocol that works for longer time scales than are often reported in previous molecular dynamics studies of N2. We demonstrated that N2-ligand complex stability and that accurate N2 150-loop dynamics, on a 350 ns time scale, are both force field-dependent (AMBER99SB-ILDN, GAFF and TIP4P water are required). Further, we showed that crystallographic waters must be retained to maintain stability of N2-ligand complexes over 350 ns. Using our modelling protocol, we were able to determine relative binding free energy values between neuraminidase ligands which correlated strongly with experimental differences in pIC50 values (R = -0.96, ρ = 0.86, N = 13, sig < 0.0001). It is anticipated that the molecular dynamics protocol and the relative binding free energy methods reported here, will both be useful in expediting the discovery of novel therapeutics for N2 and other homologous glycohydrolases.