03 – Relaxando geometrias

Em nosso primeiro exemplo, as posições dos átomos na molécula, comumente conhecida como geometria da molécula, foram já dadas em uma configuração de mínimo local de energia. Entretanto, geralmente não é o caso nos nossos dados de entradas. Geralmente é dado uma aproximação inicial da geometria da molécula, mas com distâncias e ângulos não precisos. Nessas situações, nós podemos usar simulações atomísticas para encontrar um mínimo de energia local. Esse processo, variando distâncias e ângulos das geometrias são chamados de relaxação da geometria, ou optimização da geometria. Conhecendo a energia total do sistema, nós podemos fazer um procedimento simples (e não muito inteligente) de variação manual das distâncias e ângulos entres os átomos, e calcular as energias totais em função dessas variáveis (distâncias e ângulos). Para sistemas um pouco mais complexos, esse procedimento torna-se impraticável e será preciso utilizar outro procedimento para relaxar a geometria de uma molécula.

EXEMPLO 02 (Distância H-H)

Para uma molécula de H2, há apenas uma única distância que nós podemos variar para calcular a energia total. Vamos criar várias simulações atomísticas com diferentes distância entre os átomos de hidrogênio. Para criar essa variação, vamos ajustar a posição de um dos átomos em 0.0 0.0 0.0, e o outro átomo na posição 0.0 0.0 VAR_DISTANCE, onde essa palavra-chave será substituída por valores entre 0.50 e 1.2 Å, com passos de 0.05 Å. Isso significa 15 simulações diferentes. Nós podemos criar um pequeno programa em shell (script) para fazer isso. Nós poderíamos usar outra linguagem de programação, mas o programa é tão simples que vale a pena fazer em shell mesmo. Esse programa deve criar um loop com valores de 0.50 a 1.20, criar uma pasta para cada valor de distância, copiar os arquivos de pseudopotenciais e criar um arquivo de entrada input.fdf em que a palavra-chave VAR_DISTANCE é substituída por um novo valor para a distância entre hidrogênios. Essa última tarefa é feita com o comando “sed” do linux, como mostrado abaixo.

# file: create_inputs.sh

for dist in $(seq -w 0.50 0.05 1.201)
  do
  mkdir dist_$dist
  cd dist_$dist
  cp ../H.gga.psf .
  sed "s/VAR_DISTANCE/$dist/g" ../input.fdf > input.fdf
  cd ..
done

Note que usamos o valor 1.201 para finalizar o loop, ao invés de 1.2. Isso foi feito aqui apenas para incluir o valor 1.2. Por causa da precisão das variáveis tipo float, há alguns casos em que precisamos fazer isso para incluir o último valor. Isso não é necessário se o loop fosse de valores inteiros.

Após salvar o programa, nós podemos executá-lo com o comando:

$ sh create_inputs.sh

Esse programa vai criar 15 pastas com os arquivos de entrada input.fdf e os pseudopotenciais. O programa em shell também pode ser executado com o interpretador bash, ao invés do sh. Ambos devem ter o mesmo resultado.

Para realizar as simulações, vamos criar um outro programa muito simples. Um exemplo de programa para realizar com o loop variando as distâncias é mostrado abaixo:

# file: run.sh

for dist in $(seq -w 0.50 0.05 1.201)
  do
  cd dist_$dist
  siesta-omp input.fdf > job.out
  cd ..
done

A execução de último script em shell é feita pelo comando:

$ sh run.sh &

Após a finalização de todas as 15 simulações, nós podemos analisar os dados para encontrar a distância que minimiza a energia total. Para fazer essa análise, vamos criar um script de shell apenas para extrair os dados de energias totais e salvar em um arquivo. O programa get_energy.sh é mostrado abaixo:

# file: get_energy.sh

for dist in $(seq -w 0.50 0.05 1.201)
  do
  cd dist_$dist
  energy=`grep 'Total =' job.out | cut -d = -f 2`
  echo $dist $energy
  cd ..
done

Após salvar o programa, podemos executá-lo com:

$ sh get_energy.sh > energyXdistance.dat

Agora o arquivo energyXdistance.dat contém as energias totais para distâncias entre 0.50 e 1.20 Å. Plotando esses dados com um programa de gráficos (xmgrace, gnuplot, pyplot, …), nós vamos encontrar o resultado mostrado na Figura 3.1. O mínimo de energia é em torno de d = 0.8 Å.

HH-distance
Figura 3.1. Energia total da molécula H2 em função da distância H-H.

Algoritmos para minimização de energia

Para moléculas maiores, há muito mais distâncias e ângulos que precisam ser levados em conta. Por exemplo, para uma única molécula de H2O (água), há nove graus de liberdade, três coordenadas de cada átomo. Entretanto, três desses graus de liberdade estão relacionados ao movimento de translação da molécula, e outros três graus estão associados a rotação da molécula. Dessa forma, restam apenas três graus de liberdade para movimentos internos, que podem ser identificados como duas distâncias H-O, e um ângulo H-O-H. A energia total dessa molécula é uma função dessas duas distâncias e desse ângulo. Essa função energia de três variáveis formam uma hiper-superfície no espaço formado pelos graus de liberdade da molécula. A geometria relaxada da molécula é encontrada pela minimização dessa função energia, e para encontrar esse mínimo de energia podemos empregar métodos numéricos de minimização. Para moléculas maiores, há muitos mais graus de liberdade, e podem haver também vários mínimos locais de energia. O processo de relaxação da geometria leva para um dos mínimos locais da energia, mas não garante que esse mínimo é global. Outros algoritmos devem ser empregados para esse fim.

No programa SIESTA, há pelo menos três algoritmos que podem ser usados para relaxar uma geometria molecular. Os métodos implementos no programa são: Gradiente Conjugado (CG), Broyden, e Fast Inertial Relaxation Engine (FIRE). Esse algoritmos são iterativos, sendo formado por loops externos ao loop do ciclo de auto-consistência para resolver as equação de Kohn-Sham. Para cada passo da relaxação, primeiro são solucionadas as equações de Kohn-Sham, e depois são encontras as forças atuando em cada átomo através do gradiente da energia total.

EXEMPLO 03 (Água)

Para relaxar a geometria de uma molécula de H2O, nós vamos partir de uma geometria não relaxada, com a molécula H2O com distância H-O de 1.0 Å, e ângulo H-O-H de 90°, como mostrado na Figura 3.2. Nós vamos incluir o “novo bloco” no nosso arquivo de entrada (input.fdf) com especificações de qual algoritmo será usado (CG, Broyden, FIRE) [MD.TypeOfRun], o número máximo de passos (iterações) para minimizar as forças máximas toleradas [MD.NumCGsteps], e o critério de tolerância da força máxima [MD.MaxForceTol]. A geometria molecular será dita relaxada ser todas as forças atuando em todos os átomos forem menor do que o critério de tolerância especificado no arquivo de entrada.

O bloco “molecular dynamics” que será incluído no input.fdf é:

# in file: input.fdf

    # -- MOLECULAR DYNAMICS --

MD.TypeOfRun Broyden
MD.NumCGsteps 1000
MD.MaxForceTol 0.01 eV/Ang
initial-final
Figura 3.2. Geometria da molécula H2O inicial e final.

Comparando algoritmos

Nós podemos comparar o desempenho dos três algoritmos implementos do SIESTA com a convergência das forças com o número de passos. Como pode ser visto na Figura 3.3, o algoritmo Broyden foi o mais eficiente para a molécula de água, atingindo o critério de tolerância das forças antes dos outros algoritmos. Para outras moléculas e sólidos, o algoritmo FIRE ou CG podem ser mais eficientes. O algoritmo FIRE é mais eficiente quando há muitos átomos de massas atômicas muitos próximas.

relax
Figura 3.3. Evolução da força máxima com o número de passos do processo de minimização.

Anterior: 02 – Primeiras simulações com o SIESTA

Próximo: 04 – Energia de ligação