Computation of IR spectra form NVE molecular dynamics.

[alberto_ricchebuono]

Hi,

I am having issues while post-treating data from NVE molecular dynamics with the aim of extracting the IR spectrum of a material to take into account anharmonic effects.
After running the calculation for the first 3000 steps with the INCAR pasted below, I find that the dipole moments at each step are not directly printed in the OUTCAR. Instead, I get three files (namely "fort.77", "fort.78", "fort.79") which I think they should contain the information I need. The problem is that I cannot find any description of these files on the manual or on the forum and I am having troubles at interpreting them. The file "fort.77" seem to contain time-dependent data, since the first column increases of 0.5 at each row as the time step, but then I do not understand why it stops at 750, while the whole simulation counts 3000 steps (i.e., 1500 fs).

Here's the INCAR:

SYSTEM = F4_MIL53
#Functional
GGA = PE
#Electronic relaxation
PREC = NORMAL
ALGO = FAST
ENCUT = 500.00 eV
LREAL = Auto
NELMIN = 4
NELM = 100
EDIFF = 1E-6
ISPIN = 1
ISYM = 0
SYMPREC= 1E-4
IWAVPR = 11
#Writing info
LEPSILON = .TRUE.
LCALCEPS = .TRUE.
NWRITE = 3
#File saving
LWAVE =.FALSE.
LCHARG =.FALSE.
LVTOT =.FALSE.
#Parallelization
NCORE = 1
LPLANE =.TRUE.
NSIM = 4
#Molecular Dynamics
IBRION = 0
MDALGO = 2
ISIF = 2
TEBEG = 300
NSW = 3000
POTIM = 0.5
SMASS =-3
#Dispersion
IVDW = 12
#Smearing
SIGMA = 0.01
ISMEAR = 0

Here you have the first rows of the three files as an example:

Fort.77:
0.0000000 0.0050702 0.0016978 0.0003656 0.0005242 0.0002606 0.0002770
0.5000000 0.0050276 0.0016902 0.0003619 0.0005418 0.0002656 0.0002981
1.0000000 0.0049188 0.0016659 0.0003511 0.0005561 0.0002695 0.0003164
1.5000000 0.0047461 0.0016260 0.0003342 0.0005669 0.0002723 0.0003316
2.0000000 0.0045127 0.0015720 0.0003126 0.0005741 0.0002739 0.0003431
Fort.78:
0.0000000 1.6665172 1.0816635 0.2121455 0.3730962 0.1161474 0.0971492
2.0943951 1.6665172 1.0816635 0.2121455 0.3730962 0.1161474 0.0971492
4.1887902 1.4753037 0.9609274 0.1959481 0.3009773 0.0516829 0.0932578
6.2831853 1.4134017 0.9224934 0.1902445 0.2812916 0.0399658 0.0940523
8.3775804 1.4848128 0.9740355 0.1956911 0.3130827 0.0705979 0.0971755
Fort.79:
0.0000000 7.0403036 3.4572350 0.6866060 5.2825771 6.2362398 0.9738934
2.0943951 7.0403036 3.4572350 0.6866060 5.2825771 6.2362398 0.9738934
4.1887902 3.6369920 1.8300958 0.3589595 2.6801056 3.1444588 0.4876858
6.2831853 2.3810364 1.2238514 0.2421737 1.7647925 2.0695996 0.3226839
8.3775804 1.7568648 0.9401773 0.1860434 1.3119908 1.5423208 0.2414113

PS: The MD run was pre-equilibrated in a NVT ensamble and I am planning to extend the NVE run to at least 10 ps to maximize the S/N ratio.

[fabien_tran1]

Hi,

First, it may be that the combination that you are trying to use (LEPSILON=.TRUE. with IBRION=0) is not fully supported, as the presence of the files fort.7x may indicate. Then, regarding these files, some information is written in the file linear_response.F within the block:

Code: Select all

    IF (NSTEP==NSW .AND. SQRT(DFIELD(1)**2+DFIELD(2)**2+DFIELD(3)**2)==0) THEN

which contains the WRITE statements.

Finally, the time scale of the correlation function (750 fs), which is half of 1500 fs, comes from the fact that it has to be averaged over different starting times. Thus, only half of the simulation time can be taken, otherwise differently converged averages for correlation functions starting after half the total simulation time would be obtained. You may find some explanations in https://doi.org/10.1103/PhysRevMaterials.2.073604 (Eq. (10)) and in https://doi.org/10.1103/PhysRevB.105.024302 (Eq. (18)).