へっぽこ博士のなんちゃって研究室

PythonやR、分子軌道計算等に関する記事を書きます

R : 原子量表のスクレイピング

R

0.事始め

原子量一覧が欲しかったので、日本化学会にて公開されているものからスクレイピングでデータフレームを構築しようとしたら、結構苦労したので、Rのソースコードを晒しておきます。Rにおける文字列処理方法の勉強にはなったと思います。apply関数？知らない子ですね。

1.ゴリ押し大作戦

pdftoolsパッケージを使って6～7ページを読み込むと、以下のようなデータが得られます。

"     1 水           素 H      1.008               44 ル  テ  ニ  ウ  ム Ru  101.1"                                   
"     2 ヘ   リ   ウ   ム He     4.003               45 ロ   ジ   ウ   ム Rh  102.9"                                  
"                                      3 リ   チ   ウ   ム Li     6.94                46 パ  ラ  ジ  ウ  ム Pd  106.4"
"     4 ベ  リ  リ  ウ  ム Be     9.012               47       銀       Ag  107.9"                                    
"     5 ホ     ウ     素 B     10.81                48 カ  ド  ミ  ウ  ム Cd  112.4"                                  
"     6 炭           素 C     12.01                49 イ  ン  ジ  ウ  ム In  114.8"

日本語の元素名が均等割り付けされているため、不要なスペースが紛れ込んでおり、そのままでは区切り文字で各要素を分割することができません。
こういった体裁だけのための書式は、データサイエンティストにとって、本当に面倒です。
とりあえず、以下の作戦でスクレイピングすることにしました。
・1 数字またはアルファベットの前後のスペースを、
　　別の区切り文字に置き換える。
・2 余分なスペース、重複した区切り文字を除去。
・3 データフレームに入れて整形。

2.プログラム本体

で、最終的にこうなりました。

require(pdftools)
require(tidyverse) #stringr, tidyr, dplyr

# check multiple-byte character : return bool
is_multiple <- function(char) {
  pattern <- "\\p{Block=Katakana}|\\p{Script=Han}"
  invisible(return(str_detect(char,pattern=pattern)))
}

# check digit, alphabet and period : return bool
is_alnum <- function(char) {
  if (is_multiple(char)) invisible(return(FALSE))

  pattern <- "[[:alnum:]]|\\."
  invisible(return(str_detect(char,pattern=pattern)))
}

# check blank : return bool
is_blank <- function(char) {
  if (is_multiple(char)) invisible(return(FALSE))

  pattern <- "[[:blank:]]"
  invisible(return(str_detect(char, pattern=pattern)))
}

# extract character at a specified position : return character
chr_sub <-function(string, locate) {
  invisible(return(str_sub(string, locate, locate)))
}

# replace character at a specified position : return string
chr_replace <- function(string, locate, char) {
  forward <- str_sub(string, 1, (locate-1))
  backward <- str_sub(string, (locate+1), -1)
  invisible(return(str_c(forward,char,backward)))
}

# rebuild string : return data.frame
rebuild_str <- function(string) {
  new_str <- str_trim(string)
  N <- str_length(new_str)
  if (is.na(N)|(N < 2)) invisible(return(new_str))
  
  for (i in 2:(N-1)) {
    char <- chr_sub(new_str, i)
    if (is_blank(char)) {
      fchar <- chr_sub(new_str, (i-1))
      bchar <- chr_sub(new_str, (i+1))
      if (is_alnum(fchar)|is_alnum(bchar)) {
        new_str <- chr_replace(new_str, i, ",")
      }
    }
  }
  new_str <- str_replace_all(new_str, " ","" )
  sp_str <- str_split(new_str, "[,]+") %>% as.data.frame()
  invisible(return(t(sp_str)))
}

# extended bind_rows : return data.frame
bind_rows_ex <- function(df1, df2) {
  colnames(df2) <- colnames(df1)
  bind_df <- bind_rows(df1, df2)
  invisible(return(bind_df))
}

# rebuild page : return data.frame
rebuild_page <- function(source, page, start=1L, end=-1L) {
  if ((end==-1L)|(end > length(source[[page]]))) end <- length(source[[page]])

  txt <- source[[page]][start:end]
  txt.df <- data.frame(matrix(nrow=(end-start), ncol=8))
  for (i in 1:length(txt)) {
    txt.df[i,] <- rebuild_str(txt[i])
  }
  txt.df <- bind_rows_ex(txt.df[,1:4], txt.df[,5:8])
  invisible(return(txt.df))
}

# main process
fname <- "atom_2022.pdf"
suppressMessages(pdf_file <- pdf_text(fname))

pdftxt <- str_replace_all(pdf_file, c("＊\n"="", "†\n"="", "（"="", "）"=""))
pdftxt <- str_split(pdftxt, "\n", simplify = F)
p1 <- rebuild_page(pdftxt, 6, 12, 54) 
p2 <- rebuild_page(pdftxt, 7, 3, 18)

atom_table <- bind_rows(p1, p2)
colnames(atom_table) <- c("Number","Element","Symbol","Weight")
View(atom_table)


3.ソースコードの解説

is_multiple(char)
is_alnum(char)
is_blank(char)
文字種類の判定用関数です。
最初に全角、いわゆる2バイト文字の判定処理を入れないと、1バイト目の文字コードからアルファベットとして判断されることがあるようです。
長音符号(ー)はUnicodeではKatakanaブロックにあります。 解説サイトでは、片仮名の判定をスクリプト(Script, sc, または省略)で指定されているほうが多いようで、この場合は、長音符号が含まれないので注意が必要です。
今回の場合、漢字の範囲はScriptで十分足りますが、人名や地名を判定する必要がある場合は不十分です。 詳しくはこちらのサイト様が参考になると思います。

chr_sub(string, locate)
stringからlocateで指定した位置の文字を返すだけのラッパー関数です。

chr_replace(string, locate, char)
string中のlocateで指定した位置の文字をcharに置き換える関数です。意外とありそうでない関数の１つで、 pythonでも同様の質問がされていました。 C言語だとポインタを渡して簡単にアクセスできるのに・・・

rebuild_str(string)
stringから空白を取り除き、要素に分解したデータフレームを返します。作戦1～2のメイン処理部分です。
日本語の元素名は不要という方は、下記のような実装でもいいかもしれません。

rebuild_str <- function(string) {
  new_str <- str_trim(string)
  N <- str_length(new_str)
  if (is.na(N)|(N<2)) invisible(return(new_str))

  for (i in 2:(N-1)) {
    char <- chr_sub(new_str, i)
    if (is_multiple(char)) {
      new_str <- str_remove(new_str, char)
    }
  }
  
  sp_str <- str_split(new_str, "[ ]+") %>% as.data.frame()
  invisible(return(t(sp_str)))
}


連続した区切り文字を1文字として扱うには、正規表現では"[ ]+"のように後ろに+を付けます。

bind_rows_ex(df1, df2)
2つのデータフレームを縦に結合するラッパー関数です。
dplyr::bind_rowsは、base::rbindと比較して、列数や列名が異なっていても結合できる便利な関数です。 ただ、思っていたのと挙動が違ったので、用意しました。 色々試行錯誤した結果、列名にNULLを入れて結合することに落ち着いたので、rbindで良かったような・・・

bind_rowsの挙動は、以下のような2つのデータフレームを結合すると、

AB
12

CD
34

こんな感じになります。

ABCD
12NANA
NANA34


こうなって欲しかったんですけどね・・・

A'B'
12
34


rebuid_page(source, page, start, end)
pdf_textで読み込んだデータ(source)の指定ページ(page)を処理する関数です。
startとendで読み込み開始行と終了行を指定します。 WindowsとLinuxとで改行コードの違いからか、若干ずれるので、お手持ちの環境に合わせて調整ください。
データ処理部分をfor文で回していますが、sapply関数でも代替できます。 ただ、実行時間に差が出なかった(どちらも0.5秒程度)のと、後々の処理が面倒だったので、書き慣れている方で処理しました。 必ずしもapply系関数の方が優れているわけではないという一例ですね。

メイン処理部分
pdftools::pdf_textを使うと、"PDF error: Invalid Font Weight"という大量のエラーを吐くことがあります。 これは、内部で使用しているPopplerライブラリが原因のようなので、ヘルプにある通り、suppressMessagesで出力を抑制しています。

後は、必要に応じてwrite.csv等で出力しておくのも良いかもしれません。

Quantum ESPRESSO : qe-7.0ビルド 備忘録

Quantum ESPRESSO

1. ビルド環境

Ubuntu 20.04.3LTS on VirtualBox 6.1.30
https://www.ubuntulinux.jp/ubuntu
日本語Remix、最小インストール、アップデート済み(kernel 5.13.0-28-generic)
VirtualBox Extension Pack適用済み

2. 前準備

2-1. ビルドツール類のインストール

sudo apt install build-essential gfortran
sudo apt libtool autoconf
sudo apt install cmake cmake-qt-gui
sudo apt install curl git ssh


2-2. Quantum ESPRESSO 7.0のダウンロード

https://www.quantum-espresso.org/
公式のdownloadページからアクセスするにはユーザー登録が必要になりました。利用者の職種や目的の調査用だそうです*1。
GitHub、GitLabは認証なしでアクセス可能なので、面倒な方はgitで複製するのが吉。
https://github.com/QEF/q-e/
https://gitlab.com/QEF/q-e/

git clone https://github.com/QEF/q-e.git ./qe-7.0

または、

git clone https://gitlab.com/QEF/q-e.git ./qe-7.0

ダウンロード先を変更したい場合は、末尾のディレクトリ指定を任意に変更ください。
なお、GitHubからは、.tar.gzでも一応ダウンロード可能です。

3. ビルド＆インストール

3-1. 同梱パッケージのみでビルド

とりあえず使えればいいという方向け。
README通りのやり方です。

./configure --prefix=$HOME/.local/qe
make all -j
make install

pw.xを実行して、クレジットが表示されたら、インストール成功です。

$HOME/.local/qe/bin/pw.x

Ctrl+Cで中断してください。

cmakeでビルドする場合は、少し注意が必要です。
READMEには明記されていませんが、
QE_LAPACK_INTERNAL=ON
と
QE_FFTW_VENDOR=Internal
の2つを指定しないと、configureに失敗します。

mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$HOME/.local/qe -DCMAKE_Fortran_COMPILER=mpif90 -DCMAKE_C_COMPILER=mpicc -DQE_LAPACK_INTERNAL=ON -DQE_FFTW_VENDOR=Internal ..
make -j
make install


3-2. APTパッケージによるビルド

実行速度を求めたり、並列計算や様々な汎関数を使いたい方向け。
特にこだわりがなければ、以下が一番楽です。

sudo apt install openmpi-bin libopenmpi-dev 
sudo apt install libopenblas-dev liblapack-dev libfftw3-dev libfftw3-mpi-dev libxc-dev

./configure --prefix=$HOME/.local/qe --with-libxc
make all -j
make install


こちらの場合でcmakeする場合は、README通りの方法で問題なくconfigureが通ります。

mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$HOME/.local/qe -DCMAKE_Fortran_COMPILER=mpif90 -DCMAKE_C_COMPILER=mpicc ..
make -j
make install

ただ、追加機能分のオプション指定が必要で、オプションパラメータの前に「-D」を付けます。 コマンドラインによるオプション指定は面倒なので、cmake-guiをインストールしている環境なら、CMakeCache.txtをダブルクリックすると、より簡単に編集できます。

(編集の一例)

CMAKE_INSTALL_PREFIX  $HOME/.local/qe
ENABLE_SCALAPACK_MPI  ON
ESPRESSO_PSEUDO       .
QE_ENABLE_LIBXC       ON
QE_ENABLE_SCALAPACK   ON
QE_FFTW_VENDOR        FFTW3

疑ポテンシャルファイルを自動探索するディレクトリが、デフォルトではソースディレクトリのpseudoになっているので、指定し直したほうが使い勝手が良いと思います。

Intel製CPUなら、openblasをintel-mklに変更することも可能です。ただ、configureが通らなくなるので(lapackとfftwを見つけられないようです)、cmakeでビルドします。

sudo apt install intel-mkl-full

インストール途中の選択肢はデフォルトのままOKを選択して問題ありません。

mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$HOME/.local/qe -DCMAKE_Fortran_COMPILER=mpif90 -DCMAKE_C_COMPILER=mpicc -DENABLE_SCALAPACK_MPI=ON -DESPRESSO_PSEUDO=. -DQE_ENABLE_LIBXC=ON -DQE_ENABLE_SCALAPACK=ON -DQE_FFTW_VENDOR=Intel_FFTW3 ..
make -j
make install

上述のとおり、オプション部分はCMakeCache.txtを編集してもOKです。


3-3. ライブラリのソースビルド

最新のライブラリにこだわりたい方向け。
正直、makeに失敗することが多いので、お勧めできません。
具体的には、手動インストールしたもののパスを指定しても、そんな場所にはないと怒られてconfigure自体が通らなかったり、configureが通った場合でも、FFTWに関連する部分でコンパイルに失敗したりします。
makefileに問題があるため*2、自身で編集できるレベルの知識が必要です。ちなみに、現時点の筆者の知識レベルではすべてを解決できませんでした。


○ Open MPI

https://www.open-mpi.org/
Downloadページからダウンロードします(2022/2/13時点では、4.1.2が最新の安定版)。
3種類用意されていますが、.tar.gzを選択するのが無難です。
並列処理を担うソフトとしては、他にmpich*3があります。
実行ファイル名や環境変数がかぶっていたり、BLAS類の設定方法が異なるなど面倒事が多いので、使っているソフトの仕様上、切り替えて使用する必要があるというような特殊な事情でもない限り、どちらか１つに絞ったほうが良いと思います。

tar -xvf openmpi-4.1.2.tar.gz
cd openmpi-4.1.2

./configure --prefix=$HOME/.local/openmpi
make all
make install

環境設定
.bash_aliases(デフォルトでは存在しません)を開いて、

gedit ~/.bash_aliases

以下を追記します。

export MPIROOT=$HOME/.local/openmpi
export LD_LIBRARY_PATH=$MPIROOT/lib
export PATH=$PATH:$MPIROOT/bin

新たなエイリアスやパス設定等は、.bash_aliasesに書き込むのが良いと、.bashrc内に記載があります。
不幸な事故*4を防ぐためにも、.bashrcは極力触らないほうが良いです。
編集が終わったら、ターミナルを再起動します。

exec $SHELL --login

mpiexecやmpirunなどを実行して、パスが通っていれば、インストール作業は成功です。

mpiexec --version


○ libxc

https://tddft.org/programs/libxc/
Downloadページよりダウンロードします(2022/2/13時点では5.2.2が最新版)。

tar -xvf libxc-5.2.2.tar.gz
cd libxc-5.2.2
./configure CC=gcc FC=gfortran --enable-fortran --prefix=$HOME/.local/libxc
make
make test
sudo make install

QEでは、Fortran用のライブラリが必要なので、 「--enable-fortran」を指定します。
このスイッチは、helpでは出てこないのですが(記載し忘れ？)、configureファイル内を探ると、きちんと設定があり、デフォルトでは「--disable-fortran」が設定されることがわかります。

GitHubからもダウンロード可能です。
https://github.com/ElectronicStructureLibrary/libxc

git clone https://gitlab.com/libxc/libxc.git ./libxc-5.2.2
cd libxc-5.2.2
autoreconf -i .

autoreconfでconfigureを再生成します。
親ディレクトリを示す最後の「.」を忘れないようにしましょう*5。
後の手順はtar.gzの場合と同じです。

cmakeでビルドする場合は、f90に対応するファイルが欠落しているので、pkgconfigディレクトリ内のlibxcf90.pc.inをcmakeディレクトリにコピーしておく必要があります。

QEで使用する場合は、以下の3つのスイッチ
--with-libxc
--with-libxc-prefix
--wuth-libxc-include
を指定するのですが、パッケージインストールした場合とファイル構成が異なるためか、configure側でうまく認識してくれません。
手動ビルドした後、すべてのディレクトリにファイルをコピーして同じ内容にしてみましたが、解決できませんでした。

○ LAPACK

http://www.netlib.org/lapack/
QE同梱のBLAS、LAPACK、CBLASです。
最新版は、3.10.0(2022/2/13時点)。

mkdir build
cd build
cmake -DCMAKE_INSTALL_LIBDIR=$HOME/.local/lapack ..
cmake --build . -j --target install

qe/installディレクトリ内のextlibs_makefileを編集したら、自動でインストールされないかと試したのですが、現実は甘くなかったです。
一応、編集したものを以下に記載しますが、cmakeの場合しか通りません(configureはエラーになります)。

14～16行目

LAPACK_NETLIB=lapack-3.10.0.tar.gz
LAPACK_NETLIB_NAME=lapack-3.10.0
LAPACK_NETLIB_URL=https://github.com/Reference-LAPACK/lapack/archive/refs/tags/v3.10.0.tar.gz

3.9.0以降はGitHubに移行しているため、URLの指定が少し面倒です*6。

○ FFTW

https://www.fftw.org/
最新版は、3.3.10(2022/2/13時点)。

./configure --enable-mpi
make
make install

QEのconfigureでは、インストール先をデフォルトにしても認識してくれませんでした。
cmakeだとQE_FFTW_VENDORにFFTW3を指定すると認識はしてくれますが、Wannier90のあたりでコケます。

○ Wannier90

http://www.wannier.org/
2022/2/13時点では同梱のものと同じですので、敢えて単体をダウンロードする必要はないと思います。

cp ./config/make.inc.gfort ./make.inc
make -j
sudo make install

gfortranでコンパイルする場合は、make.inc.gfortをコピーします。そのほかのコンパイラを使用する場合は、README.installを参照ください。
デフォルトでは/usr以下にインストールされます。

4. 疑ポテンシャルファイル

適切な疑ポテンシャルファイルは、目的や計算対象によって異なるため、まずは文献等の先例に従って、その結果を再現できるか確認することが重要です*7。そのうえで、変更を加えたときに、どのくらい信頼性を担保できるか検証することが大切です。
ここでは、公式のlegacy pseudopotentialsにあるPSlibraryのみを紹介します。
http://pseudopotentials.quantum-espresso.org/

○ PSLibrary
https://github.com/dalcorso/pslibrary
公式の周期表から１つずつダウンロードしてもよいのですが、このスクリプトを使用すると、必要なポテンシャルファイルを手元で一度に生成できます。xml形式で生成されるため、タグの分だけファイルサイズが大きくなるのがやや難点です。

git clone https://github.com/dalcorso/pslibrary.git ./pslibrary.1.0.1

実行にあたって、いくつかファイルの編集が必要です。
・QE_path

PWDIR=$HOME/.local/qe

QEのインストール場所、具体的にはld1.xが存在するディレクトリ(/bin)を指定します。最後のパス区切り文字(/)は入力しないよう注意。

・make_ps
11行目

element='H. C. N. O. Al Si Fe Au'

全元素(all)を生成すると時間がかかるので、必要な元素記号を指定します。
一文字の元素は2文字目に「.」を付け、元素間はスペースで区切ります。

21行目

#. $PWDIR/environment_variables

旧版(<5.0)用の設定なので、コメントアウトしておきます。
エラー表示が気にならなければ、そのままでも大丈夫です。

・make_ps_all
15行目～

cd ./pz
. ../make_ps
cd ..

必要な疑ポテンシャルファイルをディレクトリ名で指定します。
通常はpzまたはpbeのみで十分事足りると思います。

ここまで準備できたら、

./make_ps_all

で、それぞれのディレクトリのPSEUDOPOTENTIALS以下に疑ポテンシャルファイルが生成されます。
テスト用のディレクトリ(WORK)も生成されますが、こちらは必要ないので、削除してもかまいません。

5. 参考サイト

雑多な記録
第一原理計算ソフトQuantumESPRESSOの使い方
http://www2.yukawa.kyoto-u.ac.jp/~koudai.sugimoto/dokuwiki/doku.php?id=quantumespresso

株式会社クロスアビリティ
Quantum ESPRESSO入力ファイル作成手順
３．擬ポテンシャルファイルの選択方法
https://qiita.com/xa_member/items/006d9dd44c662a17903d

*1:2022/1/22確認時。 しばらく前からサーバー証明が切れているとは思っていたのですが、リニューアルと同時に変更になったようです。

*2:おそらくパッケージインストールしか想定されていないんだろうと思います。

*3:https://www.mpich.org/

*4:https://www.nemotos.net/?p=3529

*5:READMEの説明だと、文章の終わりを示すピリオドに見えるので、何も生成されない、configureできないと、しばらくハマりました。

*6:指定が間違っていると、一応警告は出るのですが、少し先まで処理が流れてから止まるので、ビルドできない原因を突き止めるまでにずいぶん時間がかかりました。

*7:先例がない場合は、実験値をどのくらい再現できるか、適用範囲はどのくらいかの検証を行います。

Quantum ESPRESSO : cif2pw

Quantum ESPRESSO Python

- 事始め


しばらくxlt2pw*1のお世話になっていたのですが、変換をコマンドラインだけで完結したい、カットオフエネルギーを自動でセットしたい等、自分用に改良していったら別物になってしまったので、晒しておきます。
一応、Quantum ESPRESSOのPW/toolsフォルダにcif2qe.shというスクリプトがあるのですが、awkプログラム言語で書かれているうえ、長らくメンテナンスされていないようなので、これを編集する気にはなれませんでした。

- コード

#!/usr/bin/python3

# cif2pw.py : modified xtl2pw.py
# version 1.0 programmed by Doktor-Heppoko

import sys
import os
import subprocess as sp
import numpy as np
from pymatgen.core import periodic_table as mg

# environment check
# check python version
if sys.version_info < (3, 6):
    print("Python 3.6 or lattar required")
    sys.exit(1)
    
# check NUMPY version
if np.version.version < "1.10":
    print("Numpy 1.10 or lattar required")
    sys.exit(1)

# check cif2cell installed
try:
    proc = sp.run(["which","cif2cell"],stdout=sp.PIPE)
    if proc.returncode != 0:
        print("cif2cell not installed")
        sys.exit(1)
except Exception as emsg:
    print(emsg)
    sys.exit(1)

# argment check
if len(sys.argv) < 3:
    print("Usage      : cif2pw.py [ CIF file ] [ method ] [ constraint ]")
    print("For example: > python3 cif2pw.py hoge.cif relax 1")
    print(" method    : scf, relax, vc-relax")
    print(" constraint: 0 fixed, 1 free[default]")
    sys.exit(1)

# initialize variable literal
# argments
cif = sys.argv[1]
method = sys.argv[2]
try:
    flag = sys.argv[3]
except:
    flag = 1
flag_str = "free" if flag == 1 else "fixed"
print(f"Reading File                 : {cif}") 
print(f"Calculation Method           : {method}") 
print(f"Constraint Flag              : {flag_str}") 

# file related
fbase = os.path.splitext(cif)[-2]
ciftmp = "cif2cell.tmp"
output = f"{fbase}.{method}.in"

# calcuration condition related
# CONTROL card
nstep = 100              # iterations for optimization
etot_conv_thr = "1.0D-4" # convergence threshold on total energy
forc_conv_thr = "1.0D-3" # convergence threshold on forces
pseudo_dir = '/home/ubuntu/Downloads/pslibrary.1.0.0/pbe/PSEUDOPOTENTIALS/'

# SYSTEM card
nat = 0                # number of atoms (auto set)
ntyp = 0               # kinds of atomic species (auto set)
nbnd = 0               # total valence electrons (auto set)
ecutwfc_def = 25       # kinetic energy cutoff for wavefunctions
ecutrho_def = 100      # kinetic energy cutoff for charge density and potential
occupations = "fixed"  # for insulators such as general organic compounds
input_dft = "vdw-df2"  # See Modules/funct.f90
vdw_corr = "Grimme-D3" # non-local functions should be specified by input_dft
use_nbnd = False       # auto-calculation flag : use = True (valence * 1.2)
use_vdw_corr = True    # use van der Waals correction flag : use = True

# ELECTRONS card
electron_maxstep = 100    # iterations for SCF
conv_thr = "1.0D-6"       # convergence threshold for SCF
mixing_mode = "local-TF"  # for insulators such as general organic compounds
mixing_beta = 0.4         # ratio of new wave functions to be mixed
diagonalization = "david" # diagonalization algorithm

# IONS card
ion_dynamics = "bfgs"       # mathematical optimization algorithm for ions
ion_temperature = "initial" # ion temperature control method
tempw = "297.16"            # 25 degrees Celsius

# CELL card
cell_dynamics = "bfgs"   # mathematical optimization algorithm for cell
press = "1.0D-3"         # in kbar : 1.0D-3kbar = 0.1MPa = ca. 1 atm
press_conv_thr = "0.5D0" # convergence threshold on pressure
cell_dofree = "all"      # moved cell parameters

# ATOMIC_SPECIES card
kinds = np.empty(0)      # element symbol
wt = np.empty(0)         # atomic mass from pymatgen
ps = np.empty(0)         # pseudo potensial file generated from element symbol
pstype_HA = ".pbe-kjpaw_psl.1.0.0.UPF"     # for H atom
pstype_NM = ".pbe-n-kjpaw_psl.1.0.0.UPF"   # for non-metal elements
pstype_ME = ".pbe-spn-kjpaw_psl.1.0.0.UPF" # for metal elements

atom_count = np.empty(0) # amount of element
ps_type = np.empty(0)    # pseudo potential type from file
ps_func = np.empty(0)    # dft functional type from file
valences = np.empty(0)   # valence electrons from file
wfc = np.empty(0)        # kinetic energy cutoff for wave function from file
rho = np.empty(0)        # kinetic energy cutoff for charge density and potential from file

# ATOMIC_POSITIONS card
atom = np.empty(0) # atomic label from cif2cell template
px = np.empty(0)   # x-coordinate from cif2cell template
py = np.empty(0)   # y-coordinate from cif2cell template
pz = np.empty(0)   # z-coordinate from cif2cell template

# K_POINTS card
nk = np.ones(3, dtype='int') # number of supplied special k-points (auto set)
k_point = 700                # maximum points by Monkhorst-Pack grids

# CELL_PARAMETERS card 
V = np.zeros((3,3)) # lattice vector [ 3 x 3 ] from cif2cell
cell = np.ones(3)   # lattice length [ a b c ]

# main processing
# creat pwscf .in template file 
try:
    proc = sp.run(["cif2cell","-q","-p","pwscf","-f",cif,"-o",ciftmp,"--pwscf-cartesian"],stdout=sp.PIPE)
    if proc.returncode != 0:
        sys.exit(1)
except Exception as emsg:
    print(emsg)
    sys.exit(1)

# read template file
try:
    with open(ciftmp,'r') as f:
        ftxt = f.readlines()
except OSError as emsg:
    print(emsg)
    sys.exit(1)
    
# delete template file
try:
    proc = sp.run(["rm",ciftmp],stdout=sp.PIPE)
except Exception as emeg:
    print(emsg)
    sys.exit(1)

# get cell parameters
for nline, sline in enumerate(ftxt):
    if "CELL_PARAMETERS" in sline:
        break
if nline == (len(ftxt)-1):
    print("Not Found CELL_PARAMETERS card")
    sys.exit(1)
for i in range(3):
    astr = ftxt[nline+i+1].strip().split()
    for j in range(3):
        V[i][j] = float(astr[j])

for i in range(3):
    cell[i] = np.sqrt(V[i] @ V[i])    
#print(f"Lattice Length               : {cell}")

# get atom positions
for nline, sline in enumerate(ftxt):
    if "ATOMIC_POSITIONS" in sline:
        break
if nline == (len(ftxt)-1):
    print("Not Found ATOMIC_POSITIONS")
    sys.exit(1)
for sline in ftxt[nline+1:]:
    if len(sline) == 0:
        break
    astr = sline.strip().split()
    atom = np.append(atom, astr[0])
    px = np.append(px, float(astr[1]))
    py = np.append(py, float(astr[2]))
    pz = np.append(pz, float(astr[3]))

# check atom kinds
nat = len(atom)
kinds, atom_count = np.unique(atom, return_counts=True)
formula = ""
for symbol, count in zip(kinds, atom_count):
    formula += f"{symbol}{count} "
print(f"Chemical Formula             : {formula}")

# set pseudo file
for symbol in kinds:
    wt = np.append(wt, mg.Element(symbol).atomic_mass)
    if symbol == 'H':
        ps = np.append(ps, symbol.rjust(2)+pstype_HA)
    else:
        if mg.Element(symbol).is_metal == True:
            ps = np.append(ps, symbol.rjust(2)+pstype_ME)
        else:
            ps = np.append(ps, symbol.rjust(2)+pstype_NM)

# get pseudo file information
# pseudo potensial file infomation class
class get_psinfo:
    
    def __init__(self, psfile):

        self.psinfo = []
        self.ps_type = ""
        self.ps_func = ""
        self.valence = 0.0
        self.wfc_cutoff = 0.0
        self.rho_cutoff = 0.0
    
        try:
            with open(psfile,'r') as f:
                self.psinfo = f.readlines()
            self.ps_type = self.get_pstype()
            self.ps_func = self.get_psfunc()
            self.valence = self.get_valence()
            self.wfc_cutoff = self.get_wfc()
            self.rho_cutoff = self.get_rho()
        except OSError as emsg:
            print(emsg)
            return
        
    def get_pstype(self):
        flag = False
        for sline in self.psinfo:
            astr = sline.strip().split()
            for s in astr:
                if "pseudo_type=" in s.lower():
                    val = s.split("=")[-1]
                    self.ps_type = val.replace('"','')
                    flag = True
                    break
            if len(self.ps_type) != 0:
                #print(self.ps_type)
                break
        return self.ps_type

    def get_psfunc(self):
        for sline in self.psinfo:
            if "functional:" in sline.lower():
                self.ps_func = sline.split(":")[-1].strip()
                break
        return self.ps_func

    def get_valence(self):
        flag = False
        for sline in self.psinfo:
            astr = sline.strip().split()
            for s in astr:
                if "z_valence=" in s.lower():
                    val = s.split("=")[-1]
                    val = val.replace('"','')
                    self.valence = float(val)
                    flag = True
                    break
            if flag == True:
                break
        return self.valence

    def get_wfc(self):
        flag = False
        for sline in self.psinfo:
            astr = sline.strip().split()
            for s in astr:
                if "wfc_cutoff=" in s.lower():
                    val = s.split("=")[-1]
                    val = val.replace('"','')
                    self.wfc_cutoff = float(val)
                    flag = True
                    break
            if flag == True:
                break
        return self.wfc_cutoff

    def get_rho(self):
        flag = False
        for sline in self.psinfo:
            astr = sline.strip().split()
            for s in astr:
                if "rho_cutoff=" in s.lower():
                    val = s.split("=")[-1]
                    val = val.replace('"','')
                    self.rho_cutoff = float(val)
                    flag = True
                    break
            if flag == True:
                break
        return self.rho_cutoff

print(f"Pseudo Potential Informations:")
for psfile in ps:
    if os.path.exists(pseudo_dir+psfile.strip()) == False:
        print(f"Not Found [{pseudo_dir+psfile.strip()}]")
        sys.exit(1)
    ps_info = get_psinfo(pseudo_dir+psfile.strip())
    ps_type = np.append(ps_type, ps_info.ps_type)
    if ps_info.ps_type == "USPP":
        rho_factor = 10
    else:
        rho_factor = 4
    ps_func = np.append(ps_func, ps_info.ps_func)
    valences = np.append(valences, ps_info.valence)
    wfc = np.append(wfc, ps_info.wfc_cutoff)
    rho = np.append(rho, ps_info.rho_cutoff)
    print(f"{psfile[0:2]:2s} [{ps_info.ps_type:4s}, {ps_info.ps_func:18s}]: valence={ps_info.valence:5.1f}, wfc={ps_info.wfc_cutoff:5.1f}, rho={ps_info.rho_cutoff:6.1f}")

# keyword setting
# set CONTROL card
CONTROL_Card = [
    "&CONTROL\n",
    f"    calculation = '{method}'\n",
    "    verbosity = 'low'\n",
    "    restart_mode = 'from_scratch'\n",
    f"    nstep = {nstep}\n",
    "    outdir = './work/'\n",
    "    wfcdir = './work/'\n",
    f"    prefix = '{fbase}'\n",
    f"    etot_conv_thr = {etot_conv_thr}\n",
    f"    forc_conv_thr = {forc_conv_thr}\n",
    "    disk_io = 'low'\n",
    f"    pseudo_dir = '{pseudo_dir}'\n",
    "/\n\n"
]

# set SYSTEM card
nbnd = np.sum(valences * atom_count)
if use_nbnd == True:
    nbnd = nbnd * 1.2
nbnd_str = " Occupied Only " if use_nbnd == False else "Unoccupied +20%"
print(f"nbnd ({nbnd_str:15s})       : {round(nbnd)}")

min_wfc = max(wfc)
ecutwfc = max(min_wfc, ecutwfc_def)
min_rho = max(rho)
calcd_rho = ecutwfc * rho_factor
ecutrho = max(min_rho, ecutrho_def, calcd_rho)
print(f"Cut Off Energy               : wfc = {ecutwfc:.1f}Ry, rho = {ecutrho:.1f}Ry")

SYSTEM_Card = [
    "&SYSTEM\n",
    "    ibrav = 0\n",
    f"    nat = {nat}\n",
    f"    ntyp = {len(kinds)}\n",
    f"!    nbnd = {round(nbnd)}\n",
    "    tot_charge = 0.0\n",
    f"    ecutwfc = {ecutwfc:.1f}\n",
    f"    ecutrho = {ecutrho:.1f}\n",
    "    nosym = .FALSE.\n",
    "    nspin = 1\n",
    f"!    input_dft = '{input_dft}'\n",
    f"    vdw_corr = '{vdw_corr}'\n",
    "/\n\n"
]

# set ELECTRONS card
ELECTRONS_Card = [
    "&ELECTRONS\n",
    f"    electron_maxstep = {electron_maxstep}\n",
    f"    conv_thr = {conv_thr}\n",
    f"    mixing_mode = '{mixing_mode}'\n",
    f"    mixing_beta = {mixing_beta}\n",
    f"    diagonalization = '{diagonalization}'\n",
    "/\n\n"
]

# set IONS card for relax, md, vc-relax or vc-md
IONS_Card = [
    "&IONS\n",
    f"    ion_dynamics = '{ion_dynamics}'\n",
    f"    ion_temperature = '{ion_temperature}'\n",
    f"    tempw = {tempw}\n",
    "/\n\n"
]

# set CELL card for vc-relax or vc-md
CELL_Card = [
    "&CELL\n",
    f"    cell_dynamics = '{cell_dynamics}'\n",
    f"    press = {press}\n",
    f"    press_conv_thr = {press_conv_thr}\n",
    f"    cell_dofree = '{cell_dofree}'\n",
    "/\n\n"
]

# set ATOMIC_SPECIES card
ATOMIC_SPECIES_Card = [
    "ATOMIC_SPECIES\n"
]
for i in range(len(kinds)):
    ATOMIC_SPECIES_Card += [f"   {kinds[i].rjust(2)} {wt[i]:9.6f} {ps[i]}\n"]
ATOMIC_SPECIES_Card += ["\n"]

# set ATOMIC_POSITIONS card
ATOMIC_POSITIONS_Card = [
    "ATOMIC_POSITIONS {angstrom}\n"
]
for i in range(nat):
    symbol = atom[i].rjust(2)
    PX = format(px[i], "10.6f")
    PY = format(py[i], "10.6f")
    PZ = format(pz[i], "10.6f")
    FLAG = f"{flag}  {flag}  {flag}"
    ATOMIC_POSITIONS_Card += [f"   {symbol}  {PX}  {PY}  {PZ}  {FLAG}\n"]
ATOMIC_POSITIONS_Card += ["\n"]

# set K_POINTS card
cell_max = np.max(cell[0:3])
for i in range(3):
    nk[i] = round(cell[i]/cell_max)
try:
    k_point = k_point / nat
except ZeroDivisionError:
    k_point = 100
factor = np.cbrt(k_point / (nk[0]*nk[1]*nk[2]))
for i in range(3):
    nk[i] = round(nk[i]*factor)
    if nk[i] ==0:
        nk[i] = 1
print(f"K Point = {round(k_point): 4d}               : calcd = {nk[0]*nk[1]*nk[2]} [{nk[0]} x {nk[1]} x {nk[2]}]")

K_POINTS_Card = [
    "K_POINTS {automatic}\n",
    f"  {nk[0]} {nk[1]} {nk[2]} 0 0 0\n",
    "\n"
]

# set CELL_PARAMETERS card
CELL_PARAMETERS_Card = [
    "CELL_PARAMETERS {angstrom}\n",
    f"  {V[0][0]:>9.6f} {V[0][1]:>9.6f} {V[0][2]:>9.6f}\n",
    f"  {V[1][0]:>9.6f} {V[1][1]:>9.6f} {V[1][2]:>9.6f}\n",
    f"  {V[2][0]:>9.6f} {V[2][1]:>9.6f} {V[2][2]:>9.6f}\n",
    "\n"
]

# writing file
print(f"Writing File                 : {output}") 
try:
    with open(output,'w') as f:
        for ftxt in CONTROL_Card:
            f.write(ftxt)
        for ftxt in SYSTEM_Card:
            if "nbnd" in ftxt:
                if use_nbnd == True:
                    ftxt = ftxt.replace('!','')
            if "input_dft" in ftxt:
                if use_vdw_corr == False:
                    ftxt = ftxt.replace('!','')
            if "vdw_corr" in ftxt:
                if use_vdw_corr == False:
                    ftxt = "!"+ftxt
            f.write(ftxt)
        for ftxt in ELECTRONS_Card:
            f.write(ftxt)
        if (method=='relax') or (method=='md') or (method=='vc-relax') or (method=='vc-md'):
            for ftxt in IONS_Card:
                f.write(ftxt)
        if (method=='vc-relax') or (method=='vc-md'):
            for ftxt in CELL_Card:
                f.write(ftxt)
        for ftxt in ATOMIC_SPECIES_Card:
            f.write(ftxt)
        for ftxt in ATOMIC_POSITIONS_Card:
            f.write(ftxt)
        for ftxt in K_POINTS_Card:
            f.write(ftxt)
        for ftxt in CELL_PARAMETERS_Card:
            f.write(ftxt)
except OSError as emsg:
    print(emsg)
    sys.exit(1)

- コードの簡易説明


必要環境
・python 3.6以上
・numpy 1.10以上
・pymatgen
・cif2cell

原子座標の読み込み
子プロセスでcif2cell*2を呼び出し、中間ファイルを作成して処理しています。本当は自分で実装したかったのですが、座標変換部分(対称操作)の勉強中なので、ひとまずあるものを利用させていただくということで・・・

余談ですが、pythonでCIFファイルを扱えるライブラリとしてはpymatgenのio.cifモジュール*3と、IUCrがサポートしているPyCifRWのCifFile*4が有名です。
前者は機能が豊富なのですが、ソースコードを見る限り独自の実装となっているようで、旧タイプのキーワードを使用しているCIFファイルに対応していないことが致命的です。一見読み込めたように見えても、インスタンスにアクセスするとエラーを吐いてしまうことがあります。
一方後者は、論文誌に投稿されたものだけあって、古い形式のCIFファイルも問題なく読み込めます。ただ、使用方法の情報が少ないのが欠点です。機能が少ない分は、自分でラッパーを用意するだけの手間で済むので、個人的にはこちらの使用をお勧めします。
そのうち、この記事も書きたいです。

疑ポテンシャルファイル情報読み込みクラス
価電子数とカットオフエネルギーの最低値を自動的に取得したかったので、クラスの実装練習がてらに用意してみました。

    #インスタンス作成
    psinfo = get_psinfo("ファイルパス")
    #アトリビュートへのアクセス
    valence = psinfo.valence     #価電子数
    wfc_min = psinfo.wfc_cutoff  #wfc
    rho_min = psinfo.rho_cutoff  #rho

wfc、rhoのカットオフエネルギーはそれぞれ、25、100Ryを初期値として設定しています。疑ポテンシャルファイルに設定されている最低値がこれよりも大きければ、こちらを採用します。ウルトラソフト型やPAW型はwfcが50Ry程度で設定されていることが多いので、デフォルトを50、200Ryにするべきかもしれません。

テンプレート群
書き込みエラーへの対応のため、一旦配列に格納することにしました。各カード毎に分けたので、多分、編集しやすくなっています。

*1:https://www.misasa.okayama-u.ac.jp/~masami/pukiwiki/index.php?Quantum-Espresso%E3%81%A8Vesta%E3%81%AE%E9%80%A3%E6%90%BA

*2:ネットで検索すると、原著優先ということでだいたいpython2用に誘導されるのですが、python3用がこちらからダウンロードできます。
https://github.com/kmu/cif2cell

*3:https://pymatgen.org/pymatgen.io.cif.html

*4:https://www.iucr.org/resources/cif/software/pycifrw

格子定数と格子ベクトルの相互変換

Quantum ESPRESSO Python

- 事始め


Quantum ESPRESSO[*1]の入力ファイルに格子定数を入力するには、以下の２通りの宣言方法があります[*2]。

・SYSTEMカード内にibravとcelldmを記述
・CELL_PARAMETERSカード内に格子ベクトルを記述

ibravの説明では、どちらもBravais格子に合わせて計算する必要があるように見えますが、後者の方法では三斜晶系の方法を覚えるだけで良いので便利です(その他は三斜晶系の特殊例にすぎません)。
また、計算途中の格子定数を求めたい場合、出力ファイルにはCELL_PARAMETERSしか出力されないため、格子ベクトルから変換する必要があります。
そこで、ここではそれらの相互変換方法を紹介します。

- 格子定数と格子ベクトルの関係


Quantum ESPRESSOにおける格子ベクトルは、a軸を特殊軸とする下記の仕様です[*3]。

\begin{bmatrix} a_x & 0 & 0 \\ b_x & b_y & 0 \\ c_x & c_y & c_z \\ \end{bmatrix}

格子定数をa, b, c, α, β, γとしたとき、格子ベクトルの各要素は下記式でそれぞれ求めることができます。

\begin{aligned} a_x &= a \\ b_x &= b\cos \gamma \\ b_y &= b\sin \gamma \\ c_x &= c\cos \beta \\ c_y &= c\frac{(\cos \alpha - \cos \beta \cos \gamma)}{\sin \gamma} \\ c_z &= \sqrt{c^2 - c_x^2 - c_y^2} \\ \end{aligned}


Pythonのコード例

格子定数は、カフェインとテオフィリンの共結晶(COD:7108131)[*4]から引用しました。

import numpy as np

# cell = [ a, b, c, α, β, γ ]
cell = np.array([7.0081, 8.7904, 13.0856, 95.648, 97.09, 91.667])
V = np.zeros((3,3))

cos_a = np.cos(np.deg2rad(cell[3]))
cos_b = np.cos(np.deg2rad(cell[4]))
cos_c = np.cos(np.deg2rad(cell[5]))
sin_c = np.sin(np.deg2rad(cell[5]))

V[0][0] = cell[0]
V[1][0] = cell[1] * cos_c
V[1][1] = cell[1] * sin_c
V[2][0] = cell[2] * cos_b
V[2][1] = cell[2] * (cos_a - cos_b * cos_c) / sin_c
V[2][2] = np.sqrt(cell[2] ** 2 - V[2][0] ** 2 - V[2][1] ** 2)


sin γで除算するので、本来ならZeroDivisionErrorをチェックしたほうが良いですが、そうなるのはγが0°、または180°の時だけで、そのような結晶は存在しえないため、わざわざ記述する必要はないように思います。


CELL_PARAMETERSを使用する場合は、SYSTEMカード内でibrav=0を指定します[*5]。

入力ファイルの一例

&SYSTEM
 ibrav = 0
/

CELL_PARAMETERS {angstrom}
  7.008100  0.000000  0.000000
 -0.255717  8.786680  0.000000
 -1.615133 -1.335391 12.916695


一方、格子ベクトルから格子定数へは下記式で変換することができます。

\begin{aligned} a &= a_x \\ b &= \sqrt{b_x^2 + b_y^2} \\ c &= \sqrt{c_x^2 + c_y^2 + c_z^2} \\ \alpha &= \cos^{-1} \left( \frac{\overrightarrow{b} \cdot \overrightarrow{c}}{bc} \right) \\ \beta &= \cos^{-1} \left( \frac{\overrightarrow{c} \cdot \overrightarrow{a}}{ca} \right) \\ \gamma &= \cos^{-1} \left( \frac{\overrightarrow{a} \cdot \overrightarrow{b}}{ab} \right) \\ \end{aligned}


Pythonのコード例2

numpy:ndarrayでは、@で内積を計算できます[*6]。

import numpy as np

V = np.array([[1.0, 0.0, 0.0],[1.0, 1.0, 0.0],[1.0, 1.0, 1.0]])

a = V[0][0]
b = np.sqrt(V[1] @ V[1])
c = np.sqrt(V[2] @ V[2])

cos_a = (V[1] @ V[2]) / (b * c)
cos_b = (V[2] @ V[0]) / (c * a)
cos_c = (V[0] @ V[1]) / (a * b)

alpha = np.rad2deg(np.arccos(cos_a))
beta = np.rad2deg(np.arccos(cos_b))
gamma = np.rad2deg(np.arccos(cos_c))


因みに、セル体積は下記式で格子ベクトルから簡単に求めることができます。

V = a_x \times b_y \times c_z


*1:https://www.quantum-espresso.org/

*2:https://www.quantum-espresso.org/Doc/INPUT_PW.html#idm199

*3:https://www.quantum-espresso.org/Doc/INPUT_PW.html#idm1361
この確認を怠ったために、c軸を特殊軸として計算していまい、結果が一致しないと数日悩むことになったのは内緒です。

*4:Crystallography Open Database (COD)
http://www.crystallography.net/cod/index.php

*5:cif2cellの出力を見るとAも指定するのが作法のようですが、なくても問題ありません。

*6:Python 3.5、Numpy 1.10以降で採用された仕様です。

RDKit：分子が描画されないエラー

Python

●事始め

久しぶりにRDKit(1)を使ったプログラムをJupyter Notebookで書いたら、
分子が描画されないエラー(2)が発生し、
色々と奔走してしまったので、その対処方法を備忘録として残しておきます。
結論から先に書いておくと、画像処理ライブラリPillowのバージョンアップが
原因でした。
(1) version 2019.09.1 : Ubuntu 20.04
(2) 以下ようなエラーです。

f:id:Doktor-Heppoko:20210505153225p:plain

●対処方法１[Pillowのバージョン下げ]

この問題は、Pillow 8.0.0以上で発生するので、
素直にそれ未満のバージョンに下げます。
8.0.0未満の最新バージョンは7.2.0です。

Pillow バージョンの確認

pip show pillow

結果例

Name: Pillow
Version: 8.2.0
Summary: Python Imaging Library (Fork)
Home-page: https://python-pillow.org
Author: Alex Clark (PIL Fork Author)
Author-email: aclark@python-pillow.org
License: HPND
Location: /home/ubuntu/.local/lib/python3.8/site-packages
Requires:
Required-by: pytesseract, pyocr, matplotlib, ipymol

Pillow ダウングレード

pip install pillow==7.2.0

結果例

Defaulting to user installation because normal site-packages is not writeable
Collecting pillow==7.2.0
Downloading Pillow-7.2.0-cp38-cp38-manylinux1_x86_64.whl (2.2 MB)
|████████████████████████████████| 2.2 MB 13.7 MB/s
Installing collected packages: pillow
Attempting uninstall: pillow
Found existing installation: Pillow 8.2.0
Uninstalling Pillow-8.2.0:
Successfully uninstalled Pillow-8.2.0
Successfully installed pillow-7.2.0

Pillow 8.0.0からim.tostring()メソッドが削除されたため(3)、
これを使用しているcairoCanvas.pyから呼び出せないことによるエラーです。
一応、version 2.0.0から既にim.tobytes()メソッドへの置き換えが
推奨されていたようです。
(3) https://pillow.readthedocs.io/en/stable/deprecations.html#removed-features
　
　
　 ●対処方法２[cairoCanvas.pyの書き換え]

PythonとLinuxの知識をお持ちの方限定。
初心者の方は環境を破壊してしまう恐れがありますので、お勧めできません。

エラーの原因となっているファイルはcairoCanvas.pyのみで、該当箇所は4つ。
これを書き換えることで対処します。

該当のディレクトリに移動(環境に合わせ、読み替えてください)

cd /usr/lib/python3/dist-pakages/rdkit/Chem/Draw

念のためバックアップを作成

sudo cp cairoCanvas.py cairoCanvas.py.old

テキストエディタで開く(管理者権限でオープンすることを忘れずに！)

sudo gedit carioCanvas.py

該当箇所(Line 147, 151, 194, 197)の修正
(以下は該当部分を含むソースコードの一部)

 class Canvas(CanvasBase):  
   
     def __init__(self,  
                  image=None,  # PIL image  
                  size=None,  
                  ctx=None,  
                  imageType=None,  # determines file type  
                  fileName=None,  # if set determines output file name  
                  ):  
         """  
         Canvas can be used in four modes:  
         1) using the supplied PIL image  
         2) using the supplied cairo context ctx  
         3) writing to a file fileName with image type imageType  
         4) creating a cairo surface and context within the constructor  
         """  
         self.image = None  
         self.imageType = imageType  
         if image is not None:  
             try:  
                 imgd = getattr(image, 'tobytes', image.tostring)("raw", "BGRA")  
             except SystemError:  
                 r, g, b, a = image.split()  
                 mrg = Image.merge("RGBA", (b, g, r, a))  
                 imgd = getattr(mrg, 'tobytes', mrg.tostring)("raw", "RGBA")  
   
             a = array.array('B', imgd)  
             stride = image.size[0] * 4  
             surface = cairo.ImageSurface.create_for_data(a, cairo.FORMAT_ARGB32, image.size[0],  
                                                          image.size[1], stride)  
             ctx = cairo.Context(surface)  
             size = image.size[0], image.size[1]  
             self.image = image  
         elif ctx is None and size is not None:  
             if hasattr(cairo, "PDFSurface") and imageType == "pdf":  
                 surface = cairo.PDFSurface(fileName, size[0], size[1])  
             elif hasattr(cairo, "SVGSurface") and imageType == "svg":  
                 surface = cairo.SVGSurface(fileName, size[0], size[1])  
             elif hasattr(cairo, "PSSurface") and imageType == "ps":  
                 surface = cairo.PSSurface(fileName, size[0], size[1])  
             elif imageType == "png":  
                 surface = cairo.ImageSurface(cairo.FORMAT_ARGB32, size[0], size[1])  
             else:  
                 raise ValueError("Unrecognized file type. Valid choices are pdf, svg, ps, and png")  
             ctx = cairo.Context(surface)  
             ctx.set_source_rgb(1, 1, 1)  
             ctx.paint()  
         else:  
             surface = ctx.get_target()  
             if size is None:  
                 try:  
                     size = surface.get_width(), surface.get_height()  
                 except AttributeError:  
                     size = None  
         self.ctx = ctx  
         self.size = size  
         self.surface = surface  
         self.fileName = fileName  
   
     def flush(self):  
         """temporary interface, must be splitted to different methods,  
         """  
         if self.fileName and self.imageType == 'png':  
             self.surface.write_to_png(self.fileName)  
         elif self.image is not None:  
             # on linux at least it seems like the PIL images are BGRA, not RGBA:  
             if hasattr(self.surface, 'get_data'):  
                 getattr(self.image, 'frombytes', self.image.fromstring)(bytes(self.surface.get_data()),  
                                                                         "raw", "BGRA", 0, 1)  
             else:  
                 getattr(self.image, 'frombytes', self.image.fromstring)(  
                   bytes(self.surface.get_data_as_rgba()), "raw", "RGBA", 0, 1)  
             self.surface.finish()  
         elif self.imageType == "png":  
             if hasattr(self.surface, 'get_data'):  
                 buffer = self.surface.get_data()  
             else:  
                 buffer = self.surface.get_data_as_rgba()  
             return buffer

Pillow 8.0.0未満の環境で実行する予定がないのであれば、
下のような感じに修正するだけでＯＫ 。

147 imgd = image.tobytes("raw","BGRA")  
151 imgd = mrg.tobytes("raw","RGBA")  
   
194 self.image.frombytes(bytes(self.surface.get_data()), "raw", "BGRA", 0, 1)  
197 self.image.frombytes(bytes(self.surface.get_data_as_rgba()), "raw", "RGBA", 0, 1)

getattrメソッドでtostring、fromstringをそれぞれtobytes、frombytesに
置き換えるようコードされているため、今回のようにアトリビュートから
削除されてしまうとエラーになってしまいます。
汎用性を求めるなら、try～except～で対処すべきと思うのですが・・・
最新版のソースコード(4)では、上記コードに置き換えられていますので、
RDKitの方を2021.03.1へバージョンアップするのが確実です。
2021/5/5時点だと、Ubuntu 20.04上のaptコマンドでインストールできる
バージョンは 2019.09.1のため、condaでインストールまたは
ソースコードからの ビルドになります。

(4)https://github.com/rdkit/rdkit/blob/master/rdkit/Chem/Draw/cairoCanvas.py