Tutorial QE: Sistemas moleculares (0D) - Parte 2

---

Nesse tutorial aprenderemos:

• Como calcular diretamente a energia de formação diretamente;
• Como calcular por DFPT os modos vibracionais diretamente;

---

Calculando a energia de formação

O átomo de nitrogênio possui configuração eletrônica \([1s^2]2s^22p^3\), e em seu estado fundamental, segundo a regra de Hund, apresenta spin S=3/2: \(2s\) \(\downarrow \uparrow\) \(2p\) \(\uparrow\) \(\uparrow\) \(\uparrow\) tendo assim termo espectroscópico \(^4S_{3/4}\) Nos cálculos que fizemos até agora, não levamos isso em consideração. Para pequenas variações da distância de equilíbrio não há uma grande influência pois a ligação está formada e os spins estão acoplados na formação da ligação. Entretanto, se quisermos calcular a energia de formação da molécula com precisão precisamos fazer um cálculo spin-polarizado.

As regras de Hund, conjunto de regras formuladas em 1927 pelo físico alemão Friedrich Hund, é uma forma de determinar o termo espectroscópico (\(^{2S+1}L_J\)) correspondente ao estado fundamental de um átomo polieletrônico. São elas:

• Para uma dada configuração eletrônica, o termo com máxima multiplicidade apresenta a menor energia. A multiplicidade é definida como \(2S+1\) onde \(S\) é o momento angular de spin total;
• Para uma dada multiplicidade, o termo com maior valor do número quântico de momento angular orbital, \(L\) apresenta menor energia;
• Para um determinado termo, em um átomo com sub-camada mais externa com metade do preenchimento ou menos, o nível com o valor mais baixo do número quântico do momento angular total \(J=L + S\), apresenta a energia mais baixa. Se a camada mais externa estiver mais da metade preenchida, o nível com o valor mais alto de \(J\) é o mais baixo em energia.

Para esse tipo de cálculo vamos modificar o script como o apresentado à seguir, adicionando as variáveis nspin=2 para informar que o cálculo é spin-polarizado, tot_magnetization = 3.0 para um spin total de 3 magnetons de Bohr.

#!/bin/bash

PREFIX="N"
PSEUDO_DIR="/home/pseudo"
CALC="scf"
CMD_PW="mpirun -np 4 /home/interlab/qe/bin/pw.x"

#Make SCF calculation
cat>$PREFIX.$CALC.in<<EOF
 &control
  calculation = '$CALC' ,
  restart_mode = 'from_scratch' ,
  wf_collect = .false. ,
  outdir = '.' ,
  pseudo_dir = '$PSEUDO_DIR' ,
  prefix = '$PREFIX' ,
/
&system
  ibrav = 1,  celldm(1) = 20
  nat =  1 ,  ntyp = 1 ,
  ecutwfc = 80.0 ,  ecutrho = 600.0 ,
  nspin = 2, tot_magnetization = 3.0, occupations = 'smearing',
  degauss = 0.001,
/
&electrons
  conv_thr =  1.0D-10 , electron_maxstep = 5000 ,
/
ATOMIC_SPECIES
  N   14.0107      N.pbesol-n-rrkjus_psl.1.0.0.UPF
ATOMIC_POSITIONS (angstrom)
  N        0.000000000   0.000000000   0.000000000
K_POINTS automatic
  1 1 1   0 0 0
EOF
$CMD_PW < $PREFIX.$CALC.in > $PREFIX.$CALC.out

Ao executarmos esse cálculo veremos que a energia total calculada para o átomo de Nitrogênio no estado fundamental \(^4S_{3/4}\) é de -20.08714260 Ry.

O valor absoluto das energias calculadas (autovalores das equações de Kohn-Sham) não possuem nenhum significado físico. Entretanto, diferenças de energias, sempre calculadas no mesmo nível de cálculo, podem ser comparadas com valores experimentais. Por esse motivo, sempre que vamos reportar dados subtraímos de todos os valores o \(U_0\) obtido do ajuste e renomeamos o eixo y para Energia de formação.

A energia de formação \(E_f\) de uma molécula homonuclear constituída por \(n\) átomos pode ser calculada como \(E_f = n \cdot E_{\text{átomo}} - E_{\text{molécula}}\)

1. Com base na distância de equilíbrio obtida anteriormente faça o cálculo para a molécula de \(N_2\) em seu estado de equilíbrio e utilize a equação acima para obter a energia de dissociação.
2. O valor experimental da energia de dissociação é de \(225.07 \pm 0.11 kcal/mol\). [2] Qual a discrepância relativa entre o valor calculado e o valor experimental?
3. Compare o valor obtido para o átomo isolado com o valor de \(U_0\) obtido do ajuste do potencial de Morse. Qual é a possível origem dessa diferença?

Calculando os modos vibracionais

Por fim, para confirmarmos que estamos mesmo com uma estrutura no estado de equilíbrio, precisamos calcular os modos vibracionais para a molécula. Se todos os modos vibracionais apresentarem frequências reais, i.e., autovalores positivos podemos assumir que a estrutura está no estado de equilíbrio. Para isso precisamos primeiro fazer um cálculo SCF com a molécula apresentando o comprimento de ligação obtido pelo ajuste do potencial de Morse. Em seguida precisamos fazer o cálculo da matriz dinâmica para a molécula, utilizando o programa ph.x, com o input:

#!/bin/bash

CMD_PH="mpirun -np 4 /home/interlab/qe-6.3/bin/ph.x"
CMD_DYNMAT="mpirun -np 4 /home/interlab/qe-6.3/bin/dynmat.x"
PREFIX="n2_eq"
PSEUDO_DIR="/home/interlab/pseudo/"
CALC='phG'
cat > $PREFIX.$CALC.in << EOF
$PREFIX - Phonons at Gamma
&inputph
  tr2_ph=1.0d-14,  prefix='$PREFIX',  epsil=.true.,
  amass(1)=14.0067,  outdir='.',  fildyn='$PREFIX.dynG',
  asr=.true. ,
/
0.0 0.0 0.0
EOF
$CMD_PH < $PREFIX.$CALC.in > $PREFIX.$CALC.out
#Make DynMat Acoustic Sum Rules calculation
cat>$PREFIX.dynmat.in<<EOF
$PREFIX - Acoustic Sum Rules - DynMat
&input
  fildyn='$PREFIX.dynG',  asr='zero-dim' ,
/
EOF
$CMD_DYNMAT < $PREFIX.dynmat.in > $PREFIX.dynmat.out

O programa ph.x utiliza on density-functional perturbation theory (DFPT) [3] para calcular as frequências vibracionais diretamente da estrutura em equilíbrio, usando teoria de perturbação.

Se olharmos para o arquivo n2_eq.phG.out veremos que existem algumas frequências negativas. Esse valores negativos são apenas uma convenção que é utilizada para indicar que a diagonalização da matriz dinâmica (veremos mais sobre isso futuramente) apresentou autovalores \(\omega^2 < 0\). Nesse caso o código apresenta a quantidade \(-\sqrt{\|\omega^2\|}\), sendo o sinal de menos somente um aviso de que algo está errado.

Como estamos lidando com uma molécula isolada esperamos obter \(\omega = 0\) para os 5 primeiros modos (3 translações e duas rotações), e um modo de alta frequência referente ao estiramento no eixo da molécula. O fato de a molécula estar em uma supercéula periódica faz com que essas rotações e translações envolvam uma pequena quantidade de energia fictícia, levando à um artefato de cálculo e ao aparecimento dessas frequências negativas.

Podemos corrigir esse artefato de cálculo impondo a chamada regra de soma acústica. Esse procedimento corresponde a modificar a matriz dinâmica de forma a descontar esses movimentos de rotação e translação. Esse cálculo de pós-processamento é feito pelo programa dynmat.x. Para olharmos o resultado basta abrir o arquivo n2_eq.dynmat.out. Nele vemos que os 5 primeiros modos agora apresentam frequência zero e a frequência de oscilação do estiramento no eixo da molécula é de \(\omega_0 = 2338.16 cm^{-1}\).

1. O valor experimental para o estado fundamental(\(X^1\Sigma_g^+\)) da molécula de \(N_2\) é \(\omega_0 = 2329.92 \pm 3 cm^{-1}\). [4] Calcule a discrepância relativa entre o valor calculado por DFPT e o experimental.
   • Use \(c = 2.998\cdot 10^{10} cm.s^{-1}\) , \(1 Ry = 2.179872 \cdot 10^{-18} kg.m^2.s^{-2}\), \(1 u.m.a = 1.660538\cdot 10^{-27} kg\), \(h = 6.62607\cdot 10^{-34} m^2.kg.s^{-1}\)
2. Calcule a ZPE com base na frequência vibracional obtida por DFPT. Compare com o valor obtido na questão anterior e com o valor “experimental” \(1175.78 \pm 5 cm^{-1}\) [5].

[1] BOWEN, Humphrey John Moule. Tables of interatonic distances and configuration in molecules and ions. 1958.

[2] VEDENEEV, Vladimir Ivanovich; VEDENEEV, Vladimir Ivanovich. Bond energies, ionization potentials, and electron affinities. London: E. Arnold, 1966.

[3] BARONI, Stefano et al. Phonons and related crystal properties from density-functional perturbation theory. Reviews of Modern Physics, v. 73, n. 2, p. 515, 2001.

[4] BENDTSEN, Jørgen. The rotational and rotation‐vibrational Raman spectra of \(^{14}N_2\), \(^{14}N^{15}N\) and \(^{15}N_2\). Journal of Raman Spectroscopy, v. 2, n. 2, p. 133-145, 1974.

[5] IRIKURA, Karl K. Experimental vibrational zero-point energies: Diatomic molecules. Journal of physical and chemical reference data, v. 36, n. 2, p. 389-397, 2007.

Respostas:

• Calculando a energia de formação 1. 2. 3.

• Calculando os modos vibracionais 1. 2.