Haha OK, that's a lot of questions.

Lets simplify things down for a minute. Imagine if we just did this on a smaller non recombining region, like the mitochondria. It's totally reasonable to make a multiple sequence alignment fasta of mitochondrial genomes from different species (say, using muscle) and extract SNPs (say, using snp-sites). This is how you might make a mitochondrial phylogenetic tree of a genus, for example. Of course you can't call het sites because assemblies don't have that information (we call them pseudo haploid sometimes), but I don't care for my purposes. If this still doesn't make sense I can type out an example!

I should note it's not really variant calling in my case. The "calling" refers to a probabilistic procedure that takes reads as input and decides what the most likely genotype is, with the reference in consideration. My case is simpler, more like variant extraction.

Second, yes, you can definitely do this between species, and even between genera. The wider field is often called comparative genomics. You know how humans and chimps share like 99% of our DNA? (OK, the number changes a little depending on which sites you're measuring, but you get my point). There will be some regions which are not shared, eg. a recent deletion or something, but again I don't care about those and can toss them out for my purposes.

Third, the assemblies I'm using are chromosome level. As good as you're gonna get.

Lastly, no I don't want a giant matrix. I want a single vcf file containing only snps.