I think you want
h = F.relu(self.fc3(z)) # KW use tanh instead of relu
mu = self.fc4(h)
sigma = self.fc5(h)
return mu, sigma

and then pass mu, sigma and x into your criterion to give the negative log likelihood of x as a vector sampled from N(mu_i, sigma_i^2) (i.e. you need to take the log PDF).
So the Gaussian at the reconstruction step has nothing to do (well, except being conditional on the latents) with the Gaussian from the latents (which is the bit where you do the reparametrization and things).
I think what might cause the confusion that the “reconstruction error” is easily taken too literally when looking at the example implementation. In the paper they only take and maximize the (log) likelihood of the input data under the output distribution rather than doing actual reconstruction (by sampling or somesuch). In evaluation mode, you would be expected to sample from the Gaussian.
The last paragraph of section 2.3 in the paper talks about how in the probabilistic VAE framework this term is a […] reconstruction error in auto-encoder parlance (notice the indefinite article) or maybe rather the analogue of the reconstruction error.
Best regards
Thomas