bash script 用於計算化學軟體間的分子軌域資料格式轉換

0317025 李仲淳





一、相關背景知識

1.

MO(Molecular Orbital)用來描述分子的能量及電子在空間中分布的情形(與薛丁格方程式相關)。分子軌域可藉由SCF(自治場方法)運 算出來,而市面上有各種不同的計算化學或分子模擬軟體善於利用不同的計算方法來算出各種不同分子的MO能量以及係數。每一套軟體都 有其各自計算的方法與擅長的部份,而我以與報告內容較有關聯性: Orca 以及 Gaussian 分別作介紹。

2.

以Gausse View 所繪出的 Mn2O2(NH3)8,後續將以此模型為例子介紹MO及相關應用。


3.

分子軌域形成時的一種特殊狀況:Broken symmetry,指的是形成分子軌域的兩個原子軌域一個都為Alpha spin另一個全為Beta spin的狀況,這種狀況讓分子能量最低、最穩定,常發生在過度金屬上,這邊以上述例子中的兩個錳金屬為例,最上面的狀態為電子都沒有成對,而且都是Alph a spin 的狀態。中間則為成對的狀況,而最下面的例子便是Broken symmetry的情形。隨著電子Alpha spin及Beta spin的不同,兩個錳 金屬間的距離也不相同。若電子全為Alpha spin,兩金屬間的排斥力相對較大,造成兩金屬的距離較遠。相反的,若Alpha 與Beta spin 成對於兩金屬之間,則造成兩金屬的距離較近。



4.ORCA & Gaussian優點

1. ORCA

1.) 此款軟體善於利用DFT(Density functional theory:利用電子密度來計算分子能量)

2.) 善於計算Broken symmetry(利用flip spin的方式)

2.Gaussian

1.) Optimization最佳化結構

2.) 找尋Transition state


因此在計算具有Broken symmetry的分子結構時,通常會以Orca這套軟體進行,再藉由其他軟體像是Gaussian來進行分子構型上的最佳 化。而這些不同的軟體,在Input資料檔以及Output結果檔的檔案格式當然各不相同,因此我這次要介紹我是如何將Orca輸出後的Output 檔格式轉換為Gaussian的Input檔案,以利於在計算完Broken symmetry後進行結構最佳化。

二、Shell Script 內容

在轉換這兩個檔案前必須先確認,兩檔案的格式。

ORCA的Output檔格式

ORCA Output檔中將屬於各個能量的所有MO係數表現在同一個直欄上,並且每次只有五個欄位(五個能量),當係數結束後,同一個直攔便會有下個能量及其所屬係數,所以第一個能量及第六個能量與係數為同一直欄,第二與第七的能量與係數在同一直欄,以次類推。


Gaussian的Input檔格式

Gaussain Input檔的格式則將屬於各個能量的所有係數,以橫行表示,並且每五個係數變換行。


Script 的內容

我的Script利用ORCA紀錄能量以及係數的特性(每五個能量與係數為一個大的矩陣),先將所有能量抓出來丟到屬於能量的陣列,再透過陣列的編號,將屬於各個能量的所有係數丟到各個係數陣列(例如:將第一個能量編號為一,並將其所對應的係數丟到一號陣列),最後再寫一個迴圈,利用陣列的index把各個能量以及各個能量的係數print出來,並可以將Orca output檔案的格式轉換成Gaussian input檔案的格式了。以下為Script的內容

#!/bin/bash
printf "(5E16.8) \n"

#check if the file is broken symetry
test=$(grep -n 'BETA' $1 | awk -F "_" '{print $2}' | head -n 1)
if [[ "$test" == "BETA " ]]; then
   t=0
else
   t=1
fi

#The line number distinguish alpha and beta energy
Beta=$(grep -n 'COEFF_BETA' $1 | awk -F ":" '{print $1}')

#The output without Broken symetry
if [[ $t -eq 1 ]]; then
  line=$(grep -n 'a1g' $1 | awk -F ":" '{printf " " $1}')
  declare -a linenumber=(`echo "$line"`)
  coe=$((${linenumber[1]}-${linenumber[0]}-2))
  Eandcoe=$(($coe+1))
  declare -a E=()

#form array of E

  for (( i=0 ; i<${#linenumber[@]} ; i++ ))
  do
     startline=$((${linenumber[$i]}+1))
     Energyline=$(tail -n +$startline $1 | head -n 1)
     E+=($Energyline)

#form coefficient correspond to E
     startitem=$(($i*5))
     for (( j=$startitem ; j<${#E[@]} ; j++ ))
     do
        counting=$(($j+1))
        startcolumn=$(($j-($i*5)+1))
        co=($(tail -n +$startline $1 | head -n $Eandcoe | tail -n $coe | sed -e 's/   / /g' -e 's/  / /g' -e 's/^ //g' |
        cut -d ' ' -f $startcolumn ))
        printf "    $counting. Alpha MO OE="
        printf "%+16.8e" ${E[$j]}
        printf "\n"
        for (( h=0 ; h<${#co[@]} ; h++ ))
        do
               I=$(($h%5))
               if [[ $h>1 && $I -eq 0 ]]; then
                   printf "\n"
               fi
               printf "%+16.8e" ${co[$h]}
        done
        printf "\n"
     done
  done
fi

#The output with broken symetry
if [[ $t -eq 0 ]]; then
  line=$(grep -n 'a1g' $1 | awk -F ":" '{printf " " $1}')
  declare -a linenumber=(`echo "$line"`)
  coe=$((${linenumber[1]}-${linenumber[0]}-2))
  Eandcoe=$(($coe+1))
#form array of E alpha

  declare -a Ea=()
  declare -a Eb=()
  for (( i=0 ; i<${#linenumber[@]} ; i++ ))
  do
     if [[ "${linenumber[$i]}" -lt "$Beta" ]]; then
        startline=$((${linenumber[$i]}+1))
        Energyline=$(tail -n +$startline $1 | head -n 1)
        Ea+=($Energyline)

#form alpha coefficient correspond to E

        startitem=$(($i*5))
        for (( j=$startitem ; j<${#Ea[@]} ; j++ ))
        do
            counting=$(($j+1))
            startcolumn=$(($j-($i*5)+1))
            co=($(tail -n +$startline $1 | head -n $Eandcoe | tail -n $coe | sed -e 's/   / /g' -e 's/  / /g' -e 's/^ //g' |
            cut -d ' ' -f $startcolumn ))
            printf "    $counting Alpha MO OE="
            printf "%+16.8e" ${Ea[$j]}
            printf "\n"
            for (( h=0 ; h<${#co[@]} ; h++ ))
            do
               I=$(($h%5))
               if [[ $h>1 && $I -eq 0 ]]; then
                   printf "\n"
               fi
               printf "%+16.8e" ${co[$h]}
            done
            printf "\n"
        done
     fi
  done
          printf "0 \n"

#form array of E beta till second last line
  last=$((${#linenumber[@]}-1))
  for (( i=0 ; i<$last ; i++ ))
  do
     if [[ "${linenumber[$i]}" -gt "$Beta" ]]; then
          startline=$((${linenumber[$i]}+1))
          Energyline=$(tail -n +$startline $1 | head -n 1)
          Eb+=($Energyline)
          startitem=$((${#Eb[@]}-5))
          for (( g=$startitem ; g<${#Eb[@]} ; g++ ))
          do
              counting=$(($g+1))
              startcolumn=$(($g-$startitem+1))
              co=($(tail -n +$startline $1 | head -n $Eandcoe | tail -n $coe | sed -e 's/   / /g' -e 's/  / /g' -e 's/^ //g' |
              cut -d ' ' -f $startcolumn ))
              printf "    $counting Beta MO OE="
              printf "%+16.8e" ${Eb[$g]}
              printf "\n"
              for (( h=0 ; h<${#co[@]} ; h++ ))
              do
                 I=$(($h%5))
                 if [[ $h>1 && $I -eq 0 ]]; then
                    printf "\n"
                 fi
                 printf "%+16.8e" ${co[$h]}
              done
              printf "\n"
          done
     fi
  done

#form array of E beta of the last line
  declare -a Elast=()
  startline=$((${linenumber[$last]}+1))
  Energyline=$(tail -n +$startline $1 | head -n 1)
  Elast+=($Energyline)
  for (( i=0 ; i<${#Elast[@]} ; i++ ))
  do
     lastN=$(($counting+$i+1))
     startcolumn=$(($i+1))
     co=($(tail -n +$startline $1 | head -n $Eandcoe | tail -n $coe | sed -e 's/   / /g' -e 's/  / /g' -e 's/^ //g' |\
           cut -d ' ' -f $startcolumn ))
     printf "    $lastN Beta MO OE="
     printf "%+16.8e" ${Elast[$i]}
     printf "\n"
     for (( h=0 ; h<${#co[@]} ; h++ ))
     do
        I=$(($h%5))
        if [[ $h>1 && $I -eq 0 ]]; then
            printf "\n"
        fi
        printf "%+16.8e" ${co[$h]}
     done
     printf "\n"
  done
fi

三、結果

將ORCA計算能量的結果以及利用Shell Script轉換成Gaussian Input檔並利用Gaussian計算能量後的結果做比較

利用ORCA跑出的結果


利用Gaussian跑出的結果


從上面的結果可以得到轉換之後跑出來的能量是相同的,但是後續再確認MO構型時,3D模擬結果卻並不相同,還有待後續的確認,預計在大專生計畫後,繼續後續的確認。