DynMat calculations

[vladimirovpv]

I've compiled vasp with VTST tools and tried to use this code for the dynamic matrix calculations. As recommended I've set ICHAIN=1, POTIM=0.0, but got an error: * 252 Floating-point zero divide PROG=step ELN=228.
After looking at the source code of step.F I supposed that it could be due to the setting POTIM to zero. Indeed, after changing POTIM=0.000001calculation continued without error.
Questions:
1. Why using recommended settings results in the error?
2. This job is still continuing (due to the small value of POTIM ion relaxations are extremely slow). Should I manually stop it and if yes can I trust the results obtained using the post processing scripts?
3. Script dymmatrix.pl produced non-zero frequencies using OUTCAR of the still continuing job. Where I could find a reference on the format of the output files (eigs.dat, freq.dat, modes.dat and freq.mat) and their interpretation?

For your convenience I reproduce the INCAR file here:

Code: Select all

# dynamical matrix calculation
ICHAIN=1
POTIM = 0.000001
NSW=40 # 3*n+1

ISMEAR = 2
SIGMA = 0.2
ENCUT = 450

NPAR = 4
GGA = 91
ISPIN = 1

IALGO = 48
LREAL = .FALSE.
NSIM = 4

NELMIN = 4
EDIFFG = 1E-4

[andri]

Did you set IBRION = 3 ? The calculations will not run correctly without that.

[vladimirovpv]

Thanks! I missed this point.
Do you know how to interpret the output files?

[andri]

I think the content of the output files is listed on the website ... http://theory.cm.utexas.edu/vtsttools/scripts/#DYM .

The unit for the eigenvalues in freq.dat is cm^{-1}. The freq.mat is the mass-scaled Hessian matrix (using eV and A as units), the modes.dat contains its eigenvectors and eigs.dat are the eigenvalues. The eigenvalues are then converted to cm^{-1} in freq.dat

[vladimirovpv]

Unfortunately, setting IBRION=3 did not help.
I've got slightly modified error message: * 252 Floating-point zero divide PROG=poscar.outpos_trail ELN=795(4000a498c)
It means that division by zero occurred now in the subroutine outpos_trail (module poscar). This is again due to setting of POTIM=0.0.
The line 795, where division by zero occurs is VTMP(1)= DYN%VEL(1,NI)/DYN%POTIM

Then the question 1 from my previous posting is still valid.

Dear developers, could you, please, help me to resolve this problem?

[graeme]

I would like to get to the bottom of this. Here are a few ideas:

Do you have velocities in your POSCAR file? If so, try removing them.

Which version of the vasp source are you using. These POTIM=0 problems used to show up many versions ago, but for 4.6.31 and 4.6.34 (well, and any release in the past year or more), these problems seem to have gone away. I wonder now, looking at these lines in the OUTPOS_TRAIL, if we are relying on getting a zero from a 0/0 calculation. In 4.6.31 and 4.6.34, there are several checks for POTIM==0 to avoid dividing by zero. I don't know the logical flow well enough to know if we should be getting into the OUTPOS_TRAIL with our POTIM=0 set.

Why don't you send us a .tar.gz file of this calculation, and let's see how it runs on our machines. If it runs fine (this is what I expect), we can try to figure out what is different on your setup, and narrow it down to the compiler, platform, or vasp version. With your input file I can add a print statement in OUTPOS_TRAIL to see if we get there, and if so, what is happening at these lines.

[andri]

I suspect we are actually getting a null 0/0 situation. I recall a while back we were getting NaN
in the velocity field of the CONTCAR, but I haven't looked at it recently. It's then up to the compiler
how floating point exceptions are handled. Judging from your error-output you're not using ifort.
Is there an option with your compiler to set how exceptions like this one are handled? Of course
a preferable solution is to avoid the situation, but this might allow you to finish your calculation. Another
options is, if your system isn't too large, to try the IBRION = 5 option (see the official VASP manual for
that one).

[vladimirovpv]

1. Do you have velocities in your POSCAR file? If so, try removing them.
Yes I have velocities in POSCAR. I've removed them and got the same error.

2. Which version of the vasp source are you using?
vasp.4.6.31
VTST: version 2.03c (08/15/07)
Is vasp 4.6.34 publically available?

3. Why don't you send us a .tar.gz file of this calculation?
I've sent the file to you.

[vladimirovpv]

I recompiled the code with nozdiv option as Andri suggested. The error messages changed as follows:
Before recompilation:
* 252 Floating-point zero divide PROG=initio ELN=237(40045d9b8)
* 252 Floating-point zero divide PROG=initio ELN=237(40045d9b8)
* 252 Floating-point zero divide PROG=initio ELN=237(40045d9b8)
* 252 Floating-point zero divide PROG=initio ELN=237(40045d9b8)
* 252 Floating-point zero divide PROG=poscar.outpos_trail ELN=795(4000a498c)
* 253 Invalid operation PROG=poscar.outpos_trail ELN=797(4000a4a18)
* 253 Invalid operation PROG=poscar.outpos_trail ELN=797(4000a4a18)
After recompilation:
* 253 Invalid operation PROG=poscar.outpos_trail ELN=797(4000a4998)
* 253 Invalid operation PROG=poscar.outpos_trail ELN=797(4000a4998)
* 253 Invalid operation PROG=poscar.outpos_trail ELN=797(4000a4998)

The task still fails to finish.
Any suggestions are appreciated.

[graeme]

Thanks for the files - I'll try running this when I get back to the office.

First though, since you have chosen 'selective dynamics' in the poscar file, try adding ' T T T' after each line with coordinates in it. There may be some problem with reading that poscar.

Also, looking at your INCAR: if you set a low force criteria, you will make sure that the dynmat iterations are not terminated early due to reaching the convergence criteria. Also, you will need to ask for precise forces for the finite difference. I suggest using:

EDIFFG = -1E-7
EDIFF=1e-8

[graeme]

Update:

When I ran your calculation directly, it hangs with the statement:
No initial positions read in

When I add those selective dynamics flags (T T T), it started properly.

[vladimirovpv]

Thanks, Graeme.
The trick was to use "Selective dynamics" and "T T T" flags.
When I removed them the same division by zero occurs.

Now the task runs, but ...
In spite of the fact that I set
EDIFFG = -1E-9
EDIFF=1e-10
calculation stops after the first ionic relaxation step ("aborting loop because EDIFF is reached"). Only one call to DYNMAT is performed:
DOING DYMAT, READING DISPLACECAR
FOR 1 OUT OF 1
HIPREC POSITION (A)
-----------------------------------------------------------------------------------
DYNMAT: Loop at J= 1
DYNMAT: Before DIRKAR
DYNMAT: After DIRKAR
-4.508072507371597 7.821228892928404 10.647821387491449
DYNMAT: Loop at J= 2
DYNMAT: Before DIRKAR
DYNMAT: After DIRKAR
-0.095760356486480 1.248947042673503 1.735996852146915
...
HIPREC TOTAL-FORCE (eV/A)
-----------------------------------------------------------------------------------
0.000424468422787 0.000336111899726 0.000369259718467
0.001755498975477 0.001013537806054 0.000363648648985
0.000000000000000 0.000175999433220 0.001016033137806
0.000000000000000 -0.000448313962350 0.001698927762991
-0.000424468422787 0.000336111899726 0.000369259718467
-0.001755498975477 0.001013537806054 0.000363648648985
...
DYMAT: ******************************
DYMAT: DISPLACEMENT 0 0
DYMAT: ENERGY -358.293162
DYMAT: ------------------------------
DYMAT: DEGREE OF FREEDOM 4 1 1
DYMAT: ------------------------------
DYMAT: DEGREE OF FREEDOM 4 2 2
DYMAT: ------------------------------
DYMAT: DEGREE OF FREEDOM 4 3 3
DYMAT: ------------------------------
...
DYMAT: DEGREE OF FREEDOM 97 2 38
DYMAT: ------------------------------
DYMAT: DEGREE OF FREEDOM 97 3 39
DYMAT: LEAVING
FORCES: max atom, RMS 0.008750 0.002370
FORCE total and by dimension 0.023343 0.008750

Here the OUTCAR stops and dymmatrix.pl produces all zero frequencies. Are the displacements of 0.005 too small?

[andri]

Is it possible that you're running into problems with using so many k-points (13x13x13) but so
few nodes (NPAR = 4 & NSIM = 4)? What happens if you reduce the number of k-points, say 5x5x5?

[graeme]

I think that the problem could be due to a symmetry setting.
For dynamical matrix calculations, you need to turn symmetry off: ISYM=0

[vladimirovpv]

Andry, I do not think that the problem is in the number of k-points, as normal VASP relaxation proceeds without problems.
Yes, I am using only 4 CPU, but on rather powerful vector parallel machine (NEC SX8R).

Graeme, I set ISYM=0 and restarted the task. The result is the same: the task stops after the first ion relaxation step. Is it normal behavior? And what behavior is anticipated?
Dymmatrix produces all zero frequencies again.

I suspect that the problem could be somewhere else: e.g. in preparation of DISPLACECAR. I used dymselsph.pl to restrict vibration mode calculation only to those atoms, which are around interstitial defect.

Any suggestions?