Way to calculate the total tetrad twist of a rubik's cube

It is rather difficult to extract this information directly from the current locations of the corner pieces. By far the easiest way is to actually try to solve those corner pieces using only half turns (while ignoring the edge pieces), and see how far you get.

For now I'll assume that the corner pieces are already located in their correct tetrad orbits {UFR, UBL, DFL, DBR} and {UFL, UBR, DFR, DBL}. You can solve the pieces of one tetrad really easily, no more than one half turn for each piece, at most 3 moves in total. For example, solve DBR using at most one of {D2, B2, R2}, then DFL using at most one of {F2, L2}, and finally UBL using {U2} if necessary, which also leaves UFR solved.

You then solve one piece of the second tetrad, for example DBL, using one of the move sequences {F2 L2 F2 U2, U2 F2 U2 L2, L2 U2 L2 F2}. These move sequences perform a double swap on the four pieces of the second tetrad, and are the only permutations possible that keep the first tetrad fixed.

This leaves three unsolved pieces, {UFL, UBR, DFR}. These can be in any one of 3!=6 permutations. These 6 possibilities represent the tetrad twist in combination with the permutation parity, so if you map this permutation to a number from 0 to 5, you have encoded both the permutation parity and the tetrad twist into a single number.

For the Thistlethwaite algorithm you'll probably want to encode an arbitrary position of the third stage of the algorithm. This has to be done in a consistent manner, by which I mean that if two different positions are brought into the fourth stage by the same move sequence (i.e. after applying the move sequence to those positions they both become solvable using only half turns) then those two positions must have the same encoding for stage 3.

Presumably your program lists the corner locations of the cube in a particular fixed order. For example, you may have an array of length 8 representing the locations in the order
UFR, UFL, UBL, UBR, DFR, DFL, DBL, DBR.
I have bolded those locations that represent one of the tetrads, those at index 0, 2, 5, 7 in the array. You may have chosen a different ordering convention in your program, but the method is the same.

Suppose you have now an arbitrary stage 3 cube position, some random permutation of those 8 corners, for example:
UBR, UBL, DBR, DFR, DFL, UFR, UFL, DBL.
A simple consistent way to separate the pieces into the two tetrads is to literally separate out the two types of pieces without changing their relative order:
UBL, DBR, DFL, UFR
UBR, DFR, UFL, DBL.
And then put them into the storage array, in order, each into their correct set of tetrad locations. So the first set goes into indices 0,2,5,7, the other at indices 1,3,4,6.
UBL, UBR, DBR, DFR, UFL, DFL, DBL, UFR.
Now you can apply the solving technique I explained at the start, to end up with a consistent encoding of the positions tetrad twist and parity.
The above assumes you use a single standardised way to represent the cube, and you apply moves to that. You may instead want to keep the two separate lists of the tetrad pieces as a simplified representation of this position, and permute those directly as you solve it to extract the twist+parity encoding.

You may take a look at some of the programs at this old cube programming contest, though I am not sure they will be terribly helpful since they are written for conciseness, not for intelligibility.