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 (n(\mathbf{r})), a energia potencial de Hartree (V_H(\mathbf{r})), e a energia potencial total (V_T(\mathbf{r})). Neste tutorial vamos aprender a como fazer simulações com o SIESTA que escrevem arquivos que podem ser representados no espaço real, e qual programa podemos para visualizar isso. Também vamos poder utilizar o arquivo da densidade eletrônica para fazer uma análise de cargas parciais de cada átomo da nossa geometria molecular.

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.

glycine
Figura 5.1. Representação esquemática de uma molécula de glicina.

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.

vesta_gui
Figura 5.2. Interface gráfica do programa VESTA.

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.

GlycineAminoAcid_rho_isosurface
Figura 5.3. Representação da densidade eletrônica da glicina por isosuperfícies (isosurfaces).

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.

GlycineAminoAcid_rho_volume_slice
Figura 5.4. Representação da densidade eletrônica da glicina por mapa de cor em fatias de volume (volume slice).

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.

GlycineAminoAcid_electrostaticmap
Figura 5.5. Mapa de potencial eletrostático da molécula de glicina.

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.

Glycine_zwitterion
Figura 5.6. Mapa de potencial eletrostático do zwitterion de glicina.

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.

Para a análise de cargas parciais com o algoritmo de Voronoi, execute o comando:

$ bader -c voronoi arquivo.cube

Nesse caso, os valores das cargas totais dentro de cada região da partição são mostradas na tela do terminal, no fim da execução do programa. Os dados de saída serão na forma:

  VORONOI ANALYSIS RESULT
    #         X           Y           Z        CHARGE      ATOMIC VOL
  ----------------------------------------------------------------------
    1      7.1203     13.1751      9.4051      4.9757        25.5680
    2      9.8848     12.4171      9.7547      5.5135        75.6987
    3     12.0614      9.6473     10.2032      2.3579     12228.8258
    4      3.5133     11.7594      9.4016      2.4258      8627.6532
    5      5.6179      9.7853      8.2176      2.4330      7487.5844
    6      6.7376     14.5999     10.9157      2.4425      7085.9796
    7      7.0944     14.2648      7.5793      2.4592      6732.9367
    8     11.6461     13.9256      9.7237      8.3980      9420.0241
    9     10.2187      9.8631     10.0469      7.7727      1397.8741
   10      5.3246     11.0946      9.6244      6.2215       904.4350
  -----------------------------------------------------------------------
           NUMBER OF ELECTRONS:       44.99985

Enquanto isso, o arquivo de saída ACF.dat continua mostrando apenas os valores calculados pelo algoritmo de Bader, não pelo algoritmo de Voronoi.

Anterior: 04 – Energia de ligação

Próximo: 06 – Simulando sólidos