How to include Berry connection in Hamiltonian? When we calculate Berry connection,
$A(R)=i<\psi(x,y)|\frac{d}{dR}|\psi(x,y)>\hat{R}$ corresponding to the Berry phase of any system, the gauge potential is related to the $R$ of the parameter space. It is not dependent on $x$ like $\psi(x)$ because in the inner product the spatial part is integrated over.
We also know that in presence of a vector potential, the Hamiltonian can be written as
$\mathcal{H}=\dfrac{(p-A)^2}{2m}+V$. Here $p=p_x \hat{i}+p_y \hat{j}$, (in 2D) and $A$ is also supposed to be $A_x \hat{i}+A_y\hat{j}$.
But in case of Berry connection, $A=A(R) \hat{R}$ corresponding to parameter space. So, how to include this Berry gauge potential in Hamiltonian? I cannot just subtract a vector in parameter coordinates from a vector in spatial coordinate, right?
 A: The Berry connection lives in the parameter space, thus it appears not in the microscopic Hamiltonian given in the question but in  the effective Hamiltonian equation of motion in the parameter space.  The aim, in the following,  is to show the details in the variational approximation.
To be precise, the bra-ket notation that I'll be using is explained in the following two equations:
$$\mathbf{A} = \langle\psi(R) |\mathbf{\nabla_R} | \psi(R) \rangle \triangleq
 \int d^3x \ \psi(x,R)^{\dagger} \mathbf{\nabla_R}( \psi(x,R)) $$
Where $\psi(x,R)$ are the scalar wave functions corresponding to the state vectors $| \psi(R) \rangle$
In components the above equation takes the form:
$$A_i = \langle\psi(R) |\frac{\partial} {\partial R^i} | \psi(R) \rangle = \int d^3x \ \psi(x,R)^{\dagger}\frac{\partial} {\partial R^i}\psi(x,R)) $$
Given a microspcopic Hamiltonian (which can also have an explicit dependence on the parameter space$H(x, R)$, (For a fixed $R$,   $H(x, R)$ can be a Schroedinger operator in the real space), the effective Hamiltonian on the parameter space is defined by:
$$\mathcal{H}(R) = \langle\psi(R) |H | \psi(R) \rangle = \int d^3x \ \psi(x,R)^{\dagger}H\psi(x,R)) $$
$H$ can have the form given in the question with $A$ being an additional external field not related to the Berry connection.
The exact time dependent shroedinger equation:
$$ \frac{\partial  | \psi \rangle }{\partial t} = H | \psi\rangle $$
can be derived from by the variation of the Lagrangian:
$$ L = \langle\psi(R) |(\frac{\partial}{\partial t}- H )| \psi(R) \rangle$$
According to the variational approximation, we seek a solution by varying the state vectors not on the whole Hilbert space but only within the parameter space. The meaning of this approximation that we are not allowing the state vectors to vary in the directions of the real space, thus we are close to the lowest excitation state of the hamiltoniam for a fixed $R$.
Thus, we vary the Lagrangian only with respect to the parameter space and find a solution to the Lagrange equation of motion:
$$\frac{d}{dt} (\frac{\partial L}{\partial \dot{R}^i}) - \frac{\partial L}{\partial R^i} = 0$$
Using the above notation we have:
$$ L = A_i \dot{R}^i - \mathcal{H}(R) $$
The Lagrange equations of motion:
$$\frac{dA_i }{dt} - \frac{\partial A_j}{\partial R^i} \dot{R}^j - \frac{\partial \mathcal{H}}{\partial R^i} = 0$$
Using:
$$\frac{dA_i }{dt} = \frac{\partial A_i}{\partial R^j} \dot{R}^j$$
We obtain:
$$( \frac{\partial A_i}{\partial R^j}  - \frac{\partial A_j}{\partial R^i} )\dot{R}^j =  \frac{\partial \mathcal{H}}{\partial R^i} $$
Recognizing the expression of the Berry curvature:
$$ F_{ij}(R) = \frac{\partial A_i}{\partial R^j}  - \frac{\partial A_j}{\partial R^i}$$
We obtain:
$$ F_{ij}(R)\dot{R}^j =  \frac{\partial \mathcal{H}}{\partial R^i} $$
If the Berry curvature is invertible, i.e., there is a matrix $\Omega^{ij}(R)$ such that:
$$\Omega^{ij}(R) F_{jk}(R) = \delta^i_k,$$
we get the equation of motion on the parameter space:
$$ \dot{R}^i =  \Omega^{ij}\frac{\partial \mathcal{H}}{\partial R^j} $$
These are classical Hamilton equation of motion, with the symplectic structure equal to the Berry curvature. Thus the quantum evolution can be approximated by a classical evolution in the parameter space. 
Update: Effect on the system spectrum
In many cases, the parameter space is compact integrable system. In this case, its quantization splits each spectral line of the system's spectrum into a finite number of energy levels. In this case there will be a finite number of periodic solutions of the equations of motion 
$$R(t) = R(t+T),$$
living on fixed energy hypersurfaces $\mathcal{H} = E = const.$. We should fix our system parameters such that the Bohr-Sommerfeld quantization condition:
$$\int_{0}^{T} A_i(R(t) )\dot{R}^i(t) dt = 2\pi n$$
Where $n$ is a fixed integer. The above conditions are implicit equations which fix the energy levels $E_k$ and the corresponding periods $T_k$.
The corresponding state is given by:
$$| \psi_k\rangle = \frac{1}{T_k}\int_0^{T_k} dt \ e^{i\int_0^{t} A_i(R(\tau) )\dot{R}^i(\tau) d\tau} | \psi(R(t))\rangle $$
A: My apologies if I've misestimated where your confusion lies, but what I answer will address some ambiguities that bothers me in some pedagogical expositions. It will also address some common misconceptions.
If you are asking what the Hamiltonian is given a parameterized vector potential then it is supposed to be an operator, but the notation about which operator can be confusing.
When there is no vector potential the Hamiltonian is
$$\hat H=\frac{1}{2m}\hat p^2+V$$ which equals $$\hat H=\frac{1}{2m}\left(
\hat p_x^2+\hat p_y^2+\hat p_z^2\right)+
V.$$
And in this case that momentum is the kinetic and the canonical momentum because they are the same. It bothers me that the notion makes it seem like there is a momentum operator that operates twice but there you go.
When there is a vector potential then there are three self adjoint operators $\hat A_x,$ $\hat A_y,$ and $\hat A_z.$ Right away that could bother you, because it makes it seem like you can observe the vector potential, which is gauge dependant. But the point is that $\hat A_i$ is a real valued scalar function of position so on the position representation they each as Hermitian as the scalar potential $V.$ And the scalar potential was gauge dependent too, even when there was no vector potential. So observable doesn't mean you can see it in the lab, it means you can project onto different eigenspaces and separate those projections by coupling them differently to other things. In the case of energy, you don't measure energy per se, you measure energy differences. Though people do sometimes try to talk about a zero of energy in quantum mechanics even though the original classical potential could be adjusted by an arbitrary constant.
The second thing is that in the Hamiltonian $$ \frac{1}{2m}\left(
\left(\hat p_x-q\hat A_x\right)^2+
\left(\hat p_y-q\hat A_y\right)^2+
\left(\hat p_z-q\hat A_z\right)^2
\right)+V,$$
the $\hat p_i$ that appear are the canonical momentum and the $p-qA$ is then the kinetic momentum where $q$ is the charge of the particle (and sometimes people even write the charge in an easily ambiguous manner).  Finally the product of Hermitian operators isn't always Hermitian, but since the products are all terms like $\left(\hat p_x-q\hat A_x\right)^2$=$\hat p_x^2$-$q\hat p_x\hat A_x$-$q\hat A_x\hat p_x$+$q^2\hat A_x^2$ where it is the same Hermitian operator squared then you are fine (i.e. the square of a Hermitian operator is Hermitian even though a product of Hermitian operators might not be Hermitian). Note also that since the canonical momentum doesn't commute with the vector potential we need all four terms in any square. Finally the result is the kinetic momentum, but remember that the different components if kinetic momentum no longer commute with each other. The kinetic momentum has components like $\hat p_x-q\hat A_x$ and while each of the components are Hermitian they do not compute with each other. This means that like spin or any other angular momentum any component can be measured and the square can be measured, but you can't measure all three components of kinetic momentum, you can't measure all three components of velocity. That's life. Back when kinetic and canonical were the same you might build up wrong ideas about kinetic momentum and now is the time to separate them if you haven't already.
So I tried to address every possible confusion and every thing that often gets glossed over or done ambiguously. So if you have an $A$ you should be able to make your Hamiltonian.
Nothing here seems specific to a Berry-Pancharatnam phase. But there is a nice article in the American Journal of Physics about the Berry phase that solves it exactly, without any appeal to an adiabatic approximation. It could be good yo look at to compare in case you do do an adiabatic approximations. Studying an approximation in a case where you have the exact solution can help to make sure you are following what is going on.
