Using a scalar relativistic all-electron density functional method, we studied uranium(VI) complexation with benzoic acid and its derivatives in aqueous solution as models of uranyl humates. We explored monodentate, bidentate, and chelate coordination of various isomers of methyl and hydroxyl substituted benzoic acid ligands. Monodentate complexes were determined to be energetically preferred as long as entropy effects were neglected. However, bidentate structures were favored at the Gibbs free energy level. Coordination of aromatic carboxylic acids tends to be weaker than that of aliphatic ones, while structural characteristics were determined to be quite similar. Optimized geometries yield uranyl bonds and U-C distances in agreement with EXAFS results for monocarboxylate of benzoate and p-hydroxy benzoate. Average uranyl-oxygen distances to equatorial ligands, U-O(eq), are shorter than in experiment, which is tentatively rationalized by variations in the coordination numbers. As for aliphatic monocarboxylate complexes studied earlier, U-O(eq) values of benzoic acid derivatives do not discriminate mono- and bidentate coordinated species. Structures and energies determined support the interpretation of uranyl humate complexes as bidentate carboxylate species with fivefold coordination of uranyl.