Dinamica molecolare - Molecular dynamics

Esempio di simulazione di dinamica molecolare in un sistema semplice: deposizione di un atomo di rame (Cu) su un cristallo freddo di rame ( superficie dell'indice di Miller (001) ). Ogni cerchio rappresenta la posizione di un atomo. L'energia cinetica dell'atomo che si avvicina dall'alto viene ridistribuita tra gli altri atomi, quindi invece di rimbalzare rimane attaccata a causa delle forze attrattive tra gli atomi.
Le simulazioni di dinamica molecolare sono spesso utilizzate per studiare i sistemi biofisici. Qui è raffigurata una simulazione dell'acqua a 100 ps.
Una descrizione semplificata dell'algoritmo di simulazione della dinamica molecolare standard, quando viene utilizzato un integratore di tipo predittore-correttore. Le forze possono provenire da potenziali interatomici classici (descritti matematicamente come ) o da metodi quantomeccanici (descritti matematicamente come ). Esistono grandi differenze tra i diversi integratori; alcuni non hanno esattamente gli stessi termini di ordine più alto come indicato nel diagramma di flusso, molti usano anche derivate temporali di ordine superiore e alcuni usano sia il passo temporale corrente che quello precedente negli schemi a passo temporale variabile.

La dinamica molecolare ( MD ) è un metodo di simulazione al computer per analizzare i movimenti fisici di atomi e molecole . Gli atomi e le molecole possono interagire per un determinato periodo di tempo, dando una visione dell' "evoluzione" dinamica del sistema. Nella versione più comune, le traiettorie di atomi e molecole sono determinate risolvendo numericamente le equazioni del moto di Newton per un sistema di particelle interagenti, dove le forze tra le particelle e le loro energie potenziali sono spesso calcolate utilizzando potenziali interatomici o campi di forza della meccanica molecolare . Il metodo è applicato principalmente in fisica chimica , scienza dei materiali e biofisica .

Poiché i sistemi molecolari sono tipicamente costituiti da un vasto numero di particelle, è impossibile determinare analiticamente le proprietà di tali sistemi complessi ; La simulazione MD aggira questo problema utilizzando metodi numerici . Tuttavia, lunghe simulazioni MD sono matematicamente mal condizionate , generando errori cumulativi nell'integrazione numerica che possono essere minimizzati con una corretta selezione di algoritmi e parametri, ma non eliminati del tutto.

Per i sistemi che obbediscono all'ipotesi ergodica , l'evoluzione di una simulazione di dinamica molecolare può essere utilizzata per determinare le proprietà termodinamiche macroscopiche del sistema: le medie temporali di un sistema ergodico corrispondono alle medie di insieme microcanoniche . La MD è stata anche definita "meccanica statistica in base ai numeri" e " visione di Laplace della meccanica newtoniana " di predire il futuro animando le forze della natura e consentendo la comprensione del movimento molecolare su scala atomica.

Storia

MD è stato originariamente sviluppato nei primi anni '50, in seguito ai precedenti successi con le simulazioni Monte Carlo , che risalgono al XVIII secolo, ad esempio nel problema dell'ago di Buffon , ma è stato reso popolare per la meccanica statistica al Los Alamos National Laboratory da Rosenbluth e Metropolis in quello che oggi è noto come algoritmo di Metropolis-Hastings . L'interesse per l'evoluzione temporale dei sistemi di N-corpi risale a molto prima del XVII secolo, a partire da Newton, e continuò nel secolo successivo in gran parte concentrandosi sulla meccanica celeste e su questioni come la stabilità del sistema solare. Molti dei metodi numerici utilizzati oggi sono stati sviluppati durante questo periodo di tempo, che precede l'uso dei computer; ad esempio, l'algoritmo di integrazione più comune utilizzato oggi, l' algoritmo di integrazione Verlet , è stato utilizzato già nel 1791 da Jean Baptiste Joseph Delambre . I calcoli numerici con questi algoritmi possono essere considerati MD "a mano".

Già nel 1941 l'integrazione delle equazioni del moto a molti corpi veniva effettuata con computer analogici. Alcuni hanno intrapreso il lavoro intensivo di modellazione del movimento atomico costruendo modelli fisici, ad esempio utilizzando sfere macroscopiche. Lo scopo era di disporli in modo tale da replicare la struttura di un liquido e utilizzarla per esaminarne il comportamento. JD Bernal ha detto, nel 1962: " ... ho preso un certo numero di palline di gomma e le ho attaccate insieme con aste di una selezione di lunghezze diverse che vanno da 2,75 a 4 pollici. Ho cercato di farlo in primo luogo il più casualmente possibile , lavorando nel mio ufficio, essendo interrotto ogni cinque minuti circa e non ricordando cosa avevo fatto prima dell'interruzione. "

In seguito alla scoperta di particelle microscopiche e allo sviluppo dei computer, l'interesse si è esteso oltre il terreno di prova dei sistemi gravitazionali alle proprietà statistiche della materia. Nel tentativo di comprendere l'origine dell'irreversibilità, Fermi propose nel 1953, e pubblicò nel 1955, l'uso di MANIAC I , sempre al Los Alamos National Laboratory , per risolvere l'evoluzione temporale delle equazioni del moto per un soggetto di sistema a molti corpi a diverse scelte di leggi di forza; oggi, questo lavoro fondamentale è conosciuto come il problema di Fermi–Pasta–Ulam–Tsingou . L'evoluzione temporale dell'energia dall'opera originale è mostrata nella figura a destra.

Una delle prime simulazioni di un sistema a N-corpi è stata effettuata sul MANIAC-I da Fermi e collaboratori per comprendere le origini dell'irreversibilità in natura. Qui viene mostrata l'energia in funzione del tempo per un sistema a 64 particelle.

Nel 1957, Alder e Wainwright usarono un computer IBM 704 per simulare collisioni perfettamente elastiche tra sfere dure . Nel 1960, forse nella prima simulazione realistica della materia, Gibson et al. danno da radiazioni simulato di rame solido utilizzando un tipo di interazione repulsiva Born-Mayer insieme a una forza superficiale coesiva. Nel 1964, Rahman pubblicò simulazioni di argon liquido che utilizzavano un potenziale di Lennard-Jones ; calcoli delle proprietà del sistema, come il coefficiente di autodiffusione , si confrontano bene con i dati sperimentali.

Aree di applicazione e limiti

Utilizzato per la prima volta nella fisica teorica , il metodo MD ha guadagnato popolarità nella scienza dei materiali poco dopo e dagli anni '70 è comune anche in biochimica e biofisica . La MD viene spesso utilizzata per affinare strutture tridimensionali di proteine e altre macromolecole sulla base di vincoli sperimentali dalla cristallografia a raggi X o dalla spettroscopia NMR . In fisica, la MD viene utilizzata per esaminare la dinamica di fenomeni a livello atomico che non possono essere osservati direttamente, come la crescita di film sottili e la subplantazione ionica, e anche per esaminare le proprietà fisiche dei dispositivi nanotecnologici che non sono stati o non possono ancora essere creati . In biofisica e biologia strutturale , il metodo è frequentemente applicato per studiare i moti di macromolecole come proteine e acidi nucleici , che possono essere utili per interpretare i risultati di alcuni esperimenti biofisici e per modellare le interazioni con altre molecole, come nel docking di ligandi . In linea di principio MD può essere utilizzato per la previsione ab initio della struttura proteica simulando il ripiegamento della catena polipeptidica da una bobina casuale .

I risultati delle simulazioni MD possono essere testati attraverso il confronto con esperimenti che misurano la dinamica molecolare, di cui un metodo popolare è la spettroscopia NMR. Le previsioni della struttura derivata dalla MD possono essere testate attraverso esperimenti a livello di comunità nella valutazione critica della previsione della struttura proteica ( CASP ), sebbene il metodo abbia storicamente avuto un successo limitato in questo settore. Michael Levitt , che ha condiviso il Premio Nobel in parte per l'applicazione della MD alle proteine, ha scritto nel 1999 che i partecipanti al CASP di solito non usavano il metodo a causa di " ... un imbarazzo centrale della meccanica molecolare, vale a dire che la minimizzazione dell'energia o la dinamica molecolare in generale porta a un modello che è meno simile alla struttura sperimentale. " I miglioramenti nelle risorse computazionali che consentono traiettorie MD sempre più lunghe, combinati con i moderni miglioramenti nella qualità dei parametri del campo di forza , hanno prodotto alcuni miglioramenti sia nella previsione della struttura che nel perfezionamento del modello di omologia , senza raggiungere il punto di utilità pratica in queste aree; molti identificano i parametri del campo di forza come un'area chiave per ulteriori sviluppi.

La simulazione MD è stata segnalata per lo sviluppo di farmacofori e la progettazione di farmaci. Ad esempio Pinto et al. ha implementato simulazioni MD di complessi Bcl-Xl per calcolare le posizioni medie degli amminoacidi critici coinvolti nel legame del ligando. Carlson et al. ha implementato la simulazione della dinamica molecolare per identificare i composti che completano il recettore causando al contempo un'interruzione minima della conformazione e della flessibilità del sito attivo. Le istantanee della proteina a intervalli di tempo costanti durante la simulazione sono state sovrapposte per identificare le regioni di legame conservate (conservate in almeno tre fotogrammi su undici) per lo sviluppo del farmacoforo. Spyrakis et al. si è basato su un flusso di lavoro di simulazioni MD, impronte digitali per ligandi e proteine ​​(FLAP) e analisi discriminante lineare per identificare le migliori conformazioni ligando-proteina per agire come modelli farmacoforici basati sull'analisi ROC retrospettiva dei farmacofori risultanti. Nel tentativo di migliorare la modellazione della scoperta di farmaci basata sulla struttura, di fronte alla necessità di molti composti modellati, Hatmal et al hanno proposto una combinazione di simulazione MD e analisi dei contatti intermolecolari ligando-recettore per discernere i contatti intermolecolari critici (interazione di legame) da quelli ridondanti in un singolo complesso ligando-proteina. I contatti critici possono quindi essere convertiti in modelli farmacoforici che possono essere utilizzati per lo screening virtuale.

I limiti del metodo sono legati agli insiemi di parametri utilizzati e ai campi di forza della meccanica molecolare sottostanti . Un'esecuzione di una simulazione MD ottimizza l' energia potenziale , piuttosto che l' energia libera della proteina, il che significa che tutti i contributi entropici alla stabilità termodinamica della struttura proteica sono trascurati, inclusa l' entropia conformazionale della catena polipeptidica (il principale fattore che destabilizza la struttura proteica ) ed effetti idrofobici (le principali forze trainanti del ripiegamento delle proteine). Un altro fattore importante sono i legami idrogeno intramolecolari , che non sono esplicitamente inclusi nei moderni campi di forza, ma descritti come interazioni coulombiane di cariche puntiformi atomiche. Questa è un'approssimazione grossolana perché i legami idrogeno hanno una natura parzialmente quantistica e chimica . Inoltre, le interazioni elettrostatiche vengono solitamente calcolate utilizzando la costante dielettrica del vuoto , sebbene la soluzione acquosa circostante abbia una costante dielettrica molto più elevata. L'uso della costante dielettrica macroscopica a brevi distanze interatomiche è discutibile. Infine, le interazioni di van der Waals nella MD sono solitamente descritte dai potenziali di Lennard-Jones basati sulla teoria di Fritz London che è applicabile solo nel vuoto. Tuttavia, tutti i tipi di forze di van der Waals sono in definitiva di origine elettrostatica e quindi dipendono dalle proprietà dielettriche dell'ambiente . La misurazione diretta delle forze di attrazione tra diversi materiali (come la costante di Hamaker ) mostra che "l'interazione tra gli idrocarburi attraverso l'acqua è circa il 10% di quella attraverso il vuoto". La dipendenza dall'ambiente delle forze di van der Waals è trascurata nelle simulazioni standard, ma può essere inclusa sviluppando campi di forza polarizzabili.

Vincoli di progettazione

La progettazione di una simulazione di dinamica molecolare dovrebbe tenere conto della potenza computazionale disponibile. La dimensione della simulazione ( n = numero di particelle), il passo temporale e la durata totale devono essere selezionati in modo che il calcolo possa terminare entro un periodo di tempo ragionevole. Tuttavia, le simulazioni dovrebbero essere sufficientemente lunghe da essere rilevanti per le scale temporali dei processi naturali studiati. Per trarre conclusioni statisticamente valide dalle simulazioni, l'intervallo di tempo simulato dovrebbe corrispondere alla cinetica del processo naturale. Altrimenti, è analogo a trarre conclusioni su come cammina un essere umano quando guarda solo meno di un passo. La maggior parte delle pubblicazioni scientifiche sulla dinamica delle proteine ​​e del DNA utilizzano dati provenienti da simulazioni che vanno da nanosecondi (10 -9 s) a microsecondi (10 -6 s). Per ottenere queste simulazioni sono necessari da diversi giorni CPU ad anni CPU. Algoritmi paralleli consentono di distribuire il carico tra le CPU; un esempio è l'algoritmo di decomposizione spaziale o forzata.

Durante una classica simulazione MD, il compito più impegnativo per la CPU è la valutazione del potenziale in funzione delle coordinate interne delle particelle. All'interno di tale valutazione energetica, quella più costosa è la parte non legata o non covalente. Nella notazione Big O , le simulazioni di dinamica molecolare comuni scalano in base al fatto che tutte le interazioni elettrostatiche e di van der Waals a coppie debbano essere considerate esplicitamente. Questo costo computazionale può essere ridotto impiegando metodi elettrostatici come la sommatoria di Ewald mesh mesh ( ), particella-particella-particella-mesh ( P3M ) o buoni metodi di taglio sferico ( ).

Un altro fattore che influisce sul tempo totale della CPU necessario per una simulazione è la dimensione del timestep di integrazione. Questo è l'intervallo di tempo tra le valutazioni del potenziale. Il timestep deve essere scelto sufficientemente piccolo da evitare errori di discretizzazione (cioè, minore del periodo relativo alla frequenza vibrazionale più veloce nel sistema). I tempi tipici per la MD classica sono dell'ordine di 1 femtosecondo (10 −15 s). Questo valore può essere esteso utilizzando algoritmi come l' algoritmo di vincolo SHAKE , che fissa le vibrazioni degli atomi più veloci (ad esempio gli idrogeni) in posizione. Sono stati inoltre sviluppati più metodi su scala temporale, che consentono tempi estesi tra gli aggiornamenti di forze a lungo raggio più lente.

Per simulare le molecole in un solvente, è necessario scegliere tra solvente esplicito e implicito . Le particelle di solvente esplicite (come i modelli dell'acqua TIP3P , SPC/E e SPC-f ) devono essere calcolate in modo costoso dal campo di forza, mentre i solventi impliciti utilizzano un approccio a campo medio. L'uso di un solvente esplicito è computazionalmente costoso, poiché richiede l'inclusione di circa dieci volte più particelle nella simulazione. Ma la granularità e la viscosità del solvente esplicito sono essenziali per riprodurre determinate proprietà delle molecole di soluto. Ciò è particolarmente importante per riprodurre la cinetica chimica .

In tutti i tipi di simulazioni di dinamica molecolare, la dimensione del riquadro di simulazione deve essere sufficientemente grande da evitare artefatti da condizioni al contorno . Le condizioni al contorno sono spesso trattate scegliendo valori fissi ai bordi (che possono causare artefatti) o impiegando condizioni al contorno periodiche in cui un lato della simulazione torna al lato opposto, imitando una fase bulk (che può causare anche artefatti) .

Rappresentazione schematica del campionamento della superficie di energia potenziale del sistema con dinamica molecolare (in rosso) rispetto ai metodi Monte Carlo (in blu).

Insieme microcanonico (NVE)

Nel microcanonico , il sistema è isolato da variazioni moli (N), il volume (V), e dell'energia (E). Corrisponde ad un processo adiabatico senza scambio termico. Una traiettoria di dinamica molecolare microcanonica può essere vista come uno scambio di energia potenziale e cinetica, con conservazione dell'energia totale. Per un sistema di N particelle con coordinate e velocità , la seguente coppia di equazioni differenziali del primo ordine può essere scritta nella notazione di Newton come

La funzione di energia potenziale del sistema è una funzione delle coordinate delle particelle . Viene indicato semplicemente come il potenziale in fisica o il campo di forza in chimica. La prima equazione deriva dalle leggi del moto di Newton ; la forza che agisce su ciascuna particella del sistema può essere calcolata come il gradiente negativo di .

Per ogni passo temporale, la posizione e la velocità di ciascuna particella possono essere integrate con un metodo di integrazione simplettico come l' integrazione di Verlet . L'evoluzione temporale di e viene chiamata traiettoria. Date le posizioni iniziali (ad esempio, dalla conoscenza teorica) e le velocità (ad esempio, gaussiana randomizzata), possiamo calcolare tutte le posizioni e le velocità future (o passate).

Una frequente fonte di confusione è il significato della temperatura in MD. Comunemente abbiamo esperienza con temperature macroscopiche, che coinvolgono un numero enorme di particelle. Ma la temperatura è una grandezza statistica. Se c'è un numero sufficiente di atomi, la temperatura statistica può essere stimata dalla temperatura istantanea , che si trova eguagliando l'energia cinetica del sistema a nk B T /2 dove n è il numero di gradi di libertà del sistema.

Un fenomeno correlato alla temperatura si verifica a causa del piccolo numero di atomi utilizzati nelle simulazioni MD. Ad esempio, si consideri la simulazione della crescita di un film di rame partendo da un substrato contenente 500 atomi e un'energia di deposizione di 100 eV. Nel mondo reale, i 100 eV dell'atomo depositato verrebbero rapidamente trasportati e condivisi tra un gran numero di atomi ( o più) senza grandi cambiamenti di temperatura. Quando ci sono solo 500 atomi, tuttavia, il substrato viene quasi immediatamente vaporizzato dalla deposizione. Qualcosa di simile accade nelle simulazioni biofisiche. La temperatura del sistema in NVE è naturalmente aumentata quando le macromolecole come le proteine ​​subiscono cambiamenti conformazionali esotermici e legame.

Ensemble canonico (NVT)

Nel insieme canonico , quantità di sostanza (N), il volume (V) e temperatura (T) sono conservati. A volte è anche chiamato dinamica molecolare a temperatura costante (CTMD). In NVT, l'energia dei processi endotermici ed esotermici viene scambiata con un termostato.

Sono disponibili una varietà di algoritmi di termostato per aggiungere e rimuovere energia dai confini di una simulazione MD in modo più o meno realistico, approssimando l' insieme canonico . I metodi popolari per controllare la temperatura includono il ridimensionamento della velocità, il termostato Nosé-Hoover, le catene Nosé-Hoover, il termostato Berendsen , il termostato Andersen e la dinamica di Langevin . Il termostato Berendsen potrebbe introdurre l' effetto cubo di ghiaccio volante , che porta a traslazioni e rotazioni non fisiche del sistema simulato.

Non è banale ottenere una distribuzione canonica di insieme di conformazioni e velocità utilizzando questi algoritmi. Come questo dipenda dalle dimensioni del sistema, dalla scelta del termostato, dai parametri del termostato, dal time step e dall'integratore è oggetto di molti articoli nel settore.

Insieme isotermico-isobarico (NPT)

Nel complesso isotermica-isobarico , quantità di sostanza (N), pressione (P) e temperatura (T) sono conservati. Oltre a un termostato, è necessario un barostato. Corrisponde maggiormente alle condizioni di laboratorio con un pallone aperto a temperatura e pressione ambiente.

Nella simulazione delle membrane biologiche, il controllo della pressione isotropa non è appropriato. Per i doppi strati lipidici, il controllo della pressione avviene sotto un'area di membrana costante (NPAT) o una tensione superficiale costante "gamma" (NPγT).

Insiemi generalizzati

Il metodo di scambio di repliche è un insieme generalizzato. È stato originariamente creato per affrontare le dinamiche lente dei sistemi di spin disordinati. È anche chiamato rinvenimento parallelo. La formulazione del replica exchange MD (REMD) cerca di superare il problema dei minimi multipli scambiando la temperatura delle repliche non interagenti del sistema in esecuzione a diverse temperature.

Potenzialità nelle simulazioni MD

Una simulazione di dinamica molecolare richiede la definizione di una funzione potenziale , o una descrizione dei termini con cui le particelle nella simulazione interagiranno. In chimica e biologia questo è solitamente indicato come un campo di forza e nella fisica dei materiali come potenziale interatomico . I potenziali possono essere definiti a molti livelli di accuratezza fisica; quelli più comunemente usati in chimica si basano sulla meccanica molecolare e incarnano un trattamento meccanico classico delle interazioni particella-particella che può riprodurre cambiamenti strutturali e conformazionali ma solitamente non può riprodurre reazioni chimiche .

La riduzione da una descrizione completamente quantistica a un potenziale classico comporta due principali approssimazioni. La prima è l' approssimazione di Born-Oppenheimer , che afferma che la dinamica degli elettroni è così veloce che si può considerare che reagiscano istantaneamente al movimento dei loro nuclei. Di conseguenza, possono essere trattati separatamente. La seconda tratta i nuclei, che sono molto più pesanti degli elettroni, come particelle puntiformi che seguono la classica dinamica newtoniana. Nella dinamica molecolare classica, l'effetto degli elettroni è approssimato come una superficie di energia potenziale, che di solito rappresenta lo stato fondamentale.

Quando sono necessari livelli di dettaglio più fini, vengono utilizzati potenziali basati sulla meccanica quantistica ; alcuni metodi tentano di creare potenziali ibridi classico/quantistico in cui la maggior parte del sistema viene trattata in modo classico ma una piccola regione viene trattata come un sistema quantistico, che di solito subisce una trasformazione chimica.

Potenziali empirici

I potenziali empirici usati in chimica sono spesso chiamati campi di forza , mentre quelli usati nella fisica dei materiali sono chiamati potenziali interatomici .

La maggior parte dei campi di forza in chimica sono empirici e consistono in una somma di forze legate associate a legami chimici , angoli di legame e diedri di legame e forze non legate associate a forze di van der Waals e carica elettrostatica . I potenziali empirici rappresentano effetti quantomeccanici in modo limitato attraverso approssimazioni funzionali ad hoc. Questi potenziali contengono parametri liberi come la carica atomica , i parametri di van der Waals che riflettono le stime del raggio atomico e la lunghezza del legame di equilibrio , l'angolo e il diedro; questi sono ottenuti confrontando calcoli elettronici dettagliati (simulazioni chimiche quantistiche) o proprietà fisiche sperimentali come costanti elastiche , parametri reticolari e misurazioni spettroscopiche .

A causa della natura non locale delle interazioni non legate, coinvolgono almeno interazioni deboli tra tutte le particelle del sistema. Il suo calcolo è normalmente il collo di bottiglia nella velocità delle simulazioni MD. Per ridurre il costo computazionale, i campi di forza utilizzano approssimazioni numeriche come raggi di taglio spostati, algoritmi del campo di reazione , sommatoria Ewald mesh di particelle o la più recente mesh particella-particella-particella ( P3M ).

I campi di forza chimici impiegano comunemente disposizioni di legame preimpostate (un'eccezione è la dinamica ab initio ) e quindi non sono in grado di modellare esplicitamente il processo di rottura del legame chimico e le reazioni. D'altra parte, molti dei potenziali utilizzati in fisica, come quelli basati sul formalismo dell'ordine di legame, possono descrivere diverse coordinazioni di un sistema e la rottura del legame. Esempi di tali potenziali includono il potenziale del Brennero per gli idrocarburi e i suoi ulteriori sviluppi per i sistemi C-Si-H e COH. Il potenziale ReaxFF può essere considerato un ibrido completamente reattivo tra potenziali dell'ordine di legame e campi di forza chimica.

Potenziali di coppia contro potenziali a molti corpi

Le funzioni potenziali che rappresentano l'energia non legata sono formulate come una somma sulle interazioni tra le particelle del sistema. La scelta più semplice, impiegata in molti campi di forza popolari , è il "potenziale di coppia", in cui l'energia potenziale totale può essere calcolata dalla somma dei contributi energetici tra coppie di atomi. Pertanto, questi campi di forza sono anche chiamati "campi di forza additivi". Un esempio di tale potenziale di coppia è il potenziale di Lennard-Jones non legato (chiamato anche potenziale 6-12), utilizzato per calcolare le forze di van der Waals.

Un altro esempio è il modello Born (ionico) del reticolo ionico. Il primo termine nella prossima equazione è la legge di Coulomb per una coppia di ioni, il secondo termine è la repulsione a corto raggio spiegata dal principio di esclusione di Pauli e il termine finale è il termine di interazione di dispersione. Di solito, una simulazione include solo il termine dipolare, sebbene a volte sia incluso anche il termine quadrupolare. Quando n l = 6, questo potenziale è anche chiamato potenziale di Coulomb-Buckingham .

Nei potenziali a molti corpi , l'energia potenziale include gli effetti di tre o più particelle che interagiscono tra loro. Nelle simulazioni con potenziali a coppie, esistono anche interazioni globali nel sistema, ma avvengono solo attraverso termini a coppie. Nei potenziali a molti corpi, l'energia potenziale non può essere trovata da una somma su coppie di atomi, poiché queste interazioni sono calcolate esplicitamente come una combinazione di termini di ordine superiore. Dal punto di vista statistico, la dipendenza tra le variabili non può in generale essere espressa utilizzando solo prodotti a coppie dei gradi di libertà. Ad esempio, il potenziale di Tersoff , originariamente utilizzato per simulare carbonio , silicio e germanio , e da allora utilizzato per un'ampia gamma di altri materiali, comporta una somma su gruppi di tre atomi, con gli angoli tra gli atomi essendo un fattore importante nel potenziale. Altri esempi sono il metodo dell'atomo incorporato (EAM), l'EDIP e i potenziali Tight-Binding Second Moment Approximation (TBSMA), in cui la densità elettronica degli stati nella regione di un atomo viene calcolata dalla somma dei contributi degli atomi circostanti , e il potenziale contributo energetico è quindi una funzione di questa somma.

Potenziali semi-empirici

I potenziali semi-empirici fanno uso della rappresentazione matriciale della meccanica quantistica. Tuttavia, i valori degli elementi della matrice si trovano attraverso formule empiriche che stimano il grado di sovrapposizione di specifici orbitali atomici. La matrice viene quindi diagonalizzata per determinare l'occupazione dei diversi orbitali atomici e vengono utilizzate nuovamente formule empiriche per determinare i contributi energetici degli orbitali.

Esiste un'ampia varietà di potenziali semi-empirici, chiamati potenziali di legame stretto , che variano a seconda degli atomi che vengono modellati.

Potenziali polarizzabili

La maggior parte dei campi di forza classici include implicitamente l'effetto della polarizzabilità , ad esempio aumentando le cariche parziali ottenute dai calcoli della chimica quantistica. Queste cariche parziali sono stazionarie rispetto alla massa dell'atomo. Ma le simulazioni di dinamica molecolare possono modellare esplicitamente la polarizzabilità con l'introduzione di dipoli indotti attraverso metodi diversi, come le particelle di Drude o le cariche fluttuanti. Ciò consente una ridistribuzione dinamica della carica tra gli atomi che risponde all'ambiente chimico locale.

Per molti anni, le simulazioni MD polarizzabili sono state propagandate come la prossima generazione. Per liquidi omogenei come l'acqua, è stata ottenuta una maggiore precisione attraverso l'inclusione della polarizzabilità. Alcuni risultati promettenti sono stati ottenuti anche per le proteine. Tuttavia, è ancora incerto come approssimare al meglio la polarizzabilità in una simulazione.

Potenzialità nei metodi ab initio

Nella dinamica molecolare classica, una superficie di energia potenziale (solitamente lo stato fondamentale) è rappresentata nel campo di forza. Questa è una conseguenza dell'approssimazione di Born-Oppenheimer . Negli stati eccitati, nelle reazioni chimiche o quando è necessaria una rappresentazione più accurata, il comportamento elettronico può essere ottenuto dai primi principi utilizzando un metodo quantomeccanico, come la teoria del funzionale della densità . Questo è chiamato Ab Initio Molecular Dynamics (AIMD). A causa del costo del trattamento dei gradi di libertà elettronici, il costo computazionale di queste simulazioni è di gran lunga superiore alla dinamica molecolare classica. Ciò implica che AIMD è limitato a sistemi più piccoli e tempi più brevi.

I metodi quantistici e chimici ab initio possono essere utilizzati per calcolareal volo l' energia potenziale di un sistema, come necessario per le conformazioni in una traiettoria. Questo calcolo viene solitamente effettuato nelle immediate vicinanze della coordinata di reazione . Sebbene possano essere utilizzate varie approssimazioni, queste si basano su considerazioni teoriche, non su adattamenti empirici. I calcoli ab initio producono una grande quantità di informazioni che non sono disponibili da metodi empirici, come la densità degli stati elettronici o altre proprietà elettroniche. Un vantaggio significativo dell'utilizzo dimetodi ab initio è la capacità di studiare reazioni che comportano la rottura o la formazione di legami covalenti, che corrispondono a più stati elettronici. Inoltre, imetodi ab initio consentono anche di recuperare effetti oltre l' approssimazione di Born-Oppenheimer utilizzando approcci come la dinamica mista quanto-classica .

QM/MM ibrida

I metodi QM (quantum-meccanici) sono molto potenti. Tuttavia, sono computazionalmente costosi, mentre i metodi MM (meccanica classica o molecolare) sono veloci ma soffrono di diversi limiti (richiedono parametrizzazioni estese; le stime di energia ottenute non sono molto accurate; non possono essere utilizzati per simulare reazioni in cui i legami covalenti sono rotti/formati ; e sono limitate nelle loro capacità di fornire dettagli accurati sull'ambiente chimico). È emersa una nuova classe di metodi che combina i vantaggi dei calcoli QM (precisione) e MM (velocità). Questi metodi sono chiamati metodi misti o ibridi quantomeccanici e di meccanica molecolare (ibridi QM/MM).

Il vantaggio più importante del metodo ibrido QM/MM è la velocità. Il costo di fare dinamica molecolare classica (MM) nel caso più semplice bilance O (n 2 ), dove n è il numero di atomi nel sistema. Ciò è dovuto principalmente al termine interazioni elettrostatiche (ogni particella interagisce con ogni altra particella). Tuttavia, l'uso di radio cutoff, aggiornamenti pair-list periodiche e più recentemente le variazioni della particella a maglie di Ewald metodo (PME) ha ridotto questo a tra O (n) a O (n 2 ). In altre parole, se si simulasse un sistema con il doppio degli atomi, sarebbe necessaria una potenza di calcolo da due a quattro volte superiore. D'altra parte, i calcoli ab initio più semplici tipicamente scalano O(n 3 ) o peggio ( è stato suggerito che i calcoli ristretti di Hartree-Fock scalano ~ O(n 2.7 )). Per superare il limite, una piccola parte del sistema viene trattata in modo quantistico (tipicamente il sito attivo di un enzima) e il sistema rimanente viene trattato in modo classico.

In implementazioni più sofisticate, esistono metodi QM/MM per trattare sia i nuclei leggeri suscettibili di effetti quantistici (come gli idrogeni) che gli stati elettronici. Ciò consente di generare funzioni d'onda dell'idrogeno (simili alle funzioni d'onda elettroniche). Questa metodologia è stata utile nello studio di fenomeni come il tunneling dell'idrogeno. Un esempio in cui i metodi QM/MM hanno fornito nuove scoperte è il calcolo del trasferimento di idruro nell'enzima alcol deidrogenasi epatica . In questo caso, il tunneling quantistico è importante per l'idrogeno, poiché determina la velocità di reazione.

Rappresentazioni a grana grossa e ridotte

All'altra estremità della scala di dettaglio ci sono modelli a grana grossa e reticolari. Invece di rappresentare esplicitamente ogni atomo del sistema, si usano "pseudo-atomi" per rappresentare gruppi di atomi. Le simulazioni MD su sistemi molto grandi possono richiedere risorse informatiche così grandi da non poter essere facilmente studiate con i metodi tradizionali basati su tutti gli atomi. Allo stesso modo, le simulazioni di processi su scale temporali lunghe (oltre 1 microsecondo circa) sono proibitive, perché richiedono molti passaggi temporali. In questi casi, a volte si può affrontare il problema utilizzando rappresentazioni ridotte, chiamate anche modelli a grana grossa .

Esempi di metodi a grana grossa (CG) sono la dinamica molecolare discontinua (CG-DMD) e i modelli Go. La grana grossa viene eseguita a volte prendendo pseudoatomi più grandi. Tali approssimazioni atomiche unite sono state utilizzate nelle simulazioni MD di membrane biologiche. L'implementazione di tale approccio su sistemi in cui le proprietà elettriche sono di interesse può essere impegnativa a causa della difficoltà di utilizzare una corretta distribuzione di carica sugli pseudo-atomi. Le code alifatiche dei lipidi sono rappresentate da pochi pseudo-atomi, raccogliendo da 2 a 4 gruppi metilenici in ogni pseudo-atomo.

La parametrizzazione di questi modelli a grana molto grossa deve essere eseguita empiricamente, abbinando il comportamento del modello a dati sperimentali appropriati o simulazioni di tutti gli atomi. Idealmente, questi parametri dovrebbero tenere conto in modo implicito dei contributi sia entalpici che entropici all'energia libera. Quando la grana grossa viene eseguita a livelli più elevati, l'accuratezza della descrizione dinamica potrebbe essere meno affidabile. Ma modelli a grana molto grossa sono stati utilizzati con successo per esaminare un'ampia gamma di questioni di biologia strutturale, organizzazione dei cristalli liquidi e vetri polimerici.

Esempi di applicazioni di grana grossa:

La forma più semplice di grana grossa è l' atomo unito (a volte chiamato atomo esteso ) ed è stato utilizzato nella maggior parte delle prime simulazioni MD di proteine, lipidi e acidi nucleici. Ad esempio, invece di trattare esplicitamente tutti e quattro gli atomi di un gruppo metilico CH 3 (o tutti e tre gli atomi del gruppo metilenico CH 2 ), si rappresenta l'intero gruppo con uno pseudo-atomo. Deve, ovviamente, essere adeguatamente parametrizzato in modo che le sue interazioni di van der Waals con altri gruppi abbiano la giusta dipendenza dalla distanza. Considerazioni simili si applicano ai legami, agli angoli e alle torsioni a cui partecipa lo pseudo-atomo. In questo tipo di rappresentazione atomica unita, in genere si eliminano tutti gli atomi di idrogeno espliciti tranne quelli che hanno la capacità di partecipare ai legami idrogeno ( idrogeni polari ). Un esempio di questo è il campo di forza CHARMM 19.

Gli idrogeni polari vengono solitamente mantenuti nel modello, poiché un trattamento adeguato dei legami idrogeno richiede una descrizione ragionevolmente accurata della direzionalità e delle interazioni elettrostatiche tra i gruppi donatore e accettore. Un gruppo ossidrile, ad esempio, può essere sia un donatore di legami idrogeno, sia un accettore di legami idrogeno, e sarebbe impossibile trattarlo con uno pseudo-atomo OH. Circa la metà degli atomi in una proteina o in un acido nucleico sono idrogeni non polari, quindi l'uso di atomi uniti può fornire un notevole risparmio di tempo per il computer.

Incorporando effetti solventi

In molte simulazioni di un sistema soluto-solvente l'attenzione principale è sul comportamento del soluto con scarso interesse del comportamento del solvente in particolare in quelle molecole di solvente che risiedono in regioni lontane dalla molecola di soluto. I solventi possono influenzare il comportamento dinamico dei soluti tramite collisioni casuali e imponendo una resistenza per attrito sul movimento del soluto attraverso il solvente. L'uso di condizioni al contorno periodiche non rettangolari, confini stocastici e gusci di solvente può aiutare a ridurre il numero di molecole di solvente necessarie e consentire di dedicare una parte maggiore del tempo di calcolo alla simulazione del soluto. È anche possibile incorporare gli effetti di un solvente senza che siano presenti molecole di solvente esplicite. Un esempio di questo approccio è l'uso di una forza media potenziale (PMF) che descrive come cambia l'energia libera al variare di una particolare coordinata. La variazione di energia libera descritta da PMF contiene gli effetti medi del solvente.

Forze a lungo raggio

Un'interazione a lungo raggio è un'interazione in cui l'interazione spaziale cade non più velocemente di dove è la dimensionalità del sistema. Gli esempi includono interazioni carica-carica tra ioni e interazioni dipolo-dipolo tra molecole. La modellazione di queste forze rappresenta una vera sfida in quanto sono significative su una distanza che può essere maggiore della metà della lunghezza della scatola con simulazioni di molte migliaia di particelle. Sebbene una soluzione sarebbe quella di aumentare significativamente la dimensione della lunghezza della scatola, questo approccio di forza bruta è tutt'altro che ideale in quanto la simulazione diventerebbe molto costosa dal punto di vista computazionale. Anche il troncamento sferico del potenziale è fuori questione poiché si può osservare un comportamento irrealistico quando la distanza è vicina alla distanza di taglio.

Dinamica molecolare guidata (SMD)

Le simulazioni di dinamica molecolare diretta (SMD), o simulazioni di sonde di forza, applicano forze a una proteina per manipolarne la struttura trascinandola lungo i gradi di libertà desiderati. Questi esperimenti possono essere usati per rivelare cambiamenti strutturali in una proteina a livello atomico. SMD viene spesso utilizzato per simulare eventi come lo spiegamento meccanico o lo stretching.

Esistono due protocolli tipici di SMD: uno in cui la velocità di trazione è mantenuta costante e uno in cui la forza applicata è costante. Tipicamente, parte del sistema studiato (es. un atomo in una proteina) è trattenuto da un potenziale armonico. Le forze vengono quindi applicate ad atomi specifici a velocità costante o forza costante. Il campionamento ad ombrello viene utilizzato per spostare il sistema lungo la coordinata di reazione desiderata variando, ad esempio, le forze, le distanze e gli angoli manipolati nella simulazione. Attraverso il campionamento ad ombrello, tutte le configurazioni del sistema, sia ad alta che a bassa energia, vengono adeguatamente campionate. Quindi, la variazione di energia libera di ciascuna configurazione può essere calcolata come il potenziale della forza media . Un metodo popolare di calcolo del PMF è attraverso il metodo di analisi dell'istogramma pesato (WHAM), che analizza una serie di simulazioni di campionamento a ombrello.

Molte importanti applicazioni di SMD sono nel campo della scoperta di farmaci e delle scienze biomolecolari. Ad esempio, l'SMD è stato utilizzato per studiare la stabilità delle protofibrille di Alzheimer, per studiare l'interazione del ligando proteico nella chinasi 5 ciclina-dipendente e persino per mostrare l'effetto del campo elettrico sul complesso trombina (proteina) e aptamero (nucleotide) tra molti altri studi interessanti .

Esempi di applicazioni

Simulazione di dinamica molecolare di un motore molecolare sintetico composto da tre molecole in un nanoporo (diametro esterno 6,7 nm) a 250 K.

La dinamica molecolare è utilizzata in molti campi della scienza.

  • La prima simulazione MD di un processo di ripiegamento biologico semplificato è stata pubblicata nel 1975. La sua simulazione pubblicata su Nature ha aperto la strada alla vasta area del moderno ripiegamento computazionale delle proteine.
  • La prima simulazione MD di un processo biologico è stata pubblicata nel 1976. La sua simulazione pubblicata su Nature ha aperto la strada alla comprensione del movimento delle proteine ​​come essenziale in funzione e non solo accessorio.
  • La MD è il metodo standard per trattare le cascate di collisione nel regime dei picchi di calore, ovvero gli effetti che i neutroni energetici e l'irradiazione ionica hanno sui solidi e sulle superfici solide.

I seguenti esempi biofisici illustrano notevoli sforzi per produrre simulazioni di sistemi di dimensioni molto grandi (un virus completo) o tempi di simulazione molto lunghi (fino a 1,112 millisecondi):

  • Simulazione MD del virus completo del mosaico del tabacco satellitare (STMV) (2006, dimensione: 1 milione di atomi, tempo di simulazione: 50 ns, programma: NAMD ) Questo virus è un piccolo virus vegetale icosaedrico che peggiora i sintomi dell'infezione da virus del mosaico del tabacco (TMV). Sono state utilizzate simulazioni di dinamica molecolare per sondare i meccanismi dell'assemblaggio virale . L'intera particella STMV consiste di 60 copie identiche di una proteina che compongono il capside virale (rivestimento) e un genoma di RNA a singolo filamento di 1063 nucleotidi . Una scoperta chiave è che il capside è molto instabile quando non c'è RNA all'interno. La simulazione richiederebbe un computer desktop del 2006 circa 35 anni per essere completata. È stato così fatto in molti processori in parallelo con una comunicazione continua tra di loro.
  • Simulazioni di piegatura del copricapo di Villin in dettaglio all-atomo (2006, dimensione: 20.000 atomi; tempo di simulazione: 500 μs = 500.000 ns, programma: Folding@home ) Questa simulazione è stata eseguita in 200.000 CPU dei personal computer partecipanti in tutto il mondo. Questi computer avevano installato il programma Folding@home, uno sforzo di calcolo distribuito su larga scala coordinato da Vijay Pande della Stanford University. Le proprietà cinetiche della proteina Villin Headpiece sono state sondate utilizzando molte traiettorie brevi e indipendenti gestite dalla CPU senza una comunicazione continua in tempo reale. Un metodo impiegato era l'analisi del valore Pfold, che misura la probabilità di piegarsi prima dello spiegamento di una specifica conformazione iniziale. Pfold fornisce informazioni sulle strutture dello stato di transizione e sull'ordinamento delle conformazioni lungo il percorso di piegatura . Ogni traiettoria in un calcolo Pfold può essere relativamente breve, ma sono necessarie molte traiettorie indipendenti.
  • Sono state eseguite lunghe simulazioni a traiettoria continua su Anton , un supercomputer a parallelismo massiccio progettato e costruito attorno a circuiti integrati specifici per le applicazioni (ASIC) e interconnessioni da DE Shaw Research . Il risultato più lungo pubblicato di una simulazione eseguita utilizzando Anton è una simulazione di 1,112 millisecondi di NTL9 a 355 K; è stata inoltre eseguita una seconda simulazione indipendente di 1,073 millisecondi di questa configurazione (e molte altre simulazioni di oltre 250 μs di tempo chimico continuo). In How Fast-Folding Proteins Fold , i ricercatori Kresten Lindorff-Larsen, Stefano Piana, Ron O. Dror e David E. Shaw discutono "i risultati delle simulazioni di dinamica molecolare a livello atomico, su periodi compresi tra 100 μs e 1 ms, che rivelano una serie di principi comuni alla base del ripiegamento di 12 proteine ​​strutturalmente diverse". L'esame di queste diverse lunghe traiettorie, rese possibili da hardware specializzato e personalizzato, consente loro di concludere che "Nella maggior parte dei casi, la piegatura segue un unico percorso dominante in cui gli elementi della struttura nativa appaiono in un ordine altamente correlato con la loro propensione a formarsi nel stato dispiegato". In uno studio separato, Anton è stato utilizzato per condurre una simulazione di 1,013 millisecondi della dinamica dello stato nativo dell'inibitore della tripsina pancreatica bovina (BPTI) a 300 K.

Un'altra importante applicazione del metodo MD beneficia della sua capacità di caratterizzazione tridimensionale e analisi dell'evoluzione microstrutturale su scala atomica.

  • Le simulazioni MD sono utilizzate nella caratterizzazione dell'evoluzione della dimensione del grano, ad esempio, quando si descrivono l'usura e l'attrito di materiali nanocristallini Al e Al(Zr). L'evoluzione delle dislocazioni e l'evoluzione della dimensione dei grani vengono analizzate durante il processo di attrito in questa simulazione. Poiché il metodo MD ha fornito tutte le informazioni sulla microstruttura, l'evoluzione della dimensione dei grani è stata calcolata in 3D utilizzando i metodi Polyhedral Template Matching, Grain Segmentation e Graph clustering. In tale simulazione, il metodo MD ha fornito una misurazione accurata della dimensione del grano. Utilizzando queste informazioni, sono state estratte, misurate e presentate le strutture dei grani reali. Rispetto al metodo tradizionale di utilizzo del SEM con una singola fetta bidimensionale del materiale, MD fornisce un modo tridimensionale e accurato per caratterizzare l'evoluzione microstrutturale su scala atomica.

Algoritmi di dinamica molecolare

Integratori

Algoritmi di interazione a corto raggio

Algoritmi di interazione a lungo raggio

Strategie di parallelizzazione

Dinamica molecolare ab-initio

Hardware specializzato per simulazioni MD

  • Anton : un supercomputer specializzato e massicciamente parallelo progettato per eseguire simulazioni MD
  • MDGRAPE – Un sistema per scopi speciali costruito per simulazioni di dinamica molecolare, in particolare per la previsione della struttura delle proteine

Scheda grafica come hardware per simulazioni MD

Guarda anche

Riferimenti

Riferimenti generali

  • MP Allen, DJ Tildesley (1989) Simulazione al computer di liquidi . La stampa dell'università di Oxford. ISBN  0-19-855645-4 .
  • JA McCammon, SC Harvey (1987) Dinamica di proteine ​​e acidi nucleici . Cambridge University Press. ISBN  0-521-30750-3 (rilegato).
  • DC Rapaport (1996) L'arte della simulazione della dinamica molecolare . ISBN  0-521-44561-2 .
  • M.Griebel ; S. Knapek; G. Zumbusch (2007). Simulazione numerica in dinamica molecolare . Berlino, Heidelberg: Springer. ISBN 978-3-540-68094-9.
  • Frenkel, Daan ; Smit, Berend (2002) [2001]. Capire la simulazione molecolare: dagli algoritmi alle applicazioni . San Diego: stampa accademica. ISBN 978-0-12-267351-1.
  • JM Haile (2001) Simulazione della dinamica molecolare: metodi elementari . ISBN  0-471-18439-X
  • RJ Sadus, Simulazione molecolare dei fluidi: teoria, algoritmi e orientamento agli oggetti , 2002, ISBN  0-444-51082-6
  • Oren M. Becker, Alexander D. Mackerell, Jr., Benoît Roux, Masakatsu Watanabe (2001) Biochimica e biofisica computazionali . Marcel Dekker. ISBN  0-8247-0455-X .
  • Andrew Leach (2001) Modellazione molecolare: principi e applicazioni . (2a edizione) Prentice Hall. ISBN  978-0-582-38210-7 .
  • Tamar Schlick (2002) Modellazione e simulazione molecolare . Springer. ISBN  0-387-95404-X .
  • William Graham Hoover (1991) Meccanica statistica computazionale , Elsevier, ISBN  0-444-88192-1 .
  • DJ Evans e GP Morriss (2008) Meccanica statistica dei liquidi non in equilibrio , seconda edizione, Cambridge University Press, ISBN  978-0-521-85791-8 .

link esterno