Fabio Schoen

Ottimizzazione Globale
Fabio Schoen
[email protected]
http://globopt.dsi.unifi.it/users/schoen
Scuola CIRO 2002
Typeset with LATEX
Problemi generici di ottimizzazione globale
min f (x)
x∈S⊆Rn
Primo problema: cosa intendiamo per ”ottimizzazione globale”: in astratto
vorremmo determinare
f ∗ = min n f (x)
x∈S⊆R
e
x∗ = arg min f (x) : f (x∗ ) ≤ f (x) ∀ x ∈ S
Questa definizione non è soddisfacente:
• il problema è ”mal posto” nelle incognite x (due funzioni la cui distanza è
piccola possono avere ottimi globali arbitrariamente distanti)
• è invece ben posto nel valore dell’ottimo: ||f − g|| ≤ δ =⇒ |f ∗ − g ∗ | ≤ ε
1
Spesso ci accontenta quindi di determinare f ∗ ed, al massimo, una o più
soluzioni ammissibili x̄ tali che
f (x̄) ≤ f (x∗ ) + ε
Più frequentemente, questo è già un obiettivo troppo ambizioso
Situazione della ricerca in ottimizzazione globale
• problema molto rilevante dal punto di vista applicativo
• molto (forse troppo) difficile da risolvere
• un enorme numero di articoli dedicati ad algoritmi di ottimizzazione
globale per classi specifiche di problemi
• un numero relativamente piccolo di articoli con contenuti teorici rilevanti
• spesso i lavori più eleganti hanno prodotto algoritmi insoddisfacenti e,
viceversa, i migliori algoritmi si appoggiano ad una teoria poco raffinata.
• le riviste di classe A non accettano più sottomissioni di articoli di
ottimizzazione globale indipendentemente dal loro valore. Mathematical
Programming non ha (e non ritiene utile avere) editors associati nel
campo dell’ottimizzazione globale
• nascita di “ghetti” di pubblicazione: Computational Optimization and
Applications, Journal of Global Optimization, JOTA, pochi altri.
2
• molti lavori di ottimizzazione globale vengono pubblicati su riviste
applicative, molto autorevoli nel loro campo, non valutate da parte della
comunità di ottimizzazione
• Bazaraa, Sherali, Shetty “Nonlinear Programming: theory and
algorithms”, 1993:
la parola “ottimo globale” viene citata per la prima volta a pag. 99, la
seconda a pag. 132, poi, a pag. 247:
“A desirable property of an algorithm for solving [an optimization]
problem is that it generates a sequence of points converging to a global
optimal solution. In many cases however we may have to be satisfied with
less favorable outcomes.”
dopodiché mai più (su 638 pagine). Il termine “Global optimization” non
viene mai nominato.
Il problema di ottimizzazione globale è ”difficile” in senso lato:
In assenza di informazioni di tipo ”globale” sulla funzione nessun algoritmo può
determinare e certificare un ottimo globale, a meno che non venga effettuato
un campionamento denso nell’insieme ammissibile. La nozione di informazione
”globale” può essere data in modo rigoroso – esempi di informazioni globali:
• numero di ottimi locali
• valore dell’ottimo globale
• valore di una limitazione superiore sulla costante di Lipschitz (nel caso
della minimizzazione di funzioni lipschitziane su insiemi ammissibili
”semplici”, come ad es. iper-rettangoli)
• concavità della funzione, convessità dell’insieme ammissibile
• possibilità di esprimere la funzione f come differenza di due funzioni
convesse delle quali sia nota l’espressione analitica.
3
Il problema è difficile anche secondo la teoria della complessità
computazionale:
problema di programmazione quadratica:
min
l≤Ax≤u
1 T
x Qx + cT x
2
È N P –hard [Sahni, 1974] e, posto come problema di decisione, è
N P -completo [Vavasis, 1990].Sono N P –hard anche molti casi particolari:
• massimizzazione della norma su un parallelotopo
• minimizzazione quadratica su iper–rettangoli (A = I ) nel caso in cui anche
un solo autovalore di Q sia negativo
• minimizzazione su un simplesso:
1
min xT Qx + cT x
x≥0 2
xj =
1
j
• E’ hard anche verificare che un punto è un ottimo locale!
Alcune applicazioni dell’ottimizzazione globale
• Problemi di minimizzazione di funzionali di costo concavi (sconti per
grandi quantità – economie di scala)
• Problemi di minimizzazione in presenza di costi fissi
• Problemi di ottimizzazione discreta (programmazione lineare intera):
min cT x + KxT (1 − x)
Ax = b
x
∈
oppure:
min cT x
Ax = b
x
∈
[0, 1]
x (1 − x) = 0
T
4
[0, 1]
• altri problemi di ottimizzazione combinatoria (es: max-clique):
Motzkin - Strauss:
1
max xT Ax
2
1T x = 1
x ≥ 0
ha ottimo globale
1−
1
α(G)
dove α(G) è la dimensione massima delle cliques in G ed A la matrice di
adiacenza del grafo
• Problemi di minimizzazione di funzionali di costo non convessi nè concavi
- ad esempio, configurazione di minima energia di molecole complesse:
micro-clusters di Lennard-Jones, folding di proteine, docking
proteine/ligandi
Lennard-Jones: modello valido per clusters di Argon, Krypton, Oro,
Nickel: l’energia di una coppia di atomi posizionati in X1 , X2 ∈ R3 , a
distanza r l’uno dall’altro, è la somma di una componente attrattiva e di
una repulsiva:
1
2
v(r) = 12 − 6
r
r
L’energia di un cluster di N atomi i cui centri sono posti in
X1 , . . . , XN ∈ R3 è definita come:
v(||Xi − Xj ||)
i=1,...,N j<i
Questa funzione ha un numero di ottimi locali (non globali) che cresce
esponenzialmente con N
5
3
attractive(x)
repulsive(x)
lennard-jones(x)
2
1
0
-1
-2
-3
0.5
1
1.5
2
2.5
3
3.5
4
4.5
Problemi di folding e di docking di proteine
Modello di energia potenziale E = El + Ea + Ed + Ev + Ee dove:
El =
1
i∈L
2
Kib (ri − ri0 )2
(contributo dovuto alla lunghezza dei legami chimici):
Ea =
1
i∈A
2
Kiθ (θi − θi0 )2
(angolo fra i legami)
Ed =
1
i∈T
2
Kiφ [1 + cos(nφi − γ)]
(angoli diedrali)
6
5
Ev =
(i,j) ∈ C
(van der Waals)
Ee =
Aij
Bij
12 − R6
Rij
ij
1 qi qj
2
εRij
(i,j) ∈ C
(interazione elettrostatica)
Problema di docking: date due macromolecole M1 , M2 , determinare la
configurazione di “accoppiamento” di minima energia.
Se gli atomi delle due molecole non si legano fra loro, il contributo energetico
da minimizzare per la determinazione del docking ottimale è dato da
qi qj
Bij
Aij
1
Ev + E e =
− 6 +
12
Rij
Rij
2
εRij
i∈M1 ,j∈M2
i∈M1 ,j∈M2
7
Principali approcci algoritmici
Due famiglie di algoritmi:
1. in assenza di informazioni globali (problemi ”non strutturati”)
2. in presenza di informazioni globali
Per problemi strutturati si utilizzano sia algoritmi stocastici che deterministici.
Per problemi non strutturati, tipicamente si utilizzano metodi stocastici.
Ogni algoritmo di ottimizzazione globale deve bilanciare lo sforzo
computazionale tra
• esplorazione della regione ammissibile
• approssimazione dell’ottimo
Un esempio: minimizzazione dell’energia potenziale nei clusters di
Lennard-Jones:
LJN = min LJ(X) = min
N
−1
i=1
N
1
2
−
12
Xi − Xj Xi − Xj 6
j=i+1
E’ un problema strutturato!
Ma conviene usare la struttura??
E come??
8
La mappa
N (N −1)/2
F1 : R3N → R+
F1 (X1 , . . . , XN ) → X1 − X2 2 , . . . , XN −1 − XN 2
è convessa, e la funzione
N (N −1)/2
F2 : R+
R
→
1
1
F2 (r12 , . . . , rN −1,N ) →
−2
6
3
rij
rij
è la differenza di due funzioni convesse. Quindi, LJ(X) è esprimibile come
differenza di due funzioni convesse (d.c. programming)
NB: ogni funzione di classe C 2 è d.c. (ma spesso non se ne conosce la
rappresentazione)
Ottimizzazione globale d.c.: molto elegante, esiste una teoria della dualità,
molte proprietà teoriche rilevanti, ma algoritmi tipicamente inefficienti.
Un algoritmo primale per problemi d.c.
metodo “cutting plane”
Il problema d.c. non vincolato si può trasformare in un equivalente problema
con obiettivo lineare, un vincolo convesso, un vincolo “reverse-convex”:
min g(x) − h(x)
con g, h convesse, è equivalente a
min z
g(x) − h(x) ≤ z
che è equivalente a
min z
g(x) ≤ w
h(x) + z
9
≥ w
Forma canonica d.c.:
min cT x
g(x) ≤ 0
h(x) ≥ 0
con h, g convesse. Siano
Ω
=
{x : g(x) ≤ 0}
C
= {x : h(x)≤0}
Hp:
0 ∈ intΩ ∩ intC, cT x > 0∀x ∈ Ω \ intC
Proprietà fondamentale: se il problema d.c. ammette ottimo, ne ammette
almeno uno appartenente a
∂Ω ∩ ∂C
4
cT x = 0
3
2
1
C
Ω
0
-1
-2
-3
-4
-9
-8
-7
-6
-5
-4
-3
10
-2
-1
0
1
2
3
Sia x̄ la miglior soluzione ammissibile conosciuta.
Sia
D(x̄) = {x ∈ Ω : cT x ≤ cT x̄}
Se D(x̄) ⊆ C allora x̄ è ottimale;
Verifica: si costruisce un politopo P (con vertici noti) che contiene D(x̄) e si
verifica se qualche vertice di P non appartiene a C . Se tutti appartengono a C
⇒ x̄ è ottimale; altrimenti sia v il miglior vertice ammissibile; l’intersezione del
segmento [0, v] con la frontiera di C (se ammissibile) fornisce un punto
migliore x. Altrimenti si inserisce un taglio nel politopo P tangente a Ω in x.
4
D(x̄) = {x ∈ Ω : cT x ≤ cT x̄}
cT x = 0
3
2
1
C
Ω
0
-1
x̄
-2
-3
-4
-9
-8
-7
-6
-5
-4
-3
11
-2
-1
0
1
2
3
4
P politopo : D(x̄) ⊆ P con vertici V1 , . . . , Vk . V := arg max h(Vj )
cT x = 0
3
2
1
C
Ω
0
-1
x̄
-2
-3
V
-4
-9
4
-8
-7
-6
-5
-4
-3
-2
-1
0
1
2
3
2
3
xk = ∂C ∩ [V , 0]
cT x = 0
3
2
1
C
Ω
0
-1
xk
x̄
-2
-3
V
-4
-9
-8
-7
-6
-5
-4
-3
12
-2
-1
0
1
4
Se xk ∈ Ω, si pone x̄ := xk
cT x = 0
3
2
1
C
Ω
0
-1
x̄
-2
-3
-4
-9
4
-8
-7
-6
-5
-4
-3
-2
-1
0
1
2
3
2
3
Se invece xk ∈ Ω, si suddivide il politopo
cT x = 0
3
2
1
C
Ω
0
-1
-2
-3
-4
-9
-8
-7
-6
-5
-4
-3
13
-2
-1
0
1
Dualità per problemi d.c.
Si consideri
min g(x) − h(x)
x∈S
con f, g : convesse. Siano
h (u) := sup{uT x − h(x) : x ∈ Rn }
g (u) := sup{uT x − g(x) : x ∈ Rn }
le funzioni coniugate di h e g . Allora il problema
inf{h (u) − g (u) : u : h (u) < +∞}
è detto duale di Fenchel-Rockafellar . Se il problema min g(x) − h(x) ammette
ottimo, il duale di Fenchel è un duale forte.
Se x ∈ arg min g(x) − h(x) allora
u ∈ ∂h(x )
è ottima per il duale, e se u ∈ arg min h (u) − g (u) allora
x ∈ ∂g (u )
è soluzione ottima per il primale.
14
Algoritmo primale/duale
Pk : min g(x) − (h(xk ) + (x − xk )T yk )
e
Dk : min h (y) − (g (yk−1 ) + xTk (y − yk−1 )
α-BB
Schema Branch-and-Bound. Upper bound: miglior soluzione ammissibile
trovata (nel caso dei clusters di LJ, miglior configurazione molecolare nota).
Lower bound?
Rilassamento convesso - sottostima convessa della funzione obiettivo.
Idealmente, dato il problema
min f (x)
con f non convessa, se fosse nota
φ ∈ arg Sup{g convessa : g(x) ≤ f (x)
allora l’ottimo globale si otterrebbe come
min φ(x)
x
15
∀ x}
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0
2
4
6
8
10
Metodo α-BB: sottostima convessa di f (x) su [a, b]:
f (x) − α(x − a)(b − x)
o, in generale,
f (x) −
n
U
αj (xj − xL
j )(xj − xj )
j=1
Se α sufficientemente grande ⇒ sottostima convessa di f .
Dopo la valutazione del lower bound (minimizzazione sottostima convessa), si
sceglie la regione con il minimo lower bound e la si suddivide (branching).
16
2
0
-2
-4
-6
-8
-10
-12
-14
-16
0
0.5
1
1.5
2
1.5
2
1
0
-1
-2
-3
-4
-5
-6
0
0.5
1
17
1
0.5
0
-0.5
-1
-1.5
-2
-2.5
-3
-3.5
0
0.5
1
1.5
2
Il metodo è molto simile al classico metodo di Shubert per funzioni
Lipschitziane: se
|f (x) − f (y)| ≤ Lx − y ∀x, y ∈ S
ed L è una costante nota, si possono facilmente costruire minoranti lineari a
pezzi di f e basare su queste un algoritmo di ricerca:
f (y) − Lx − y ≤ f (x) ≤ f (y) + Lx − y
18
Se la funzione f è stata osservata in y1 , . . . , yn , vale
max f (yi ) − Lx − yi ≤ f (x) ≤ min f (yi ) + Lx − yi i=1,n
i=1,n
2
0
-2
-4
-6
-8
-10
0
0.5
1
19
1.5
2
2
0
-2
-4
-6
-8
-10
0
0.5
1
1.5
2
0
0.5
1
1.5
2
2
0
-2
-4
-6
-8
-10
20
2
1
0
-1
-2
-3
-4
-5
0
0.5
1
1.5
2
1
0.5
0
-0.5
-1
-1.5
-2
0
0.5
1
21
1.5
2
Metodi basati su “merit functions”
Metodo Bayesiano: si considera la funzione obiettivo come una realizzazione
di un processo stocastico:
f (x) = F (x; ω)
Si introduce una funzione di perdita, ad esempio
L(x1 , ..., xn ; ω) = min F (xi ; ω) − min F (x; ω)
x
i=1,n
e si sceglie il prossimo punto da campionare in modo da minimizzare la perdita
attesa:
xn+1
=
arg min E (L(x1 , ..., xn , xn+1 ) | x1 , ..., xn )
=
arg min E (min(F (xn+1 ; ω), F (x; ω)) | x1 , ..., xn )
Metodo radial basis
Date n osservazioni (x1 , f1 ), . . . , (xn , fn ), la prossima si sceglie nel minimo di
un’interpolante:
xn+1 ∈ arg min s(x) =
n
λi Φ(x − xi ) + p(x)
i=1
p è un polinomio di grado massimo prefissato m. Φ è una funzione radiale. Ad
esempio:
Φ(r)
Φ(r)
Φ(r)
Φ(r)
= r
= r
lineare
3
cubica
2
= r log r
= e−γr
2
thin plate spline
gaussiana
Il polinomio p compensa i casi in cui l’interpolazione radiale “pura” è
impossibile (per singolarità della matrice {Φij = Φ(xi − xj )})
22
Metodi stocastici
Pure Random Search - campionamento casuale (uniforme) sull’insieme
ammissibile
3
2
+
+
+
+
rs
rs
+
+
+
sr
rs
rs
sr
rs
0
rs
+
rs
sr
1
+
-1
+
-2
-3
0
1
2
3
4
5
Best Start - come Pure Random Search, ma con l’esecuzione di una ricerca
locale inizializzata dal ”record”
3
2
+
+
+
+
rs
rs
rs
rs
+
rs
+
rs
rs
rs
0
rs
+
rs
+
1
+
-1
+
-2
-3
0
1
2
3
23
4
5
Multistart - discese locali inizializzate da tutti i punti del un campione
3
2
+
+
+
+
rs
rs
rs
rs
+
rs
+
rs
rs
rs
0
rs
+
rs
+
1
+
-1
+
-2
-3
0
1
2
3
4
5
3
2
+
+
+
+
rs
rs
rs
rs
+
rs
+
rs
rs
rs
0
rs
+
rs
+
1
+
-1
+
-2
-3
0
1
2
3
24
4
5
Metodi di clustering
Campione uniforme
Trasformazione (concentrazione del campione) ad es., eliminando una frazione
dei punti ”peggiori” o eseguendo qualche passo di discesa
i punti rimanenti vengono raggruppati in clusters
da ciascun cluster viene fatta partire una ricerca locale
5
rs
sr
rs
rs
Campione uniforme:
rs
rs
0
1
rs
0
1
3
2
rs
sr srrs
rs rs sr
rs rs
rs
rs
rs
rs
rs
sr
rsrs
−1
4
rs
−5
rs
rs rs
rs
rs
rs
−3
0
2
3
4
5
25
5
rs
sr
rs
rs
Concentrazione del campione:
rs rs
rs
rs
rs
rs
−3
rs
rs
rs
−5
0
+
0
1
2
3
3
2
+
+ ++
++ +
++
+
++ +
rs
+ rs
+
−1
4
4
1
0
5
r
Clustering:
rr
r
r
r
r
rr
−5
r
5
−3
4
r
u
u3
r
r
0
2
−1
1
0
0
1
2
3
4
5
26
r
Ottimizzazione locale:
rr
r
r
r
r
rr
−5
r
5
−3
4
r
u
u3
r
r
0
2
−1
1
0
0
1
2
3
4
5
27
Clustering: MLSL
Dati i punti X1 , . . . , XN ∈ [0, 1]n , etichettare Xj come “clustered” se e solo se
∃ Y ∈ X1 , . . . , XN tale che
||Xj − Y || ≤ ∆N
1
:= √
2π
log N
n
σΓ 1 +
N
2
n1
e
f (Y ) ≤ f (Xj )
Esempio (N = 10, σ = 2, ∆N ≈ 0.8):
3
2
r
+
+r
r
r
r
+
r
0
r
+
+
r
+
+
r
+
1
r
+
-1
+
-2
-3
0
1
2
3
28
4
5
Simple Linkage
Viene generato un campione sequenziale e, da ciascun punto campionato,
viene iniziata una discesa locale a meno che non esista un altro punto del
campione ”sufficientemente vicino” e con valore funzionale ”strettamente
migliore”
3
2
+
+r
r
+
+
+
r
r
r
r
r
+
r
0
r
+
r
+
1
+
-1
+
-2
-3
0
1
2
3
29
4
5
Simulated annealing
Problema: trovare il punto più basso sulla terra.
Soluzione: ? Fossa delle Marianne? No! perforazione effettuata dai russi in
Siberia
Non è importante (solo) il valore dell’ottimo globale, ma anche la sua
raggiungibilità, o la sua regione di attrazione
Una molecola in un bagno di calore a temperatura T (Kelvin), fluttua tra varie
configurazioni. Occuperà la configurazione R con probabilità
p(R) = e−f (R)/kB T
e−f (x)/kB T dx
dove kB : costante di Boltzmann, T : temperatura, f (R) energia potenziale.
La quantità
β −1 = kB T
si dice energia termica.
30
Le molecole fluttuano, ad una data temperatura, entro regioni connesse dello
spazio degli stati corrispondenti a bacini metastabili a temperatura T .
Si può pensare a questi come a zone del grafico della funzione riempite ad una
profondità ≈ kB T .
Data una regione α, la quantità
−f (R)/k T
B
e
dR
pα = α −f (x)/k T
B
e
dx
rappresenta la probabilità (Gibbs/Boltzmann) di occupare il metastato α.
La quantità Fα tale che
pα = e−βFα
e−f (x)/kB T dx
si dice energia libera ed è legata all’energia media
f (R)e−f (R)/kB T dR
U α = α −f (x)/k T
B
e
dx
dall’equazione fondamentale
Fα = U α − T Sα
con S : entropia (logaritmo del volume del metastato)
31
Ma, tutto questo, con l’ottimizzazione globale cosa ha a che fare ... ?
Se T è sufficientemente vicina a zero,
p(R) ≈ e−(R−R
/kB T )2
Si possono costruire metodi di ottimizzazione globale che anzichè misurare il
valore della funzione (e minimizzarlo) cerchino di valutare la probabilità
associata ai bacini che si formano a temperatura T , dipendenti sia dai valori
funzionali Uα che dalla dimensione della regione di attrazione Sα
Ciò che è rilevante, in fisica, è la minimizzazione dell’energia libera, non
dell’energia
MonteCarlo Simulated Annealing (Metropolis)
1. k := 0; sia T la temperatura iniziale;
2. x0 : configurazione iniziale;
3. finché la temperatura non è prossima a zero
(a) finché non si è raggiunto l’equilibrio a temperatura T :
i. Sia yk : configurazione generata a caso in un intorno di xk ;
ii. con probabilità
f (yk )−f (xk )
−
kB T
min 1, e
si pone
xk+1 := yk , k := k + 1
(b) Aggiornamento temperatura (ad es.: T := γT con γ ∈ (0, 1))
32
• Come scegliere ed aggiornare la temperatura?
• Come generare soluzioni nell’intorno della soluzione corrente?
• Quando si raggiunge l’equilibrio?
• Quando fermarsi?
Sia la temperatura che la forma dell’intorno possono essere scelti in funzione
della differenza
f (xk ) − fk
tra il valore funzionale corrente e la stima corrente dell’ottimo globale.
Se la differenza è grande ⇒ occorre esplorare la regione ⇒ temperatura alta,
generazione in un intorno “grande” ;
se la differenza è piccola ⇒ raffinamento di un ottimo locale ⇒ si abbassa la
temperatura, si diminuisce il raggio dell’intorno
Ma tutto questo funziona davvero??
Varianti di MonteCarlo simulated annealing
idea: usare scale diverse nella valutazione del panorama energetico: all’inizio
scala grande ⇒ trascurare i dettagli.
Metodologie: smoothing della superficie di energia potenziale.
Smoothing gaussiano:
T −2
1
1
f (x, Λ) ∝ e− 2 βf (y) e− 2 (x−y) Λ (x−y) dy
Packet annealing:
pB (x) ≈
Φα (Vα , xα , Λα ; x)
con
1
1
Φα (Vα , xα , Λα ; x) = e− 2 (βVα + 2 (x−xα )
33
T
Λ−2 (x−xα ))
Potenziali trasformati
Idea elementare: l’ottimizzazione locale trasforma il potenziale eliminando
molte oscillazioni ad “alta frequenza”!
10
9
8
7
6
5
4
3
2
1
0
0
1
2
3
4
34
5
6
7
8
9
10
10
9
8
7
6
5
4
3
2
1
0
0
1
2
3
4
5
6
7
8
9
10
0
1
2
3
4
5
6
7
8
9
10
10
9
8
7
6
5
4
3
2
1
0
35
Basin-hopping monotono
k := 0; E := +∞;
while k < M axIter do
Xk : configurazione random
Xk = arg min E(·; Xk );
Ek = E(Xk );
se Ek < E =⇒ E := Ek
N oImprove := 0;
while N oImprove < M axImprove do
X = perturbazione random di Xk
Y = arg minE(·; X) ;
se E(Y ) < E =⇒ Xk := Y ; N oImprove := 0; E := E(Y )
altrimenti N oImprove + +
end while
end while
10
9
8
7
6
5
4
3
2
1
0
0
1
2
3
4
36
5
6
7
8
9
10
Sia per Multistart che per Basin-Hopping, viene usata una
Ricerca Locale a due fasi:
Prima fase: ottimizzazione di un potenziale modificato (artificiale), costruito
in modo da favorire configurazioni a basso potenziale
Seconda fase: ottimizzazione del potenziale reale, usando l’ottimo della prima
fase come punto di partenza
clusters di LJ
Risultati sperimentali migliori di due ordini di grandezza rispetto alla
letteratura in particolare per clusters difficili (38, 75, 98, 102 atomi)
Prima fase: potenziale modificato con l’aggiunta di penalità lineare tra ogni
coppia di particelle e di una penalità sul diametro.
vM (r) =
2
1
2
− 6 + µr + β max 0, r2 − D2
12
r
r
Seconda fase: potenziale di LJ
Altri esperimenti:
Cluster di Morse: v(r) = (exp(ρ(1 − r)) − 1)2 − 1
Cluster di molecole di acqua.
37
docking di bio-molecole
Potenziale modificato (prima fase):
v(r) =
√
Aij
Bij
−
+αkr
6
12
r
(r/β)
38