セミナーI

この講義について

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

講義目標

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

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

講義予定

スケジュール

第01回 <2024-10-01 火> 環境構築
第02回 <2024-10-08 火> シェル操作とEmacs
第03回 <2024-10-22 火> Python演習基本文法
第04回 <2024-10-29 火> Python演習 class
第05回 <2024-11-05 火> Python演習 numpy
第06回 <2024-11-12 火> Python演習 matplotlib
第07回 <2024-11-19 火> Python演習 scipy
第08回 <2024-12-03 火> Python演習 re
第09回 <2024-12-10 火> Python演習 subprocess
第10回 <2024-12-17 火> Python演習 pandas
第11回 <2024-12-24 火> Python演習データ解析
第12回 <2025-01-07 火> C言語演習基本構文1
第13回 <2025-01-14 火> C言語演習基本構文2
第14回 <2025-01-21 火> C言語演習構造化

成績評価

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

セットアップ

はじめに

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

Windows

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

パフォーマンスを優先する
Windows標準の視覚効果を使わないよう設定します。以下のようにたどってパフォーマンス優先に設定してください。
```
コントロールパネル→システム→詳細情報→システムの詳細設定→パフォーマンス
```

タスクバーを整理する

タスクバーによくわからないアイコンが並んでいる場合は、整理したほうが良さそうです。右クリックして、不要なものは表示しないようにしましょう。

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

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

WSL (WindowsでLinux環境を作る)

学生はデータ処理や計算機シミュレーションのためにUNIXのコマンドライン操作に慣れる必要があります。Windows環境を全部消してLinux等のUNIXを導入しても良いですが、いきなりWindows環境を消してしまうとハードルが高くなって挫折してしまう可能性があるので、まずはWindows上でLinux環境を構築します。WSLとLinuxのセットアップは以下のURLを参照してください。
WSLを利用したLinux環境の構築

MacOS

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

macports https://www.macports.org/
brew https://brew.sh/index_ja

macportsしか大窪は使ったことがありません。macportsでは、/opt以下にファイルが配置されます。まずcommand line toolをインストールします。

xcode-select --install

うまくいかない場合は、appleのサイトから、Command Line Toolsを入手して、インストールしておく。

macportsのページからinstallerをダウンロードしてインストールします。これで、パーッケージを管理するportコマンドが使えるようになります。

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

sudo cat /opt/local/etc/macports/sources.conf | sed 's|^rsync://rsync.macports.org/macports/release/tarballs/ports.tar \[default\]|# &|' > /opt/local/etc/macports/sources.conf 
sudo echo "https://distfiles.macports.org/ports.tar.gz [default]" >> /opt/local/etc/macports/sources.conf

編集したら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\
 ispell\
 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\
 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 pandas -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/memo.files/wsl/files/data.tgz | tar xzf - -C ~/

Google Drive

講義で扱うファイルや課題は、Google Driveで配布します。各自、クライアントをインストールして大学アカウントでログインし、PCからGoogle Drive上のファイルにアクセスできるようにします。
https://www.google.com/intl/ja_jp/drive/download/

/mnt/gのマウント

Windows側からGoogle DriveのG:が見えるにもかかわらず、Debian側から /mnt/gでファイルが見えない場合、DebianでG:をマウントするよう/etc/fstabを作成します。

curl -O -u glass https://amorphous.tf.chiba-u.jp/memo.files/wsl/files/fstab
sudo mv fstab /etc

PowerShellを管理者権限で実行して、WSLをshutdownする。

wsl --shutdown

シェルと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を身につければCやFortranを理解するのも非常に早くなるかと思います。世の中のpythonプログラミングは、webアプリケーションやデータベース、それからGUIの開発としてのニーズが多いかと思います。そのためpythonの書籍は、そのニーズに関連した内容を中心にまとめたものが多いです。これは分子シミュレーションではあまり必要ない知識です。そこで研究のためのプログラム技法を効率よく身につけるため、分子シミュレーションに関連したpythonのサンプルプログラムを以下にまとめました。処理結果をvesta等の分子描画ソフトウェアで確認しながらプログラムの学習を進めれば良いかと思います。サンプルプログラムを参考にして箇条書きの処理を自分で作成すれば理解の助けになるでしょう。

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

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

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

引数解析を行う。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足して表示する。

H₂Oの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ファイルを作る

H₂Oの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))

H₂O原子を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]))

指定した原子に一番近い元素だけを表示する。

構造を調べる

H₂Oの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]])

周期境界を跨ぐH₂Oだけを出力する。

H₂Oの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動径分布関数を計算する。

H₂Oの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))

三斜晶系でも行えるように変更する。

CH₄を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()

ファイルを読んでデータをプロットするオプションを追加する。
データの最大値で規格化するオプションを追加する。

ほかに必要な勉強

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

量子力学量子論の基礎の1,2,3章 (スペードマークは必要ありません)
分子動力計算 HOW TO分子シミュレーションの1,2,3章