a few questions on "Phonons.py" run

[liuzf07]

I have a few questions regarding the use of inelastica for phonon and e-ph coupling calculations (I am using siesta-trunk-433 and the newest inelastica version that I cloned from github - let me know I should use a newer version of siesta):

(1) I recognize that in OSrun, for the AuH2 example given in the distribution of the code, the 6 structures are different only in x, y, and z coordinates for all the atoms (not only the FC active ones). What is the consideration behind this run? Why do we require multiple overlap matrices for the calculation of e-ph coupling? And do we always need these 6 structures with the above displacement pattern (with the corresponding MD.FCdispl) for every system?

(2) It looks like in the example, both TS.HS.Save and TS.SaveHS are set to true. What is the difference between the two keywords? Also, how many .TSHS files (and what are they, as I noticed there are many of these files generated) should be generated in the FCrun directory and what is the motivation behind that?

(3) Does inelastica generate the same vibrational frequencies as the vibrator utility shipped with the siesta distribution, found in Util/Vibra/Src/vibrator? I seem to find slight discrepancies. If there are indeed discrepancies, are they somehow related to the “eliminating interactions through PBCs along the z-direction”?

(4) Generally speaking, is it typically good enough to only include the molecular atoms in the FC runs? Or a few electrode atoms (such as an adatom) need to be included?

(5) If I want to use inelastica to calculate phonon properties of systems other than molecular junctions, how do I set the “device first”, “device last”, “PBC first”, and “PBC last”? Related to this question, in order to reduce computational cost, can I remove those “electrode atoms” (in transiesta sense) and do FCruns only for the device region (so that devicefirst=1 and devicelast=last atom in the coordinates) for a molecular junction? I thought if the device region is long enough, those electrode atoms outside the device region should hardly affect the phonon frequencies and e-ph couplings. Is this correct?

(6) In the output, there are a few files, such as output.mol.real-displ.asxf, output.mol.charlength-displ.asxf, etc, Could someone explain the meaning and the format (also units used) of each file?

(7) Does inelastica enforce the use of “Bohr” in all units (and the output files)? Or it will differentiate Angstrom and Bohr automatically based on the corresponding keywords in the .fdf? I guess it's the latter, but I just want to confirm.

A few other suggestions:

A) I think it's better to write in the instructions that different runs, such as FCrun, OSrun, etc. must reside in different directories starting with the name "FC", "OS", etc, and the “Phonons” run needs to be outside these directories. Also, the default file name is FCrun/RUN.fdf. I believe it will benefit the users if the author could explicitly list the required file names to run in each directory;

B) It would be helpful to tell the user up front that they should use transiesta instead of siesta to run the FCrun and OSrun (if this is correct), and the transiesta has better be compiled with NetCDF4 support;

C) In Phonons.py, the code can only read coordinates that are explicitly included in the RUN.fdf. If the coordinates are given by, e.g., %block .... < coordinate.file.name, then the inelastica does not work. I think this can be improved;

D) The setupOSrun only reads .XV files to generate the RUN_(1-6).fdf - if I do not start with a CG run (e.g., I get the coordinates from someone else), then I do not have the .XV file. I think the SetupOSrun should be able to figure out the coordinates based on the RUN.fdf in the FCrun, instead of relying on a .XV file from geometry optimization. In other words, I feel that the requirement for the file names and files present in each directory can be loosened a lit, to give the users more flexibilities.

[zerothi]

Let me answer the transiesta related questions:

(2) They are equivalent, TS.HS.Save is the keyword that is meant to be there in the future. It is more consistent with the rest of the TS keywords. So prefer that keyword. There will be 1 + N * 3 * 2 TSHS files where N is the number of displaced atoms. Each atom is displaced +-x, +-y and +-z, i.e. 6 times per atom.

B) Since Siesta 4.1 this requirement is lifted. So it is only in older versions of Siesta this holds.

D) Simply use geom2geom which can translate a fdf to XV file format.

The rest I'll leave standing as open questions. :)

[tfrederiksen]

Hi @liuzf07,

Here some answers to your other questions:

(1) The overlaps of basis orbitals between the initial and displaced geometries are required to take into account that the Hamiltonian matrices in the FCruns are given in different bases (basis orbitals follow the displaced atom). See Eqs. (14)-(16) in PRB 75, 205413 (2007). The OSrun is only necessary if the user wants to compute electron-phonon couplings (and not just the vibrational modes and frequencies).

(3) Inelastica performs various symmetry operations on the force constants obtained from SIESTA. I suspect the differences you have seen originate here. Please post more details if you believe there is an issue.

(4) This depends on the system. For a typical molecular junction, with an organic molecule of light atoms suspended between heavy gold electrodes, many molecular vibrations are above the phonon bands of the electrodes and these will be well described by just moving the molecular atoms. However, if you are interested in low-energy modes resonant with bulk phonons or stretch modes against the electrode atoms, it may be necessary to include some of the electrode atoms too. Generally speaking, it is a good idea to check how the IETS spectrum is modified when you enlarge the region of the vibrating atoms.

(5) DeviceFirst=1 and DeviceLast=<last atom N> is allowed. However, the motion of atom 1 may induce a significant change on atom N due to periodic boundary conditions. The PBC keywords are intended to allow the user to remove such couplings across cells (in the z direction only).

(6) These output files are used to visualize the vibration modes in, e.g., XCrysDen. The filename strings have the following meaning:

• .mol. (.per) excludes (includes) the cell vectors in the file,
• .real-displ. indicates that the mode vectors are transformed into real-space displacements, i.e., transformed back from the mass-scaled coordinates (as found in, e.g., the *.mol.axsf and .per.axsf files without this string),
• .charlength-displ. is an experimental development involving the characteristic length of the normal mode.

(7) Internally, Inelastica works with Ångström and eV as units of length and energy. If a user specifies coordinates in Bohr in the fdf, it will be converted upon reading the input.

(A) Some of these are just defaults. The input fdf file can be passed to the scripts with the argument -f <filename> or --fdf <filename>. Similarly, the defaults can be changed with the options --OSdir <OSdir> and --FCwildcard <FC*>. Please try the help functionality, e.g., Phonons -h.

(C) Yes, this kind of input is currently not supported.

[liuzf07]

Thanks a lot, @zerothi and @tfrederiksen . They are very helpful and answered most of my questions. Two remaining ones are the following:

Regarding (1), now I see that they are only needed to compute the derivatives of the overlap matrix elements (thanks for referring me back to that part of the paper which I certainly overlooked the details when I first read it). With that said, does it mean that this displacement we put in OSrun files does not need to be identical to those in FCruns (MD.FCdispl)?

Regarding (6), how about the e-ph couplings? Are they only included in the .nc file? If so, what are the variables and dimensions that I should be looking for (if I use, e.g., netcdf4 python package to read it - or do you recommend other tools?) It would be certainly great to have a document listing the data structures of this .nc file.

Thanks!

[tfrederiksen]

(1) Yes, in principle the derivatives of the overlap can be computed with a different amplitude than the one used in FCrun. However, I see no particular reason to do so.

(6) Yes, the e-ph couplings are only provided in the .nc file. You can readily interrogate the contents of NetCDF files with tools like ncdump.

[liuzf07]

Thanks a lot! Problems all solved for now. I would close this issue.