Fysisk Organisk Kemi

Ufuldstændig vejledning i brug af Gaussian programkomplekset


Gaussian har nok været det hyppigst anvendte kvantekemiske redskab gennem årene. Det har udviklet sig meget undervejs, og der er en del mindre forskelle mellem de forskellige Gaussian-versioner. Den nyeste hedder Gaussian 03 (g03); den umiddelbart foregående hed g98. Den officielle vejledning (User's Guide) til g98 findes her, mens vejledningen til g03 findes på Gaussians hjemmeside.


Begrundelse
Inddata
%-Angivelser
Route card
Z-Matrix
De første trin
Eksempler på inddatafiler (.com filer)
Litteratur
Lokale retningslinier

Når ikke-teoretikere griber til kvantekemiske metoder som beregninger med Gaussian programkomplekset skyldes det ofte et ønske om at undersøge hvordan et måske utilgængeligt molekyle ser ud i sin grundtilstand, eller at bestemme hvilken energi det har (herunder nulpunktsenergibidrag). Mere ambitiøst kunne ønsket være at undersøge overgangstilstandes egenskaber eller at beregne aktiveringsenergier.

Det første ønske, en nogenlunde pålidelig bestemmelse af geometrien, lader sig normalt opfylde, uden alt for megen møje. Det er ligeledes realistisk at sigte på at bestemme relative energier, fx foretage en sammenligning af isomere molekylers dannelsesvarme. Med moderne sammensatte beregningsmetoder kan man desuden påregne at kunne bestemme dannelsesvarmer for molekyler med op til 8-9 tunge atomer (dvs atomer tungere end hydrogen) med god nøjagtighed.

Beregninger der vedrører overgangstilstandes egenskaber er ligeledes ret pålidelige, omend resultaternes relevans (dvs deres relation til virkeligheden) afhænger kritisk af hvorvidt beregningens forudsætninger (fx om geometri eller reaktionskoordinat) stemmer overens med molekylets faktiske opførsel.

Opsætning af inddata

En beregning med Gaussian 03 kan under Unix afvikles så simpelt som

  g03 < navn-på-inddatafil > navn-på-uddatafil
forudsat at programmet g03 er på ens PATH, og at de forskellige environment variable som g03 bruger er defineret på passende vis (dette sker automatisk på vore egne maskiner). Inddatafilens indhold skal dog være struktureret på ganske bestemt måde, som kan synes knudret. Det er opdelt i fire afsnit:

(betegnelsen route card er fra yngre hulkorttid). Disse udgør hver sin adskilte del af inddata, og med undtagelse af eventuelle %-angivelser skal de alle være tilstede og alle afsluttes med en tom (blank) linie. Formatet bærer præg af Fortran-omgivelser og af fortidige inddatamedier. Det er reelt en forudsætning for brug af Gaussian 03 her på stedet at man besidder et i det mindste rudimentært kendskab til Unix; der er nødvendigt at kunne skrive inddatafilerne, at kunne rette fejl deri, at kunne flytte dem og slette dem.

%-angivelser. Før route card afsnittet kan inddata indeholde særlige instruktioner. Det drejer sig om angivelse af et eventuelt unormalt stort brug af kernelager (ram) og om navnet på en eventuel checkpointfil (se nedenfor). Disse instruktioner skal hver stå på en linie for sig og skal indledes med et %-tegn. At de står først og for sig skyldes at de oprindelig var bestemt for operativsystemet og skulle angives før det egentlige program kunne afvikles.

  %Mem=30000000
  %Chk=iprop.631ss.chk
Her instrueres programmet om højst at anvende 30 mega-ord (millioner ord, hvor et ord er 8 bytes) kernelager under beregningerne, dvs 240 Mb, og om at anvende en checkpointfil benævnt iprop.631ss.chk (filnavne må kun indeholde hvad der i USA anses for normale tal og bogstaver; ingen sære tegn og heller ikke æøå. Se i øvrigt afsnittet om hensigtsmæssig navngivning af Gauss-filer).

Checkpointfiler benyttes både til at indlæse resultater fra tidligere beregninger og til at gemme resultater fra den igangsatte; der læses kun i checkpointfilen i det omfang det angives explicit i route card afsnittet, hvorimod der altid skrives i den. Checkpointfiler er binære (prøv ikke på at læse dem eller skrive dem ud!)

Route card. Denne del af inddata indledes med et obligatorisk #-tegn, efterfulgt af en specifikation af beregningens art og metode. Dette gøres gennem angivelse af keywords; betydningen af de fleste keywords kan yderligere specificeres ved hjælp af forskellige options. Det er obligatorisk at anføre de keywords der angiver hvordan de elektroniske energier skal beregnes (HF, UMP2, CIS, osv), hvilket basissæt der skal benyttes (som fx 3-21G, 6-31G**), og hvilken beregning der skal udføres (Sp, Opt, Freq, osv) (disse besværgelser betegner hhv single point, dvs energien ved en bestemt geometri; geometrioptimering og beregning af energien for den optimerede struktur; beregning af vibrationsfrekvenser).

Herudover kan route card afsnittet indeholde angivelser af i hvilket omfang beregningen skal bygge på tidligere resultater, fx om geometri, kraftkonstanter eller første gæt på bølgefunktionen skal læses ind fra checkpointfilen. Beregningsmæssige detaljer kan angives, fx hvordan visse trin udføres eller hvor mange iterationstrin der maksimalt må bruges, og endelig kan udøves en vis (omend for ringe) kontrol over mængden af uddata.

De enkelte keywords kan adskilles med mellemrum, komma eller skråstreg, eller de kan stå på hver sin linie. Det sidste er det mest overskuelige, omend beregningsmetode og det benyttede basissæt oftest angives sammen. Options til keywords knyttes til dette med et lighedstegn; er der flere skal de stå indenfor parentes, adskilt af kommaer.

Et par eksempler:

  #N UHF/6-31G**
Route card begynder (#), uddata skal ikke omfatte visse beregningsmæssige detaljer (N), beregningen skal være en unrestricted (U) Hartree-Fock (HF) beregning med et udvidet (**) 6-31G basissæt.

  Opt
(eller OPT eller opt; kun i filnavne skelnes mellem store og små bogstaver) Beregningen skal involvere en geometrioptimering (opt).

  Opt=ReadFC	  Opt=(TS,ReadFC)
Geometrioptimering og beregning af den elektroniske energi ved et minimum på potentialfladen, med udgangspunkt (ReadFC) i et tidligere bestemt sæt kraftkonstanter (læses i checkpointfilen); ditto ved et førsteordens saddelpunkt (TS, transition state).
  Freq		  Freq=NoRaman		Freq=(ReadFC,ReadIsotopes)
Beregning af vibrationsfrekvenser (altid efter en forudgående optimering med samme basissæt, typisk med udgangspunkt i dennes checkpointfil); beregningen afvikles hurtigere hvis raman-intensiteter ikke beregnes (NoRaman). I det tredie eksempel beregnes vibrationsfrekvenser for et isotopsubstitueret molekyle (ReadIsotopes, hvor isotopsubstitutionen anføres i inddatafilens sidste afsnit), med udgangspunkt (ReadFC) i en checkpointfil der er fremkommet ved beregning af vibrationsfrekvenserne for det usubstituerede molekyle.
  Geom=Check	  Geom=Modify
I begge tilfælde læses molekylets geometri i checkpointfilen; med modify angives at checkpointfilens geometrioplysninger skal modificeres på enkelte punkter (som oftest ved at en fastholdt variabel, fx en bindingslængde, gives en anden værdi; de konkrete ændringer anføres i sidste inddataafsnit).
  Guess=Check	  Guess=Read
(disse har samme betydning) Det første gæt på bølgefunktionen skal tages fra checkpointfilen.
  Maxdisk=100200300
Visse beregninger har særligt stort behov for diskplads til mellemregninger, og det kan derfor angives separat. Det er dog kun sjældent nødvendigt hos os, da default er sat meget højt.

Tekstfelt. Inddata skal altid indeholde et kommentarfelt (title card), adskilt fra det foranstående route card afsnit og fra den efterfølgende geometrispecifikation med en tom linie. Det er nyttigt her at angive hvilket molekyle beregningen vedrører og hvilken beregning en eventuel checkpointfil stammer fra. Tekstfeltet må højst være på 5 linier med hver max 80 tegn.

Specifikation af molekyle, Z-matrix.

Her skal altid være mindst én linie hvor molekylets ladning og spinmultiplicitet anføres (0 1 for et neutralt molekyle, 1 2 for et kationradikal, osv). Hvis molekylets geometri indlæses fra checkpointfilen (Geom=Check eller Geom=Modify) afsluttes inddata på dette sted med en tom linie. I modsat fald følger (uden tom linie imellem) en såkaldt Z-matrix, hvor molekylet beskrives ved angivelse af bindinger, bindingsvinkler og toplansvinkler. Det er hensigtsmæssigt at dele denne i to afsnit (der skal være adskilt af en tom linie), hvor molekylet i det første specificeres symbolsk, og de symbolske variable i det andet tillægges passende startværdier. I den symbolske beskrivelse er der en linie for hvert atom; først angives atomets art (grundstof), dernæst angives dets binding til et andet atom, så en bindingsvinkel (der involverer et tredie atom), og sidst en toplansvinkel (hvorved et fjerde atom inddrages). Det kan være nyttigt at nummerere atomerne; efter grundstofangivelsen skrives (uden mellemrum) et tal; nummereringen bør af hensyn til uddatas læsbarhed simpelt hen være fortløbende.

  N1
  C2 N1 BND-C2
  C3 C2 BND-C3 N1 ANG-C3
  H4 C3 BND-H4 C2 ANG-H4 N1 DIH-H4
  ...
Her står: C2 er bundet til N1 med en bindingslængde (i Ångström) på BND-C2; C3 er bundet til C2 med en binding der er BND-C3 lang, og de danner en vinkel (i grader) på ANG-C3 til N1; H4 er bundet til C3 med en binding på BND-H4, H4-C3-C2 vinklen er på ANG-H4, og HCCN toplansvinklen er DIH-H4 (også i grader).

I det efterfølgende afsnit (glem nu ikke den mellemliggende tomme linie) tillægges de forskellige bindinger og vinkler startværdier:

  BND-C2=1.25
  BND-C3=1.45
  ANG-C3=120.0
  BND-H4=1.07
  ANG-H4=109.5
  DIH-H4=60.0
  ...
Bemærk at alle tal skal være kommatal (dvs med punktummer...). Og husk endelig at der afsluttes med en tom linie. Ønskes en variabel fastholdt på en bestemt værdi indføjes et F sidst i den pågældende linie i Z-matricens variabelangivelser. Indlæses geometrien fra en checkpointfil med keyword Geom=Modify anføres den frosne variable og dens ny værdi på en linie for sig umiddelbart efter angivelsen af ladning og multiplicitet (dog adskilt fra denne linie med en tom linie).

  BND-C5=1.55 F

Oplysningen om at en variabel er frossen kommer til at indgå i checkpointfilen og overføres til efterfølgende beregninger. Skal en hidtil fastholdt variabel frigøres (inddrages i optimeringen) erstattes F med A (activate), eller Opt=NoFreeze indføjes i route card afsnittet.

Det sværeste ved at komme i gang med Gaussian beregninger er faktisk at få opstillet fornuftige Z-matricer. De variable skal helst være lineært uafhængige, ellers bliver nogle egenskaber specificeret flere gange og andre slet ikke. Det bedste er at tegne molekylerne, med projektionsformler, således at det bliver rimeligt nemt at overskue; derved bliver det lettere at give toplansvinklerne brugbare startværdier.

De første trin

Først opstilles en Z-matrix. Sammen med route card og kommentarafsnit bruges denne til en beregning på relativt simpelt niveau, fx HF/3-21G (UHF erstatter HF dersom molekylet har uparrede elektroner). Inddatafilen kunne gives navnet ethylamin.321.com (se afsnittet om lokale navnekonventioner) og have følgende udseende:

  %Chk=ethylamin.321.chk
  #N HF/3-21G
  Test
  Opt

Indledende beregning; ethylamin på 3-21G niveau

0 1 N1 C2 N1 BND-C2 C3 C2 BND-C3 N1 ANG-C3 H4 C3 BND-H4 C2 ANG-H4 N1 DIH-H4 ..resten af Z-matricen..

(glem ikke de tomme linier, herunder den afsluttende). Keyword Test angiver blot at der ikke automatisk skal produceres en arkiv-version af uddata. Hvis man har datamaten for sig selv og har til hensigt at vente indtil kørslen bliver færdig kan den afvikles fra kommandolinien som følger:
  g03 < ethylamin.321.com > ethylamin.321.out
Dersom inddata var uden fejl og molekylet ikke for stort kan man efter ikke så lang tid læse beregningsresultaterne i ethylamin.321.out; de findes tillige i maskinlæsbar form i ethylamin.321.chk. Går noget galt kommer en stribe aldeles uforståelige fejlmeddelelser (et såkaldt trace-back). I denne sammenhæng kan hvad der er ulæseligt ignoreres; uddatafilen indeholder som oftest de nødvendige oplysninger om hvad problemet var.

Næste trin kunne være en beregning på 6-31G** niveau (i stedet for de to stjerner kan angives d,p i parentes, 6-31G(d,p); dette er mere generelt og foretrækkes). Checkpointfilen fra 3-21G beregningen kopieres (for at bevare den oprindelige version intakt dersom næste beregning går galt; checkpointfiler modificeres undervejs), og gives et nyt navn: ethylamin.631ss.chk (ss da filnavne ikke bør indeholde stjerner). Den ny inddatafil, ved navn ethylamin.631ss.com, kunne have følgende indhold:

  %Chk=ethylamin.631ss.chk
  #N HF/6-31G**
  Test
  Opt
  Geom=Check
  Guess=Check

Ethylamin, HF/6-31G** beregning; ud fra 3-21G checkpoint

0 1

(glem ikke den afsluttende tomme linie). Bemærk at der her ikke anføres nogen Z-matrix, men i stedet Geom=Check. Beregning på 6-31G** niveau tager normalt længere tid end man har lyst at vente, så denne bør afvikles ved hjælp af batch-systemet (se afsnittet om lokal afvikling af Gauss-beregninger).

Herefter kunne man fx gå videre med en beregning af vibrationsfrekvenser, som regel ikke for at korrelere de beregnede frekvenser med et konkret IR-spektrum, men for at bestemme molekylets nulpunktsvibrationsenergi: checkpointfilen kopieres og kopien bruges ved den ny beregning, i inddata erstattes Opt med Freq=NoRaman, og i kommentarfeltet noteres checkpointfilens oprindelse.

Alternativt kunne næste trin, dersom molekylet ikke er for stort, være en beregning med elektronkorrelation, som oftest på MP2-niveau (Møller-Plesset perturbationsregning til anden orden). I inddata erstattes HF/6-31G** med MP2/6-31G**, og hvis Opt linien bevares vil geometrien blive optimeret på MP2-niveau (og det er da bedst med Opt=ReadFC at gå ud fra en checkpointfil fra en forudgående frekvensberegning med samme basissæt). Fjernes Opt udføres beregningen som single-point beregning, dvs med molekylet i den geometri som det forudgående trin (hvorfra checkpointfilen kommer) nåede frem til. Bemærk at 6-31G** beregninger kan tage forholdsvis lang tid selv på moderne arbejdsstationer hvis molekylet har flere end 4-5 tunge atomer, og at frekvensberegninger og MP2-beregninger med geometrioptimering tager endnu længere. Ambitionsniveauet må holdes i ave hvis man deler datamat med andre.

Eksempler på Gaussian inddatafiler:

Nyttig litteratur

hvad angår den konkrete udførelse af beregninger er der ikke meget af. Gaussian 03 - User's Guide (papirversion eller den der er tilgængelig på nettet, se øverst) indeholder blandt andet beskrivelser af alle de mulige keywords og deres options (ind imellem mange), men er ikke velegnet som novicens introduktion til emnet. Bedre er Exploring Chemistry with Electronic Structure Methods: A Guide to Using Gaussian af Foresman og Frisch; vi har en slidt kopi (men den er sikkert udlånt).

Lokale regler, retningslinier og tips ved brug af Gaussian

Vi har Gaussian 94, 98 og 03 kørende på vore maskiner; g94 og g98 på Electron, og g98 og g03 på HP-maskinerne. Der er forskel på hvor meget kernelager disse maskiner har, ligesom der er forskel på deres regnehastighed. De bruges derfor lidt forskelligt: Electron som legeplads og øvelsesområde, og HPerne til mere alvorlige ting (sådan da).

Lager. Der må helst ikke allokeres (med %Mem) mere end 120 Mb kernelager (15 Mord) pr beregning på Electron og ikke mere end 320 Mb (40 Mord) på HPerne, med mindre man træffer særlig aftale.

På Electron er der plads til ca 16 GB diskfiler ialt til mellemregninger (det der specificeres med Maxdisk), så det kan være nødvendigt at fare med lempe. På HPerne er der i almindelighed plads nok (de fleste har 100 GB scratch-plads). Forulykkede beregninger kan efterlade store mellemregningsfiler; på vore maskiner slettes sådanne automatisk (en gang i timen), men dette er normalt ikke tilfældet (bør erindres hvis man afvikler jobs andetsteds).

Afvikling. Som hovedregel benyttes altid et af kø-systemerne. På Electron sker dette med batch:

  batch  gauss  navn-på-batchfil
hvor batchfilen normalt blot indeholder en enkelt linie, g98 < inddatafil > uddatafil; den kan dog indeholde alle former for shell-kommandoer (som fx filkopieringer). Brugeren informeres pr email når beregningen er afsluttet. Batchkøernes længde kan konstateres med kommandoen batchqueue. Vil man slette et (eget) job i køen benyttes kommandoen batchcancel, med angivelse af jobnummer (som findes i batchqueues output; se også man batchcancel),

  batchcancel  jobnummer
Jobs der skal afvikles på HPerne kan kun startes fra Neon. Dette sker med submit:
  submit  navn-på-fil
hvor denne fil har ganske samme indhold som under batch (se ovenfor). En oversigt over aktive jobs kan ses med kommandoen qstat, og jobs kan anulleres med qdel der kræver et jobnummer; dette fremgår af uddata fra qstat.

Filnavne. Ved brug af Gaussian frembringer de fleste på meget kort tid et meget stort antal filer; for at kunne holde rede på hvad der er hvad, og for at tillade vore resultat-ekstraherende scripts at finde de relevante oplysninger bør uddatafilernes navne bygges op efter bestemte konventioner. Det er hensigtsmæssigt at disse retningslinier følges generelt, også for inddatafiler og checkpointfiler (men konventionerne håndhæves ikke af programmellet). Filnavnene bør have tre led, et førsteled der betegner molekylet, et (evt todelt) mellemled der angiver beregningstypen, og et suffix der angiver filtypen. Inddatafiler gives suffix .com, uddatafiler gives suffix .out, og checkpointfiler gives suffix .chk. Beregningsart og -metode angives ligeud ad landevejen: 321 for HF/3-21G, 631ss for HF/6-31G**, MP2 for MP2/6-31G**, freq for frekvensberegning. På denne måde fås navne som

  ipropoh.321.chk	et2me.MP2.out	   meetim.freq.631ss.com
Dersom man regner på flere forskellige udgaver af samme molekyle, fx grundtilstand og overgangstilstand, kan dette bygges ind i forbindelsesnavnet.
  dietS0.MP2.chk	me2iprTfreq.631ss.out	   ethylX3.631.com
Filnavne må indeholde tal og bogstaverne a-z (både store og små), men udover punktum, tankestreg og understreg ingen sære tegn (i særdeleshed ingen skråstreger og ingen mellemrum); de må højst være 255 tegn lange. Sædvanligvis benytter vore scripts tilstedeværelsen af bestemte elementer i filnavnet (såsom out, freq, MP2) til at lokalisere oplysninger, uden hensyn til navnets øvrige udformning, som derfor kan vælges frit.

Katalogstruktur. Det store antal Gauss-filer gør en strukturering hensigtsmæssig. En fornuftig fremgangsmåde består i at benytte et enkelt hoved-katalog, fx gauss, og herunder have kataloger for hvert molekyle (og eventuelt herunder igen underkataloger).

Uddata. Uddata fra Gauss-beregninger er normalt overordentlig omfangsrigt, og det må ikke skrives ud i sin helhed. Det er ganske enkelt for dyrt i papir og toner.

Uanset om det haves på skærmen eller på papir er uddata vanskeligt at overskue. For at finde og udtrække alle de ønskede resultater er det yderst nyttigt at have noget kendskab til Unix-redskaber som grep, sed og (især) awk (eller perl, for de viderekomne). Til de mest almindelige opgaver har vi nogle forholdsvis enkle scripts til hjælp (se nedenfor).

De vigtigste resultater fra energiberegninger står til slut i uddatafilen (den meget tættrykte del). De enkelte felter er her adskilt af \-tegn. Ved frekvensberegninger hentes de fleste resultater inde midt i filen, de vigtigste i afsnittet Thermochemistry.

Scripts. Det kan være nyttigt, og i hvert fald tildels tilfredsstille ens nysgerrighed, at følge aktive beregninger. Her kan yesno, scfdone og eump være nyttige (eump er kun relevant for MP2-beregninger). De afvikles fra det katalog hvor beregningens uddatafil skrives; scfdone og eump viser hvordan den beregnede energi ændrer sig efterhånden som iterationerne skrider frem; yesno viser hvorledes det går med konvergensen.

Resultater der er frembragt med henblik på efterfølgende RRKM-beregninger kan med lidt held udtrækkes med getfreq.

Beregninger med sammensatte metoder kan sættes op med makeall. Denne fordrer at der findes en checkpointfil fra en geometrioptimering på det pågældende molekyle. Scriptet be'r selv om de nødvendige oplysninger og afleverer den fil der efterfølgende skal bruges når jobbet startes med batch eller submit

Resultaterne af en beregning med fx den sammensatte G3-metode ekstraheres fra uddatafilerne med scriptet getg3. Resultater fra andre beregninger med sammensatte metoder kan høstes med get32 (G3MP2-metoden), getg3b (G3//B3LYP), get32b (G3MP2//B3LYP), getcbs (CBS-Q metoden), getqb3 (CBS-QB3 metoden), og enkelte mindre vigtige lignende scripts. Masochister er velkomne til at studere disse scripts; de befinder sig i /usr/local/bin.

Visualisering. Ønsker man at se beregningsresultaterne på grafisk maner benyttes programmet molden. Det kan vise den endelige struktur (og illustrere den måde molekylet vred og vendte sig på under optimeringen), det kan vise hvilke vibrationer frekvensberegningen kom frem til (og sågar et simuleret IR-spektrum), og ret mange andre ting. Det kan kun anvendes fra en af vore Unix-maskiner (kræver X-windows).

  molden  navn-på-uddatafil
Diskplads. Der er altid for lidt, især da Gauss-filer fylder meget. For at undgå vrøvl bør udtjente filer slettes, og filer man måske får brug for senere bør komprimeres (gælder især .chk og .out filer), ved hjælp af programmet gzip:
  gzip  molekyle.type.chk
Herved fremkommer en fil ved navn molekyle.type.chk.gz (dvs filnavnet udstyres med et ekstra suffix, .gz); gzippede filer fylder betydeligt mindre end originalerne. Den oprindelige fil kan gendannes med gunzip:
  gunzip  molekyle.type.chk.gz
(.gz suffixet kan udelades efter behag når gunzip påkaldes).


Artikler; nogle af vore egne bidrag.
Beslægtede emner: G3, CBS-Q og andre sammensatte metoder og Termokemiske beregninger med Gaussian
Oversigt over forskningsområder.
Ordliste; ordforklaringer og terminologi.


12/94; rev 4/05 steen@kiku.dk