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
_______________________________________________
Quantum ESPRESSO is supported by MaX (www.max-centre.eu/quantum-espresso)
users mailing list [email protected]
https://lists.quantum-espresso.org/mailman/listinfo/users
