Nova versão do programa para extrair dados do arquivo de log do LAMMPS

Estou disponibilizando uma nova versão do programa em Python para extrair informações impressas com o comando thermo dos arquivos de log produzidos nas simulações de dinâmica molecular feitas com o LAMMPS. Para baixar essa versão, clique no link abaixo:

https://drive.google.com/file/d/1Ac7zsS8RPGQYqPw89H0VoDfbILUDP3Wy/view?usp=sharing

Além de extrair a evolução de cada variável definida no comando thermo_style em arquivos distintos e calcular as médias, desvios e médias em bloco como na versão anterior, a nova versão inclui duas novas funcionalidades:

1- Permite plotar diretamente as variáveis em qualquer unidade de tempo desejada ao invés de plotar em função do número de passos como na primeira versão. Para isso, basta informar o passo de integração usado na simulação na variável dt no programa. (o default é dt = 1, que pode ser visto como plotar em função do número de passos que era o caso da versão anterior).

2- Faz regressão linear de cada propriedade calculada em função do tempo usando os dados após o mesmo número de passos em que se inicia o cálculo de médias e desvios (definido pela variável avstart). Os coeficientes angular e linear de cada propriedade em cada run são dados no arquivo linear_regression.txt . Essas regressões tem duas utilidades principais. Primeiramente, algumas variáveis podem apresentar realmente um comportamento linear, como é o caso do desvio médio quadrático das partículas em um líquido, e o coeficiente angular pode se relacionar a propriedades interessantes (no caso do desvio médio quadrático, relaciona-se com o coeficiente de difusão). Para outras propriedades, as quais espera-se que apenas flutuem em torno da média após atingir o equilíbrio, um coeficiente angular não-desprezível pode indicar que seu sistema ainda não relaxou, havendo alguma mudança lenta em andamento, ou que o número de passos excluído do início da simulação não foi grande o bastante. Caso o valor de avstart seja maior que o número de passos da respectiva simulação, as regressões não são calculadas e uma mensagem é impressa no arquivo linear_regression.txt.

Para cálculo das regressões, o programa necessita do módulo “scipy” instalado, algo que não era necessário na primeira versão. Caso prefira usar a versão anterior, ela ainda está disponível no link https://drive.google.com/file/d/1zH7i8ouuDwEtLxmc9pUUbc87Zshz0t46/view?usp=sharing .

Para mais informações sobre a realização e análise de simulações no LAMMPS, cheque nosso tutorial.

Script em python para extrair dados dos arquivos de log do LAMMPS

Olá,

Estou disponibilizando o programa extract_lammps.py, escrito por mim, para extrair propriedades calculadas e impressas no arquivo de log do LAMMPS por meio do comando thermo ao longo de uma trajetória. Para baixá-lo, basta acessar o link abaixo:

https://drive.google.com/open?id=1zH7i8ouuDwEtLxmc9pUUbc87Zshz0t46

O programa irá gerar arquivos com a evolução temporal de cada variável impressa com o comando thermo em função do número de passos além de calcular as respectivas médias e desvios após um determinado número de passos determinado pelo usuário. Em caso de simulações sequenciais com um mesmo arquivo de log, o programa analisa cada simulação separadamente e identifica os arquivos e dados de média e desvio por meio do prefixo “run0_” para a primeira simulação, “run1_” para a segunda, etc.

O único requisito do programa é que a primeira variável a ser impressa com o comando thermo seja o passo de simulação (“Step”), assim, o input do LAMMPS deve conter o comando

thermo_style custom step (outras variáveis que se deseja imprimir)

Para mais detalhes sobre simulações no LAMMPS, veja o tutorial na página https://kalilbn.wordpress.com/tutorial-lammps-mistura-de-agua-e-etanol-parte-1/

 

pmx – pacote para análises de trajetórias usando python

Olá,

Recentemente comecei a trabalhar com o pacote pmx, que possibilita desenvolver ferramentas de análise para extrair informações de trajetórias no formato xtc, que é o formato empregado pelo Gromacs. Ele pode ser baixado da página https://code.google.com/p/pmx/ e necessita ter o numpy e o scipy instalados, esses dois últimos presentes nos repositórios de muitas distribuições Linux. Eu achei bem simples o desenvolvimento de ferramentas usando-o e recomendo para quem está trabalhando com dinâmica molecular e muitas vezes precisa de um dado relativamente simples mas que não possui nenhuma ferramenta padrão para obtê-lo.

Como exemplo de aplicação, fiz um programa para calcular o número de moléculas de água na primeira camada de coordenação de um soluto em função do tempo. Ele basicamente calcula a distância entre o oxigênio da água e o soluto e armazena quantas moléculas estão a uma distância menor que a especificada. Eu o apliquei para analisar o número de moléculas na primeira camada de coordenação de um íon sódio na simulação feita no tutorial de íons em água (https://kalilbn.wordpress.com/dinamica-molecular-conceitos-fundamentais/371-2/), sendo a distância de corte tomada como a posição do primeiro mínimo no gráfico de distribuição radial. Tal programa encontra-se disponível em https://drive.google.com/file/d/0BxbfUlGvt1wJT0FNLWNlRzdPLWM/edit?usp=sharing e está totalmente comentado, de modo a auxiliar quem quiser começar a desenvolver suas próprias ferramentas de análise (note que ele só funcionará se estiver com o pmx devidamente instalado).

O único problema que tive com o pmx após a instalação foi que ele não estava reconhecendo automaticamente a biblioteca libgmx, o que resolvi copiando o arquivo como root:

cp /usr/lib/libgmx.so.8 /usr/lib/libgmx.so

(observação: copie a biblioteca ao invés de apenas renomear para evitar problemas com outros softwares como o gromacs)

e definindo o caminho:

GMX_DLL=/usr/lib/

export GMX_DLL

Para não ter que usar os dois últimos comandos toda vez que reiniciar a máquina, pode acrescentar essas duas linhas ao final do arquivo .profile que está em seu diretório home. Após essas pequenas modificações, o pmx funcionou perfeitamente. Outras instruções mais detalhadas podem ser achadas na página (em inglês): https://code.google.com/p/pmx/ . Ele possui outras funções além da análise de xtc, porém eu mesmo ainda não testei.

Referência do pmx: D. Seeliger and Bert L. de Groot, Biophys. J. 98(10):2309-2316 (2010)  

Página sobre cinética química em sistemas não-isotérmicos e nuvem de tags

Em mais uma página sobre cinética química, serão abordados os efeitos da temperatura sobre a velocidade das reações químicas e também o efeito da variação de energia interna resultante da reação sobre a temperatura do sistema. É apresentada nessa página uma nova versão do programa kinetica, chamada de kineticaT, onde, além dos parâmetros cinéticos das reações químicas e concentrações iniciais, é especificada também a temperatura do sistema, a variação de energia interna da reação química, a capacidade calorífica do sistema e sua taxa de troca de calor com as vizinhanças.

Além de uma apresentação dos conceitos fundamentais referentes à dependência da velocidade das reações com a temperatura, são tratados 4 casos distintos como exemplos de aplicações do programa kineticaT:

1. Efeito da temperatura sobre a velocidade da reação em um sistema isotérmico
2. Efeitos da taxa de transferência de calor entre o sistema e as vizinhanças em um sistema não-isotérmico para uma reação exotérmica e uma endotérmica
3. Relaxação de um sistema em equilíbrio como resultado de uma variação brusca na temperatura das vizinhanças
4. Histerese sobre uma reação reversível em um sistema submetido a oscilações de temperatura.

Além de gerar os mesmos gráficos do programa kinetica, o kineticaT gera ainda gráficos de taxas de troca de calor pelo sistema, temperatura e pressão em função do tempo, concentrações em função da temperatura entre outros.

Link para a página: https://kalilbn.wordpress.com/cinetica-quimica-i-introducao-e-reacoes-de-primeira-ordem/cinetica-quimica-com-temperatura-variavel/

Além dessa página e do programa kineticaT, adicionei também ao diretório do blog no Google Drive os programa kinetica e kinetica_oscilatorio. Clique no link abaixo para acessá-lo:

https://drive.google.com/folderview?id=0BxbfUlGvt1wJOTczaGZuVEZ6cEU&usp=sharing

Outra adição que acabo de fazer foi a introdução de uma nuvem com as 40 tags mais usadas no blog, acessível por meio do menu lateral na página principal, sendo essa mais uma forma de localizar conteúdos dentro do blog.

Pasta de arquivos do blog no Google Drive

Olá,

Dado ao WordPress não permitir armazenar arquivos de texto simples para download, impedindo colocar para download direto arquivos dos programas em python e de input dos programas utilizados nos tutoriais,  eu estava até o momento incluindo esses ou como página web ou como arquivos .doc para download. Isso gerava alguns problemas de perda de formatação dependendo inclusive de qual navegador ou editor de texto era usado para abrí-lo, o que impedia em alguns casos o uso dos programas ou execução apropriada de tutoriais.

Para resolver isso de modo mais prático, estou hospedando primeiro os programas em python e depois hospedarei todos os arquivos de input para os tutoriais do Gromacs em uma pasta no Google Drive, acessível pelo endereço:

https://drive.google.com/folderview?id=0BxbfUlGvt1wJOTczaGZuVEZ6cEU&usp=sharing

Vou terminar de colocar todos os arquivos nos próximos dias e criar uma página detalhando cada programa. Estou aproveitando também para revisar o conteúdo das páginas em que esses são utilizados e já alterar o link da página para o download no Google Drive.

Para executar qualquer um dos programas em python, basta digitar no terminal “python” seguido pelo nome do programa, no mesmo diretório em que o programa estiver salvo.

Curvas de titulação geradas pelo programa pyH, um dos programas já hospedados no Google Drive do blog. pKa(HX) = 2,0 e pKa(HY) = 5,0 titulados com base forte 0,1 mol/L. Mais detalhes na página:https://kalilbn.wordpress.com/acidos-e-bases-definicoes/acidos-e-bases-tratamento-exato/

Programa para gerar topologias para o Gromacs

Olá,

Após um longo período sem atualizações no blog dada minha defesa de mestrado, estou postando a primeira versão do programa Gentop, que é um programa em python para gerar arquivos de topologia molecular formatados para o Gromacs dada uma estrutura em formato xyz. O programa utiliza um critério de distâncias para determinar quais átomos formam ligações entre si na molécula e, a partir do conjunto determinado de ligações químicas, determina os ângulos e diedros presentes na molécula e estima os tipos atômicos no campo de força OPLS-AA com base em algumas funções orgânicas simples, como alcano, álcool e amina primária. Tais parâmetros devem ser tomados apenas como guia para diferenciar entre alguns dos átomos presentes na molécula, visto que o algoritmo para determiná-los é no momento bastante limitado, sendo que o usuário deve buscar no campo de força quais são os parâmetros mais apropriados para a molécula. Além disso, como o programa se baseia unicamente em critérios de distâncias para determinar as ligações químicas, é necessário que a estrutura da molécula esteja minimamente otimizada para que essas sejam determinadas corretamente.

Essa versão ainda está sendo testada e pretendo em breve melhorá-la com um algoritmo que permita identificar de modo correto uma maior variedade de funções orgânicas. Segue abaixo o link para download:

https://drive.google.com/file/d/0BxbfUlGvt1wJT2xuRDF2LTRfOEE/edit?usp=sharing

Para executar, basta salvar o programa no mesmo diretório contendo a estrutura em formato xyz e digitar python gentop.py

Exemplo de programa simples em python comentado

Para aqueles interessados em aprender um pouco de programação em python, estou colocando em anexo nesse post um arquivo .doc com um programa simples para o cálculo numérico da integral definida, da derivada primeira e da derivada segunda de uma função de uma variável. São fornecidos com outputs quatro gráficos (arquivos .dat que podem ser abertos no Grace, gnuplot ou mesmo importados no Excel para plotagem dos gráficos) com a função propriamente dita, sua integral definida entre o ponto inicial ini e cada ponto x posterior e as derivadas primeira e segunda, além de imprimir na tela do terminal os pontos de máximo e mínimo da função (determinados com base nos sinais das derivadas primeira e segunda)

A integral é calculada usando a idéia da soma de Riemann, ou seja, a área abaixo da curva f(x) é dividida em retângulos de altura f(x) e largura dx e, sendo dx um número muito pequeno, a soma das áreas desses retângulos, feita através de um loop no programa, fornece o valor da integral I. Para o cálculo da derivada primeira, toma-se a diferença entre a função f(x) e a função f(x+dx) dividida por dx. Para o cálculo da derivada segunda, que é a “derivada da derivada”, poderíamos ter feito da mesma maneira, mas tomamos a diferença entre f'(x-dx) e f'(x) dividida por dx, onde f'(x) é a derivada primeira de f(x), escolher entre fazer f(x) – f(x-dx) ou f(x+dx) – f(x) para a derivada no ponto x gera sim uma pequena diferença numérica, mas essa torna-se desprezível ao usar valores muito pequenos de dx.

Ao baixar o arquivo, verá que ele tem várias linhas iniciadas com # e com a fonte em azul que são comentários, ou seja, informações que não serão executadas pelo programa, mas são informações sobre o que cada linha do programa faz para que o leitor possa estudá-lo passo a passo e aprender um pouco de programação na prática. As estruturas usadas nesse programa permitem resolver diversos outros problemas matemáticos além dos apresentados sem grandes modificações.

Para executar o programa, simplesmente copie o texto do arquivo .doc e cole num arquivo de texto simples (usando o Bloco de Notas, gedit ou o próprio vi no terminal) e execute em um terminal digitando python seguido pelo nome do arquivo. Eu não anexo já como arquivo de texto simples pois o Worpress não me permite, e colar diretamente no post desformata completamente o texto.

Clique para fazer o download.

Claro que, como exemplo de função, selecionei uma que não possui integral analítica (pelo menos não em um intervalo arbitrário) e que é de interesse para a química:

f(y) = (2y² – 1)².exp[-x²]

Essa função é essencialmente a densidade de probabilidade, ou seja, o quadrado da função de onda, do segundo estado excitado do oscilador harmônico quântico não-normalizada e com α = 1. A integral calculada pelo programa, em um intevalo de menos infinito a mais infinito seria, portanto, a integral de normalização da função, é claro que não é possível integrar por soma de Riemann em um intervalo infinito, mas antes e após o intervalo especificado no programa o valor da função é tão pequeno que o erro do truncamento nessa integral torna-se desprezível. Você pode (e recomendo que faça) modificar o programa para calcular a integral de qualquer outra função matemática em qualquer intervalo que desejar. Para se convencer que funciona, teste com algumas funções que saiba calcular a integral e as derivadas analiticamente e compare os resultados com os fornecidos pelo programa. Experimente também variar o número de pontos usados no cálculo e ver como isso afeta o resultado.

Espero que isso seja útil para quem tem interesse em começar a trabalhar com python para aplicações científicas ou matemáticas. No mais, existem muitos tutoriais de python na internet e diversos livros (muitos de domínio público) voltados para públicos de diversos níveis.

Corrigido link quebrado para o programa kinetica

Ao acrescentar a página sobre campos vetoriais, percebi que o link para baixar o programa kinetica, na página introdutória sobre cinética química estava quebrado. Não sei se alguma atualização do WordPress quebrou o link, já que era para um documento .doc e não para uma página, mas agora resolvi o problema e coloquei o código programa em python em uma página convencional ao invés de em um arquivo para download.

Clique aqui para baixar o código do programa. Para uma descrição detalhada, veja essa página,

Página sobre Reações Oscilantes

Acabo de adicionar uma página ao blog sobre cinética de reações químicas oscilantes, onde a concentração de uma ou mais espécies químicas apresenta variações periódicas ao longo do processo, e reações químicas auto-catalíticas, dois temas que estão fortemente relacionados entre si, onde é apresentada uma nova versão do programa kinetica feita para o estudo de reações oscilante utilizando o mecanismo de Lotka-Volterra, dado pela sequência de reações:

A + X → 2 X

X + Y → 2 Y

Y → B

São discutidas a variação das concentrações e das velocidades ao longo do tempo em diversas situações envolvendo esse mecanismo, sendo que em cada ponto discutido acrescentamos um pouco mais de complexidade ao modelo.

Link para a página:

https://kalilbn.wordpress.com/cinetica-quimica-i-introducao-e-reacoes-de-primeira-ordem/reacoes-oscilantes/

Link para a página introdutória de cinética química:

https://kalilbn.wordpress.com/cinetica-quimica-i-introducao-e-reacoes-de-primeira-ordem/

Adicionada página sobre titulações ácido base

Acabo de adicionar mais uma página sobre equilíbrio ácido-base, onde são apresentados resultados obtidos com o pyH para titulações ácido-base simples e de misturas de ácidos.

Link para a página: https://kalilbn.wordpress.com/acidos-e-bases-definicoes/acidos-e-bases-titulacoes/