The thermodynamic coefficients of a free standing infinite graphene monolayer are calculated using the quasi-classical unsymmetrized self-consistent field method (USF). The basic nonlinear integral equations of this theory are solved numerically in the strong anharmonic approximation. The isothermal and adiabatic elastic bulk moduli, the isochoric and isobaric heat capacities, the thermal expansion, thermal pressure coefficient, and the macroscopic Gruneisen parameter are calculated in terms of the derivatives of a specifically chosen interatomic potential function for different values of stress and for temperatures ranging from below room temperature up to the point of loss of thermodynamic stability. The nearest-neighbor distances vary from 1.4 Å to 1.8 Å for zero stress. Under stress, these distances decrease. At room temperature the molar heat capacities are ∼5.0Jmol −1 K −1 . The elasticity moduli vary from 15.0eV 2 -Å up to zero at the temperature of loss of stability and are increased by stress. The thermal expansion coefficient has a strong dependence on the temperature and is negative for temperatures lower than ∼340K. For high temperatures it monotonically increases and decreases with stress. The macroscopic Gruneisen parameter has a strong nonlinear dependence with temperature and is estimated in about 3.0 3.7;at ∼340K its value decreases to ∼1.0K and for even lower temperature it shows a peak and deep structure similar to what has been earlier reported for fullerene C 60 .The thermodynamic properties of crystals are related to atomic vibrations around their equilibrium positions. The relative displacements are small in comparison with the interatomic distance. Even in quantum crystals, where the zero point vibrations are large, the mean square deviations from the sites are about 30% of interatomic distances while for most of the solids they do not exceed about 10% up to the melting temperatures. This allows for the expansion of the potential energy in a power series of the nuclei displacements in the Born-Oppenheimer approximation [22]. For crystals with strong anharmonicity, perturbative techniques based on harmonic or quasiharmonic approximations do not work. The pseudoharmonic approximation [19][20][21] reduces the study of strong anharmonicity to the analysis of weak interacting phonons and because of that it is more suited to the investigation of the phonon spectrum and related phenomena.The USF has been applied over the years to several types of crystals but mostly with Bravais lattices, and for non Bravais ionic crystals [23,24]. It has been also applied to noncrystalline solids [25]. Recently [26], we have reported on a generalization of the USF to the case of layered structured crystals. In that work though, the thermodynamics of a free standing graphene monolayer was investigated in the weak anharmonic approximation. The anharmonicity has a close relation to the thermal properties of graphene and on-going debates about the contribution of different phonon modes to heat conduction. The...