Retro Programming Nostalgia: quando l’Amiga 500 calcolava le cariche atomiche

/ Danilo Roccatano

Un programma in GFA BASIC, il metodo di Gasteiger–Marsili e un ricordo della mia tesi di laurea in Chimica.

All’epoca della ricerca per il mio progetto di laurea in chimica, ero ancora un entusiasta possessore dell’Amiga 500. Il kernel del Linux OS era stato appena ufficialmente rilasciato nella versione 0.01 da Linus Torvalds, e la distribuzione Slackware sarebbe stata rilasciata nel luglio del 1993. Come descritto nel mio articolo sulla mia tesi di laurea, all’epoca l’Apple Macintosh per scrivere la tesi e il mainframe IBM 3090 per i calcoli. Il mio primo PC lo acquistai durante il mio dottorato di ricerca con l’idea di installarlo con Linux in dual boot con MS Windows. All’epoca conoscevo bene il BASIC disponibile sui vari microcomputer che ho avuto (vedi i miei articoli sul retrocomputing) e l’Amiga Basic del mio Amiga 500 che acquistai durante i miei primi anni universitari. Quando iniziai la mia tesi di laurea, cominciai a imparare a programmare in FORTRAN 77, linguaggio usato per scrivere i codici per le simulazioni di dinamica molecolare. Continui limitatamente ad usare anche il mio Amiga 500 per la grafica e la programmazione.

Nella mia tesi di laurea ho effettuato degli studi di dinamica molecolare di un composto cristallino del rame. Per le simulazioni era necessario ottenere delle informazioni sulla distribuzione di equilibrio delle cariche atomiche parziali del composto. Quando si descrive una molecola mediante un modello computazionale, uno degli aspetti più importanti è la distribuzione della carica elettrica tra i suoi atomi. Sebbene una molecola possa essere complessivamente neutra, gli elettroni non sono necessariamente distribuiti in maniera uniforme tra i suoi costituenti.Questa distribuzione può essere rappresentata mediante cariche atomiche parziali, che costituiscono una componente fondamentale di molti modelli di meccanica molecolare e permettono di descrivere, almeno approssimativamente, le interazioni elettrostatiche.

Il calcolo quantomeccanico delle cariche parziali per composti con un elevato numero di atomi, con in aggiunta un elemento di transizione, era proibitivo; pertanto, fu necessario usare dei metodi empirici. Tra i metodi più semplici scelsi di usare metodi basati sull’equilibrazione dell’elettronegatività. Uno di essi è il metodo empirico proposto da Johann Gasteiger e Mario Marsili nel 1980, basato sul principio dell’equalizzazione dell’elettronegatività attraverso il trasferimento iterativo di carica tra atomi legati.

Il metodo di Gasteiger-Marsili

Il metodo di Gasteiger–Marsili, noto anche come Partial Equalization of Orbital Electronegativity (PEOE), è un metodo empirico per il calcolo delle cariche atomiche parziali basato sul principio dell’equalizzazione dell’elettronegatività. L’idea fondamentale è che, quando due atomi sono legati chimicamente, la differenza tra le loro elettronegatività determina una redistribuzione della carica elettronica.

Il metodo descrive l’elettronegatività di ciascun atomo mediante una funzione quadratica della sua carica parziale:

\chi_i(q_i)=a_i+b_iq_i+c_iq_i^2

dove (q_i) rappresenta la carica parziale dell’atomo (i), mentre (a_i), (b_i) e (c_i) sono parametri empirici che dipendono dal tipo di atomo e dal suo ambiente chimico.

Il calcolo procede attraverso una serie di iterazioni.

1. Assegnazione delle cariche iniziali

A ciascun atomo viene assegnata una carica iniziale (q_i^{(0)}), generalmente nulla per gli atomi di una molecola neutra, salvo quando siano necessarie cariche iniziali differenti.

2. Calcolo delle elettronegatività

Utilizzando i parametri empirici, il programma calcola l’elettronegatività di ogni atomo in funzione della sua carica corrente:


\chi_i^{(n)}=a_i+b_iq_i^{(n)}+c_i(q_i^{(n)})^2

3. Trasferimento di carica tra atomi legati

Per ogni coppia di atomi direttamente legati, il programma calcola una correzione della carica proporzionale alla differenza di elettronegatività:

\frac{\chi_j^{(n)}-\chi_i^{(n)}}{D_{ij}^{(n)}}

dove (D_{ij}^{(n)}) è un parametro di normalizzazione determinato dalle caratteristiche degli atomi coinvolti.

Nell’implementazione PEOEP, questo parametro è ricavato dai coefficienti empirici dell’atomo con elettronegatività minore, utilizzando una costante specifica per l’idrogeno.

4. Aggiornamento iterativo delle cariche

La carica di ciascun atomo viene aggiornata sommando i contributi degli atomi direttamente legati e applicando un fattore di smorzamento:

q_i^{(n)} + f_n\sum_{j\in N(i)}\Delta q_{ij}^{(n)}

dove (N(i)) rappresenta l’insieme degli atomi legati all’atomo (i).

Nel programma originale il fattore di smorzamento è:

f_n=(0.5)^n

con le iterazioni numerate a partire da (n=1). Il contributo dei trasferimenti di carica diminuisce quindi progressivamente durante il calcolo.

5. Calcolo delle cariche finali

Al termine di ogni iterazione vengono ricalcolate le elettronegatività utilizzando le nuove cariche atomiche. Il procedimento viene ripetuto per un numero prefissato di iterazioni, pari a dieci nel programma PEOEP originale.

Al termine del calcolo, il programma restituisce le cariche parziali di ciascun atomo e confronta la carica totale iniziale con quella finale, permettendo di verificare la conservazione della carica elettrica della molecola.

L’implementazione in GFA BASIC

Non disponendo del programma originale di Gasteiger–Marsili, decidemmo, insieme al mio supervisore di tesi, di implementare autonomamente l’algoritmo in FORTRAN 77, il linguaggio che stavo imparando proprio in quel periodo e che utilizzavamo per i programmi di simulazione di dinamica molecolare.

Tuttavia, non disponendo di un compilatore FORTRAN sul mio Amiga 500, decisi di sviluppare anche una versione dell’algoritmo in BASIC, un linguaggio con cui avevo già acquisito una buona familiarità grazie ai microcomputer utilizzati negli anni precedenti.

Per questo progetto volli sperimentare il GFA BASIC, un dialetto BASIC sviluppato da Frank Ostrowski nel 1986, inizialmente per i computer Atari ST e successivamente disponibile anche per Amiga. Rispetto all’Amiga BASIC, il GFA BASIC offriva numerose funzionalità avanzate, sia per il calcolo numerico sia per la gestione dei file e della grafica, oltre alla possibilità di compilare i programmi. Una caratteristica che trovai particolarmente utile fu l’istruzione FILESELECT, che permetteva di selezionare il file di input attraverso una finestra di dialogo, conferendo al programma un’interfaccia più pratica e professionale.

Nacque così il piccolo programma scientifico che scrissi nel 1992 e il cui listato originale è riportato nell’appendice di questo articolo.

Questo programma può essere adattato rapidamente al QB64PE BASIC (che ho usato negli altri articoli di questa serie) sostituendo le linee:

FILESELECT "PEOEP","LOAD MT","DF0:PEOEP/",nmol$
OPEN "i",#1,nmol$

con 

nmol$ = _OPENFILEDIALOG$("Seleziona file molecolare", "", "*.*", "File molecolari", 0)
IF nmol$ = "" THEN END

Il nome file di input si deve trovare nello stesso direttorio in cui si trova l’eseguibile el programma, altrimenti occorre fornire il percorso completo della posizione sul disco. Il file di input contiene la lista degli atomi e dei legami atomici della molecola.

Struttura generale del programma

Definizione dei parametri atomici:

Il programma include un set predefinito di parametri per 15 tipi di atomi e gruppi atomici, come idrogeno (H), carbonio (C), azoto (N), ossigeno (O), e altri. Ogni atomo è descritto tramite tre parametri numerici:

• Valore base dell’elettronegatività.

• Coefficiente lineare per il calcolo della carica.

• Coefficiente quadratico per il calcolo della carica.

Esempio:

DATA H , 7.17, 6.24,-0.56

DATA C , 7.98, 9.18, 1.88

Input della topologia molecolare:

Gli utenti forniscono un file con la descrizione della molecola, inclusi:

• Numero di atomi.

• Tipo di ogni atomo.

• Connettività tra gli atomi.

Calcolo iniziale dell’elettronegatività

Il programma utilizza una formula quadratica per stimare l’elettronegatività iniziale di ogni atomo:

chi(i) = par(t%,1) + par(t%,2)*q(i) + par(t%,3)*q(i)*q(i)

Routine principale

Il calcolo dell’elettronegatività viene raffinato iterativamente. In ogni ciclo:

• Gli atomi “donano” o “ricevono” carica in base alla differenza di elettronegatività con gli atomi vicini.

• L’aggiornamento delle cariche tiene conto di un fattore di smorzamento (damp) per garantire la stabilità della convergenza.

q(i) = q(i) + (somh + somk) * (damp ^ alpha)

Output

L’output include:

• Cariche iniziali e finali per ogni atomo.

• Cariche totali iniziali e finali della molecola.

• Numero di iterazioni e fattore di smorzamento utilizzati.

Esempio pratico

Analizziamo la molecola del bromuro di etano (CH_3CH_2-Br). Fornendo al programma la topologia della molecola

AT
1-BROMOETANO
8
1,2,4,2,3,4,5,0,0,0,0
2,1,1,1,0,0,0,0,0,0,0
3,1,1,1,0,0,0,0,0,0,0
4,1,1,1,0,0,0,0,0,0,0
5,2,4,1,6,7,8,0,0,0,0
6,1,1,5,0,0,0,0,0,0,0
7,1,1,5,0,0,0,0,0,0,0
8,12,1,5,0,0,0,0,0,0,0
0
0
0
0
0
0
0
0

Utilizzando i parametri atomici predefiniti nel programma, PEOEP calcola le cariche parziali di ogni atomo e l’elettronegatività risultante.

Conclusioni

Rileggere oggi il listato originale di PEOEP, scritto nel 1992 durante la preparazione della mia tesi di laurea in Chimica, è stato un piacevole viaggio nel passato. All’epoca il mio Amiga 500 rappresentava soprattutto uno strumento per sperimentare con la programmazione e la grafica, ma con questo progetto divenne anche un piccolo laboratorio di chimica computazionale.

Il metodo di Gasteiger–Marsili si prestava particolarmente bene a un’implementazione su un microcomputer: non richiedeva risorse di calcolo particolarmente elevate e permetteva di affrontare un problema scientifico concreto attraverso un algoritmo relativamente semplice.

A distanza di oltre trent’anni, il programma può essere recuperato e adattato a un moderno ambiente BASIC come QB64PE, dimostrando quanto alcuni algoritmi scientifici possano essere longevi, anche quando le piattaforme hardware e i linguaggi di programmazione cambiano radicalmente.

Questo piccolo programma rappresenta quindi non soltanto un ricordo dei miei primi esperimenti di programmazione scientifica, ma anche una testimonianza di un periodo in cui i confini tra informatica domestica e calcolo scientifico cominciavano a diventare sempre meno netti.

BIBLIOGRAFIA

  1. J. Gasteiger e M. Marsili, Iterative Partial Equalization of Orbital Electronegativity—A Rapid Access to Atomic Charges, Tetrahedron, 1980, 36, 3219–3228.

APPENDICE

Il programma originale scritto in GFA BASIC per microprocessore Amiga 

'
' PEOEP: PARTIAL EQUALIZATION of ORBITAL
' ELECTRONEGATIVITY PROGRAM
'
' Copyright (c): 1992 Danilo Roccatano
'

DIM smatr(300,12),q(300),par(20,4),chi(300)
DIM naat$(20),les(10),hig(10),qi(300)
' DEFINT "SMART,NATOM,NAP,NITER"
' DEFFLT "Q,PAR,CHI,DAMP,SOMK,SOMH,eka"
damp=0.5
niter=10
'
' INPUT PARAMETRI
'
FOR t%=1 TO 15
  READ naat$(t%)
  FOR k=1 TO 3
    READ par(t%,k)
  NEXT k
NEXT t%
DATA H , 7.17, 6.24,-0.56
DATA C , 7.98, 9.18, 1.88
DATA C1, 8.79, 9.32, 1.51
DATA C2,10.39, 9.45, 0.73
DATA N ,11.54,10.82, 1.36
DATA N1,12.87,11.15, 0.85
DATA N2,15.68,11.70,-0.27
DATA O ,14.18,12.92, 0.47
DATA O1,17.07,13.79, 0.47
DATA F ,14.66,13.85, 2.31
DATA CL,11.00, 9.69, 1.35
DATA BR,10.08, 8.47, 1.16
DATA I , 9.90, 7.96, 0.96
DATA S ,10.14, 9.13, 1.38
DATA Cu, 4.48, 5.07, 1.21
'
' INPUT MOLECULAR TOPOLOGY
'
FILESELECT "PEOEP","LOAD MT","DF0:PEOEP/",nmol$
OPEN "i",#1,nmol$
INPUT #1,title$
INPUT #1,natom
FOR t=1 TO natom
  FOR y=1 TO 11
    INPUT #1,smatr(t,y)
  NEXT y
NEXT t
FOR t=1 TO natom
  INPUT #1,q(t)
  qi=qi+q(t)
  qi(t)=q(t)
NEXT t
CLOSE #1
'
' CALCULATION OF INITIAL ELECTRONEGATIVITIES
'
FOR i=1 TO natom
  t%=smatr(i,2)
  chi(i)=par(t%,1)+par(t%,2)*q(i)+par(t%,3)*q(i)*q(i)
NEXT i
PRINT
'
' MAIN ROUTINE
'
FOR alpha=1 TO niter
  FOR i=1 TO natom
    g=0
    l=0
    somh=0
    somk=0
    FOR k=1 TO smatr(i,3)
      IF chi(smatr(i,k+3))<=chi(i) THEN
        g=g+1
        les(g)=smatr(i,k+3)
      ELSE
        l=l+1
        hig(l)=smatr(i,k+3)
      ENDIF
    NEXT k
    FOR h=1 TO l
      co=smatr(i,2)
      IF co=1 THEN
        eka=20.02
      ELSE
        eka=par(co,1)+par(co,2)+par(co,3)
      ENDIF
      somh=somh+(chi(hig(h))-chi(i))/eka
    NEXT h
    FOR h=1 TO g
      co=smatr(les(h),2)
      IF co=1 THEN
        eka=20.02
      ELSE
        eka=par(co,1)+par(co,2)+par(co,3)
      ENDIF
      somk=somk+(chi(les(h))-chi(i))/eka
    NEXT h
    q(i)=q(i)+(somh+somk)*(damp^alpha)
  NEXT i
  FOR i=1 TO natom
    t%=smatr(i,2)
    chi(i)=par(t%,1)+par(t%,2)*q(i)+par(t%,3)*q(i)*q(i)
  NEXT i
NEXT alpha
'
' OUTPUT CARICHE
'
OPEN "O",#1,nmol$+".OUT"
PRINT #1,title$
PRINT #1,"OUTPUT CARICHE"
PRINT #1,"DUMPING Factor:";damp
PRINT #1,"Number of iterations:";niter
PRINT "S.N.";"  TIPO","C.INIZIALE";"      C.FINALE"
PRINT #1,"S.N.";"  TIPO","C.INIZIALE";"      C.FINALE"
FOR t%=1 TO natom
  PRINT t%;"       ";naat$(smatr(t%,2));" ",qi(t%)," ";q(t%)
  PRINT #1,t%;"       ";naat$(smatr(t%,2));" ",qi(t%)," ";q(t%)
  qt=qt+q(t%)
NEXT t%
PRINT "Carica totale iniziale: ";qi
PRINT "CARICA TOTALE FINALE  : ";qt
PRINT #1,"Carica totale iniziale: ";qi
PRINT #1,"CARICA TOTALE FINALE  : ";qt
CLOSE #1
END

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.