UT theoretical chemistry code forum

I do dynamical matrix calculation and don’t sure my procedure is correct. Please help me check.
1. Initial data: I take initial and transition state got by previously NEB result to be my initial and finial state of dynamical matrix calculation.
2. Producing file: I make 4 image file between initial and transition state by nebmake.pl, so the NN = 4.
3. Making DISPLACECAR: Using dymseldsp.pl to produce the DISPLACECAR. The usage is dymseldsp.pl (POSCAR 1) (POSCAR 2) (atoms to include) (displacement). I take initial and Transition state as POSCAR1 and POSCAR 2, respectively. The atoms to include set total number atoms. The displacement takes maximizing displacement between initial state and transition.
4. NSW setting: My total number of atoms is 39, As a result of fixing z- direction of 10 atoms, the degrees of freedom is 39x3-10 = 107.the NSW is 107/4+1 = 27.25, which one I should choose, 27 or 28?

If my procedure is correct, how cloud I do to get prefactor. should I make use of dymprefactor.pl? The usage of dymprefactor.pl is dymprefactor.pl (freq.dat of minimum) (freq.dat of transition state). which image file I could find the freq.dat of minimum and freq.dat of transition state?

By the way,could you take me more detail of dymseldsp.pl. Does the atoms to include mean total number of atoms? Does the displacement mean the max displacment? should i fix some atoms when calculating dynamical matrix calculation or NEB or Dimer method?
thanks for you reply

Your method does not sound quite right. The first step is correct - finding the minimum and saddle point. Now, do 2 separate (single image) dynamical matrix calculations, one for each of these states. You can use the dymseldsp.pl script to generate a DISPLACECAR, but since you do not have too many atoms in the cell, it would be safer to make it by hand. Try modifying your DISPLACECAR by putting a displacement of about 0.005 for each non-frozen degree of freedom (107). Now, set NSW=108 in your INCAR. Also, make sure to set ediff=1e-6 or 1e-7.

For both states (initial and saddle), use the dymmatrix.pl script to extract the dynamical matrix from your calculations. Then you can use the dymprefactor.pl to find the ratio of frequencies between these states, and the classical TST prefactor.

Dear Graeme:
thanks for your response. because i don't quiet understand the structure of DISPLACECAR. so I don't understand you say '"displacement of about 0.005 for each non-frozen degree of freedom '" is to edit which part of DISPLACECAR. the content of my DISPLACECAR is:

1 2 3 4 5
0.87283 0.87283 0.87283 22 0.339738797066248
0.87283 0.87283 0.87283 23 0.202063631301066
0.87283 0.87283 0.87283 24 0.206673667532054
0.87283 0.87283 0.87283 25 0.0842033300284822
0.87283 0.87283 0.87283 26 0.0649211115820173

you mean i can edit the DISPLACECAR by hand on fifth column ?
could you explain other column meaning roughly? I deeply appreciate your help.

The DISPLACECAR file requires only a list of 3xN numbers, indicating how much to displace each atom along each degree of freedom. Columns 4 and 5 are just for reference - they indicate the atom number and some distance, for example how much that atom has moved in a reaction. These two columns are not required.

Finite difference displacements of 0.87 (in your example) are much too large. Start with a file containing a 3xN list of displacements all set to 0.005. Then, set the displacements along all the frozen directions to 0. The total number of non-zero values should then be 107. This should be enough to do a dynamical matrix calculation.

Dear Graeme:
thanks for your help, i have finished the dynamical matrix calculation and have produced the freq.dat of minimum and freq.dat of transition state by dymmatrix.pl.but The message appear "Initial state freq.dat has imaginary frequencies" when i use dymprefactor.pl. my content of initial state freq.dat is :

941.341742 cm^{-1} ... 1 337.796850 cm^{-1} ... 0 508.951826 cm^{-1} ... 0
148.560811 cm^{-1} ... 1 340.269190 cm^{-1} ... 0 509.881012 cm^{-1} ... 0
110.037403 cm^{-1} ... 1 341.960290 cm^{-1} ... 0 510.510215 cm^{-1} ... 0
101.811833 cm^{-1} ... 1 352.042069 cm^{-1} ... 0 518.323032 cm^{-1} ... 0
47.278062 cm^{-1} ... 1 355.501714 cm^{-1} ... 0 532.113501 cm^{-1} ... 0
65.144570 cm^{-1} ... 0 359.681937 cm^{-1} ... 0 535.750753 cm^{-1} ... 0
89.273846 cm^{-1} ... 0 362.608675 cm^{-1} ... 0 537.159422 cm^{-1} ... 0
114.066732 cm^{-1} ... 0 366.363141 cm^{-1} ... 0 549.426563 cm^{-1} ... 0
129.746334 cm^{-1} ... 0 369.726918 cm^{-1} ... 0 550.733288 cm^{-1} ... 0
149.837520 cm^{-1} ... 0 372.665459 cm^{-1} ... 0 553.186060 cm^{-1} ... 0
161.264727 cm^{-1} ... 0 376.901599 cm^{-1} ... 0 558.218072 cm^{-1} ... 0
176.045045 cm^{-1} ... 0 384.248972 cm^{-1} ... 0 567.707400 cm^{-1} ... 0
190.255533 cm^{-1} ... 0 387.272226 cm^{-1} ... 0 581.024128 cm^{-1} ... 0
195.993649 cm^{-1} ... 0 394.152708 cm^{-1} ... 0 592.190292 cm^{-1} ... 0
208.587467 cm^{-1} ... 0 397.143846 cm^{-1} ... 0 594.623996 cm^{-1} ... 0
214.608097 cm^{-1} ... 0 405.598287 cm^{-1} ... 0 602.290251 cm^{-1} ... 0
222.628267 cm^{-1} ... 0 408.235954 cm^{-1} ... 0 608.988422 cm^{-1} ... 0
224.438223 cm^{-1} ... 0 420.612303 cm^{-1} ... 0 621.808843 cm^{-1} ... 0
234.068159 cm^{-1} ... 0 431.989191 cm^{-1} ... 0 628.179834 cm^{-1} ... 0
237.258076 cm^{-1} ... 0 436.832023 cm^{-1} ... 0 635.743849 cm^{-1} ... 0
246.918589 cm^{-1} ... 0 440.332087 cm^{-1} ... 0 670.825354 cm^{-1} ... 0
251.920487 cm^{-1} ... 0 443.474118 cm^{-1} ... 0 688.857654 cm^{-1} ... 0
257.715576 cm^{-1} ... 0 447.569250 cm^{-1} ... 0 697.604660 cm^{-1} ... 0
269.260617 cm^{-1} ... 0 452.681178 cm^{-1} ... 0 720.076755 cm^{-1} ... 0
273.126583 cm^{-1} ... 0 456.422861 cm^{-1} ... 0 720.557545 cm^{-1} ... 0
279.054076 cm^{-1} ... 0 461.810284 cm^{-1} ... 0 722.443945 cm^{-1} ... 0
283.306721 cm^{-1} ... 0 466.924467 cm^{-1} ... 0 750.633956 cm^{-1} ... 0
292.490982 cm^{-1} ... 0 469.385716 cm^{-1} ... 0 766.912826 cm^{-1} ... 0
294.711081 cm^{-1} ... 0 473.217741 cm^{-1} ... 0 769.112977 cm^{-1} ... 0
302.023312 cm^{-1} ... 0 480.833441 cm^{-1} ... 0 772.988666 cm^{-1} ... 0
308.267294 cm^{-1} ... 0 486.574312 cm^{-1} ... 0 787.342636 cm^{-1} ... 0
314.384194 cm^{-1} ... 0 489.989054 cm^{-1} ... 0 841.286620 cm^{-1} ... 0
320.492063 cm^{-1} ... 0 494.573925 cm^{-1} ... 0 883.249955 cm^{-1} ... 0
324.405949 cm^{-1} ... 0 495.401596 cm^{-1} ... 0 1325.276492 cm^{-1} ... 0
329.099517 cm^{-1} ... 0 499.690424 cm^{-1} ... 0 2801.389197 cm^{-1} ... 0

Does the ... 1 after figure stand for imaginary frequencies or minus prior to figure? Could you tell me what situation will produce this error ? due to change the DISPLACECAR by hand? thanks for your reply.

Your initial state should not have negative modes (imaginary frequencies, marked with the 1 tag).

Make sure that your initial state in fully optimized to a minimum. Check that your displacements are small enough that your calculated frequencies are insensitive to the displacement magnitude (typically 0.01-0.001 Ang).

The imaginary frequencies are likely due to the fact that you have only frozen the z direction of some atoms. If the system is able to translate, these modes will show up as small (real or imaginary) frequencies. If you can calculate the frequencies for a subset of atoms that move the most in your transition (leaving the rest fixed), you should be able to find a converged prefactor.

Dear Graeme:
the Initial state have negative modes is my human error indeed, thank you for help. when i do this calculation, i find a very important problem, this is, how to decide which atom i should fix? for any kind of calculation(NEB、dimer、dynamical matrix), is it have the same standard? my model will fix the lowest layer to do adsorption calculation. but when i calculate Hessian matrix(IBRION = 5), i be taught to fix the all surface atom. when i do calculation on dimer, i be taught to free all atom include lowest layer? different standard make me confuse, please told me the standard roughly, thanks.

if i fix the surface atom, the negative modes of my initial state will disappear, and my Transition state will produce only one negative modes. it seems work, Should i still do on dynamical matrix calculation?
i also calculation the dimer method, when i do this, i will free all atom. but the calculations will cost lots of time(NSW greater than 1500) , is it caused by my free atoms?
thanks for your kind help.

For slab calculations, we freeze one or two layers for all calculations (minimization, dimer, dynamical matrix, neb ... ). With frozen atoms, the slab is held in place, an you will avoid the low modes associated with rotation and translation - even when you include all non-frozen atoms in your dynamical matrix calculation.

It is possible that you can get away with a subset of non-frozen atoms for your dynamical matrix calculations, particularly if you have a large system and some non-frozen atoms do not move significantly in the processes. Our dymseldsp.pl and dymselsph.pl scripts can be used to help generate such a subset. Then, like all calculations, you must check for convergence with respect the system size and the subset of atoms in your dynamical matrix.

Dear Graeme:
After taking your advice, I start frozen lowest layers to calculate the dimer method in my system. It still can’t converge after running 2000 step. My content of INCAR:
System = dimer
ISMEAR = -1
SIGMA = 0.1
NSW = 3000
ENCUT = 600
ALGO = Fast
ISIF = 2
GGA = 91
LREAL = Auto
NPAR = 16
NELM = 100
LCHARG = F
LWAVE = F
#dimer
ICHAIN = 2
IBRION = 3
POTIM = 0.0
EDIFF = 1E-7
EDIFFG = -0.01
DdR = 5E-3
DRotMax = 4
DFNMax = 1.0
DFNMin = 0.01
#opt
IOPT = 2
My content of DIMER:
217 0.11210 20.38554 -326.47049 -8.06343 22.68353
217 0.11210 19.09383 -326.47049 -6.28473 18.69901
217 0.11210 13.12400 -326.47049 -7.46599 10.25663
217 0.11210 7.20344 -326.47049 -10.50839 5.79657
218 0.12190 16.18353 -326.47048 -9.58247 19.61230
218 0.12190 12.30862 -326.47048 -9.74127 13.56443
218 0.12190 6.36175 -326.47048 -11.40274 7.97270
218 0.12190 7.51736 -326.47048 -11.55634 9.95456
219 0.11471 18.01268 -326.47046 -8.18103 16.23952
219 0.11471 12.17351 -326.47046 -10.00418 17.34036
219 0.11471 11.07192 -326.47046 -10.08963 13.56821
219 0.11471 12.12466 -326.47046 -9.85816 13.68469
220 0.10188 22.88539 -326.47049 -7.28138 29.29719
220 0.10188 14.43288 -326.47049 -7.66799 17.94874
220 0.10188 11.34474 -326.47049 -10.72440 12.65176
220 0.10188 11.64271 -326.47049 -9.68169 11.14884
221 0.10002 16.53122 -326.47051 -8.94181 17.78788
221 0.10002 13.08230 -326.47051 -9.08386 17.66327
221 0.10002 11.62297 -326.47051 -9.82039 11.28711
221 0.10002 13.53880 -326.47051 -9.88836 14.24793
222 0.09513 14.82995 -326.47049 -8.77070 14.95829
222 0.09513 9.97007 -326.47049 -10.79701 9.23518
222 0.09513 5.91326 -326.47049 -11.60459 7.64933
222 0.09513 6.63784 -326.47049 -11.10644 5.47970
223 0.07647 9.76276 -326.47051 -10.68836 13.38485
223 0.07647 15.31228 -326.47051 -9.72642 9.90551
223 0.07647 8.46597 -326.47051 -10.45832 7.57404
223 0.07647 4.95181 -326.47051 -11.26126 4.17546
224 0.10713 12.55186 -326.47050 -10.49312 16.70653
224 0.10713 11.49414 -326.47050 -9.63934 12.77147
224 0.10713 8.81184 -326.47050 -10.35836 8.78363
224 0.10713 8.33831 -326.47050 -11.08773 11.32726
225 0.12906 10.00408 -326.47047 -9.97633 13.16871
225 0.12906 9.59266 -326.47047 -10.06810 8.15130
225 0.12906 4.03706 -326.47047 -11.22104 3.20346
Is it normal situation? Can my system converge? Please give me some suggestion, thanks.

This is certainly slow convergence, but the numbers at this point look reasonable. As long as the 5th column remains negative, and forces (2nd column) remain small, it is likely near a saddle. Take a look at the XDATCAR (you can use xdat2xyz.pl to see a movie, or using vmd) to make sure that the geometry makes sense.