The third-order elastic constants (TOECs) are fundamental to describe crystal’s nonlinear response to stress, and can be applied to explore anharmonic properties of crystals such as Gruneisen parameters, thermal expansion coefficient, and the effect of pressure on second-order elastic constants (SOECs). Here, we report an open-source python package, Elastic3rd, which is able to calculate the SOECs and TOECs using the strain–energy method for crystals with any symmetry from first-principles calculations. An algorithm to generate necessary strain modes and the corresponding coefficients for a given symmetry is proposed. These strain modes are then applied to the fully relaxed structure to generate the deformed structures. The total energies of the strained structures are calculated by a chosen first-principles code, and the SOECs and TOECs are determined by fitting the resulted strain–energy data. The present code has been validated by several case studies of C, Si and Mg, and the case of MnP4 shows the ability for low-symmetry crystals. Program summary Program title: Elastic3rd CPC Library link to program files: https://doi.org/10.17632/n54vr2kwx8.1 Developer’s repository link: https://github.com/hitliaomq/ELASTIC3rd Licensing provisions: GNU General Public License 3 Programming language: Python 2.7.X and Python 3.X. External routines/libraries: Numpy [1], Scipy [2], and Matplotlib [3]. Nature of problem: To automatically calculate the third-order elastic constants for crystals with any symmetry from first-principles calculations. Solution method: Firstly, the deformed patterns and the corresponding coefficients for the second- and third-order elastic constants are automatically generated according to the symmetry. Secondly, the deformed patterns are applied to the initial lattice and the corresponding input files for first-principles calculations are generated. Thirdly, the energies of the deformed structures are calculated by first-principles software such as CASTEP [4] and VASP [5]. Then, the parameters of strain–energy functions are obtained by fitting the strain and energy relationships. Finally, the third-order elastic constants are determined by solving the linear-independent equations, and the second-order elastic constants are calculated by solving the overdetermined equation using the least squares method. Addition comments including restrictions and unusual features: This code is well modularized. Thus, the energy can be calculated using energy-generated first-principles or other packages only by writing few functions. [1] S. Van Der Walt, S.C. Colbert, G. Varoquaux, Comput. Sci. Eng. 13(2011):22–30 [2] P. Virtanen, R. Gommers, T.E. Oliphant, et al., ArXiv:1907.10121 (2019). [3] J.D. Hunter, Comput. Sci. Eng. 9(2007):99–104. [4] S.J. Clark, M. D. Segall, C.J. Pickard, et al. Z. Kristallogr. 220(2005):567-570. [5] G. Kresse, J. Furthmuller, Phys. Rev. B 54 (1996) 11169.