UT theoretical chemistry code forum

Hello everybody!

I've searched the forum for answers, but answers to similar problems did not help. Excuse me in advance if this has been posed before.

After performing a dynamical matrix calculation I obtain two freq.dat files corresponding to the equilibrium and the transition state of a particular migration. When trying to obtain the prefactor, I get:

> dymprefactor.pl freq1.dat freq2.dat
Illegal division by zero at dymprefactor.pl line 54.

At line 54, the script shows

$RatioEig *= ($line1[0]/$line2[0]);

Looking at the script, I thought it could be that a frequency that is not zero in the first file is divided by a zero in the second. So I tried to delete all lines in both files that happened to be in that situation. But the error message persists.

I would be very grateful for any suggestion. Please feel free to request more info if necessary.

Thanks in advance!

Roberto

If you post your two frequency files, I'll see what's going wrong.

Hi!

Thanks a lot for your quick reply! Please find below the two frequency files, being freq1.dat the one corresponding to the equilibrium state and freq2.dat the one for the transition state.

freq1.dat:

46.628538 cm^{-1} ... 1
44.812202 cm^{-1} ... 1
0.000002 cm^{-1} ... 1
0.000001 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000001 cm^{-1} ... 0
0.000003 cm^{-1} ... 0
0.000003 cm^{-1} ... 0
140.497444 cm^{-1} ... 0
160.059569 cm^{-1} ... 0

freq2.dat:

46.271306 cm^{-1} ... 1
44.204973 cm^{-1} ... 1
0.000003 cm^{-1} ... 1
0.000002 cm^{-1} ... 1
0.000002 cm^{-1} ... 1
0.000001 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 1
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000000 cm^{-1} ... 0
0.000001 cm^{-1} ... 0
143.097903 cm^{-1} ... 0
162.122181 cm^{-1} ... 0

Thanks a lot for your help!

Kind regards

Roberto

Ah, there is something very wrong with your frequency calculations. There are far too many degrees of freedom with zero modes.

Check first that the vasp calculation that was used to generate the Hessian has completed the correct number of force/energy evaluations. If the calculation ended early, there could be many zeros in the Hessian matrix. If you are using our dynamical matrix code, check to see if the calculation terminated due to reaching the specified convergence setting. If that is the case, set ediffg=-1e-10 so that vasp will complete all finite difference calculations even if the force is low.

But again, the problem is related to the Hessian matrix calculation and the modes, rather than the prefactor calculation.

Thank you very much, Graeme! We'll try and do as you suggest, since there are a number of zero rows in the dynamical matrix as seen in modes.dat.
Just one more thing: I think you mean setting EDIFFG=1E-10, with a "plus", not a "minus", is this correct?

Cheers!
Roberto

Hello again!

After setting EDIFFG=1E-10 according to your advice, there are no zeros in the Hessian matrix, so that all finite difference calculations are complete. Now we are left with two freq.dat files, one for the equilibrium and the other for the saddle point, both having three imaginary frequencies. Should we simply delete them since they should be related to translational modes? We performed the calculations using a DISPLACECAR with finite small displacements for all the atomic coordinates. If we do that, the saddle point would have no imaginary frequencies at all and if I remember correctly, it should have at least one. I paste the content of the two frequency files below, in case they can be helpful:

freq1.dat (equilibrium)

434.162238 cm^{-1} ... 1
30.645143 cm^{-1} ... 1
20.783570 cm^{-1} ... 1
80.692775 cm^{-1} ... 0
95.776080 cm^{-1} ... 0
104.551240 cm^{-1} ... 0
117.347212 cm^{-1} ... 0
119.678558 cm^{-1} ... 0
125.420911 cm^{-1} ... 0
131.149286 cm^{-1} ... 0
131.466655 cm^{-1} ... 0
132.945374 cm^{-1} ... 0
139.542806 cm^{-1} ... 0
141.787311 cm^{-1} ... 0
143.058077 cm^{-1} ... 0
146.186204 cm^{-1} ... 0
151.021756 cm^{-1} ... 0
153.067789 cm^{-1} ... 0
158.128922 cm^{-1} ... 0
167.180574 cm^{-1} ... 0
169.934848 cm^{-1} ... 0
172.322479 cm^{-1} ... 0
174.041873 cm^{-1} ... 0
178.991087 cm^{-1} ... 0
180.917631 cm^{-1} ... 0
185.456012 cm^{-1} ... 0
186.984125 cm^{-1} ... 0
189.042042 cm^{-1} ... 0
190.173936 cm^{-1} ... 0
195.480411 cm^{-1} ... 0
199.643105 cm^{-1} ... 0
202.415537 cm^{-1} ... 0
206.265648 cm^{-1} ... 0
208.756349 cm^{-1} ... 0
212.402052 cm^{-1} ... 0
217.196335 cm^{-1} ... 0
221.519066 cm^{-1} ... 0
222.866453 cm^{-1} ... 0
223.604743 cm^{-1} ... 0
225.843226 cm^{-1} ... 0
234.205940 cm^{-1} ... 0
242.488028 cm^{-1} ... 0
259.020825 cm^{-1} ... 0
276.078625 cm^{-1} ... 0
481.929513 cm^{-1} ... 0

freq2.dat (saddle)

421.647527 cm^{-1} ... 1
39.595782 cm^{-1} ... 1
20.197448 cm^{-1} ... 1
57.787814 cm^{-1} ... 0
81.863918 cm^{-1} ... 0
86.813537 cm^{-1} ... 0
107.543002 cm^{-1} ... 0
112.683635 cm^{-1} ... 0
114.385697 cm^{-1} ... 0
122.366075 cm^{-1} ... 0
125.665857 cm^{-1} ... 0
127.476201 cm^{-1} ... 0
132.407272 cm^{-1} ... 0
133.868668 cm^{-1} ... 0
134.658598 cm^{-1} ... 0
139.845632 cm^{-1} ... 0
143.341992 cm^{-1} ... 0
150.248728 cm^{-1} ... 0
154.478528 cm^{-1} ... 0
160.918225 cm^{-1} ... 0
164.413568 cm^{-1} ... 0
167.080485 cm^{-1} ... 0
170.047050 cm^{-1} ... 0
175.107874 cm^{-1} ... 0
177.978101 cm^{-1} ... 0
180.695766 cm^{-1} ... 0
181.361999 cm^{-1} ... 0
184.908142 cm^{-1} ... 0
185.997860 cm^{-1} ... 0
190.097774 cm^{-1} ... 0
194.982785 cm^{-1} ... 0
196.564978 cm^{-1} ... 0
200.790759 cm^{-1} ... 0
203.680527 cm^{-1} ... 0
209.157029 cm^{-1} ... 0
213.572175 cm^{-1} ... 0
218.843799 cm^{-1} ... 0
220.215700 cm^{-1} ... 0
222.909768 cm^{-1} ... 0
224.679642 cm^{-1} ... 0
230.636342 cm^{-1} ... 0
244.305616 cm^{-1} ... 0
264.669228 cm^{-1} ... 0
288.379642 cm^{-1} ... 0
453.762407 cm^{-1} ... 0

Thank you again!

Roberto

Unfortunately, the -400 mode at the minimum can not be ignored or deleted. The most likely explanation is that your minimum is not fully converged. Try setting ediff=1e-6 and try to get convergence at the ediffg=-0.005 or even 0.001 level. Make sure that ISYM=0. Then use ediff=1e-8 and ediffg=1e-10 for the finite difference calculation. The three translational modes at the minimum should have absolute values < 100 cm^-1.

You will probably have to do the same at the saddle to get three small modes and one more clearly negative mode. Only when you have this will the prefactor calculation make sense.

One last thought: I think the fact that both calculations have a ~ -400 mode is indicating something, but I don't know what. Perhaps check that the Hessian is symmetric -- at least roughly by inspection.

Hi!

Thanks a lot again for your input, Graeme! I have to say that I didn't find any obvious symmetry in the Hessian matrix, at least by inspection. So I tried to set EDIFF=1E-6 and and EDIFFG=-0.001, with ISYM=0. After convergence, I have used EDIFF=1E-8 and EDIFFG=1E-10 for the finite difference calculation, both for the equilibrium and the saddle points. These are now the frequencies I found, well below 100 cm^-1:

freq1.dat (equilibrium):

37.562250 cm^{-1} ... 1
36.309112 cm^{-1} ... 1
33.286362 cm^{-1} ... 1
24.363433 cm^{-1} ... 1
18.214586 cm^{-1} ... 1
17.372552 cm^{-1} ... 1
8.943577 cm^{-1} ... 1
3.868448 cm^{-1} ... 1
2.942990 cm^{-1} ... 1
81.490834 cm^{-1} ... 0
86.305779 cm^{-1} ... 0
86.397974 cm^{-1} ... 0
87.140818 cm^{-1} ... 0
89.667391 cm^{-1} ... 0
90.228575 cm^{-1} ... 0
90.522423 cm^{-1} ... 0
122.082209 cm^{-1} ... 0
122.348162 cm^{-1} ... 0
122.654460 cm^{-1} ... 0
124.049703 cm^{-1} ... 0
124.492376 cm^{-1} ... 0
124.632002 cm^{-1} ... 0
132.706376 cm^{-1} ... 0
145.459785 cm^{-1} ... 0
146.762930 cm^{-1} ... 0
148.363984 cm^{-1} ... 0
158.605797 cm^{-1} ... 0
159.157169 cm^{-1} ... 0
159.227567 cm^{-1} ... 0
159.379854 cm^{-1} ... 0
159.699933 cm^{-1} ... 0
166.889195 cm^{-1} ... 0
167.017335 cm^{-1} ... 0
167.399744 cm^{-1} ... 0
176.642032 cm^{-1} ... 0
176.766234 cm^{-1} ... 0
176.819091 cm^{-1} ... 0
178.071498 cm^{-1} ... 0
178.622169 cm^{-1} ... 0
190.945463 cm^{-1} ... 0
191.949127 cm^{-1} ... 0
192.358611 cm^{-1} ... 0
211.680582 cm^{-1} ... 0
211.816572 cm^{-1} ... 0
211.912224 cm^{-1} ... 0

freq2.dat (saddle):

6.863244 cm^{-1} ... 1
2.461273 cm^{-1} ... 1
1.479681 cm^{-1} ... 1
53.276502 cm^{-1} ... 0
55.755411 cm^{-1} ... 0
57.630405 cm^{-1} ... 0
62.257022 cm^{-1} ... 0
63.010138 cm^{-1} ... 0
63.734515 cm^{-1} ... 0
102.451306 cm^{-1} ... 0
104.165971 cm^{-1} ... 0
104.200718 cm^{-1} ... 0
107.316190 cm^{-1} ... 0
110.511269 cm^{-1} ... 0
110.886137 cm^{-1} ... 0
111.643278 cm^{-1} ... 0
136.969733 cm^{-1} ... 0
137.351782 cm^{-1} ... 0
138.422457 cm^{-1} ... 0
139.295848 cm^{-1} ... 0
140.865138 cm^{-1} ... 0
141.011164 cm^{-1} ... 0
146.618687 cm^{-1} ... 0
160.194901 cm^{-1} ... 0
161.420275 cm^{-1} ... 0
162.729913 cm^{-1} ... 0
169.823629 cm^{-1} ... 0
170.367050 cm^{-1} ... 0
171.194343 cm^{-1} ... 0
171.545273 cm^{-1} ... 0
173.484833 cm^{-1} ... 0
179.323106 cm^{-1} ... 0
179.683457 cm^{-1} ... 0
181.693145 cm^{-1} ... 0
188.189234 cm^{-1} ... 0
188.369030 cm^{-1} ... 0
188.771148 cm^{-1} ... 0
189.786261 cm^{-1} ... 0
190.407225 cm^{-1} ... 0
202.412371 cm^{-1} ... 0
203.131103 cm^{-1} ... 0
203.636639 cm^{-1} ... 0
221.019022 cm^{-1} ... 0
221.287542 cm^{-1} ... 0
221.369403 cm^{-1} ... 0

Do they appear as corrrect? I don't know if thhee first modes for the saddle point are small enough and the first one is the one you said that should be clearly more negative. On top of that, the 9 imaginary modes at the equilibrium position mistify me quite a lot, since without freezing any atom at all, I would expect to have just three of them.

Thanks again for any suggestion!

Roberto

This does not look good. If you upload your min and saddle calculations I might be able to check for problems.

Dear Graeme,

please find attached two tar.gz files with the main files from the equilibrium and saddle dynamical matrix calculations. Please let me know if you need something else.

Thanks a lot for your keen interest in my problem!

Roberto[attachment=0]sad.tar.gz[/attachment][attachment=1]eq.tar.gz[/attachment]

The calculation in the saddle directory looks great. The equilibrium calculation is messed up, but I have not yet figured out what's wrong -- everything looks good and the settings are the same as the saddle calculation. Just one additional general suggestion is to set LREAL=.FALSE. for more accurate forces.

More importantly, however, is that your saddle geometry appears to be a minor translation of your equilibrium structure. Thus, I expect that the saddle is really a minimum and equivalent to your equilibrium geometry.

Thank you very much for your quick reply!

Looking closely at the saddle, you are right, apparently the NEB calculation found another version of the equilibrium configuration, obtained by a translation of about 0.0066 Ang. However, why is it then that the equilibrium dynamical matrix calculation is not correct? This is indeed puzzling!

Roberto

Dear Graeme,

Much to my amazement, after including LREAL=.FALSE. in the INCAR file for the equilibrium finite difference calculation, I obtained just three imaginary frequencies, that if I understand correctly correspond to the translations. If I delete those three modes, and since both freq.dat files should have the same number of modes to correctly apply the dymprefactor.pl script, I guess in principle I should delete the corresponding first three modes at the saddle point as well. However, the saddle point should have one imaginary mode left, and there are only precisely those three on top of the frquency file. So my question would be, which modes should I delete at the saddle point? I paste below both frequency files for convenience:

freq1.dat (equilibrium point):

2.211535 cm^{-1} ... 1
1.504660 cm^{-1} ... 1
0.875579 cm^{-1} ... 1
61.207587 cm^{-1} ... 0
61.722839 cm^{-1} ... 0
63.346729 cm^{-1} ... 0
63.443165 cm^{-1} ... 0
63.470131 cm^{-1} ... 0
63.844567 cm^{-1} ... 0
103.303453 cm^{-1} ... 0
106.247530 cm^{-1} ... 0
106.567140 cm^{-1} ... 0
106.670935 cm^{-1} ... 0
114.782854 cm^{-1} ... 0
115.159988 cm^{-1} ... 0
115.383027 cm^{-1} ... 0
139.640850 cm^{-1} ... 0
139.871726 cm^{-1} ... 0
140.150218 cm^{-1} ... 0
142.238428 cm^{-1} ... 0
142.648165 cm^{-1} ... 0
142.684831 cm^{-1} ... 0
145.972095 cm^{-1} ... 0
162.169478 cm^{-1} ... 0
163.057311 cm^{-1} ... 0
164.046054 cm^{-1} ... 0
171.652296 cm^{-1} ... 0
172.319739 cm^{-1} ... 0
173.828317 cm^{-1} ... 0
173.881831 cm^{-1} ... 0
174.332567 cm^{-1} ... 0
181.466990 cm^{-1} ... 0
181.613144 cm^{-1} ... 0
181.948467 cm^{-1} ... 0
190.458632 cm^{-1} ... 0
190.574435 cm^{-1} ... 0
190.628240 cm^{-1} ... 0
191.605674 cm^{-1} ... 0
192.179256 cm^{-1} ... 0
204.069180 cm^{-1} ... 0
204.944209 cm^{-1} ... 0
205.341459 cm^{-1} ... 0
222.918778 cm^{-1} ... 0
223.085943 cm^{-1} ... 0
223.177381 cm^{-1} ... 0

freq2.dat (saddle point, as re-calculated as well with LREAL=.FALSE.):

2.193931 cm^{-1} ... 1
1.630131 cm^{-1} ... 1
0.853173 cm^{-1} ... 1
61.285624 cm^{-1} ... 0
61.741450 cm^{-1} ... 0
63.380363 cm^{-1} ... 0
63.419169 cm^{-1} ... 0
63.519849 cm^{-1} ... 0
63.905098 cm^{-1} ... 0
103.368422 cm^{-1} ... 0
106.271225 cm^{-1} ... 0
106.564234 cm^{-1} ... 0
106.721074 cm^{-1} ... 0
114.660668 cm^{-1} ... 0
115.174286 cm^{-1} ... 0
115.395858 cm^{-1} ... 0
139.536576 cm^{-1} ... 0
139.818656 cm^{-1} ... 0
140.101946 cm^{-1} ... 0
142.241464 cm^{-1} ... 0
142.640068 cm^{-1} ... 0
142.747110 cm^{-1} ... 0
146.069431 cm^{-1} ... 0
162.269625 cm^{-1} ... 0
163.054037 cm^{-1} ... 0
163.996223 cm^{-1} ... 0
171.676552 cm^{-1} ... 0
172.325848 cm^{-1} ... 0
173.821065 cm^{-1} ... 0
173.990820 cm^{-1} ... 0
174.349602 cm^{-1} ... 0
181.405004 cm^{-1} ... 0
181.601569 cm^{-1} ... 0
181.977833 cm^{-1} ... 0
190.482601 cm^{-1} ... 0
190.545142 cm^{-1} ... 0
190.601962 cm^{-1} ... 0
191.613240 cm^{-1} ... 0
192.250427 cm^{-1} ... 0
204.624934 cm^{-1} ... 0
204.813618 cm^{-1} ... 0
205.395036 cm^{-1} ... 0
222.976702 cm^{-1} ... 0
223.071737 cm^{-1} ... 0
223.113123 cm^{-1} ... 0

Thank you again!

Roberto