Dimerization of formic acid has been simulated using ab initio molecular dynamics at conditions mimicking rare gas matrix isolation experiments. Aggregation product distributions and the corresponding reaction pathways have been studied as a function of temperature. At higher temperatures, the cyclic, C-2h symmetric, global minimum structure A with two O-H center dot center dot center dot O=C hydrogen bonds predominates over the second most stable, acyclic, local minimum isomer B with one O-H center dot center dot center dot O=C and one C-H center dot center dot center dot O=C hydrogen bond, whereas the latter is the main species at low temperature. Significant concentrations of two additional, less stable, local minimum species, C exhibiting an O-H center dot center dot center dot O=C and an O-H center dot center dot center dot O(H)-C hydrogen bond, and D with an O-H center dot center dot center dot O(H)-C and a C-H center dot center dot center dot O=C hydrogen bond, which should be detectable in experiment, are predicted at low temperature. Theoretical vibrational spectra are provided to guide the experimental search for these species. Furthermore, free-energy barriers for thermal interconversion of different dimer species have been calculated using targeted molecular dynamics in conjunction with thermodynamic integration. The small barrier for the C -> A reaction of 7.0 kJ/mol indicates that the C species can only be stabilized at ultracold conditions. The data presented here thus hold important clues for the execution and interpretation of low-temperature vibrational spectroscopy at matrix isolation or ultracold helium droplet conditions.