I think I have arrived to something useful, and it has to do with my attempt to ‘fix’ the function free_energy(). The problem is that, when ‘wx_b’ overflows, the computetion of hidden_term involves a non-continuous (and hence non-differentiable) transformation: the substitution of the cell in hidden_term that evaluates to inf by the corresponding value in wx_b. In an attempt to solve it, I have turned to the log-sum-exp trick, having
def free_energy(self,v):
    vbias_term = v.mv(self.v_bias)
    wx_b = F.linear(v,self.W,self.h_bias)
    zr = Variable(torch.zeros(wx_b.size()))
    mask = torch.max(zr, wx_b)
    hidden_term = (((wx_b - mask).exp() + (-mask).exp()).log() + (mask)).sum(1)
    return (-hidden_term - vbias_term).mean()

This code, for now, is giving a finite value for a gradient that was NaN in an explicit example that I have been able to isolate. I will report if any other problems appear.
EDIT: This indeed does the work. I have trained the model for 50 epochs, while before the NaNs would appear no later than in the 20th epoch. Thank you all for the tips!