Il calcolo delle cariche parziali atomiche

Le cariche parziali da usare nel campo di forze si ottengono generalmente da calcoli quantomeccanici. Nel caso di molecole rigide, il calcolo delle cariche è abbastanza semplice. Nel caso di molecole flessibili, occorre valutare quanto le diverse conformazioni influenzano la distribuzione di carica e, quindi, stimare la carica parziale come media pesata tra i vari conformeri. Il calcolo QM fornisce le cosidette cariche di Mulliken. Questo tipo di cariche possono portare ad una elevata inaccuratezza nel riprodurre proprietà chimico-fisiche di piccole molecole. Per evitare questo inconveniente sono state introdotte varie procedure per ottenere delle cariche parziali che tengano conto della diversa capacitá dei singoli atomi di accomodare una diversa distribuzione di carica. Queste procedure vanno sotto il nome di metodi di Electrostatic potential fitting tra cui i più usati solo il RESP e il CHELPG. Vediamo come questi metodi funzionano.

Il potenziale elettrostatico {\phi(\vec{r})} nel punto {\vec{r}} può essere determinato dalla densità elettronica attraverso la relazione:

\displaystyle \phi(\vec{r}) = \phi_N(\vec{r}) + \phi_e(\vec{r}) =\sum_{i=1}^N \frac{Z_i}{\vec{r} - \vec{R_i}} -\int\frac{dr' \rho(\vec{r})}{|\vec{r'}-\vec{r}|} \ \ \ \ \ (22)

dove {\phi_N(\vec{r})} e {\phi_e(\vec{r})} sono i contributi nucleare ed elettrronico alla carica totale, N è il numero di atomi, {\vec{R_i}} è al posizione dell’atomo i-esimo. Per ottenere una descrizione a cariche atomiche puntiformi di questo potenziale, è necessario effettuare una procedura di modellazione dei valori di queste cariche in modo da minimizzare la funzione:

\displaystyle R=\sum_{i=1}^M w_i (\phi_(\vec{r_i})^{QM} -\phi_(\vec{r_i})^{Coul}) \ \ \ \ \ (23)

ovvero, la somma dei minimi quadrati della differenza tra il potenziale quantomeccanico e quello calcolato da queste cariche in un certo numero di punti {M} opportunamente scelti. {w_i} sono degli opportuni pesi da attribuire ai vari punti in base alla loro importanza. La carica sull’atomo N-esimo ({q_N}) è legata alle altre attraverso la relazione:

\displaystyle q_N = Z - \sum_{j=1}^{N-1} q_i \ \ \ \ \ (24)

dove {Z} è la carica molcolare totale, pertanto il potenziale coulombiano descritto dalle N cariche si può esprimere come:

\displaystyle \phi_(\vec{r_i})^{Coul}= \sum_{j=1}^{N-1} \frac{q_j}{4 \pi \epsilon_0 r_{ij}} + \frac{Z-\sum_{j=1}^{N-1} q_i}{4 \pi \epsilon_0 r_{iN}} \ \ \ \ \ (25)

derivando {R} in {q_k} si ottiene

\displaystyle \frac{\partial R}{\partial q_k} =-2 \sum_{i=1}^M w_i (\phi (\vec{r_i})^{QM} -\phi(\vec{r_i})^{Coul})\left( \frac{\partial \phi (\vec{r_i})^{Coul}}{\partial q_k}\right) =0 \ \ \ \ \ (26)

che può essere scritta come un set di M equazioni

\displaystyle \sum_{i=1}^{M} w_i\left(\phi (\vec{r_i})^{QM}-\frac{Z}{r_{iN}}\right)\left(\frac{1}{r_{ik}}-\frac{1}{r_{iN}}\right) = \sum_{j=1}^{N-1}\left[\sum_{i=1}^{M} w_i \left(\frac{1}{r_{ik}}-\frac{1}{r_{iN}}\right) \left(\frac{1}{r_{ij}}-\frac{1}{r_{iN}}\right)\right] \frac{q_j}{4\pi \epsilon_0} \ \ \ \ \ (27)

o, in forma matriciale {\vec{a}=\vec{A}\vec{q}}. Usando i metodi dell’algebra lineare, la precedente equazione può essere risolta ottenendo le cariche {\vec{q} = \vec{A^{-1}}\vec{a}}. La scelta degli M punti viene fatta in differenti modi, ma devono essere presi in una regione di spazio che sia significativa per definire un corretto modello di interazioni molecolari. Due sono le procedure più utilizzate per questo calcolo:

  • CHELPG: sviluppato da Breneman e Wiberg nel 1990. In questo metodo, la molecola viene posta in una reticolo cubico di punti (con spaziatura compresa tra 0.3 e 0.8 {\AA}). Tutti i punti che si trovano entro il raggio di van der waals e {2.8 \AA} al di fuori di questo vengono scartati, e i punti rimanenti usati per il calcolo.
  • RESP: sviluppato dal gruppo di P. Kolmann per calcolare le cariche parziali da usare nel campo di forza AMBER. In questo caso, i punti vengono presi da superfici molecolari costruite usando raggi di van der waals progressivamente crescenti, e imponendo dei vincoli su particolari atomi (per esempio carboni localizzati all’interno della molecola) che potrebbero altrimenti assumere cariche eccessivamente alte. A differenza del CHELPG, le cariche ottenute con questo metodo sono meno sensibili alla conformazione molecolare.

Per effettuare il calcolo delle cariche parziali si può ricorrere ai programmi per calcoli di QM, GAUSSIAN o GAMESS/US (quest’ultimo è un programma liberamente distribuito). Per il calcolo è necessario disporre di coordinate iniziali della molecola. Queste possono essere ottenute dal data base di strutture (come il Cambridge Structural Database) oppure costruite usando programmi di modellistica molecolare. Un programma molto diffuso di questo tipo è MOLDEN che può essere liberamento scaricato dal sito: http://www.cmbi.ru.nl/molden

dicloroetano

Figure 1: Molecola di 1,2-dicloroetano in conformazione anti

Usando questo programma è possibile costruire la struttura iniziale della molecola e salvare su file le coordinate nel formato GAUSSIAN o GAMESS/US, per il calcolo delle cariche. Come esempio di calcolo delle cariche viene mostrato il contenuto del file di input per GAUSSIAN per al calcolo a punto singolo delle conformazioni syn, gauge e anti del 1,2,dicloroetano.

%Chk=DCLETANO
# RHF/6-31G(d) FOpt Test Pop=(CHelpG,dipole)
Calcolo delle cariche con il metodo di CHELPG

0 1
 c
cl   1 clc2
 c   1 cc3        2 cccl3
 h   1 hc4        2 hccl4        3 dih4
 h   1 hc5        2 hccl5        3 dih5
cl   3 clc6       1 clcc6        2 60.
 h   3 hc7        1 hcc7         6 dih7
 h   3 hc8        1 hcc8         6 dih8
clc2        1.750000
cc3         1.540000
cccl3       109.471
hc4         1.089000
hccl4       109.471
dih4        120.000
hc5         1.089000
hccl5       109.471
dih5       -120.000
clc6        1.750000
clcc6       109.471
hc7         1.089000
hcc7        109.471
dih7        120.000
hc8         1.089000
hcc8        109.471
dih8        240.000
dih6          0.0

Come si può vedere dalla prima linea, i programmi effettuano dapprima una minimizzazione delle coordinate, e una volta che si è raggiunta la convergenza viene effettuato un calcolo ESP usando il metodo CHELPG.

Nella tabella I, sono riportati i valori delle cariche ottenute per le varie configurazioni e il relativo momento di dipolo.

6.2. Calcolo di un profilo di energia potenziale per rotazione di un legame.

Quando in una molecola sono presenti vari angoli diedri, può essere utile calcolare il profilo dell’energia conformazionale ottenuta dalla rotazione intorno a questi. Questo profilo può essere usato per modellare l’energia potenziale classica dovuta alla rotazione intorno al suddetto diedro, usando il potenziale di Rycheman-Belleman ed escudendo le interazioni di non-legame intramolecolari tra gli atomi coinvolti.

In GAUSSIAN è possibile calcolare automaticamente un profilo di energia conformazionale. Nell’esempio che segue viene mostrato come ottenere, il profilo di energia conformazionale per la molecola di 1,2-dicloroetano. Di seguito viene mostrato il file di input per GAUSSIAN per il calcolo del profilo di energia.

%Chk=DCLETANO
# RHF/6-311G(D) Opt=Z-matrix NoSymm Test

Relaxed PES scan
0 1 
c 
cl 1 clc2 
. 
. 
. 

dih6 0. 18 10.

Nell’ultima riga viene indicato l’angolo diedro Cl-C-C-Cl (dih6) interessato alla rotazione. La scansione viene effettuata ruotando 18 volte di 10 gradi l’angolo diedro. Per estrarre dal file di uscita il profilo di energia in funzione dell’angolo si può usare lo script in awk riportato alla fine di questo articolo.

Nella Figura 2, viene riportato il grafico del profilo di energia ottenuto dalla scansione dell’angolo diedro.

PES_DCLET.jpg

Figure 2. Andamento della profilo di energia potenziale ottenuto dalla rotazione intorno a un diedro Cl-C-C-Cl del 1,2-dicloroetano.

 

#! /usr/bin/awk -f
BEGIN {
# Incremental step of the scansion
        step = 10.
# Conversion factor Hartree --> Kcal/mol
        conv = 627.5095
        cc = 0
        count = 0
     }
$1 == "EIGENVALUES" { for (i=3;i<=NF;i++) {
                     k=split($i,a,"-")
                      for (ii=1;ii<=k;ii++) {
                         if (a[ii] != "") {
                           val = -a[ii]
                           en[cc] = val
                           if (val < min) {min = val}
                           count+=step
                           cc++
                         }
                        }
                       }
		    }
END {
      count = 0
      for (ii=0;ii<cc;ii++) {
         printf "% 4.1f   %8.3f\n", count, (en[ii]-min)*conv
         count+=step
      }
    }

 

About Danilo Roccatano

I have a Doctorate in chemistry at the University of Roma “La Sapienza”. I led educational and research activities at different universities in Italy, The Netherlands, Germany and now in the UK. I am fascinated by the study of nature with theoretical models and computational. For years, my scientific research is focused on the study of molecular systems of biological interest using the technique of Molecular Dynamics simulation. I have developed a server (the link is in one of my post) for statistical analysis at the amino acid level of the effect of random mutations induced by random mutagenesis methods. I am also very active in the didactic activity in physical chemistry, computational chemistry, and molecular modeling. I have several other interests and hobbies as video/photography, robotics, computer vision, electronics, programming, microscopy, entomology, recreational mathematics and computational linguistics.
This entry was posted in Research, Science Topics. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s