I'm very sorry had mistyped the gnm function, now it is corrected:

gnm[x_, y_, n_, m_] := (2*Lx/(n*Pi))*Sin[n*Pi*(x + Lx)/(2*Lx)]*Sin[m*Pi*(y + Ly)/(2*Ly)]

Adding the Exclusions option, the Integration becomes a little slower, but the result is the same and the error message is still there:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 36 recursive bisections in x near {y,x,y1,x1} = {-0.0115588,-0.0025,-0.0115588,-0.000779382}. NIntegrate obtained -4.16285*10^-6 and 0.000030333190157773638` for the integral and error estimates.

Then, I impose WorkingPrecision -> 30, MaxRecursion -> 100 and the message becomes

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.00052147542424724789220173235532717301517707858296659716371526180090667256502648835`80. and 0.0017515012502079903997882341525304894398325903605063888758977587791494059315495883`80. for the integral and error estimates.

After this, I increase the suggested (by the error message) MaxErrorIncreases option: Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 3000}, and the error message stays the same, just now complaining that the error increased 3000 times, and my NIntegrate command already looks like this:

atonce[i_, j_, n_, m_] :=
NIntegrate[current[x, y, i, j]*current[x1, y1, n, m]/Sqrt[(x - x1)^2 + (y - y1)^2],{y, -Lx, Lx}, {x, -Lx, Lx}, {y1, -Lx, Lx}, {x1, -Lx, Lx},
Exclusions -> (Sqrt[(x - x1)^2 + (y - y1)^2] == 0), Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 3000}, WorkingPrecision -> 30, MaxRecursion -> 100]
atonce1[1, 2, 1, 1] // AbsoluteTiming

Thanks for the suggestions anyway! I have finally arrived at something promising! Basically I changed the method of handling the singularities, and the way I declare the exclusions. Any further suggestions that may speed things up are more than welcome! Here's the code:

Clear[Lx, Ly, gnm, current, try1]
{Lx, Ly} = {1/50, 1/50};

gnm[x_, y_, n_, m_] := (2*Lx/(n*\[Pi]))*Sin[n*\[Pi]*(x + Lx)/(2*Lx)]*Sin[m*\[Pi]*(y + Ly)/(2*Ly)];
current[x_, y_, n_, m_] = D[gnm[x, y, n, m], y];

try[i_, j_, n_, m_] := NIntegrate[jxnm[x, y, i, j]*jxnm[x1, y1, n, m]/Sqrt[(x - x1)^2 + (y - y1)^2],
{y, -Lx, Lx}, {x, -Lx, Lx}, {y1, -Ly, y, Ly}, {x1, -Lx, x, Lx},Method -> {"GlobalAdaptive", "SingularityHandler" -> "DuffyCoordinates"}];

try[1, 1, 1, 1] // AbsoluteTiming

Edit I think it's especially slow when the value tends to 0, how can one deal with this?