GDAL: sum of all bands in a raster

I’m familiar with using gdal_calc.py to sum individual raster bands. However, I am curious whether there is an easy way to programmatically sum all bands of a raster into a single value. I have a list of multi-band .vrts, each with a variable number of bands (FYI, each .vrt represents harvested area of a number of different crops in different groups such as cereals, fruits, etc. So each .vrt may have differing numbers of bands depending on how many crops are in that group.) For each individual .vrt, I want to sum up all the bands into a single value and write it to another .vrt. What is the recommended syntax for this?

This question is sort of a duplicate of this. But the linked question is >6 years old so I was hoping to give some new attention to the topic, in case there is a new and better answer.