GAMESSで電子密度分布を出したい

はじめに

GAMESSで電子密度をcube形式で欲しくなったのに苦戦したので共有を
（Gaussianならcubegenで一発なんですけどね……）

バージョン違いの影響は確認していないので、悪しからず

インプット

$ELDENSを指定すれば電子密度を出してくれる。

ただ、電子密度を計算する位置の指定の仕方に一癖あって、
計算したい位置(複数可)を直接指定するWHERE=POINTSと
計算したい面を指定するWHERE=GRIDの2種類しか実質ない。
（なので、直接3次元格子点を指定する方法はない様子）
今回はcubeのグリッド座標を$POINTSに羅列することにした。

$POINTSの一行目には座標の単位(BOHR)とグリッドの点数(N, 勿論整数)を指定する。
（最初公式ドキュメントを読んでもよく意味が分からずここで地味に苦労した）

 $ELDENS
IEDEN=1,  ! 1: compute the electron density
MORB=0,   ! 0: total electron density
WHERE=POINTS, ! calculate on the grids specified by $POINTS
OUTPUT=PUNCH  ! output only to PUNCH
 $END
 $POINTS
BOHR N
-9.666705 -8.941050 -6.512752
-9.666705 -8.941050 -6.304549
-9.666705 -8.941050 -6.096346
(以降続く)
 $END

これで動くかと思いきや、

     $SYSTEM OPTIONS
     ---------------
  REPLICATED MEMORY=     3000000 WORDS (ON EVERY NODE).
 DISTRIBUTED MEMDDI=           0 MILLION WORDS IN AGGREGATE,
 MEMDDI DISTRIBUTED OVER   4 PROCESSORS IS           0 WORDS/PROCESSOR.
 TOTAL MEMORY REQUESTED ON EACH PROCESSOR=     3000000 WORDS.
 TIMLIM=    80000000.00 MINUTES, OR   55555.6 DAYS.
 PARALL= T  BALTYP=  DLB     KDIAG=    0  COREFL= F
 MXSEQ2=     300 MXSEQ3=     150  mem10=         0
 **** ERROR, MAXIMUM NO. OF POINTS IN $POINTS IS  100
 EXECUTION OF GAMESS TERMINATED -ABNORMALLY- AT Sat Jul 12 04:32:10 2025

とエラー落ちしてしまう。
どうもデフォルトでは$POINTSには100点までしか指定できないらしく、
これを弄るオプションも見つからなかったため、ソースコードのパラメータを弄った。

cat /path/to/gamess/source/prplib.src | nl -b a | grep MXPTPT
    74        PARAMETER (MXPTPT=100)
    83        COMMON /POINTS/ NPOINT,IPUNIT,XPOINT(MXPTPT),YPOINT(MXPTPT),
    84       *                ZPOINT(MXPTPT)
  1119        PARAMETER (MXPTPT=100)
  1137        COMMON /POINTS/ NPOINT,IPUNIT,XPOINT(MXPTPT),YPOINT(MXPTPT),
  1138       *                ZPOINT(MXPTPT)
  1169        CALL VCLR(XPOINT,1,MXPTPT)
  1170        CALL VCLR(YPOINT,1,MXPTPT)
  1171        CALL VCLR(ZPOINT,1,MXPTPT)
  1187        IF(NPOINT.GE.MXPTPT) THEN
  1188           IF(MASWRK) WRITE(IW,920) MXPTPT
  1589        PARAMETER (MXPTPT=100)
  1593        COMMON /POINTS/ NPOINT,IPUNIT,XPOINT(MXPTPT),YPOINT(MXPTPT),
  1594       *                ZPOINT(MXPTPT)

ここのMXPTPT=100を適当に大きい数に変えてコンパイルし直せば動く。

アウトプット

OUTPUT=PUNCHにした場合電子密度がdatファイルに書き出される。

 ELECTRON DENSITY, IPOINT,X,Y,Z,EDENS
      1      -9.66671    -8.94105    -6.51275     0.156157E-08
      2      -9.66671    -8.94105    -6.30455     0.188021E-08
      3      -9.66671    -8.94105    -6.09635     0.224508E-08
      4      -9.66671    -8.94105    -5.88814     0.265850E-08
(以降続く)

後はawkなどで加工すればcubeファイルとしてデータを作れますね。
例えば

cat datfile | 
    tac | # SCF計算の度に出力されるので最後のモノをもってくる
    sed -n '1,/ELECTRON DENSITY/p' | 
    tac | 
    sed -e 1d | # ELECTRON DENSITY, IPOINT,X,Y,Z,EDENSの部分を消す
    awk 'NR==1{n=NF}{if(n==NF) print $5; else exit}' | # 電子密度が書き出されている箇所を抜き出す
    tr '\n' ' ' | 
    sed -r 's/(([^ ]+ ){N3})/\1\n/g' | # N3にはcubeの3軸目のグリッド点数を入れる
    sed -r 's/(([^ ]+ ){6})/\1\n/g' > cubfile