セミナーI

目次

講義詳細情報

  • 開講時間: 3年後期火曜4限
  • 単位数: 1.0 単位
  • 場所: 工学部2号棟201

講義目標

4年生で研究をスタートさせるために必須となる「コンピュータ」と「無機材料化学」の知識を身につける。予習復習を必ず行って講義に参加すること。主体的な態度で講義に参加することを求めます。疑問に思ったことは自主的に調べ、教員や先輩に積極的に質問することを心がけてください。 なお大窪担当のコンピュータ処理は必須ですので、履修していない(または受講しても理解していない)場合は、受講してください。

課題をこなせないなど計算機に向いていない学生は、別途実験することにします。

講義予定

スケジュール

  • 第1回 セットアップ <2023-10-03 火 14:30>–<2023-10-03 火 16:00>
  • 第2回 シェル操作 <2023-10-10 火 14:30>–<2023-10-10 火 16:00>
  • 第3回 Emacs <2023-10-17 火 14:30>–<2023-10-17 火 16:00>
  • 第4回 Python演習 <2023-10-24 火 14:30>–<2023-10-24 火 16:00>
  • 第5回 Python演習 <2023-11-07 火 14:30>–<2023-11-07 火 16:00>
  • 第6回 Python演習 <2023-11-14 火 14:30>–<2023-11-14 火 16:00>
  • 第7回 Python演習 <2023-11-21 火 14:30>–<2023-11-21 火 16:00>
  • 第8回 Python演習 <2023-12-05 火 14:30>–<2023-12-05 火 16:00>
  • 第9回 Python演習 <2023-12-12 火 14:30>–<2023-12-12 火 16:00>
  • 第10回 Python演習 <2023-12-19 火 14:30>–<2023-12-19 火 16:00>
  • 第11回 Python演習 <2023-12-26 火 14:30>–<2023-12-26 火 16:00>
  • 第12回 C言語演習 <2024-01-09 火 14:30>–<2024-01-09 火 16:00>
  • 第13回 C言語演習 <2024-01-16 火 14:30>–<2024-01-16 火 16:00>
  • 第14回 C言語演習 <2024-01-23 火 14:30>–<2024-01-23 火 16:00>

成績評価

  • 公欠を除いて、3回講義を休むと自動的に不可になります。
  • タイピング速度が規定のレベルに達しない場合は不可になります。
  • 講義前に少なくとも90分以上の予習を行ってください。予習を行っている前提で講義を進めます。
  • 講義にのぞむ主体性、講義中の態度で成績を加点減点します。
  • 予習していなかったり遅刻欠席の場合は成績を減点します。
  • 計算機の演習をクリアしないと、計算機の研究テーマを担当できません。
  • セミナーIの成績は課題で評価します。

セットアップ

はじめに

本文章は新しく配属された学生向けに実験やシミュレーションで得たデータを 整理してグラフを作成し、論文や発表資料で利用するための必要最低限なコン ピュータの知識についてまとめたものです。研究室では実験やシミュレーショ ンから得られる大量のデータを解析する必要があります。Excel等の一般的な 表計算ソフトでは不十分なため、UNIXを使った解析プログラムの使い方やプロ グラミングを習得する必要があります。コマンド操作などで最初は苦労するか と思いますが、使っているうちにだんだんと慣れてきます。計算機を学習する コツは、「楽する方法を常に探求」することです。研究で必要な単純作業はス クリプト化すれば単純ミスなしに行えますし、毎回行うデータ処理は、プログ ラムを書けばコマンド1行で終わります。せっかく10研に配属されたので無機 化学や物理化学だけではなく計算機を使える人間になって修了してください。

Windowsの設定

Microsoft社製OSのWindowsシリーズはWindowsXP以降、見栄えを良くするため にユーザの操作に追随して視覚効果を派手に行うようになっています。図を 作ったり論文を書いたりする際には、無駄な視覚効果は動作が遅くなるだけ なので、パフォーマンスを優先して作業を行うよう設定します。

パフォーマンスを優先する

Windows標準の視覚効果を使わないよう設定します。以下のようにたどってパフォーマンス優先に設定してください。

コントロールパネル→システム→詳細情報→システムの詳細設定→パフォーマンス

WSL (WindowsでLinux環境を作る)

学生はデータ処理や計算機シミュレーションのためにUNIXを使いこなす必要が あります。Windows環境を全部消してLinux等のUNIXを導入しても良いですが、 いきなりWindows環境を消してしまうとハードルが高くなって挫折してしまう 可能性があるので、まずはWindows上でLinux環境を構築します。

wslとlinuxのセットアップは以下のURLを参照してください。
http://chem.tf.chiba-u.jp/gacb10/memo.files/wsl/wsl_linux.html

remapkeyによるキー配列変更 (CtrlとCapsLocksのswap)

unixによるコマンド入力はctrlキーを多用します。ノートパソコンのキーボー ドは、Ctrlキーが左下にあり押し辛いです。インタネット上からrktools.exe に含まれるremapkeyで左小指の左にあるCaps Lockと左下にあるCtrlキーを入 れ替えます。Altキーも多用するのでスペースキーのすぐ左の無変換とAltを 交換した方が使いやすいかと思います。remapkeyを管理者権限で実行して元のキー配列(ベース)から変更 したいキー配列にドラックするとキー配列の変更ができます。変更が完了し たら保存して再起動するとキー配列が入れ替わっています。

キーボードリマップ用プログラム

MacOS

MacOSはベースがUNIXなので、そのままUNIX環境が使えます。各種パッケージは、macportsかbrewを利用します。

  • macports https://www.macports.org/
  • brew https://brew.sh/index_ja macportsしか大窪は使ったことがありません。macportsでは、/opt以下にファイルが配置されます。またコンパイルのためにapple storeからXcodeをインストールしておく必要があります。macportsのページからinstallerをダウンロードしてインストールします。これで、パーッケージを管理するportコマンドが使えるようになります。

macportはrsyncプロトコル(873ポート)rsyncプロトコル(873ポート)でパッケージファイル取得します。学内では、学外への873ポートが閉じられているので、tarballをhttps経由で取得して利用します。sources.confを修正します。

vi /opt/local/etc/macports/sources.conf
#rsync://rsync.macports.org/macports/release/tarballs/ports.tar [default]
https://distfiles.macports.org/ports.tar.gz [default]

編集したらupdateします。

sudo port -d sync 

必要なものをインストールします。参考に大窪のインストールしているものを 晒しておきます。大学院時代から、windowmakerで過ごしております。

sudo port -N install\
 rxvt-unicode\
 rsync\
 windowmaker\
 nkf\
 screen\
 emacs +x11 +dbus\
 gnuplot +x11 -aquaterm -luaterm -pangocairo -wxwidgets\
 ffmpeg\
 ImageMagick\
 texlive-latex-recommended\
 texlive-latex-extra\
 texlive-latex\
 texlive-bin-extra\
 texlive-math-science\
 texlive-publishers\
 texlive-lang-japanese\
 texlive-lang-cjk\
 texlive-plain-generic\
 texlive-fonts-extra\
 latexmk\
 ispell\
 inetutils\
 lookup\
 eblook

MacOS用に大窪の設定ファイルを配布します。

curl -u glass https://amorphous.tf.chiba-u.jp/memo.files/macos/files/dot.tgz | tar xzf - -C ~/

Pythonをインストールします。

curl -O https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-arm64.sh
bash Miniconda3-latest-MacOSX-arm64.sh -b -p $HOME/.miniconda3
conda install numpy scipy matplotlib lmfit -y

MacOS用のセミナーで使うバイナリ他です。

curl -u glass https://amorphous.tf.chiba-u.jp/memo.files/macos/files/local.tgz | tar xzf - -C ~/

セミナー用のデータを配布しています。次のコマンドでファイルを展開してください。

curl -u glass https://amorphous.tf.chiba-u.jp/lecture.files/seminar-I/install/data.tgz | tar xzf - -C ~/

シェルとEmacs演習

UNIX操作

UNIXはターミナルからコマンドを入力してコンピュータを操作(CUI)します。 GUIと比べて最初は難しいですが、慣れると全ての作業が早くできますし、 コンピュータ操作を簡単にスクリプト化できるので、同じような作業をする 場合に早く行えるようなります。コマンドやファイル名は正確に覚えておく 必要がありません。TABによるコマンド補完や履歴機能を使うと最小で入力 ができます。代表的なコマンド編集の操作は以下のとおりです。emacsとほ とんど同じなので丸覚えしましょう。なお矢印キーでも履歴を辿れますが効 率がわるいので全てCtrlを使うようにしてください。"C-"はCtrlキーを押し ながらという意味です。

TAB   コマンドやファイル名の補完。連打で補完候補がでてくる。
C-h   カーソルの左を消す
C-a   カーソルを行頭に移動
C-e   カーソルを行末に移動
C-f   カーソルを右に移動
C-b   カーソルを左に移動
C-k   カーソルから行末までカット
C-p   履歴を一つ前にたどる
C-n   履歴を一つ先にたどる
C-y   カットしたコマンドの挿入

UNIXコマンド

標準的に利用するコマンドをまとめています。 01unix-command.pdf

UNIX文化は小さいコマンドをパイプ等で組み合わせて利用します。例えばls した結果からhogeを含むファイルを探す場合、以下のようになります。

ls | grep hoge

またディレクトリにいくつファイルがあるか数える場合は、以下のようになります。

ls | wc -l

標準的なコマンドは全て覚える必要がありませんが、manコマンドで必要に応 じて調べる癖をつけてください。manコマンドを利用してもよいですが、大抵 の情報はネット上から得られます。大窪に聞いても教えてあげます。コマンド を知らないばっかりに、無駄な単純作業の繰り返しで時間を浪費することだけ は止めてください。

Emacs (高機能エディタ)

テキストファイルを編集するためのソフトウェアをテキストエディタと呼びま す。Windowsのメモ帳やword pad、有償版では秀丸と使用目的は同じです。 UNIXではviエディタが標準で使えます。Emacsは、高機能過ぎてやや重い面も ありますが、プログラムを書いたり論文を作成したり、メールの読み書きまで 全てのテキスト処理に利用できます。カーソルの移動やコピーペースト等、基 本的なことすべてをキーボードで行えます。キーの操作が独特で最初はとっつ きにくいが、使っているうちに慣れるので挫折しないこと。テキストの種類に よって様々なモードが存在します。

起動

emacs

.bashrcでemacsは以下のような.bashrcで以下のようにalias(別名)をつけています。

alias emacs="emacsclient -nc -a \"\""

emacsデーモンが起動していなかったらemacsデーモンを立ち上げ、すでに emacsデーモンが立ち上がっていたらemacsclientがデーモンに接続します。 タスクマネージャーで確認すれば、emacsデーモンが動いているかどうか確 認できるでしょう。

コマンド

emacsを利用するための必要最低限なコマンドを書いています。emacsは何で もできるエディタで世界中で使われていますが、何でもできる反面、操作が 複雑すぎて難しいと感じる人が多くないようです。ただし慣れるとテキスト 処理のすべてを高速で行うことができます。基本コマンドを]]してください。 emacs_reference.pdf

Emacsの生活

Emacs上ですべての作業が完結します。よく使う便利な機能は以下のとおり。

  • eperiodic (周期表)
  • yalatex (latexによる文章作成)
  • lookup (辞書)
  • wanderlust (メール)
  • text-translate (翻訳)
  • ispell (スペルチェック)
  • org-mode (スケジュール、プロジェクト管理)
  • python-mode (pythonのプログラム作成)
  • c-mode (Cのプログラム作成)

Python演習

研究室の方針として実験とシミュレーションの両方を使いこなして研究をすす めてもらいたいと考えています。就職活動や大学院試験等、学生の勉強に費や せる時間は非常に限られていますので、効率よくプログラム言語を身につけて 作成できるようになることが必要です。pythonは非常に可読性が高く容易に習 得できるプログラム言語です。pythonを身につければCやFortranを理解するの も非常に早くなるかと思います。世の中のpythonプログラミングは、webアプ リケーションやデータベース、それからGUIの開発としてのニーズが多いかと 思います。そのためpythonの書籍は、そのニーズに関連した内容を中心にまと めたものが多いです。これは分子シミュレーションではあまり必要ない知識で す。そこで研究のためのプログラム技法を効率よく身につけるため、分子シミュ レーションに関連したpythonのサンプルプログラムを以下にまとめました。処 理結果をvesta等の分子描画ソフトウェアで確認しながらプログラムの学習を 進めれば良いかと思います。サンプルプログラムを参考にして箇条書きの処理 を自分で作成すれば理解の助けになるでしょう。

サンプルプログラムは、クラスと関数を使っていません。実際にプログラムす る時は、クラスと関数(メソッド)を使って読みやすいよう作成してください。

ファイル一式: python-lesson.tgz ※pythonはお手軽script言語です。シミュレーションのpost/pre処理には有用 かもしれませんが、どうしてもスピードに限界があります。計算スピードが必 要な場合は、Cやfortran等のコンパイル言語で作成します。

pythonのインストール

minicondaを使ってインストールしてください。

pythonの基本とnumpy, scipy, matplotlib

コンピュータ処理の資料を参照してください。

http://chem.tf.chiba-u.jp/gacb10/lecture.files/chem_computer/11/11.html

ファイルの読み込みと出力

引数解析を行う。1-1.py

-p で与えた数字に10をかけた値を表示させる。

#!/usr/bin/env python
import argparse

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-i', '--infile', default="", required=True,
                 help='input file')
par.add_argument('-o', '--outfile', default="", required=True,
                 help='output file')
par.add_argument('-f', '--flag', default=False, action='store_true',
                 help='flag')
par.add_argument('-c', '--choices', default="hoge", choices=['hoge', 'huga'],
                 help='hoge or huga')
par.add_argument('-v', '--value', default=0.8, type=float,
                 help='value')

args = par.parse_args()

print("infile:", args.infile)
print("outfile:", args.outfile)
print("frag:", args.flag)
print("choices:", args.choices)
print("value:", args.value)

ファイル名から拡張子を調べる。1-2.py

  • 引数で与えたファイルの拡張子を.hogeに変更(mv)する。
#!/usr/bin/env python
import argparse
import os

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('infile', help='input file')

args = par.parse_args()

basename, ext = os.path.splitext(args.infile)
print("basename:", basename)
print("extension:", ext)

remove symmetryしたcifファイル(P)から格子定数と分率座標を読む。

#!/usr/bin/env python
import argparse
import re
import numpy as np


description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('infile',  help='input cif file')
args = par.parse_args()

L = np.empty(3, dtype=float)
A = np.empty(3, dtype=float)
body = open(args.infile, "r").read()

L[0] = re.search("_cell_length_a +([0-9\.]+)", body).group(1)
L[1] = re.search("_cell_length_b +([0-9\.]+)", body).group(1)
L[2] = re.search("_cell_length_c +([0-9\.]+)", body).group(1)
A[0] = re.search("_cell_angle_alpha +([0-9\.]+)", body).group(1)
A[1] = re.search("_cell_angle_beta +([0-9\.]+)", body).group(1)
A[2] = re.search("_cell_angle_gamma +([0-9\.]+)", body).group(1)

print("a, b, c:", L)
print("alpha, beta, gamma:", A)
  • _pd_phase_nameを抜き出す。
  • alphaに10足して表示する。

H2Oのxyzファイルを出力する

#!/usr/bin/env python
import argparse
import numpy as np


description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-o', '--outfile', default="", required=True,
                 help='output xyz file')

args = par.parse_args()

data = np.array([["O", 0.,  0., 0.],
                 ["H", 1.,  1., 0],
                 ["H", -1., 1., 0]], dtype=object)
data[:, 1:4] = data[:, 1:4] * 1/np.sqrt(3)

o = open(args.outfile, "w")
o.write("%d\n\n" % (len(data)))
for d in data:
    o.write("{:s} {:f} {:f} {:f}\n".format(d[0], d[1], d[2], d[3]))
print("{} was created.".format(args.outfile))
  • H原子だけのxyzファイルを作る

H2Oのxyzファイルからcifファイルをつくる

#!/usr/bin/env python
import argparse
import numpy as np


description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-o', '--outfile', default="", required=True,
                 help='output xyz file')
par.add_argument('-a', '--a', default=10.0, type=float,
                 help='lattice parameter, a')
par.add_argument('-b', '--b', default=10.0, type=float,
                 help='lattice parameter, b')
par.add_argument('-c', '--c', default=10.0, type=float,
                 help='lattice parameter, c')
par.add_argument('-alpha', '--alpha', default=90, type=float,
                 help='lattice parameter, alpha')
par.add_argument('-beta', '--beta', default=90, type=float,
                 help='lattice parameter, beta')
par.add_argument('-gamma', '--gamma', default=90, type=float,
                 help='lattice parameter, gamma')

args = par.parse_args()

CIF = """
_pd_phase_name                         'test'
_cell_length_a                         %f
_cell_length_b                         %f
_cell_length_c                         %f
_cell_angle_alpha                      %f
_cell_angle_beta                       %f
_cell_angle_gamma                      %f
_symmetry_space_group_name_H-M         'P 1'
_symmetry_Int_Tables_number            1

loop_
_symmetry_equiv_pos_as_xyz
   'x, y, z'

loop_
   _atom_site_type_symbol
   _atom_site_fract_x
   _atom_site_fract_y
   _atom_site_fract_z

"""

data = np.array([["O", 0.,  0., 0.],
                 ["H", 1.,  1., 0],
                 ["H", -1., 1., 0]], dtype=object)
data[:, 1:4] = data[:, 1:4]/np.sqrt(2)
a, b, c = args.a, args.b, args.c
alpha = args.alpha*np.pi/180
beta = args.beta*np.pi/180
gamma = args.gamma*np.pi/180
v1 = [a, 0, 0]
v2 = [b*np.cos(gamma),  b*np.sin(gamma), 0]
v3 = [c*np.cos(beta),
      c*(np.cos(alpha)-np.cos(beta)*np.cos(gamma))/np.sin(gamma),
      c*np.sqrt(1+2*np.cos(alpha)*np.cos(beta)*np.cos(gamma)
                - np.cos(alpha)**2-np.cos(beta)**2-np.cos(gamma)**2 /
                np.sin(gamma))]
M = np.array([v1, v2, v3])
M_ = np.linalg.inv(M)
data[:, 1:4] = np.dot(data[:, 1:4], M_)
data[:, 1:4] = data[:, 1:4] - np.min(data[:, 1:4], axis=0)
data[:, 1:4] = data[:, 1:4] + 0.5

o = open(args.outfile, "w")
o.write(CIF % (a, b, c, args.alpha, args.beta, args.gamma))
for d in data:
    o.write("{:s} {:f} {:f} {:f}\n".format(d[0], d[1], d[2], d[3]))
print("{} was created.".format(args.outfile))
  • H2O原子を2つ置いてcifファイルをつくる。

lammpstrjから座標データを読んでxyzファイルを生成する

#!/usr/bin/env python
import argparse
import numpy as np


description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-i', '--infile', default="", required=True,
                 help='input lammpstrj file')
par.add_argument('-o', '--outfile', default="", required=True,
                 help='output xyz file')

args = par.parse_args()

data = np.loadtxt(args.infile, skiprows=9, dtype=object)[:,[2,3,4,5]]

o = open(args.outfile, "w")
o.write("%d\n\n" % (len(data)))
for d in data:
    o.write("{:s} {:s} {:s} {:s}\n".format(d[0], d[1], d[2], d[3]))
print("{} was created.".format(args.outfile))
  • idが3番目から10番の範囲でxyzファイルを作る。
  • lammpstrjをcifファイルに変換する。

lammpstrjから特定の元素を抜き出してxyzファイルをつくる

#!/usr/bin/env python
import argparse
import numpy as np


description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-i', '--infile', default="", required=True,
                 help='input lammpstrj file')
par.add_argument('-o', '--outfile', default="", required=True,
                 help='output xyz file')
par.add_argument('-e', '--element', default="H",
                 help='element')

args = par.parse_args()

data = np.loadtxt(args.infile, skiprows=9, dtype=object)
data[:, 3:6] = data[:, 3:6].astype(float)
ids = np.where(data[:, 2] == args.element)[0]
data = data[ids, :]

o = open(args.outfile, "w")
o.write("{}\n\n".format(len(data)))
for d in data:
    o.write("{:s} {:f} {:f} {:f}\n".format(d[2], d[3], d[4], d[5]))
print("{} was created.".format(args.outfile))
  • x座標が5.0以下のデータだけでxyzファイルをつくる。

結晶構造の操作

cifファイルからxyz座標を求める(対称操作なし)

#!/usr/bin/env python
import argparse
import re
import numpy as np

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-i', '--infile', default="", required=True,
                 help='input cif file')
par.add_argument('-o', '--outfile', default="", required=True,
                 help='output xyz file')

args = par.parse_args()

L = np.empty(3, dtype=float)
A = np.empty(3, dtype=float)
body = open(args.infile, "r").read()
L[0] = re.search("_cell_length_a\s*(\S+)", body).group(1)
L[1] = re.search("_cell_length_b\s*(\S+)", body).group(1)
L[2] = re.search("_cell_length_c\s*(\S+)", body).group(1)
A[0] = re.search("_cell_angle_alpha\s*(\S+)", body).group(1)
A[1] = re.search("_cell_angle_beta\s*(\S+)", body).group(1)
A[2] = re.search("_cell_angle_gamma\s*(\S+)", body).group(1)

data = np.loadtxt(args.infile, skiprows=32, dtype=object)
data = data[:, [7, 2, 3, 4]]
data[:, 1:4] = data[:, 1:4].astype(float)

A = A * np.pi/180.0
a, b, c = L[0], L[1], L[2]
alpha, beta, gamma = A[0], A[1], A[2]
v1 = [a, 0, 0]
v2 = [b*np.cos(gamma),  b*np.sin(gamma), 0]
v3 = [c*np.cos(beta),
      c*(np.cos(alpha)-np.cos(beta)*np.cos(gamma))/np.sin(gamma),
      c*np.sqrt(1+2*np.cos(alpha)*np.cos(beta)*np.cos(gamma)
                - np.cos(alpha)**2-np.cos(beta)**2-np.cos(gamma)**2 /
                np.sin(gamma))]
M = np.array([v1, v2, v3])
data[:, 1:4] = np.dot(data[:, 1:4], M)

o = open(args.outfile, "w")
o.write("{}\n\n".format(len(data)))
for d in data:
    o.write("{:s} {:f} {:f} {:f}\n".format(d[0], d[1], d[2], d[3]))
print("{} was created.".format(args.outfile))
  • 格子定数を変えてxyzファイルをつくる。

cifファイルから2x2x2のスーパーセルを作成する

#!/usr/bin/env python
import argparse
import re
import numpy as np

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-i', '--infile', default="", required=True,
                 help='input cif file')
par.add_argument('-o', '--outfile', default="", required=True,
                 help='output xyz file')
par.add_argument('-nx', '--nx', default=2, type=int,
                 help='nx')
par.add_argument('-ny', '--ny', default=2, type=int,
                 help='ny')
par.add_argument('-nz', '--nz', default=2, type=int,
                 help='nz')

args = par.parse_args()

L = np.empty(3, dtype=float)
A = np.empty(3, dtype=float)
body = open(args.infile, "r").read()
L[0] = re.search("_cell_length_a\s*(\S+)", body).group(1)
L[1] = re.search("_cell_length_b\s*(\S+)", body).group(1)
L[2] = re.search("_cell_length_c\s*(\S+)", body).group(1)
A[0] = re.search("_cell_angle_alpha\s*(\S+)", body).group(1)
A[1] = re.search("_cell_angle_beta\s*(\S+)", body).group(1)
A[2] = re.search("_cell_angle_gamma\s*(\S+)", body).group(1)

data = np.empty((0, 4), dtype=object)
base = np.loadtxt(args.infile, skiprows=32, dtype=object)
base = base[:, [7, 2, 3, 4]]
base[:, 1:4] = base[:, 1:4].astype(float)
data = np.vstack((data, base))
for i in range(args.nx):
    for j in range(args.ny):
        for k in range(args.nz):
            data_ = np.copy(base)
            data_[:, 1:4] = data_[:, 1:4] + [i, j, k]
            data = np.vstack((data, data_))

A = A * np.pi/180.0
a, b, c = L[0], L[1], L[2]
alpha, beta, gamma = A[0], A[1], A[2]
v1 = [a, 0, 0]
v2 = [b*np.cos(gamma),  b*np.sin(gamma), 0]
v3 = [c*np.cos(beta),
      c*(np.cos(alpha)-np.cos(beta)*np.cos(gamma))/np.sin(gamma),
      c*np.sqrt(1+2*np.cos(alpha)*np.cos(beta)*np.cos(gamma)
                - np.cos(alpha)**2-np.cos(beta)**2-np.cos(gamma)**2 /
                np.sin(gamma))]
M = np.array([v1, v2, v3])
data[:, 1:4] = np.dot(data[:, 1:4], M)

o = open(args.outfile, "w")
o.write("{}\n\n".format(len(data)))
for d in data:
    o.write("{:s} {:f} {:f} {:f}\n".format(d[0], d[1], d[2], d[3]))
print("{} was created.".format(args.outfile))
  • [-1,-1,-1]から[2,2,2]のスーパーセルを作成する。

cifファイルから全ての原子ペアの原子間距離リストを計算する

#!/usr/bin/env python
import argparse
import re
import numpy as np

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('infile',  help='input cif file')

args = par.parse_args()

L = np.empty(3, dtype=float)
A = np.empty(3, dtype=float)
body = open(args.infile, "r").read()
L[0] = re.search("_cell_length_a\s*(\S+)", body).group(1)
L[1] = re.search("_cell_length_b\s*(\S+)", body).group(1)
L[2] = re.search("_cell_length_c\s*(\S+)", body).group(1)
A[0] = re.search("_cell_angle_alpha\s*(\S+)", body).group(1)
A[1] = re.search("_cell_angle_beta\s*(\S+)", body).group(1)
A[2] = re.search("_cell_angle_gamma\s*(\S+)", body).group(1)

data = np.loadtxt(args.infile, skiprows=32, dtype=object)
data = data[:, [0, 2, 3, 4]]
data[:, 1:4] = data[:, 1:4].astype(float)

A = A * np.pi/180.0
a, b, c = L[0], L[1], L[2]
alpha, beta, gamma = A[0], A[1], A[2]
v1 = [a, 0, 0]
v2 = [b*np.cos(gamma),  b*np.sin(gamma), 0]
v3 = [c*np.cos(beta),
      c*(np.cos(alpha)-np.cos(beta)*np.cos(gamma))/np.sin(gamma),
      c*np.sqrt(1+2*np.cos(alpha)*np.cos(beta)*np.cos(gamma)
                - np.cos(alpha)**2-np.cos(beta)**2-np.cos(gamma)**2 /
                np.sin(gamma))]
M = np.array([v1, v2, v3])

distance = np.zeros((len(data), len(data)))
for i, d in enumerate(data):
    r = (d[1:4] - data[:, 1:4]).astype(float)
    r = r - np.round(r/L)*L
    r = np.dot(r, M)
    r = np.sqrt(np.sum(r**2, axis=1))
    distance[:, i] = r

s = "({:d}-{:d}) {}-{}"
for i in range(len(distance)):
    for j in range(len(distance)):
        print(s.format(i+1, j+1, data[i, 0], data[j, 0]), end=' ')
        print("{}".format(distance[i, j]))
  • 指定した原子に一番近い元素だけを表示する。

構造を調べる

H2Oのlammpstrj中のp番目の原子からn番目まで近い原子を出力する

#!/usr/bin/env python
import argparse
import numpy as np

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('infile',  help='input lammpstrj file')
par.add_argument('-p', '--p', default=0, type=int)
par.add_argument('-n', '--nth', default=10, type=int)

args = par.parse_args()
lines = open(args.infile, "r").readlines()
L = lines[5:8]
L = np.array([l.split() for l in L], dtype=float)
L = L[:, 1] - L[:, 0]
data = [d.split() for d in lines[9:]]
data = np.array(data, dtype=object)[:, 2:6]
data[:, 1:4] = data[:, 1:4].astype(float)

rlist = (data[args.p, 1:4] - data[:, 1:4]).astype(float)
rlist = rlist - np.round(rlist/L)*L
rlist = np.sum(rlist**2, axis=1)
order = np.argsort(rlist)

for i in range(args.nth):
    print("{}({}):".format(data[order[i], 0], order[i]), end=' ')
    print(rlist[order[i]])
  • 周期境界を跨ぐH2Oだけを出力する。

H2OのlammpstrjからO-O動径分布関数を計算する

#!/usr/bin/env python
import argparse
import numpy as np
import pylab as plt

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('infile',  help='input lammpstrj file')
par.add_argument('-rmin', '--rmin', default=0, type=float)
par.add_argument('-rmax', '--rmax', default=8, type=float)
par.add_argument('-dr', '--dr', default=0.01, type=float)

args = par.parse_args()
lines = open(args.infile, "r").readlines()
L = lines[5:8]
L = np.array([l.split() for l in L], dtype=float)
L = L[:, 1] - L[:, 0]
data = [d.split() for d in lines[9:]]
data = np.array(data, dtype=object)[:, 2:6]
data[:, 1:4] = data[:, 1:4].astype(float)

Odata = data[np.where(data[:, 0] == "O")[0], :]

Orlist = np.zeros((len(Odata), len(Odata)))
for i, o in enumerate(Odata):
    rlist = (o[1:4] - Odata[:, 1:4]).astype(float)
    rlist = rlist - np.round(rlist/L)*L
    rlist = np.sqrt(np.sum(rlist**2, axis=1))
    Orlist[i, :] = rlist

rho = len(Odata)/(L[0]*L[1]*L[2])
Orlist = Orlist.reshape(-1)
bins = np.arange(args.rmin, args.rmax, args.dr)
Ogr = np.zeros((len(bins)-1, 2))
Ogr[:, 0] = bins[:-1] + args.dr*0.5
Ogr[:, 1] = np.histogram(Orlist, bins=bins)[0]
Ogr[:, 1] = Ogr[:, 1]/(4*np.pi*Ogr[:, 0]**2*args.dr)
Ogr[:, 1] = Ogr[:, 1]/rho/len(Odata)

plt.plot(Ogr[1:, 0], Ogr[1:, 1])
plt.show()
  • O-Oのrunning coordination numberを計算する。
  • O-H動径分布関数を計算する。

H2OのlammpstrjからH-O-H角度分布を計算する

#!/usr/bin/env python
import argparse
import numpy as np
import pylab as plt

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('infile', help='input lammpstrj file')
par.add_argument('-amin', '--amin', default=80, type=float)
par.add_argument('-amax', '--amax', default=120, type=float)
par.add_argument('-ar', '--ar', default=1.0, type=float)

args = par.parse_args()
lines = open(args.infile, "r").readlines()
L = lines[5:8]
L = np.array([l.split() for l in L], dtype=float)
L = L[:, 1] - L[:, 0]
data = [d.split() for d in lines[9:]]
data = np.array(data, dtype=object)[:, 2:6]
data[:, 1:4] = data[:, 1:4].astype(float)

Odata = data[np.where(data[:, 0] == "O")[0], :]
vec = np.zeros((len(Odata), 6))
for i, o in enumerate(Odata):
    rlist = (o[1:4] - data[:, 1:4]).astype(float)
    rlist = rlist - np.round(rlist/L)*L
    rlist = np.sum(rlist**2, axis=1)
    order = np.argsort(rlist)
    vec[i, 0:3] = data[order[1], 1:4] - o[1:4]
    vec[i, 3:6] = data[order[2], 1:4] - o[1:4]

a = np.linalg.norm(vec[:, 0:3], axis=1)
b = np.linalg.norm(vec[:, 3:6], axis=1)
ab = np.sum(vec[:, 0:3]*vec[:, 3:6], axis=1)
angles = np.arccos(ab/(a*b))*180/np.pi

bins = np.arange(args.amin, args.amax, args.ar)
hist = np.zeros((len(bins)-1, 2))
hist[:, 0] = bins[:-1] + args.ar*0.5
hist[:, 1] = np.histogram(angles, bins=bins)[0]

plt.plot(hist[:, 0], hist[:, 1])
plt.show()
  • O-O-O角度分布を計算する。

分子構造の作成

周期境界条件(12x12x12)の下で原子間を少なくとも2 Ang.離してランダムに置く

#!/usr/bin/env python
import argparse
import numpy as np

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-o', '--outfile', default="", required=True,
                 help='output xyz file')
par.add_argument('-Lx', '--lx', default=12.0, type=float,
                 help='Lx')
par.add_argument('-Ly', '--ly', default=12.0, type=float,
                 help='Ly')
par.add_argument('-Lz', '--lz', default=12.0, type=float,
                 help='Lz')
par.add_argument('-c', '--cutoff', default=2.0, type=float,
                 help='cutoff')
par.add_argument('-n', '--num', default=100, type=int,
                 help='Lz')

args = par.parse_args()

L = np.array([args.lx, args.ly, args.lz], dtype=float)
base = np.array([["Ar", 0, 0, 0]], dtype=object)
data = base
while len(data) < args.num:
    r = np.random.random(3)*L
    rlist = (r - data[:, 1:4]).astype(float)
    rlist = rlist - np.round(rlist/L)*L
    d = np.sqrt(np.sum(rlist**2, axis=1))
    if np.any(d < args.cutoff) == False:
        data_ = np.copy(base)
        data_[:, 1:4] = r
        data = np.vstack((data, data_))
        print(len(data), "putting ...")

o = open(args.outfile, "w")
o.write("%d\n\n" % (len(data)))
s = "{:s} {:f} {:f} {:f}\n"
for d in data:
    o.write(s.format(d[0], d[1], d[2], d[3]))
print("{} was created.".format(args.outfile))
  • 三斜晶系でも行えるように変更する。

CH4をxyz軸に自由に回転させる

#!/usr/bin/env python
import argparse
import numpy as np

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-o', '--outfile', default="", required=True,
                 help='output xyz file')
par.add_argument('-ch', '--distance', default=1.5, type=float,
                 help='output xyz file')
par.add_argument('-rx', '--rx', default=30, type=float,
                 help='rotate x')
par.add_argument('-ry', '--ry', default=30, type=float,
                 help='rotate x')
par.add_argument('-rz', '--rz', default=30, type=float,
                 help='rotate x')

args = par.parse_args()

data = np.array([["C", 0, 0, 0],
                 ["H", 1, -1, -1],
                 ["H", -1,  1, -1],
                 ["H", -1, -1,  1],
                 ["H", 1,  1,  1]], dtype=object)

data[:, 1:4] = data[:, 1:4] * args.distance * np.sqrt(1.0/3.0)


theta = np.array([args.rx, args.ry, args.rz])*np.pi/180
Rx = np.array([[1, 0, 0],
               [0, np.cos(theta[0]), -np.sin(theta[0])],
               [0, np.sin(theta[0]), np.cos(theta[0])]])

Ry = np.array([[np.cos(theta[1]), 0, np.sin(theta[1])],
               [0, 1, 0],
               [-np.sin(theta[1]), 0,  np.cos(theta[1])]])
Rz = np.array([[np.cos(theta[2]), -np.sin(theta[2]), 0],
               [np.sin(theta[2]), np.cos(theta[2]), 0],
               [0, 0, 1]])

data[:, 1:4] = np.dot(data[:, 1:4], Rx)
data[:, 1:4] = np.dot(data[:, 1:4], Ry)
data[:, 1:4] = np.dot(data[:, 1:4], Rz)

o = open(args.outfile, "w")
o.write("%d\n\n" % (len(data)))
s = "{:8s} {:18.10f} {:18.10f} {:18.10f}\n"
for d in data:
    o.write(s.format(d[0], d[1], d[2], d[3]))
print("{} was created.".format(args.outfile))
  • 3.0xを軸として30°回転させる。

頂点原子から1 Ang.離して四面体構造の中心から頂点方向に原子を置く(Hのtermination)

#!/usr/bin/env python
import argparse
import numpy as np

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-o', '--outfile', default="", required=True,
                 help='output xyz file')
par.add_argument('-ch', '--ch_distance', default=1.5, type=float,
                 help='C-H distance')
par.add_argument('-ca', '--ca_distance', default=0.8, type=float,
                 help='H-H distance')
par.add_argument('-rx', '--rx', default=0, type=float,
                 help='rotate x')
par.add_argument('-ry', '--ry', default=0, type=float,
                 help='rotate x')
par.add_argument('-rz', '--rz', default=0, type=float,
                 help='rotate x')

args = par.parse_args()

data = np.array([["C", 0, 0, 0],
                 ["H", 1, -1, -1],
                 ["H", -1,  1, -1],
                 ["H", -1, -1,  1],
                 ["H", 1,  1,  1]], dtype=object)

data[:, 1:4] = data[:, 1:4] * args.ch_distance * np.sqrt(1.0/3.0)

# https://en.wikipedia.org/wiki/Rotation_matrix
theta = np.array([args.rx, args.ry, args.rz])*np.pi/180
Rx = np.array([[1, 0, 0],
               [0, np.cos(theta[0]), -np.sin(theta[0])],
               [0, np.sin(theta[0]), np.cos(theta[0])]])

Ry = np.array([[np.cos(theta[1]), 0, np.sin(theta[1])],
               [0, 1, 0],
               [-np.sin(theta[1]), 0,  np.cos(theta[1])]])
Rz = np.array([[np.cos(theta[2]), -np.sin(theta[2]), 0],
               [np.sin(theta[2]), np.cos(theta[2]), 0],
               [0, 0, 1]])

data[:, 1:4] = np.dot(data[:, 1:4], Rx)
data[:, 1:4] = np.dot(data[:, 1:4], Ry)
data[:, 1:4] = np.dot(data[:, 1:4], Rz)

c_vec = (data[1:5, 1:4] - data[0, 1:4]).astype(float)
c_vec = c_vec/np.linalg.norm(c_vec)

Hadd = np.empty((4, 4), dtype=object)
Hadd[:, 0] = "H"
Hadd[:, 1:4] = c_vec * args.ca_distance*2 + data[1:5, 1:4]
data = np.vstack((data, Hadd))

o = open(args.outfile, "w")
o.write("{}\n\n".format(len(data)))
s = "{:8s} {:18.10f} {:18.10f} {:18.10f}\n"
for d in data:
    o.write(s.format(d[0], d[1], d[2], d[3]))
print("{} was created.".format(args.outfile))
  • C-H-Hの角度が10°になるようにHを置く

周期境界条件に基づいてメタンをランダムに回転させ1 Ang.以上離してcifファイルを作る

#!/usr/bin/env python
import argparse,re
import numpy as np

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-o', '--outfile', default="", required=True,
                 help='output xyz file')
par.add_argument('-ch', '--distance', default=1.5, type=float,
                 help='output xyz file')
par.add_argument('-s', '--seed', default=777, type=int,
                 help='seed for random number generator')
par.add_argument('-Lx', '--lx', default=12.0, type=float,
                 help='Lx')
par.add_argument('-Ly', '--ly', default=12.0, type=float,
                 help='Ly')
par.add_argument('-Lz', '--lz', default=12.0, type=float,
                 help='Lz')
par.add_argument('-n', '--num', default=10, type=int,
                 help='number of methan')
par.add_argument('-c', '--cutoff', default=2.0, type=float,
                 help='cutoff')

args = par.parse_args()

methan = np.array([["C", 0, 0, 0],
                  ["H", 1, -1, -1],                   
                  ["H",-1,  1, -1],
                  ["H",-1, -1,  1],
                  ["H", 1,  1,  1]], dtype=object)

methan[:,1:4] = methan[:,1:4] * args.distance * np.sqrt(1.0/3.0)
L = np.array([args.lx, args.ly, args.lz])
M = np.array([[L[0], 0, 0],
              [0, L[1], 0],
              [0, 0, L[2]]])
M_ = np.linalg.inv(M)
methan[:,1:4] = np.dot(methan[:,1:4], M_)

data = np.copy(methan)
while len(data)/len(methan) < args.num:
    methan_ = np.copy(methan)
    theta = np.random.random(3)
    theta = np.random.random(3)
    Rx = np.array([[1, 0, 0],
                   [0,np.cos(theta[0]),-np.sin(theta[0])],
                   [0,np.sin(theta[0]), np.cos(theta[0])]])
    Ry = np.array([[np.cos(theta[1]), 0, np.sin(theta[1])],
                   [0, 1, 0],
                   [-np.sin(theta[1]), 0,  np.cos(theta[1])]])
    Rz = np.array([[np.cos(theta[2]),-np.sin(theta[2]), 0],
                   [np.sin(theta[2]), np.cos(theta[2]), 0],
                   [0,0,1]])
    methan_[:,1:4] = np.dot(methan_[:,1:4], Rx)
    methan_[:,1:4] = np.dot(methan_[:,1:4], Ry)
    methan_[:,1:4] = np.dot(methan_[:,1:4], Rz)
    methan_[:,1:4] = methan_[:,1:4] + np.random.random(3)
    flag = 0
    for a in methan_:
        rlist = (a[1:4] - data[:,1:4]).astype(float)
        rlist = rlist - np.round(rlist)
        rlist = np.dot(rlist, M)
        rlist = np.sqrt(np.sum(rlist ** 2, axis=1))
        if np.all(rlist >  args.cutoff) == False:
            flag = 1
            break
    if flag == 0:
        data = np.vstack((data, methan_))
        print("Putting methan:", len(data)/len(methan))

data[:,1:4] = np.dot(data[:,1:4], M)

o = open(args.outfile, "w")
o.write("%d\n\n" % (len(data)))
s = "{:8s} {:18.10f} {:18.10f} {:18.10f}\n"
for d in data:
    o.write(s.format(d[0], d[1], d[2], d[3]))
print("{} was created.".format(args.outfile))
  • 三斜晶系でも行えるように変更する。

最適化

多項式でフィッテングする

#!/usr/bin/env python
import argparse
import numpy as np
import pylab as plt

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-i', '--datfile', default="polynominal.dat",
                 help='input data file')
par.add_argument('-n', '--polyn', default=3, type=int)

args = par.parse_args()

exp = np.loadtxt(args.datfile)
p = np.polyfit(exp[:,0], exp[:,1], args.polyn)
func =  np.poly1d(p)
fit_x = np.linspace(np.min(exp[:,0]), np.max(exp[:,0]), 1000)
fit_y = np.polyval(p, fit_x)

plt.plot(exp[:,0], exp[:,1], "o")
plt.plot(fit_x, fit_y, "-")
plt.show()
  • 3次と2次でフィットした時の残差を計算する。

ガウス関数でピーク分離を行う

#!/usr/bin/env python
import argparse
import numpy as np
import pylab as plt
from scipy.optimize import least_squares

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-i', '--datfile', default="spectrum.dat",
                 help='input data file')

args = par.parse_args()

def func(p, x):
    y = np.zeros(len(x))
    y = y + p[0]*np.exp(-(x-p[1])**2/p[2])
    y = y + p[3]*np.exp(-(x-p[4])**2/p[5])
    return y

def error(param):
    return data[:,1] - func(param, data[:,0])

data = np.loadtxt(args.datfile)

# I, x0, w    
p0 = np.array([[0.5, 400, 6000],
               [0.2, 600, 6000]]) 
p0 = p0.reshape(-1)
out = least_squares(error, p0).x

fit_x = np.linspace(np.min(data[:,0]), np.max(data[:,0]), 1000)
fit_y = func(out, fit_x)
print("I1, u1, w1 =", out[0], out[1], out[2])
print("I2, u2, w2 =", out[3], out[4], out[5])

plt.plot(data[:,0], data[:,1], "o")
plt.plot(fit_x, fit_y, "-")
plt.show()
  • ベースラインをパラメータとして追加しフィッテングを行う。
  • 半値幅を固定してフィッテングを行う。
  • 各ピークの面積比を求める。
  • 強度が負にならないように制限をかける。

おまけ

データをプロットする

#!/usr/bin/env python
import argparse
import numpy as np
import pylab as plt
from matplotlib.ticker import AutoMinorLocator

description = """This is a test program"""
par = argparse.ArgumentParser(description=description)

par.add_argument('-o', '--outfile', default=None, type=str)

par.add_argument('-fs', '--figsize', nargs=2, default=[6, 5.5], type=float)
par.add_argument('-yl', '--ylabel', default="ylabel", type=str)
par.add_argument('-xl', '--xlabel', default="xlabel", type=str)

par.add_argument('-mn', '--mtick', default=None, type=int)
par.add_argument('-xr', '--xrange', nargs=2, default=None, type=float)
par.add_argument('-yr', '--yrange', nargs=2, default=None, type=float)
par.add_argument('-tt', '--legends', nargs="+", default=None, type=str)
par.add_argument('-np', '--noplot', default=False, action='store_true')
args = par.parse_args()

x = np.linspace(-10, 10, 1000)
ya = np.exp(-x**2)+np.random.random(len(x))*0.1
yb = 0.2*np.exp(-(x-5)**2/0.8)+np.random.random(len(x))*0.1
yc = ya + yb + 0.04

fig, ax = plt.subplots(figsize=args.figsize)
ax.plot(x, ya)
ax.plot(x, yb)
ax.plot(x, yc)
ax.set_xlabel(args.xlabel)
ax.set_ylabel(args.ylabel)

if args.mtick != None:
    minor_locator = AutoMinorLocator(args.mtick)
    ax.yaxis.set_minor_locator(minor_locator)
    minor_locator = AutoMinorLocator(args.mtick)
    ax.xaxis.set_minor_locator(minor_locator)

if args.xrange != None: ax.set_xlim(args.xrange)
if args.yrange != None: ax.set_ylim(args.yrange)
if args.legends != None: ax.legend(args.legends, loc="best")
if args.outfile != None: fig.savefig(args.outfile)
if args.noplot == False: plt.show()
  • ファイルを読んでデータをプロットするオプションを追加する。
  • データの最大値で規格化するオプションを追加する。

計算機テスト

時間内にプログラムが作成できるテストします。

無機化学基礎

テキスト

入門固体化学 の1,2章と英語テキスト1,2章を予習しておくこと

無機化学テスト

  • 講義でとりあつかった無機材料化学の基本的な内容を問います。
  • 講義の前提にある大学低学年の化学の知識も問います。
  • 90分間の筆記試験です。
  • A4用紙1枚だけ持ち込み可です。

ほかに必要な勉強

線形代数と関連づけて計算科学と量子力学に関連する考え方を身につけてください。第一原理計算やNMRを理解する上で必要になります。以下は私の独断でわかりやすいと思った参考図書です。研究室に他にもたくさん本があるので、自分にあったものを選択して勉強することをおすすめします。簡単なプログラミングとシェルスクリプトも覚えてください。初めて学ぶプログラミング言語はpythonを身につけてください。

著者: 大窪 貴洋 (千葉大院工) ohkubo.takahiro@faculty.chiba-u.jp

Created: 2023-10-11 水 17:22

Validate