Cálculos de Dinâmica Molecular
HistoryEdited Nov 3, 2014 10:52 AM by Gustavo Seabra…
Os cálculos de dinâmica molecular normalmente incluem 4 etapas:
Além dos arquivos com informações sobre o sistema (coordenadas iniciais e topologia), precisamos definir os procedimentos envolvidos em cada uma dessas etapas, o que é feito com um terceiro arquivo de controle. Vamos analisar cada uma dessas etapas agora.
1. Minimização de Energia.
Antes dos cálculos de dinâmica molecular propriamente dita, é aconselhável fazer uma rápida minimização de energia, para eliminar possíveis contatos entre os resíduos. Para isso, vamos preparar um arquivo de controle, min.mdin, com as especificações do cálculo:
imin=1, maxcyc=10000, ncyc=1000,
ntpr=250, ntwx=250, ntwr=1000,
Uma vez preparado esse arquivo, podemos executar o programa de dinâmica molecular (em nosso caso, o sander), com o comando
$ sander -O -i min.mdin -o min.mdout -p ns3.prmtop -c ns3.inpcrd -inf min.mdinfo -r min.restrt –x min.mdcrd &
Você pode acompanhar o progresso da simulação com o comando "tail", que já usamos antes:
$ tail -f min.mdout
Quando a minimização acabar, usaremos a estrutura final escrita no arquivo min.restrt como entrada para os cálculos de dinâmica.
A estrutura que temos está a uma temperatura de 0K. Precisamos trazer a estrutura para a temperatura de 310K, aproximadamente a temperatura do corpo humano. Para isso, usaremos um esquema de aquecimento gradual, conforme o arquivo heat.mdin, abaixo. O aquecimento é feito a volume constante, de forma a permitir os átomos do sistema se ajustarem ao volume disponível.
NVT Heating to 310K in 2ns, with 1.0fs timestep
1.5 ns heating + 0.5ns stabilizing
ntb=1, cut=10.0,
tempi=0.0, temp0=310.0, nmropt=1,
nstlim=2000000, dt=0.001,
ntpr=5000, ntwx=5000, ntwr=50000,
type='TEMP0', istep1=0, istep2=1500000,
value1=0.0, value2=310.0,
type='TEMP0', istep1=1500000, istep2=2000000,
value1=310.0, value2=310.0,
O comando a ser dado agora é:
$ sander -O -i heat.mdin -o heat.mdout -p ns3.prmtop -c min.restrt -inf heat.mdinfo -r heat.restrt –x heat.mdcrd &
Note a mudança nos nomes dos arquivos. Esse cálculo levará o sistema lentamente de 0 a 310K.
3. Equilibração
Agora que o sistema está na temperatura desejada, vamos ligar o controle de pressão e permitir que ele relaxe um pouco mais, agora sob pressão e temperatura constante. Esse último passo vai deixar o sistema nas condições reais que precisamos para fazer a simulação propriamente dita, para coleta de dados. O arquivo npt_eq.mdin abaixo controla esse processo:
NPT Equilibration to 310K and 1 atm (1.01325 bar) pressure
in 2ns, with 1.0fs timestep, using berendsen barostat
tempi=310.0, temp0=310.0,
ntb=2, cut=10.0,
ntp=1, barostat=1, pres0=1.01325, taup=1.0,
ntxo=2, ioutfm=1,
O novo comando para iniciar o cálculo é:
$ sander -O -i npt_eq.mdin -o npt_eq.mdout -p ns3.prmtop -c heat.restrt -inf npt_eq.mdinfo -r npt_eq.restrt –x npt_eq.mdcrd &
Ao término desse cálculo, teremos o sistema a 310K e 1 atm de pressão. Podemos agora prosseguir para a etapa de obtenção de dados, ou "produção".
A etapa de produção envolve simplesmente simulações mais longas, em nosso caso a temperatura e pressão constantes, para gerar uma grande quantidade de estruturas compatíveis com as condições especificadas. Para isso, vamos apenas modificar o arquivo mdin anterior para um tempo maior de simulação. O arquivo prod_0-15ns.mdin, abaixo, solicita o sander a calcular 15 ns de dinâmica molecular a temperatura de 310 K e pressão 1 atm:
NPT production at 310K and 1 atm (1.01325 bar) pressure
15ns, with 2.0fs timestep
nstlim=7500000, dt=0.002,
ntpr=50000, ntwx=50000, ntwr=500000,
Ese cálculo pode ser repetido quantas vezes forem necessárias até que se acumule a amostragem desejada, sempre se iniciando a partir do último conjunto de coordenadas e velocidades. No final, teremos um conjunto de arquivos de trajetórias (mdcrd), a ser utilizado para a próxima etapa de tratamento de dados.