molecular-dynamics

Dinâmica Molecular

Visão Geral

A simulação de dinâmica molecular (MD) modela computacionalmente a evolução temporal de sistemas moleculares integrando as equações de movimento de Newton. Esta skill abrange duas ferramentas complementares:

• OpenMM (https://openmm.org/): Mecanismo de simulação de MD de alto desempenho com suporte a GPU, API Python e suporte flexível a campos de força
• MDAnalysis (https://mdanalysis.org/): Biblioteca Python para ler, escrever e analisar trajetórias de MD de todos os principais pacotes de simulação

Instalação:

conda install -c conda-forge openmm mdanalysis nglview
# ou
uv pip install openmm mdanalysis

Quando Usar Esta Skill

Use dinâmica molecular quando:

• Análise de estabilidade de proteínas: Como uma mutação afeta a dinâmica da proteína?
• Simulações de ligação de fármacos: Caracterizar o modo de ligação e o tempo de residência de um ligante
• Amostragem conformacional: Explorar a flexibilidade de proteínas e mudanças conformacionais
• Interação proteína-proteína: Modelar a dinâmica da interface e a energética de ligação
• Análise de RMSD/RMSF: Quantificar flutuações estruturais a partir de uma estrutura de referência
• Estimativa de energia livre: Calcular energia livre de ligação ou energia livre conformacional
• Simulações de membranas: Modelar proteínas em bicamadas lipídicas
• Proteínas intrinsecamente desordenadas: Estudar conjuntos conformacionais de IDR

Fluxo de Trabalho Principal: Simulação OpenMM

1. Preparação do Sistema

from openmm.app import *
from openmm import *
from openmm.unit import *
import sys

def prepare_system_from_pdb(pdb_file, forcefield_name="amber14-all.xml",
                              water_model="amber14/tip3pfb.xml"):
    """
    Prepare um sistema OpenMM a partir de um arquivo PDB.

    Args:
        pdb_file: Caminho para o arquivo PDB limpo (use PDBFixer para arquivos PDB brutos)
        forcefield_name: Arquivo XML do campo de força
        water_model: Arquivo XML do modelo de água

    Returns:
        pdb, forcefield, system, topology
    """
    # Carregar PDB
    pdb = PDBFile(pdb_file)

    # Carregar campo de força
    forcefield = ForceField(forcefield_name, water_model)

    # Adicionar hidrogênios e solvatar
    modeller = Modeller(pdb.topology, pdb.positions)
    modeller.addHydrogens(forcefield)

    # Adicionar caixa de solvente (10 Å de preenchimento, NaCl 150 mM)
    modeller.addSolvent(
        forcefield,
        model='tip3p',
        padding=10*angstroms,
        ionicStrength=0.15*molar
    )

    print(f"Sistema: {modeller.topology.getNumAtoms()} átomos, "
          f"{modeller.topology.getNumResidues()} resíduos")

    # Criar sistema
    system = forcefield.createSystem(
        modeller.topology,
        nonbondedMethod=PME,         # Malha de Partículas Ewald para eletrostática de longo alcance
        nonbondedCutoff=1.0*nanometer,
        constraints=HBonds,           # Restringir ligações de hidrogênio (permite passo de tempo de 2 fs)
        rigidWater=True,
        ewaldErrorTolerance=0.0005
    )

    return modeller, system

2. Minimização de Energia

from openmm.app import *
from openmm import *
from openmm.unit import *

def minimize_energy(modeller, system, output_pdb="minimized.pdb",
                     max_iterations=1000, tolerance=10.0):
    """
    Minimizar a energia do sistema para remover choques estéricos.

    Args:
        modeller: Objeto Modeller com topologia e posições
        system: Sistema OpenMM
        output_pdb: Caminho para salvar a estrutura minimizada
        max_iterations: Número máximo de passos de minimização
        tolerance: Critério de convergência em kJ/mol/nm

    Returns:
        objeto de simulação com posições minimizadas
    """
    # Configurar integrador (não importa para minimização)
    integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)

    # Criar simulação
    # Usar GPU se disponível (CUDA ou OpenCL), caso contrário CPU
    try:
        platform = Platform.getPlatformByName('CUDA')
        properties = {'DeviceIndex': '0', 'Precision': 'mixed'}
    except Exception:
        try:
            platform = Platform.getPlatformByName('OpenCL')
            properties = {}
        except Exception:
            platform = Platform.getPlatformByName('CPU')
            properties = {}

    simulation = Simulation(
        modeller.topology, system, integrator,
        platform, properties
    )
    simulation.context.setPositions(modeller.positions)

    # Verificar energia inicial
    state = simulation.context.getState(getEnergy=True)
    print(f"Energia inicial: {state.getPotentialEnergy()}")

    # Minimizar
    simulation.minimizeEnergy(
        tolerance=tolerance*kilojoules_per_mole/nanometer,
        maxIterations=max_iterations
    )

    state = simulation.context.getState(getEnergy=True, getPositions=True)
    print(f"Energia minimizada: {state.getPotentialEnergy()}")

    # Salvar estrutura minimizada
    with open(output_pdb, 'w') as f:
        PDBFile.writeFile(simulation.topology, state.getPositions(), f)

    return simulation

3. Equilibração NVT

from openmm.app import *
from openmm import *
from openmm.unit import *

def run_nvt_equilibration(simulation, n_steps=50000, temperature=300,
                            report_interval=1000, output_prefix="nvt"):
    """
    Equilibração NVT: N, V, T constantes.
    Equilibrar velocidades para a temperatura alvo.

    Args:
        simulation: Simulação OpenMM (após minimização)
        n_steps: Número de passos MD (50000 × 2fs = 100 ps)
        temperature: Temperatura em Kelvin
        report_interval: Passos entre relatórios de dados
        output_prefix: Prefixo de arquivo para trajetória e log
    """
    # Adicionar restrições de posição para a espinha dorsal durante NVT
    # (Opcional: restringir átomos pesados)

    # Definir temperatura
    simulation.context.setVelocitiesToTemperature(temperature*kelvin)

    # Adicionar repórteres
    simulation.reporters = []

    # Arquivo de log
    simulation.reporters.append(
        StateDataReporter(
            f"{output_prefix}_log.txt",
            report_interval,
            step=True,
            potentialEnergy=True,
            kineticEnergy=True,
            temperature=True,
            volume=True,
            speed=True
        )
    )

    # Trajetória DCD (formato binário compacto)
    simulation.reporters.append(
        DCDReporter(f"{output_prefix}_traj.dcd", report_interval)
    )

    print(f"Executando equilibração NVT: {n_steps} passos ({n_steps*2/1000:.1f} ps)")
    simulation.step(n_steps)
    print("Equilibração NVT concluída")

    return simulation

4. Equilibração NPT e Produção

def run_npt_production(simulation, n_steps=500000, temperature=300, pressure=1.0,
                        report_interval=5000, output_prefix="npt"):
    """
    Execução de produção NPT: N, P, T constantes.

    Args:
        n_steps: Passos de produção (500000 × 2fs = 1 ns)
        temperature: Temperatura em Kelvin
        pressure: Pressão em bar
        report_interval: Passos entre relatórios
    """
    # Adicionar barostato Monte Carlo para controle de pressão
    system = simulation.context.getSystem()
    system.addForce(MonteCarloBarostat(pressure*bar, temperature*kelvin, 25))
    simulation.context.reinitialize(preserveState=True)

    # Atualizar repórteres
    simulation.reporters = []
    simulation.reporters.append(
        StateDataReporter(
            f"{output_prefix}_log.txt",
            report_interval,
            step=True,
            potentialEnergy=True,
            temperature=True,
            density=True,
            speed=True
        )
    )
    simulation.reporters.append(
        DCDReporter(f"{output_prefix}_traj.dcd", report_interval)
    )

    # Salvar pontos de verificação
    simulation.reporters.append(
        CheckpointReporter(f"{output_prefix}_checkpoint.chk", 50000)
    )

    print(f"Executando produção NPT: {n_steps} passos ({n_steps*2/1000000:.2f} ns)")
    simulation.step(n_steps)
    print("MD de produção concluída")
    return simulation

Análise de Trajetória com MDAnalysis

1. Carregar Trajetória

import MDAnalysis as mda
from MDAnalysis.analysis import rms, align, contacts
import numpy as np
import matplotlib.pyplot as plt

def load_trajectory(topology_file, trajectory_file):
    """
    Carregar uma trajetória MD com MDAnalysis.

    Args:
        topology_file: Arquivo de topologia PDB, PSF ou outro
        trajectory_file: Trajetória DCD, XTC, TRR ou outra
    """
    u = mda.Universe(topology_file, trajectory_file)
    print(f"Universo: {u.atoms.n_atoms} átomos, {u.trajectory.n_frames} quadros")
    print(f"Intervalo de tempo: 0 a {u.trajectory.totaltime:.0f} ps")
    return u

2. Análise RMSD

def compute_rmsd(u, selection="backbone", reference_frame=0):
    """
    Calcular RMSD de átomos selecionados em relação ao quadro de referência.

    Args:
        u: Universo MDAnalysis
        selection: String de seleção de átomos (sintaxe MDAnalysis)
        reference_frame: Índice do quadro para estrutura de referência

    Returns:
        array numpy de valores (tempo, rmsd)
    """
    # Alinhar trajetória para minimizar RMSD
    aligner = align.AlignTraj(u, u, select=selection, in_memory=True)
    aligner.run()

    # Calcular RMSD
    R = rms.RMSD(u, select=selection, ref_frame=reference_frame)
    R.run()

    rmsd_data = R.results.rmsd  # colunas: quadro, tempo, RMSD
    return rmsd_data

def plot_rmsd(rmsd_data, title="RMSD ao longo do tempo", output_file="rmsd.png"):
    """Plotar RMSD ao longo do tempo de simulação."""
    fig, ax = plt.subplots(figsize=(10, 4))
    ax.plot(rmsd_data[:, 1] / 1000, rmsd_data[:, 2], 'b-', linewidth=0.5)
    ax.set_xlabel("Tempo (ns)")
    ax.set_ylabel("RMSD (Å)")
    ax.set_title(title)
    ax.axhline(rmsd_data[:, 2].mean(), color='r', linestyle='--',
               label=f'Média: {rmsd_data[:, 2].mean():.2f} Å')
    ax.legend()
    plt.tight_layout()
    plt.savefig(output_file, dpi=150)
    return fig

3. Análise RMSF (Flexibilidade por Resíduo)

def compute_rmsf(u, selection="backbone", start_frame=0):
    """
    Calcular RMSF por resíduo (flexibilidade).

    Returns:
        arrays de resids, rmsf_values
    """
    # Selecionar átomos
    atoms = u.select_atoms(selection)

    # Calcular RMSF
    R = rms.RMSF(atoms)
    R.run(start=start_frame)

    # Média por resíduo
    resids = []
    rmsf_per_res = []
    for res in u.select_atoms(selection).residues:
        res_atoms = res.atoms.intersection(atoms)
        if len(res_atoms) > 0:
            resids.append(res.resid)
            rmsf_per_res.append(R.results.rmsf[res_atoms.indices].mean())

    return np.array(resids), np.array(rmsf_per_res)

4. Contatos Proteína-Ligante

def analyze_contacts(u, protein_sel="protein", ligand_sel="resname LIG",
                      radius=4.5, start_frame=0):
    """
    Rastrear contatos proteína-ligante ao longo da trajetória.

    Args:
        radius: Distância de corte para contato em Angstroms
    """
    protein = u.select_atoms(protein_sel)
    ligand = u.select_atoms(ligand_sel)

    contact_frames = []
    for ts in u.trajectory[start_frame:]:
        # Encontrar átomos da proteína dentro do raio do ligante
        distances = contacts.contact_matrix(
            protein.positions, ligand.positions, radius
        )
        contact_residues = set()
        for i in range(distances.shape[0]):
            if distances[i].any():
                contact_residues.add(protein.atoms[i].resid)
        contact_frames.append(contact_residues)

    return contact_frames

Guia de Seleção de Campo de Força

Sistema	Campo de Força Recomendado	Modelo de Água
Proteínas padrão	AMBER14 (amber14-all.xml)	TIP3P-FB
Proteínas + pequenas moléculas	AMBER14 + GAFF2	TIP3P-FB
Proteínas de membrana	CHARMM36m	TIP3P
Ácidos nucleicos	AMBER99-bsc1 ou AMBER14	TIP3P
Proteínas desordenadas	ff19SB ou CHARMM36m	TIP3P

Ferramentas de Preparação do Sistema

PDBFixer (para arquivos PDB brutos)

from pdbfixer import PDBFixer
from openmm.app import PDBFile

def fix_pdb(input_pdb, output_pdb, ph=7.0):
    """Corrigir problemas comuns de PDB: resíduos faltantes, átomos, adicionar H, padronizar."""
    fixer = PDBFixer(filename=input_pdb)
    fixer.findMissingResidues()
    fixer.findNonstandardResidues()
    fixer.replaceNonstandardResidues()
    fixer.removeHeterogens(True)    # Remover água/ligantes
    fixer.findMissingAtoms()
    fixer.addMissingAtoms()
    fixer.addMissingHydrogens(ph)

    with open(output_pdb, 'w') as f:
        PDBFile.writeFile(fixer.topology, fixer.positions, f)

    return output_pdb

GAFF2 para Pequenas Moléculas (via OpenFF Toolkit)

# Para parametrização de ligantes, use o OpenFF toolkit ou ACPYPE
# uv pip install openff-toolkit
from openff.toolkit import Molecule, ForceField as OFFForceField
from openff.interchange import Interchange

def parameterize_ligand(smiles, ff_name="openff-2.0.0.offxml"):
    """Gerar parâmetros GAFF2/OpenFF para uma pequena molécula."""
    mol = Molecule.from_smiles(smiles)
    mol.generate_conformers(n_conformers=1)

    off_ff = OFFForceField(ff_name)
    interchange = off_ff.create_interchange(mol.to_topology())
    return interchange

Melhores Práticas

• Sempre minimizar antes do MD: Estruturas PDB brutas têm choques estéricos
• Equilibrar antes da produção: NVT (50–100 ps) → NPT (100–500 ps) → Produção
• Usar GPU: Simulações são 10–100× mais rápidas em GPU (CUDA/OpenCL)
• Passo de tempo de 2 fs com restrições HBonds: Padrão; use 4 fs com HMR (repartição de massa de hidrogênio)
• Analisar apenas trajetória equilibrada: Descartar os primeiros 20–50% como equilibração
• Salvar pontos de verificação: Execuções de MD podem falhar; pontos de verificação permitem reinício
• Condições de contorno periódicas: Necessário para sistemas solvatados
• PME para eletrostática: Mais preciso que métodos de corte para sistemas carregados

Recursos Adicionais

• Documentação do OpenMM: https://openmm.org/documentation.html
• Guia do usuário MDAnalysis: https://docs.mdanalysis.org/
• GROMACS (mecanismo MD alternativo): https://manual.gromacs.org/
• NAMD (alternativo): https://www.ks.uiuc.edu/Research/namd/
• CHARMM-GUI (gerador de sistemas): https://www.charmm-gui.org/