02 – Primeiras simulações com o SIESTA

Exemplo 01 (Água)

Para fazer as simulações atomísticas com o SIESTA, é preciso dos arquivos de pseudopotenciais para cada elemento do material da simulação e do arquivo input.fdf. Alguns exemplos de arquivos input.fdf podem ser encontrados no nosso repositório dos tutoriais do siesta. Copie esse tutorial com os comandos:

cd ~/Simulations/materials
git clone https://github.com/macksim-mackenzie/siesta-tutorials.git

Entre no primeiro exemplo (example01) e execute a simulação com o run_siesta:

cd siesta-tutorials/example01
run_siesta

Se tudo funcionar bem, você verá uma mensagem “Job completed” na última linha do arquivo de saída (job.out). O comando “run_siesta” executa o programa siesta lendo o arquivo de entrada (input.fdf) e escreve a saída no arquivo “job.out“.

Entendendo o seu arquivo de input

O arquivo de entrada (input.fdf) é um arquivo de texto simples com algumas chaves e valores. Para simplificar o nosso arquivo de input, nós vamos agrupar algumas chaves em alguns “blocos”. Para uma melhor compreensão, esses blocos são separados por comentários. Comentários são textos no arquivo que seguem o caractere hash (#). Os comentários são ignorados pelo programa SIESTA, e são introduzidos no arquivo apenas para facilitar o entendimento do arquivo de entrada. Por exemplo, no nosso arquivo de entrada do EXEMPLO 01, nós usamos chaves e valores agrupadas em dois blocos: “Self-consistent field” e “Material“.

O primeiro bloco é dado por:

# in file: input.fdf

         # -- SELF-CONSISTENT FIELD --

MeshCutoff                   300 Ry
PAO.EnergyShift              30 meV
XC.Functional                GGA
XC.Authors                   PBE
MaxSCFIterations             1000

A chave MeshCutoff é usada para controlar a grade do espaço real para a densidade de carga. O valor de energia é um corte para a energia cinética quando calculada usando ondas planas. A partir da versão 4.1, o valor padrão é 300 Ry. Aumente esse valor para ter uma grade mais fina, mas lembre que isso aumenta o custo computacional da simulação consideravelmente. Nós vamos optimizar esse MeshCutoff em uma seção seguinte do tutorial;
A chave PAO.EnergyShift é usada para controlar o alcance das funções bases localizadas. O programa SIESTA é baseado em combinações lineares de orbitais atômicos (LCAO). Entretanto, esses orbitais atômicos são na verdade orbitais numéricos confinado em uma região espacial. Esse confinamento espacial desloca os níveis de energia como em um poço quadrado infinito (esférico). Valores menores do PAO.EnergyShift indicam funções bases com alcances mais longos. Diferentes funções bases serão utilizadas em na seção de optimização de simulação.
XC.Functional e XC.Authors são duas chaves para escolher o funcional de troca-correlação (XC) que serão utilizados nas equações de Kohn-Sham. O funcional mais popular para simulação de sólidos é o da aproximação do gradiente generalizado (GGA) parametrizado por Perdew, Burke e Ernzerhof (PBE) [Phys. Rev. Lett. 77, 3865 (2016)];
A chave MaxSCFIterations é apenas o número máximo de iterações que poderão ser utilizadas no ciclo de auto-consistência de Kohn-Sham. O valor padrão do SIESTA é 50, mas geralmente isso é um valor muito baixo. Vamos aumentar para um valor mais alto, como 300 ou até mesmo 1000 em nossas simulações.

O segundo bloco do nosso arquivo de entrada é dado por:

# file: input.fdf
 
        # -- MATERIAL --

SystemLabel                  H2O
WriteCoorStep                true
WriteCoorXmol                true

%block ChemicalSpeciesLabel
 1   8  O
 2   1  H
%endblock ChemicalSpeciesLabel

AtomicCoordinatesFormat      Ang

LatticeConstant              1.000 Ang

%block LatticeVectors
       20.0000000000         0.0000000000         0.0000000000
        0.0000000000        20.0000000000         0.0000000000
        0.0000000000         0.0000000000        20.0000000000
%endblock LatticeVectors

%block AtomicCoordinatesAndAtomicSpecies
   10.05053122   10.05053123    9.99999999   1
   11.01902700    9.93038215   10.00000000   2
    9.93038215   11.01902699   10.00000000   2
%endblock AtomicCoordinatesAndAtomicSpecies

SystemLabel é a chave que define o nome de alguns arquivos de saída do programa. Use algo simples, sem espaços ou caracteres especiais;
WriteCoorStep e WriteCoorXmol são duas chaves para escrever as coordenadas atômicas no arquivo principal de saída (WriteCoorStep) e em um arquivo xyz (WriteCoorXmol);
ChemicalSpeciesLabel é uma chave em forma de bloco. Cada linha desse bloco recebe um índice para uma espécie química, um número atômico, e um nome de um arquivo de pseudopotencial (.psf) para essa espécie química;
AtomicCoordinatesFormat é a chave que define o formato das coordenadas atômicas. Usualmente é utilizamos Ang or ScaledByLatticeVectors. No formato Ang, as coordenadas atômicas são valores das coordenadas em ångström (Å). No formato ScaledByLatticeVectors, cada valor de coordenada multiplica um vetor da rede de Bravais para obter o valor real da coordenada. O formato Ang é mais utilizado para moléculas, enquanto que o formato ScaledByLatticeVectors é mais apropriado para cristais;
LatticeConstant é uma chave que multiplica os vetores da rede de Bravais;
LatticeVectors são chaves que definem os vetores da célula unitária que se repetem no espaço real. Para moléculas, esses vetores da rede são apenas vetores com normas muito grandes para evitar interações com as células vizinhas;
AtomicCoordinatesAndAtomicSpecies é a chave no formato de blocos que definem as posições dos átomos e a espécie química dentro da célula unitária. Cada linha no bloco é um átomo na célula unitária.

Mais informações de chaves e valores que podem ser utilizados nos arquivos de entrada do SIESTA são encontradas no manual: SIESTA manual.

Analisando seu arquivo de output

Para mostrar todo o conteúdo do arquivo de saída, use o comando cat:

cat job.out

O arquivo de saída é usualmente muito longo, com muita informação sobre a simulação. Por enquanto, vamos nos focar na parte sobre ciclos de auto-consistência (scf). Para cada iteração do ciclo, a energia de Kohn-Sham, a energia livre, e o nível de Fermi são mostrados. Outras informações também são mostradas para cada iteração (passo) do ciclo, como a diferença na matriz densidade (dDmax) com relação ao passo anterior, e a diferença no Hamiltoniano (dHmax). Essas diferenças (dDmax e dHmax) são usadas para controlar a convergência dos ciclos de auto-consistência.

        iscf     Eharris(eV)        E_KS(eV)     FreeEng(eV)     dDmax    Ef(eV) dHmax(eV)
   scf:    1     -474.831923     -476.892759     -476.892759  0.902382 -2.734301 11.509096
timer: Routine,Calls,Time,% = IterSCF        1       3.753  53.65
   scf:    2     -480.167452     -478.826398     -478.826398  0.176264 -7.267793  4.413552
   scf:    3     -479.128249     -479.090327     -479.090327  0.075113 -5.951575  0.356110
   scf:    4     -479.100644     -479.096260     -479.096260  0.014051 -5.916251  0.292365
   scf:    5     -479.097704     -479.097035     -479.097035  0.004317 -5.896236  0.268409
   scf:    6     -479.098074     -479.097594     -479.097594  0.006155 -5.839471  0.213176
   scf:    7     -479.097937     -479.097845     -479.097845  0.017635 -5.612172  0.050964
   scf:    8     -479.098079     -479.097983     -479.097983  0.003525 -5.638938  0.018375
   scf:    9     -479.097998     -479.097994     -479.097994  0.001679 -5.655289  0.002709
   scf:   10     -479.097994     -479.097994     -479.097994  0.000099 -5.655716  0.000802

Note que o ciclo scf é considerado convergido se dDmax < 1.0E-4 e dHmax < 1.0E-3 eV. Esses critérios de convergências podem ser mudados no arquivo de entrada principal. No nosso exemplo, o ciclo scf converge em 10 iterações.

Outra parte importante é a decomposição da energia final, como mostrada abaixo:

siesta: Final energy (eV):
siesta:  Band Struct. =    -109.265444
siesta:       Kinetic =     329.160653
siesta:       Hartree =     512.246332
siesta:       Eldau   =       0.000000
siesta:       Eso     =       0.000000
siesta:    Ext. field =       0.000000
siesta:   Exch.-corr. =    -128.943101
siesta:  Ion-electron =   -1312.809437
siesta:       Ion-ion =     121.247559
siesta:       Ekinion =       0.000000
siesta:         Total =    -479.097994
siesta:         Fermi =      -5.655716

Essa parte inclui a energia final total. Essa energia é de grande importância para calcular várias propriedades físicas e químicas do sistema. Aqui, a energia é -479.097994 eV. Entretanto, nós devemos lembrar que esse valor não é absoluto, já que estamos usando uma abordagem de pseudopotenciais. Apenas diferenças de energias totais podem ser comparadas com os dados experimentais.

Uma forma rápida de ver a energia total de um arquivo de saída é usando o comando siesta_getenergy, criado pelo script functions.sh da instalação:

siesta_getenergy

Anterior: 01 – Instalando o SIESTA	—	Próximo: 03 – Relaxando geometrias