I have been trying to get eigenvectors from Atomic Simulation Enviornment (ASE) which matches exactly with the eigenvectors coming from GULP(lattice dynamics software,FORTRAN).
GULP produces eigenvectors(it uses LAPACK pacakge).I need to get the same eigenvectors but with ASE(it does not have scipy at its backend,and so no LAPACK). I am using the same k-points in GULP as well as ASE softwares)
But when I include scipy.eigh for solving the eigenvectors in ASE, the eigenvectors do not match with the GULP’s eigenvectors (although the dynamical matrix is matching for both).
Question:
(1)ASE uses “finite-displacement method where the derivatives of the total energy and effective potential are obtained from finite-difference approximations”. Does GULP also uses the same technique?
(a)If yes,then could you tell me the:
(i) Series of **source files** name used to reach the eigenvector and eigenvalues (where LAPACK and EISPACK libraries are called in order to calculate eigenvalue and eigenvectors)?
(ii) The *finite displacement value* to be input (eg. *delta = 0.05 Angs* in ASE)
(b) If no, then the approach used?
(FORTRAN is new to me and so was difficult to traceback inside GULP software as to which LAPACK function is being called).
Here are the answers to your questions:
(1) By default GULP doesn’t use finite differences except for a few models like ReaxFF for example. Use of finite differences can be forced by specifying a finite difference interval with the “finite” option and adding the “numerical” keyword.
NB: No files are written out since everything is done in memory and this would be for gamma only.
(b)GULP uses analytic second derivatives wherever possible, which is why there is no need for finite differences.
As to why eigenvectors might be different, there are different conventions in reporting eigenvectors depending as to whether you use the direct solutions from the eigensolver or change the mass weighting/phase factors. All depends on what you plan to do with them afterwards.
I am attaching a comparative file of Dynamical matrix produced from GULP and ASE. Please find the attached compare_Dyn_mat.dat file.
The eigenvetors coming from the GULP is what I need for my future calculations of thermal conductivity.
The Dynamical matrix from which the eigenvectors are coming from the GULP for a wave-vector say for eg. k = [0.5,0.5,0.25] (with the use of keyword “pfinite” for finite differnce technique as you had suggested earlier) is slightly different than what ASE is producing.
ASE uses the same finite difference technique in order to calculate the Dynamical matrix but the off-diagonal terms sign and small imaginary values appear.
Question 1: Why does GULP do not have any imaginary part in the dynamical matrix by the finite differncing technique?
Question 2:Could you suggest me a possible change that may help me in producing the exact dynamical matrix from ASE as is coming from the GULP?
To address the two questions:
Question 1: Finite differences are used to generate the phonons at the gamma point & so there is no imaginary part to print out.
Question 2: As I don’t know what ASE is doing it’s hard for me to tell you how to reproduce it since this is a question for the authors of that code. As mentioned before the eigenvectors can differ by phase factors depending on the purpose. In GULP the eigenvectors are those that diagonalise the (phased) second derivative matrix divided by the product of the inverse square root masses of the appropriate atoms. You can check the eigenvectors by using them to transform the original dynamical matrix to prove that you obtain the eigenvalues.
The finite difference generation of the second derivatives for the dynamical matrix is performed in dynamicn.F90, whereas the analytic version is dynamic.F90. These are just the high level routines and the derivatives are evaluated in all the subroutines called from these routines. I’m assuming you want the dynamical matrix rather than the Hessian (since the post relates to phonons). If you want the Hessian (for optimisation) then this is generated elsewhere since it is restricted to the variables only and is adapted for symmetry and constraints.