[QE-users] Questions related to vibrational analysis (for both molecule and adsorbate)
Shen, Ziheng
zshen83 at gatech.edu
Tue Mar 3 04:01:40 CET 2020
Dear QE users,
Greetings! I’m a beginner of this program and I expected to use QE to find kinetic information of reactions in adsorption processes. Specifically, I’m trying to reproduce the PES of a published case
CH* + * -> C* + H*
where * refers to active site on the Ni(111) surface
So far, I’m able to find the initial state and final state structures. I also got one transition state structure by neb.x. Now I wanted to perform vibrational analysis to confirm the T.S. strucuture (expecting to see one negative frequency). However, when I did tests on ph.x, I met some troubles. These problems has been asked before in the forum, but I still cannot figure them out by following those given solutions. I appreciate quite a lot if anyone could provide some hints or advice for my specific case.
First, I did a simple test of ph.x by calculating the frequencies of a single methane molecule. From my previous experience with ORCA, the first 6 frequencies (which represent translational and rotational modes) should be zero for a molecule. However, the program always returned non-zero values. I followed some previously posted solutions and tried approaches including:
1) try more tight forc_thr in geometry relaxation => I decreased forc_thr from default to 1D-5
2) try more tight conv_thr in scf calculation => I decreased conv_thr from 1D-8 to 1D-12
3) try bigger ecutwfc/ecutrho in scf calculation => I increased ecutwfc from 65Ry to 85Ry
4) try to increase the size of unit cell => I increased the length of cubic cell from 20A to 30A
5) try more tight tr2_ph and smaller alpha_mix in ph.x => I decreased tr2_ph from default to 1D-16
However, none of them worked. The input files with highest thresholds for pw.x & ph.x are copied here.
=======================================scf input=======================================
&CONTROL
Calculation='scf',
restart_mode='from_scratch',
prefix = "ch4_cell2"
outdir = "./ts/tmp",
pseudo_dir = "./pseudo",
tstress = .true.
verbosity = 'high'
etot_conv_thr = 1.0D-5
forc_conv_thr = 1.0D-5
/
&SYSTEM
ibrav = 0,
nat = 5,
ntyp = 2,
ecutwfc = 85,
ecutrho = 850,
/
&ELECTRONS
conv_thr = 1.D-12,
mixing_beta = 0.1,
electron_maxstep=250,
/
&IONS
ion_dynamics='bfgs'
/
ATOMIC_SPECIES
C 12 C.pbe-n-kjpaw_psl.1.0.0.UPF
H 1 H.pbe-kjpaw_psl.1.0.0.UPF
CELL_PARAMETERS { angstrom }
30 0.0000000000 0.0000000000
0.0 30 0.0000000000
0.0000000000 0.0000000000 30
ATOMIC_POSITIONS { angstrom }
C 12.998505897 12.669983024 10.000014581
H 14.094202415 12.669376542 10.000895687
H 12.632663345 13.058381752 10.957013019
H 12.632991209 11.646964614 9.857276033
H 12.634287134 13.305064066 9.184810679
K_POINTS { automatic }
1 1 1 0 0 0
=======================================phonon input=======================================
phonons of CH4
&inputph
tr2_ph=1.0d-16,
prefix='ch4_cell2',
amass(1)=12.011,
amass(2)=1.0,
alpha_mix(1)=0.1,
outdir='./ts/tmp/',
fildyn='CH4_cell2.dynG'
/
0.0 0.0 0.0
The calculated frequencies are:
=======================================phonon output=======================================
freq ( 1) = 1.772990 [THz] = 59.140581 [cm-1]
freq ( 2) = 4.134960 [THz] = 137.927411 [cm-1]
freq ( 3) = 4.480988 [THz] = 149.469657 [cm-1]
freq ( 4) = 10.709853 [THz] = 357.242257 [cm-1]
freq ( 5) = 11.230300 [THz] = 374.602485 [cm-1]
freq ( 6) = 12.538150 [THz] = 418.227663 [cm-1]
freq ( 7) = 38.905159 [THz] = 1297.736407 [cm-1]
freq ( 8) = 39.196874 [THz] = 1307.466991 [cm-1]
freq ( 9) = 39.785039 [THz] = 1327.086048 [cm-1]
freq ( 10) = 46.262462 [THz] = 1543.149645 [cm-1]
freq ( 11) = 46.466078 [THz] = 1549.941531 [cm-1]
freq ( 12) = 89.472569 [THz] = 2984.483660 [cm-1]
freq ( 13) = 92.974263 [THz] = 3101.287588 [cm-1]
freq ( 14) = 93.136064 [THz] = 3106.684690 [cm-1]
freq ( 15) = 93.343894 [THz] = 3113.617156 [cm-1]
Following the advice from FAQ (“if the nonzero frequencies are small, you can impose the ASR to the dynamical matrix, usually with excellent results”) in QE official website, I also tried to apply asr=true but the calculation didn’t converge. I also saw some users claiming that if the translational and rotational frequencies are small, it could be due to numerical error and thus they can be negligible. However, I don’t think a frequency that’s around 300 cm-1 can be considered “small”. Could someone please take a look at my input files and point out some possible errors for me? And in which cases can we say that the frequency is “small”?
Second, my ultimate goal is to perform frequency analysis for adsorbate so that I can both determine transition state structures and apply ZPE corrections. In the maillist archive, some users suggested set “nat_todo” to only calculate adsorbates. However, I also got some crazy results by doing so.
=========================scf input, structure obtained from neb.x========================
&CONTROL
Calculation='scf',
restart_mode='from_scratch',
prefix = "Ni_ch_ts"
outdir = "./ts/tmp",
pseudo_dir = "./pseudo",
tstress = .true.
verbosity = 'high'
tefield = .true.
dipfield = .true.
/
&SYSTEM
ibrav = 0,
nat = 14,
ntyp = 3,
ecutwfc = 65,
ecutrho = 650,
Occupations='smearing',
smearing='mp',
degauss=0.01,
nspin=2,
starting_magnetization(1)=0.2,
eamp = 0.0
edir = 3
emaxpos = 0.95
eopreg = 0.05
/
&ELECTRONS
electron_maxstep=250,
conv_thr = 1.D-10,
mixing_beta = 0.1,
/
ATOMIC_SPECIES
Ni 58.69 ni_pbe_v1.4.uspp.F.UPF
C 12 C.pbe-n-kjpaw_psl.1.0.0.UPF
H 1 H.pbe-kjpaw_psl.1.0.0.UPF
CELL_PARAMETERS { angstrom }
4.9667177200 0.0000000000 0.0000000000
2.4833588600 4.3013037190 0.0000000000
0.0000000000 0.0000000000 20.000000000
ATOMIC_POSITIONS { angstrom }
Ni 0.0000000000 0.0000000000 7.9723500000
Ni 1.2416800000 2.1506500000 7.9723500000
Ni 2.4833600000 0.0000000000 7.9723500000
Ni 3.7250400000 2.1506500000 7.9723500000
Ni 2.4833600000 1.4337700000 10.0000000000
Ni 3.7250400000 3.5844200000 10.0000000000
Ni 4.9667200000 1.4337700000 10.0000000000
Ni 6.2084000000 3.5844200000 10.0000000000
Ni 1.2281125212 0.7413038538 12.1124602290
Ni 2.4833613443 2.9495126082 12.0722664849
Ni 3.7386101535 0.7413040204 12.1124604793
Ni 4.9667205066 2.8946406480 11.9560276983
C 2.4833610952 1.4972020489 13.1166864267
H 2.4833559794 3.1940064219 13.5465576823
K_POINTS { automatic }
6 6 1 0 0 0
=======================================ph.x input=======================================
phonons of CH on metal Ni at Gamma
&inputph
tr2_ph=1.0d-12,
prefix='Ni_ch_ts',
epsil=.false.,
amass(1)=58.69,
amass(2)=12.011,
amass(3)=1.0,
alpha_mix(1)=0.1,
outdir='./ts/tmp/',
fildyn='CH_2.dynG',
nat_todo=2,
/
0.0 0.0 0.0
13 14
==================================phonon output=======================================
freq ( 1 - 1) = -1514178.1 [cm-1] --> A I+R
freq ( 2 - 2) = -204586.6 [cm-1] --> A I+R
freq ( 3 - 3) = -8039.2 [cm-1] --> A I+R
freq ( 4 - 4) = -3220.6 [cm-1] --> A I+R
freq ( 5 - 5) = -2283.0 [cm-1] --> A I+R
freq ( 6 - 6) = -205.2 [cm-1] --> A I+R
freq ( 7 - 7) = -0.0 [cm-1] --> A I+R
freq ( 37 - 37) = 1712.6 [cm-1] --> A I+R
freq ( 38 - 38) = 2144.7 [cm-1] --> A I+R
freq ( 39 - 39) = 3316.3 [cm-1] --> A I+R
freq ( 40 - 40) = 7979.5 [cm-1] --> A I+R
freq ( 41 - 41) = 186893.8 [cm-1] --> A I+R
freq ( 42 - 42) = 1489603.3 [cm-1] --> A I+R
I know that tr2_ph here might not be strict enough. However, it still doesn’t make sense to me that the first 6 frequencies have such large values and are negative. Does that imply this structure is most likely not a transition state? In the manual it says “this is an approximation and may not work”, and I also saw some experienced users in the forum claimed that using “nat_todo” might have some risk. So is it possible that this function just doesn’t work for my system? How reliable is this function?
Thank you for reading this long question! Thanks in advance for anyone that could give suggestions to me!
Best regards
Ziheng Shen
PhD student @ Georgia Institute of Technology
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20200303/b0aa5533/attachment.html>
More information about the users
mailing list