05 – Análise de dados volumétricos
Alguns dados que podem ser obtidos nas simulações atomísticas com o SIESTA são dados volumétricos, ou seja, uma grandeza dada em pontos do espaço real tridimensional. Entre alguns exemplo de dados volumétricos está a própria densidade eletrônica (
Obtenção de dados volumétricos
Para obter os arquivos de dados volumétricos, vamos fazer um exemplo de um molécula de glicina, como mostrada na Figura 5.1.
Para isso, vamos incluir no input.fdf, as seguintes chaves e valores:
SaveRho true SaveElectrostaticPotential true SaveTotalPotential true
A opção SaveRho faz com que a simulação crie um arquivo *.RHO, SaveElectrostaticPotential cria um arquivo *.VH, e SaveTotalPotential cria um arquivo *.VT. Esses arquivos servem de dados de entrada para um outro programa, o grid2cube, que converte esses arquivos em um arquivo *.cube. Arquivos no formato cube podem ser abertos por programas de visualização de dados volumétricos como o VESTA.
Uma versão do grid2cube é fornecida com o próprio programa SIESTA, no diretório Util/Grid do diretório principal do SIESTA. Neste diretório, execute o comando:
make
para criar um arquivo executável grid2cube. Uma versão modificada desse programa pode ser encontrada no diretório util do repositório do tutorial do SIESTA: https://github.com/SeixasResearch/siesta-tutorials/tree/master/util. Essa versão modificada (grid2cube_diff_sum.f) será útil para os tutoriais que serão apresentados mais para frente, com a inclusão de polarização de spin. A compilação desse código em fortran deve ser feita com um compilador, como o gfortran ou ifort. Com o gfortran, o comando para compilação é:
gfortran -o grid2cube_diff_sum grid2cube_diff_sum.f
Para executar o programa grid2cube, além do arquivo *.RHO, também precisamos de um arquivo de entrada do próprio grid2cube, no formato a seguir:
GlycineAminoAcid rho 0.0 0.0 0.0 1 unformatted
A primeira linha desse arquivo deve conter o SystemLabel do arquivo *.RHO. Como o arquivo do exemplo06 do tutorial é GlycineAminoAcid.RHO, o SystemLabel é GlycineAminoAcid. Na segunda linha desse arquivo, incluímos o formato do arquivo (com todas as letras em minúsculo). Na terceira linha, são informados três números que servem para fazer um deslocamento rígido nos dados volumétricos. Esses valores para o deslocamento rígido são dados na unidade de medida bohr (1 bohr = 0.529 Å). Nas últimos duas linhas, mantenha os dados como no exemplo acima. O programa grid2cube também precisa de um arquivo *.XV que deve ter o mesmo SystemLabel do arquivo *.RHO. Salvando esse pequeno arquivo como um arquivo cube.in, podemos executar o programa grid2cube como:
grid2cube < cube.in
O resultado dessa execução deve ser algo do tipo:
Reading grid data from file GlycineAminoAcid.RHO Cell vectors 37.794537555487111 0.0000000000000000 0.0000000000000000 0.0000000000000000 37.794537555487111 0.0000000000000000 0.0000000000000000 0.0000000000000000 37.794537555487111 Grid mesh: 250 x 250 x 250 nspin = 1 Writing CUBE file GlycineAminoAcid.RHO.cube
Para conversão dos arquivos *.VH e *.VT, precisamos mudar a segunda linha do arquivo cube.in para vh ou vt, e executar o grid2cube novamente. É preciso ter um pouco de cuidado com esses arquivos *.cube, já que eles podem ser muito grandes mesmo para moléculas pequenas. Nesse nosso exemplo, o tamanho dos arquivos *.cube da molécula glicina tem em torno de 196 MB cada arquivo. Para sistemas maiores, esses arquivos podem atingir alguns gigabytes.
Visualização de dados volumétricos
Há vários programas que podem ser utilizados para visualização de arquivos no formato cube. Um dos programas que eu mais recomendo é o VESTA. Esse programa pode ser baixado no site: https://jp-minerals.org/vesta/en/. Ele possui um interface gráfica como mostrada na Figura 5.2 abaixo.
Ao abrir um arquivo *.cube, o programa já vai mostrar automaticamente a geometria da molécula e uma isosuperfície dos dados da densidade eletrônica, como na Figura 5.3.
Modo de isosuperfícies e modo de fatias de volumes
Há dois modo de visualização de dados volumétricos como o da densidade eletrônica. O primeiro modo, de isosuperfície, é mostrado automaticamente pelo VESTA ao abrir o arquivo *.cube. Nesse modo, é mostrado uma superfície em torno da geometria que possui o mesmo valor da grandeza do arquivo *.cube. Esse valor constante ao longo da superfície pode ser ajustado no menu Properties, na aba Isosurfaces. Também podemos ajustar aspectos estéticos das isosuperfícies.
Um segundo modo de visualizar os dados volumétricos é através de mapas de cores em fatias de volumes (volume slices). Neste modo, primeiro é preciso definir um plano que atravessa a geometria da molécula. Isso é feito no menu Edit > Lattice Plane. Depois, o mapa de cor dos valores da densidade eletrônica no plano definido podem ser ajustadas no menu Properties, na aba Sections. Tanto o esquema de cores, como os valores máximos e mínimos para saturação de cores podem ser definidos nesse menu. Um exemplo de volume slice é mostrado na Figura 5.4. Valores maiores da densidade eletrônica são mostradas em vermelho.
Mapa de potencial eletrostático
Há uma forma de visualizar dois dados volumétrico (densidade eletrônica e potencial de Hartree) na mesma representação gráfica chamada de mapa de potencial eletrostático. Nessa representação, há uma isosuperfície para uma das grandezas (densidade eletrônica), e um mapa de cor para a outra grandeza (potencial de Hartree) representada na mesma isosuperfície da densidade eletrônica. Esse tipo de representação gráfica são mais elucidativos do que os modos mencionados anteriormente. No exemplo mostrado abaixo, regiões da molécula em que há um acúmulo de elétrons são representadas em vermelho; enquanto que regiões em que há falta de elétrons são representadas em azul. Dessa forma, podemos visualizar a transferência de elétrons das ligações químicas da molécula pela regiões azuis e vermelhas. Regiões verdes são neutras.
Como pode ser visto na Figura 5.5, os átomos de H vão possui cargas parciais positivas, enquanto que os átomos de oxigênio vão possui cargas parciais negativas. A região negativa em torno do átomo de nitrogênio é devido ao par solitário de elétrons (lone pairs).
Para criar esse mapa de potencial eletrostático, primeiro precisamos abrir o arquivo *.cube da densidade eletrônica no modo de isosuperfícies. A seguir, devemos ir no menu Edit > Edit Data > Volumetric Data. E finalmente, na seção surface coloring, devemos importar o arquivo *.cube dos dados de potential de Hartree. Nos plots mostrados nesse tutorial não foram feitas conversões de unidades dos dados volumétricos do potencial de Hartree.
A molécula de glicina também pode apresentar uma forma de zwitterion. Nessa forma, há uma transferência de próton do grupo da carboxila para o nitrogênio, formado um grupo negativo (-COO–) e outro grupo positivo (-NH3+) na molécula. Relaxando a geometria da molécula zwitterion de glicina, encontramos que a energia total é mais alta do que a primeira geometria mostrada da glicina. A diferença da energia total é 0.88 eV. Criando os arquivos *.cube das densidade eletrônicas e potencial de Hartree do zwitterion, podemos criar um mapa de potencial eletrostático como mostrado na Figura 5.6. Nessa molécula, é bastante evidente a parte positiva da molécula (em azul), e parte negativa (em vermelho). Os dados de entrada para essa simulação do zwitterion também estão no exemplo06 do repositório do tutorial no github.
Análise de cargas parciais
Para estudar as propriedades químicas de algumas moléculas e sólidos, podemos calcular as cargas parciais de cada átomo da geometria a partir de uma partição da densidade eletrônica em pequenas regiões em torno de cada átomo. Essa partição da densidade eletrônica pode ser realizada com o programa Bader Charge Analysis [http://theory.cm.utexas.edu/henkelman/code/bader/]. Para esse programa, há dois algoritmos que podem ser usados na partição, e que podem ser passados como argumentos na execução do programa. O primeiro algoritmo é o de Bader, em que as fronteiras das regiões de um átomo para outro são definidas pelo gradiente da densidade eletrônica. A fronteira é estabelecida nas posições em que o gradiente da densidade é nulo, isto é, fluxo de carga nulo. O segundo algoritmo é o de Voronoi, em que as fronteiras são definidas pela distância média entre posições dos núcleos atômicos.
Para a análise de cargas parciais de Bader do SIESTA, primeiro precisamos incluir a chave e valor:
SaveBader true
no arquivo de entrada (input.fdf). De forma semelhante ao arquivo *.RHO, o programa cria um arquivo *.BADER de dados volumétricos que podem ser convertidos em *.cube. Para utilizar o programa bader, é preciso converter primeiro o arquivo *.BADER em *.cube. Isso é feito com o grid2cube, com o arquivo de entrada cube.in com a segunda linha com dado bader.
A execução do programa bader é feita da seguinte forma:
bader -c bader arquivo.cube
para análise de cargas parciais com o algoritmo de Bader, e na seguinte forma:
bader -c voronoi arquivo.cube
para análise de cargas parciais com o algoritmo de Voronoi. Mais opções de uso do programa bader podem ser vista com as instruções de ajuda obtidas com o comando:
bader -h
Ao executar o programa bader para a análise de cargas parciais, vamos obter um arquivo de saída ACF.dat, com conteúdo na forma:
# X Y Z CHARGE MIN DIST ATOMIC VOL -------------------------------------------------------------------------------- 1 7.120276 13.175146 9.405064 5.685763 0.994361 95.462093 2 9.884840 12.417059 9.754700 4.809201 0.739818 79.613361 3 12.061409 9.647268 10.203180 1.555516 0.394478 11941.738112 4 3.513334 11.759382 9.401606 1.712936 0.399826 8470.112541 5 5.617948 9.785329 8.217584 1.690692 0.397897 7340.087916 6 6.737633 14.599870 10.915695 1.943144 0.516477 6999.628678 7 7.094371 14.264789 7.579310 1.966490 0.561024 6669.030413 8 11.646065 13.925610 9.723651 8.917104 1.196901 9525.754889 9 10.218732 9.863150 10.046859 8.849334 1.122345 1654.739627 10 5.324648 11.094569 9.624450 7.869670 1.119712 1210.411938 -------------------------------------------------------------------------------- VACUUM CHARGE: 0.0000 VACUUM VOLUME: 0.0000 NUMBER OF ELECTRONS: 44.9999
Na quinta coluna desse arquivo, são dadas as cargas totais calculadas nas regiões definidas pelo algoritmo de Bader. É preciso fazer duas observações com esses valores de cargas. A primeira observação é que essas são cargas totais, e precisam ser comparadas com as cargas totais dos átomos isolados (ou números atômicos de cada átomo). A segunda observação é específica para o átomo de H. Para esse átomo, é preciso descontar não apenas 1 carga positiva, mas 2 cargas positivas. A análise de cargas parciais de Bader adiciona uma carga positiva para cada átomo de H da molécula. Note que para a molécula de glicina, a soma dos números atômicos da molécula é 40 (2*6+5*1+2*8+7), mas o número de elétrons totais, mostrada na última linha do arquivo ACF.dat é 45.
Para obter a carga parcial de um átomo, descontamos o número atômico da carga mostrada no arquivo ACF.dat. Por exemplo, para o N, o valor mostrado é de Q = 7.87, e o número atômico do N é 7. Assim, há 0.87 elétrons em excesso na região definida em torno do átomo de N pelo algoritmo de Bader. A carga elétrica será q = -0.87e causa do sinal negativo da carga do elétron.
MackSim 