How GMRES method finds smallest singular value and the corresponding singular vectors of a matrix?

I kinda see what you are asking. There is a relationship between singular values of a matrix $A$ and the matrix $A^HA$: $sigma_i^2(A)=lambda_i(A^HA)$. So in theory, you can use inverse power method to find $i$-th smallest eigenvalue of the matrix $A^HA$ and the corresponding eigenvector. To 'invert' the matrix $A^HA$ (during inverse power method algorithm) you can use GMRES + some preconditioner. Similarly, you can use LOBPCG, which is a CG method designed to find generalized eigenvalues and eigenvectors. Check the book Numerical Methods for Large Eigenvalue Problems - 2nd Edition by Yousef Saad for more information.

There is also an indirect relationship between singular values of a matrix and GMRES. Turns out that the eigenvalues of a matrix is not a good predictor of convergence of GMRES, it actually might be even deceptive. There is a huge body of work surrounding this idea. Some alternatives suggested are pseudospectra and field of values, both are better predictors. Latter is a set that includes singular values and here is some reading about field of values (or numerical range):

1. Horn, R.; Johnson, C.R. Matrix Analysis; Cambridge University Press: New York, NY, USA, 1985
2. Zachlin, P.F. On the Field of Values of the Inverse of a Matrix. Ph.D. Thesis, Case Western Reserve University, Cleveland, OH, USA, 2007.

If anyone else has any references to suggest, you can comment below and I will add to this list.

How GMRES method finds smallest singular value and the corresponding singular vectors of a matrix?

It doesn't. GMRES solves linear systems. Your citation probably refers to other Krylov methods: restarted Arnoldi, Golub-Kahan bidiagonalization.