(1) Which subroutine of LAPACK package is GULP using in order to calculate the eigenvectors and eigenvalues ?

(2) As I am not familiar enough with fortran so for the same subroutine (having the same dynamical matrix ) can I use Python(scipy) in order to calculate exactly the same set of eigenvectors as we are getting from GULP ?
I believe I need to call some specific subroutine inorder to calculate the eigenvectors for my Hermitian dynamical matrix.

(3) If not python then could you suggest me other ways by which I could get the same set of eigenvectors?

(1) Multiple routines are used for eigenvector/eigenvalue calculation depending on the algorithm and whether youâ€™re dealing with gamma or not. The standard routines are dsyevd or zheevd in serial.
(2) Iâ€™m sure that scipy will provide eigensolver routines, but youâ€™ll have to ask a Python person for more details.
(3) There are plenty of other maths libraries available to compute eigenvectors - just go to netlib.org

I am dealing with a reciprocal space point other than Gamma so the dynamical matrix is a Hermitian matrix. In order to get the same set of eigen values and eigenvectors I did the following (just for testing):
I took the same dynamical matrix given by GULP and took it as an input for my own program which produces eigenvalue and eigenvector using the same â€śzheevdâ€ť(used the same parameters as used in GULP) subroutine as suggested by you earlier. The eigenvalues I am getting is the same but the eigenvetors are coming out to be different.
My aim is to get the same set of eigenvectors as produced by GULP. I skimmed through the â€śpdiagl.F90â€ť subroutine which calls â€śzheevdâ€ť subroutine which I hope is the one which is being used as one of the subroutines to get the eigenvales and eigenvectors by GULP.
Am I missing something which GULP does after getting the eigenvectors from the â€śzheevdâ€ť subroutine and may be thatâ€™s why my eigenvectors are not matching?

You may want to look through the other posts on this forum since itâ€™s a pretty common question. One possibility is if you have degeneracy then the eigenvectors can mix freely and so the answer you get depends arbitrarily on the numerics of the computer you are using. A linear combination is still a valid set of eigenvectors with the same eigenvalue, so no point trying to match, though you can always mix your eigenvectors to get the other combination.