|
|
|
In cele ce urmeaza se vor scoate in evidenta trasaturile generale caracteristice celor patru campuri de forta utilizate de pachetul de programe HyperChem: MM+, AMBER, BIO+ si OPLS.
Ecuatiile celor patru campuri de forta sunt similare in ceea ce priveste tipurile de termeni pe care-i contin si anume termeni care se refera la legaturi, unghiuri de valenta, unghiuri diedre si interactii intre atomi nelegati. Exista cateva diferente in forma ecuatiilor si de care se va tine cont in alegerea unui camp sau a altuia.
In alegerea unui camp trebuie, in primul rand, avut in vedere ca rezultatele cele mai bune se obtin pentru molecule similare, de acelasi tip, cu cele care au fost folosite in parametrizarea metodei.
Reprezinta cea mai generala metoda a mecanicii moleculare, dezvoltata in primul rand pentru moleculele organice, fiind o extensie a campului de forta MM2 conceput initial de Allinger si colaboratorii (1977) si distribuit prin Quantum Chemistry Program Exchange (QCPE). Prezinta urmatoarele trasaturi particulare:
MM+ este un camp de forte care foloseste toti atomii si in acelasi timp este unic in ceea ce priveste tratarea legaturilor si unghiurilor. Atat functia care exprima energia de intindere a legaturilor chimice cat si cea referitoare la deformarea unghiurilor de valenta contin pe langa termeni patratici si termeni superiori. Aceste potentiale sunt considerate ca exprima mai bine miscarea armonica decat potentialul unui oscilator armonic.
(II.8)
(II.9)
Termenii patratici ai ecuatiilor (II.8) si (II.9)
difera de ecuatiile (II.1) si (II.3) doar prin factorul 
; coeficientii
numerici din fata sumelor sunt factori de transformare ai
unitatilor, rezultatele fiind exprimate in Kcal/mol.
MM+ contine de asemenea un termen, numit termen Urey-Bradley, care exprima cuplarea vibratiei de valenta cu deformarea unghiului de valenta.
Spre deosebire de alte campuri care evalueaza interactiile intre atomii separati prin cel putin trei legaturi. Campul MM+ include, prin acest termen interactiile 1 - 3 , interactii ce nu pot fi neglijate intr-o modelare moleculara de finete. Astfel cand unghiul de valenta se micsoreaza, cele doua legaturi care formeaza unghiul se alungesc pentru a micsora constrangerea, termenul Urey-Bradley (ecuatia II.10) luand in considerare aceste modificari structurale.
(II.10)
unde unghiul este definit de atomii i,j si k, atomul central fiind atomul k, iar ik si jk definesc cele doua legaturi.
In cazul gruparilor XH2, deci cand atomii i si j sunt atomi de hidrogen, cuplarea se neglijeaza.
MM+ introduce in locul modificarii unghiului de torsiune impropriu o schema de deformare in afara planului bazata pe aplicarea ecuatiei de deformare a unghiului de valenta (II.3) pentru calcularea componentelor acestei deformarii.
Pentru calculul interactiilor van der Waals, MM+ foloseste in locul potentialului Lennard-Jones, o combinatie dintre o repulsie exponentiala si o interatie atractiva de dispersie 1/R6; parametrii de baza sunt suma razelor atomice van der Waals rij si un parametru energetic eij care determina adancimea minimului energetic:
(II.11)
In definirea tipurilor de atomi MM+ nu atribuie sarcini atomilor si ca urmare nu calculeaza, in mod obisnuit, interactii electrostatice. Acestea sunt inlocuite prin interactii dipol-dipol intre momentele de dipol asociate legaturilor polare. Distanta intre momentele de dipol de legatura mi si mj este distanta rij dintre mijlocul legaturilor dupa cum este aratat in Figura II.7; tot in aceasta figura sunt indicate unghiurile ai aj si c care definesc pozitia relativa a celor doua momente.
Figura II.7 Definirea pozitiei relative a momentelor de dipol de legatura
Energia de interactie dipol-dipol este data de relatia:
(II.12)
unde e este constanta dielectrica presupusa a fi mai mare ca unitatea; chiar si in faza gazoasa exista o ecranare a interactiei dipolilor de catre restul moleculelor. Constanta 14,39418 este factorul de conversie in Kcal/mol.
AMBER se bazeaza pe un camp de forta conceput si dezvoltat pentru proteine si acizi nucleici de catre grupul de cercetare al lui Peter Kollman de la Universitatea din California, San Francisco.
AMBER foloseste ambele reprezentari: cea a tuturor atomilor si cea a atomului unit. Sunt disponibili o varietate mare de seturi de parametrii de la Amber 2 la Amber 96.
Termenii campului de forta AMBER sunt cei generali prezentati in paragraful II.1 cu functiile corespunzatoare. Particularitatile acestui camp au fost deja prezentate in partea generala:
introducerea in mod explicit a interactiilor prin legaturi dehidrogen (ecuatia II.6);
folosirea unor definitii standard AMBER a unghiurilor diedre improprii pentru proteine si acizi nucleici;
folosirea unui factor de scalare de 0,5 pentru interactiile 1 - 4 de tip van der Waals sau electrostatice (Tabelul II.1);
folosirea unei constante dielectrice efective constante cu factorul de scalare D = 1 cand moleculele de solvent sunt incluse in mod explicit sau o constanta dielectrica dependenta de distanta cand se incearca simularea efectului de solvent;
Sarcinile atomice sunt plasate pe atomi in functie de tipul atomului la construirea sistemului prin folosirea bazei de date sau din fisierele de date de tip PDB, sau se adauga manual folosind rezultatele unui calcul cuantic;
AMBER adauga atomilor de sulf perechile de electroni neparticipanti si calculeaza interactiile acestora cu alti atomi ca si cum acestea ar fi un tip specific de atomi.
OPLS este un camp de forta dezvoltat de grupul de cercetare al lui Bill Jorgensen, Universitatea Yale. OPLS a fost conceput si dezvoltat, intocmai ca si AMBER, pentru a corespunde cerintelor impuse de calculele efectuate asupra proteinelor si acizilor nucleici. Acest camp trateaza interactiilor intre atomii nelegati cu o mai mare acuratete, printr-o parametrizare aleasa in urma unor studii extensive prin simulari Monte Carlo asupra lichidelor moleculelor mici. Prin adaugarea acestor interactii la interactiile AMBER asupra atomilor legati s-a urmarit obtinerea unui camp de forta care sa conduca la rezultate mult mai bune privind simularile in care moleculele de solvent sunt incluse in mod explicit si unde interactiile intre atomii nelegati joaca un rol important.
Camp de forta OPLS utilizeaza conceptul de atom unit pentru majoritatea atomilor de hidrogen atasati de atomul de carbon pentru a permite marirea vitezei de calcul asupra sistemelor macromoleculare. Desi initial OPLS facea apel doar la aceasta reprezentare, ulterior au fost adaugati si parametrii tuturor atomilor.
Termenii energetici corespunzatori vibratiei de valenta, deformarii unghiului de legatura, modificarii unghiului diedru si unghiului diedru impropriu sunt cei utilizati in AMBER.
Spre deosebire de AMBER in OPLS nu sunt incluse legaturile de hidrogen si perechile de electroni neparticipanti. Regulile de combinare pentru obtinerea parametrilor necesari interactiilor van der Waals sunt diferite.
Interactiile 1 - 4 van der Waals sunt inmultite cu un factor de scalare de 0,125
Factorul de scalare pentru interactiile 1 - 4 electrostatice este de 0,5; deoarece interactiile intre atomii nelegati au fost parametrizate din studiile in lichid cu molecule de apa incluse in mod explicit, dependenta constantei dielectrice de distanta este exclusa; factorul de scalare D al constantei dielectrice este egal cu unitatea.
Campul de forta BIO+ este o implementare a campului CHARMM (Chemistry at HARvard Macromolecular Mechanics) in pachetul de programe HyperChem , dezvoltat de grupul de cercetare al lui Martin Karplus de la Universitatea Harvard.
Asemanator cu AMBER si OPLS, CHARMM a fost conceput initial in scopul efectuarii de calcule asupra macromoleculelor de interes biologic, in special asupra proteinelor, fiind un camp de forta ce face apel la reprezentarea atomului unit. Dezvoltarile recente ale parametrizarilor nu au mai fost publicate, fiind accesibile doar prin cumpararea programului CHARMM;
Parametrizarea grupului Karplus foloseste reprezentarea tuturor atomilor. Deoarece au fost publicate doar seturi limitate de parametrii, sunt disponibili parametrii doar pentru tipurile de atomi ale unui subset de aminoacizi, astfel ca multe calcule asupra unor astfel de sisteme folosind BIO+ nu pot fi efectuate si esueaza din lipsa parametrilor;
Termenii energetici corespunzatori vibratiei de valenta, deformarii unghiului de legatura, modificarii unghiului diedru si unghiului diedru impropriu sunt cei prezentati in partea generala, avand insa o parametrizare si reguli de compunere proprii;
BIO+ nu include legaturile de hidrogen si perechile de electroni neparticipanti;
Pentru interactiile 1 - 4 van der Waals nu este introdus nici un factor de scalare, dar pentru astfel de interactii sunt folositi parametrii specifici interactiilor 1 - 4;
Pentru interactiile electrostatice sarcinile atomice sunt plasate pe atomi in functie de tipul atomului la construirea sistemului prin folosirea fie a bazei de date fie a fisierelor de date de tip PDB, sau se adauga manual folosind rezultatele unui calcul cuantic. Pentru constanta dielectrica se foloseste fie un factor D = 2,5 sa D = 1 asociat cu o forma constanta pentru constanta dielectrica, atunci cand se iau explicit in calcul molecule de solvent, fie o forma dependenta de distanta in cazul simularii efectului de solvent;
Interactiile 1 - 4 electrostatice sunt factorizate cu factorul 0.5;
Calculul proprietatilor moleculare necesita, in primul rand, o structura moleculara bine definita. Pentru starea fundamentala este necesara o geometrie corespunzatoarea minimului de energie.
Obtinerea unei geometrii optimizate conduce la o noua structura stabila, aceasta reprezentand punctul de plecare in calculele cuantice care permit obtinerea unui set larg de proprietati structurale si electronice. Totodata obtinerea unor astfel de structuri reprezinta o prima etapa, obligatorie, in vederea unei simulari prin metodele dinamicii moleculare.
Determinarea minimului energetic ar fi o problema simpla daca s-ar cunoaste suprafata de energie potentiala. Pentru o molecula neliniara cu N atomi acest lucru este practic imposibil fiind vorba de o hipersuprafata n dimensionala, n fiind egal cu 3N-6 (pentru N atomi numarul gradelor de libertate este 3N din care se scad trei grade de libertate pentru translatie si alte trei pentru rotatie). Aceasta hipersuprafata corespunde unui hiperspatiu cu n+1 dimensiuni.
Metodele folosite pentru optimizarea geometriei moleculare vor fi metode de cautare a unui minim de energie pe o suprafata de energie potentiala necunoscuta. Aceasta inseamna gasirea acelor coordonate carteziene pentru care toate componentele vectorului gradient se anuleaza:
(II.13)
Pentru o suprafata monodimensionala (o curba, Figura II.8), cum ar fi cea corespunzatoare unei molecule diatomice ce poseda 3N-5 = 1 grade de libertate, gasirea configuratiei stabile, de echilibru, corespunzatoare minimului de energie, este simpla. Se pleaca de la o structura initiala, corespunzatoare de exemplu punctului a sau b si in urma unui proces iterativ de cautare a minimului se gaseste punctul c corespunzator minimului de energie.
Figura II.8 Suprafata monodimensionala de energie potentiala
corespunzatoare moleculei biatomice
In cazul moleculelor poliatomice suprafata de energie potentiala nu mai este atat de simpla si, in general, va prezenta mai multe minime, dupa cum se poate observa din sectiunea printr-o asemenea suprafata prezentata in Figura II.9.
Figura II.9.O sectiune printr-o suprafata multidimensionala ce
prezinta mai multe minime.
Aceasta suprafata de energie potentiala prezinta cinci minime dintre care minimul D este minimul global.
Metodele obisnuite de cautare pe o suprafata nu permit escaladarea sau penetrarea barierelor de potential si, din aceasta cauza, in cazul in care se pleaca initial din zonele a,b sau c, este posibila depistarea doar a unor minime locale A,B sau C. Minimul global D poate fi obtinut numai daca se pleaca din zona d sau d
O alta problema particulara de care trebuie avut grija in aplicarea procedurile de optimizare este aceea a unor structuri initiale, de plecare, riguros plane. In astfel de cazuri, metodele de optimizare, facand apel la directiile fortelor care actioneaza asupra atomilor si neexistand initial forte care sa actioneze in afara planului, vor conduce intotdeauna la structuri moleculare plane. Aceste structuri pot fi irelevante. De aceea pentru a evita astfel de situatii se scoate oricare dintre atomi in afara planului. Numai in acest fel se poate obtine, daca este cazul, un minim energetic corespunzator unei structuri neplanare.
Dintre metodele de cautare a minimului energetic cele mai cunoscute sunt: cea care foloseste algoritmul care urmeaza calea corespunzatoare celei mai abrupte pante (steepest descent method), metoda gradientului conjugat (conjugate gradient method) si metodele Newton-Raphson care fac apel si la derivatele de ordinul doi ale functiei obiectiv.
1. Metoda steepest descent
Aceasta metoda este o metoda de cautare de ordinul intai folosind derivatele de ordinul intai ale energiei potentiale in raport cu coordonatele carteziene. Metoda urmeaza calea de coborare corespunzatoare celei mai pronuntate pante, deci dupa directia vectorului gradient cel mai negativ. Coborarea este acompaniata de adaugarea la coordonate a unor incremente, in directia data de cel mai negativ gradient al energiei potentiale.
Metoda steepest descent atenueaza rapid fortele care actioneaza asupra atomilor. De aceea este folositoare pentru eliminarea interactiilor pronuntate intre atomii nelegati care apar, de regula, in structurile de la care se pleaca initial.
Un avantaj al metodei este timpul minim de calcul necesar fiecarui pas, dar, in acelasi timp, marele dezavantaj este acela al unei convergente lente catre minimul energetic.
2. Metoda gradientului conjugat
Este de asemenea o metoda de ordinul intai, dar spre deosebire de precedenta, algoritmul de cautare al minimului energetic foloseste atat gradientul curent cat si informatia obtinuta privind directia de cercetare anterioara. Cu alte cuvinte, in timp ce metoda steepest descent calculeaza de fiecare data derivatele de ordinul intai fara a acorda o atentie speciala modului in care aceste derivate se modifica, metoda gradientului conjugat utilizeaza informatia privind evolutia anterioara a derivatelor de ordinul intai in scopul de a obtine implicit informatii privind derivatele de ordinul doi. Aceste informatii ar putea fi obtinute explicit numai prin stocarea derivatelor de ordinul al doilea in raport cu coordonatele, si care sunt elemente de matrice ale asa numitului Hessian.
Metoda gradientului conjugat ofera o tehnica mai buna si mai rapida de cautare.
Exista doua variante ale metodei: Fletcher-Reeves si Polak-Ribiere
Metoda Polak-Ribiere este o metoda mai rafinata si din aceasta cauza este aleasa ca metoda implicita de cautare de catre HyperChem.
Optimizarea structurilor moleculelor mari poate dura un timp indelungat; trebuie avut in vedere ca, in timp ce numarul de cicluri cerute de catre metoda gradientului conjugat este aproximativ egal cu numarul N de atomi, timpul necesar unui ciclu este proportional cu patratul acestui numar.
Nu este recomandata folosirea acestei metode in situatia in care, initial, fortele interatomice sunt mari, deci se pleaca cu o structura moleculara foarte indepartata de minim.
3. Metodele Newton-Raphson
Aceste metode fac apel atat la derivatele de ordinul intai, cat si la cele de ordinul al doilea in raport cu coordonatele carteziene si care formeaza, dupa cum am amintit, Hessianul. In felul acesta se obtin informatii atat asupra pantei cat si al curburii suprafetei de energie potentiala
Metoda generala Newton-Raphson calculeaza integral Hessianul si in functie de acesta calculeaza incrementele coordonatelor in directia de cercetare. Diferitele variante ale acestei metode difera prin modul de aproximare al hesianului.
In cadrul metodei MM+ este folosita o varianta simplificata a metodei Newton-Raphson (Block Diagonal Newton-Raphson method); Aceasta varianta retine din Hessian doar blocurile diagonale asociate atomilor individuali si neglijeaza blocurile nediagonale care ar corespunde derivatelor de ordinul doi in raport cu coordonatele a doi atomi. Practic procedura considera la un moment dat doar un singur atom, calculeaza Hessianul (o matrice 3 x 3 ) si cele trei componente ale gradientului pentru acest atom; sunt calculate noile coordonate ale atomului procedura repetandu-se pe rand pentru toti atomii.
Metoda utilizeaza informatiile oferite de derivatele de ordinul al doilea si prin aceasta poate fi mai eficienta decat metoda gradientului conjugat, fiind in acelasi timp si mai rapida. Totusi, datorita neglijarii blocurilor nediagonale, in unele cazuri metoda poate deveni oscilanta si convergenta nu poate fi atinsa. In aceste situatii desi mai lenta este recomandata metoda gradientului conjugat, care implicit implica Hessianul integral.
La fel ca si in cazul metodei gradientului conjugat, nu se recomanda folosirea acestei metode cand se pleaca de la structuri indepartate de minimul energetic.
4. Criterii de convergenta
In optimizarea structurilor moleculare se face apel in general la doua criterii de convergenta: radacina medie patratica (root mean square - RMS) a gradientului si numarul de cicluri folosite in procesul de optimizare.
RMS - ul gradienului este definit de relatia (II.14):
RMS-gradient = 
(II.14)
Este indicata terminarea procesului de optimizare pe baza RMS - lui gradientului, care trebuie sa fie cat mai mic posibil, deoarece numarul de cicluri necesare minimizarii depind pe de o parte de structura de la care se pleaca, dar si de metoda optimizare. Nu exista nici un criteriu de a prevedea numarul de cicluri necesar minimizarii. HyperChem calculeaza automat, dar arbitrar, numarul de cicluri in functie de numarul de atomi ca fiind 15 N. Este recomandata folosirea acestui criteriu numai cand se folosesc succesiv diferiti algoritmi de optimizare, de exemplu initial se utilizeaza metoda steepest descent, urmata apoi de metoda gradientului conjugat.
