My reasoning is as follows (using Gaussian units):
Start from the second law:
$$dU=TdS+dW,$$
where $dW$ is the work done by the magnetic field. To derive $dW$, we consider a solenoid with current $I$ and voltage $V$:
$$dW=I\cdot V\cdot dt,$$
with
$$\begin{align} V&=\frac{N}{c}\cdot\frac{d(B\cdot A)}{dt} \qquad &&\text{(Faraday's law)} \\ I&=\frac{c}{4\pi}\cdot\frac{H}{n}. \qquad &&\text{(Ampere's law)} \end{align}$$
Now plug in $I$ and $V$ in $dW$. We get
$$dW=\frac{1}{4\pi}\cdot HdB.$$
Since $B=H+4 \pi M$, if we exclude the energy of the magnetic field itself, we find
$$dW=HdM.$$
Thus the second law is
$$dU=TdS+HdM.$$
For the Helmholtz free energy we have
$$F=U-TS \Rightarrow dF=-SdT+HdM.$$
Similarly, for the Gibbs free energy:
$$G=F-HM \Rightarrow dG=-SdT-MdH.$$
In experiment, we control $T$ and $H$, therefore at any given $T$ and $H$, the system should minimize its Gibbs free energy $G$. Ergo the energy we talk about for superconductor in magnetic field should be Gibbs free energy.
The questions is, why the famous superconducting textbooks (like Tinkham, Schmidt and de Gennes) use Helmholtz free energy $F$ instead of $G$? It doesn't make sense to me to minimize $F$ instead of $G$ (e.g. in deriving the GL equations using variation method).