kivantium活動日記

プログラムを使っていろいろやります

TwitterのStreaming APIについて

2018年8月にTwitterのUser Streams APIが廃止されてしまい、タイムラインのツイートをリアルタイムで取得することができなくなってしまいました。 forest.watch.impress.co.jp

それに伴い、User Streams APIを使っていたリプライによるアイコン変更スクリプトなどが動かなくなっていました。

昨日から開発を始めた画像収集スクリプトを作るためにTweepyのドキュメントを眺めていたら、Streaming APIというものが現在も存在していることに気がついたのでこれを紹介します。

(2015年からTweepyに実装されていたのに今日まで気づかなかった……。ドキュメントはちゃんと読むべき。)

statuses/filterの仕様

2020年4月15日現在、Twitter APIの紹介ページでタイトルにrealtimeとつくものにはFilter realtime TweetsSample realtime Tweetsの2つがあります。このうち、無料ユーザーが使えるのはstatuses/filterだけっぽいです。

developer.twitter.com

このAPIを使うと特定の条件を満たすツイートをリアルタイムで取得できます。パラメータは5つあります。

  • follow: ユーザーIDのリスト(上限5000個)
  • track: 検索キーワード(上限400個)
  • locations: ツイートされた場所(上限25個)
  • delimited: メッセージを長さで区切るかどうか
  • stall_warnings: stall warningsの有無

follow, track, locationsから複数を指定した場合、ORでつないだものとみなされます。

followで指定したユーザーについて、以下を満たすツイートが取得できます。(カッコ内は原文)

  • 指定したユーザーによるツイート (Tweets created by the user.)
  • 指定したユーザーがリツイートしたツイート (Tweets which are retweeted by the user.)
  • 指定したユーザーの任意のツイートに対するリプライ (Replies to any Tweet created by the user.)
  • 指定したユーザーの任意のツイートのリツイート (Retweets of any Tweet created by the user.)
  • リプライボタンを押さずに行われた手動リプライ (Manual replies, created without pressing a reply button (e.g. “@twitterapi I agree”).)

しかし、以下のツイートは含まれません。

  • 指定したユーザーにメンションしているツイート (Tweets mentioning the user (e.g. “Hello @twitterapi!”).)
  • リツイートボタンを押さずに行われた手動リツイート (Manual Retweets created without pressing a Retweet button (e.g. “RT @twitterapi The API is great”).)
  • 非公開アカウントによるツイート (Tweets by protected users.)

特に最後の「非公開アカウントによるツイート」の制約は厳しく、実験した限りではフォローしている鍵アカウントのツイートであっても取得することができません

Tweepyでの使い方

Streaming With Tweepyを読んでください。

followに自分がフォローしているアカウントと自分を指定して、取得できたツイートに含まれる画像を保存するスクリプトの例を示します。

# -*- coding: utf_8 -*-

import os
import sys
import urllib.request
from urllib.parse import urlparse
import tweepy


class MyStreamListener(tweepy.StreamListener):
    def on_status(self, status):
        if 'media' in status.entities:
            for media in status.extended_entities['media']:
                media_url = media['media_url']
                filename = os.path.basename(urlparse(media_url).path)
                if not os.path.exists(filename):
                    try:
                        urllib.request.urlretrieve(media_url, filename)
                        print("Saved :", filename)
                    except IOError:
                        print("Failed to save the image from", media_url, file=sys.stderr)


consumer_key = 'xxxxxxxxxxxxxxxxxxxxxxxxx'
consumer_secret = 'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
access_token = 'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
access_secret = 'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'

auth = tweepy.OAuthHandler(consumer_key, consumer_secret)
auth.set_access_token(access_token, access_secret)
api = tweepy.API(auth)

followee_ids = api.friends_ids(screen_name=api.me().screen_name)
watch_list = [str(user_id) for user_id in followee_ids]
watch_list.append(str(api.me().id))

# followには5000人までしか指定できない
# https://developer.twitter.com/en/docs/tweets/filter-realtime/api-reference/post-statuses-filter
assert(len(watch_list) <= 5000)

myStreamListener = MyStreamListener()
myStream = tweepy.Stream(auth=api.auth, listener=myStreamListener)
myStream.filter(follow=watch_list)

Tweepyで画像を収集する (WIP)

機械学習をするためにはデータが必要なので、頑張って集める必要があります。 とりあえずTweepyで取れるだけ取る方法をいろいろ考えていきます。

Cursorを使う

Tweepyが標準で用意しているCursorというクラスを使ってタイムラインからRate Limitの許す限り画像を収集するコードは次のようになりました。

# -*- coding: utf_8 -*-

import os
import tweepy
import urllib.request
from urllib.parse import urlparse
from pprint import pprint


consumer_key = 'xxxxxxxxxxxxxxxxxxxxxxxxx'
consumer_secret = 'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
access_token = 'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
access_secret = 'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'

auth = tweepy.OAuthHandler(consumer_key, consumer_secret)
auth.set_access_token(access_token, access_secret)
api = tweepy.API(auth, wait_on_rate_limit=True, wait_on_rate_limit_notify=True)

for tweet in tweepy.Cursor(api.home_timeline).items():
    if 'media' in tweet.entities:
        for media in tweet.extended_entities['media']:
            media_url = media['media_url']
            filename = os.path.basename(urlparse(media_url).path)
            try:
                urllib.request.urlretrieve(media_url, filename)
                print("Saved :", filename)
            except IOError:
                print("Failed to save the image from", media_url)

しかし、GET statuses/home_timelineは15分に15回しか呼び出すことができません(ドキュメント)。1回あたり20ツイート取得できるので15分で300件のツイートしか取得できないことになります。(countというパラメータがありますが、実験してみるとcountに大きい値を指定するとその分呼べる回数が減るようで、止まるまでに取得できるツイート数はあまり変わりませんでした。

タイムラインに常駐してツイートを収集しつづけるタイプのプログラムを書く必要がありますが、それはCursorでは実現できなさそうなので真面目に考える必要がありそうです。

友利奈緒判定botのときもRate Limitに悩まされたし、一生Rate Limitと戦っている……

明日以降のためのメモ

Illustration2VecをONNX経由で使う

趣味プロジェクトでIllustration2Vecを使いたくなったのですが、これは2015年の論文なのでモデルをCaffeかChainerで使うことになっています。 github.com

残念ながらCaffeもChainerも既に開発が終了しているため、Illustration2VecのモデルをONNXという共通フォーマットに変換して今後も使えるようにしました。 利用方法だけ知りたい人は「モデルの変換」を飛ばして「使い方」を見てください。

モデルの変換

まずはオリジナルのIllustration2Vecのモデルをダウンロードします。以下を実行するとCaffeのモデルがダウンロードできます。

git clone https://github.com/rezoo/illustration2vec.git
cd illustration2vec
./get_models.sh

このスクリプトでは特徴抽出モデルのprototxtがダウンロードできなかったので、Illustration2VecのInternet ArchiveからNetwork configuration file (feature vectors) illust2vec.prototxtを追加でダウンロードしました。

必要なライブラリをインストールします。

pip install onnx coremltools onnxmltools

以下のPythonスクリプトを実行すると、タグ予測モデルのillust2vec_tag_ver200.onnxと特徴ベクトル抽出モデルのillust2vec_ver200.onnxが作成されます。

import os
import onnx
import coremltools
import onnxmltools

# CaffeモデルをONNX形式で読み込む関数
# https://github.com/onnx/onnx-docker/blob/master/onnx-ecosystem/converter_scripts/caffe_coreml_onnx.ipynb


def caffe_to_onnx(proto_file, input_caffe_path):
    output_coreml_model = 'model.mlmodel'  # 中間ファイル名
    # 中間ファイルが既に存在したら例外を送出
    if os.path.exists(output_coreml_model):
        raise FileExistsError('model.mlmodel already exists')

    # CaffeのモデルをCore MLに変換
    coreml_model = coremltools.converters.caffe.convert(
        (input_caffe_path, proto_file))
    # Core MLモデルを保存
    coreml_model.save(output_coreml_model)
    # Core MLモデルを読み込む
    coreml_model = coremltools.utils.load_spec(output_coreml_model)
    # Core MLモデルをONNXに変換
    onnx_model = onnxmltools.convert_coreml(coreml_model)
    # Core MLモデルを削除
    os.remove(output_coreml_model)
    return onnx_model


# タグ予測モデルの変換・保存
onnx_tag_model = caffe_to_onnx(
    'illust2vec_tag.prototxt', 'illust2vec_tag_ver200.caffemodel')
onnxmltools.utils.save_model(onnx_tag_model, 'illust2vec_tag_ver200.onnx')

# 特徴ベクトル抽出モデルの変換・保存
onnx_model = caffe_to_onnx('illust2vec.prototxt',
                           'illust2vec_ver200.caffemodel')
# encode1レイヤーをONNXから利用できるようにする
# https://github.com/microsoft/onnxruntime/issues/2119
intermediate_tensor_name = 'encode1'
intermediate_layer_value_info = onnx.helper.ValueInfoProto()
intermediate_layer_value_info.name = intermediate_tensor_name
onnx_model.graph.output.extend([intermediate_layer_value_info])
onnx.save(onnx_model, 'illust2vec_ver200.onnx')

使い方

上のようにしてONNX形式に変換したモデルと、それを利用するためのコードを用意しました。 github.com

ONNX形式を他のフレームワークで読み込んで実行してもいいのですが、ONNX RuntimeというMicrosoft製のパフォーマンスを重視した推論専用のライブラリがあったのでこれを使うことにしました。 UbuntuでONNX RuntimeをCPU向けにインストールするコマンドは以下の通りです。

sudo apt install libgomp1
pip install onnxruntime

例を実行する前にコードとpre-trainedモデルのダウンロードを行ってください。

git clone https://github.com/kivantium/illustration2vec.git
cd illustration2vec
./get_onnx_models.sh

タグ予測

コード

import i2v
from PIL import Image
from pprint import pprint

illust2vec = i2v.make_i2v_with_onnx(
    "illust2vec_tag_ver200.onnx", "tag_list.json")

img = Image.open("images/miku.jpg")
pprint(illust2vec.estimate_plausible_tags([img], threshold=0.5))

入力

https://github.com/kivantium/illustration2vec/blob/master/images/miku.jpg?raw=true

Hatsune Miku (初音ミク), © Crypton Future Media, INC., http://piapro.net/en_for_creators.html. This image is licensed under the Creative Commons - Attribution-NonCommercial, 3.0 Unported (CC BY-NC).

出力

[{'character': [('hatsune miku', 0.9999994039535522)],
  'copyright': [('vocaloid', 0.9999999403953552)],
  'general': [('thighhighs', 0.9956372976303101),
              ('1girl', 0.9873461723327637),
              ('twintails', 0.9812833666801453),
              ('solo', 0.9632900953292847),
              ('aqua hair', 0.9167952537536621),
              ('long hair', 0.8817101716995239),
              ('very long hair', 0.8326565027236938),
              ('detached sleeves', 0.7448851466178894),
              ('skirt', 0.6780778169631958),
              ('necktie', 0.560835063457489),
              ('aqua eyes', 0.5527758598327637)],
  'rating': [('safe', 0.9785730242729187),
             ('questionable', 0.02053523063659668),
             ('explicit', 0.0006299614906311035)]}]

Chainer版とほとんど同じ結果が出力されました。Chainerではこの処理に6秒かかっていましたが、onnx-runtimeだと2秒で実行できたのでたしかにパフォーマンスにも優れているようです(ChainerではCaffeのモデルを変換する手間が掛かっているので1枚を処理する時間で比較するのは公平ではないですが)。

特徴ベクトルの抽出

コード

import i2v
from PIL import Image

illust2vec = i2v.make_i2v_with_onnx("illust2vec_ver200.onnx")

img = Image.open("images/miku.jpg")

# extract a 4,096-dimensional feature vector
result_real = illust2vec.extract_feature([img])
print("shape: {}, dtype: {}".format(result_real.shape, result_real.dtype))
print(result_real)

# i2v also supports a 4,096-bit binary feature vector
result_binary = illust2vec.extract_binary_feature([img])
print("shape: {}, dtype: {}".format(result_binary.shape, result_binary.dtype))
print(result_binary)

先ほどと同じ入力に対する出力

shape: (1, 4096), dtype: float32
[[ 7.474596    3.6860986   0.537967   ... -0.14563629  2.7182112
   7.3140917 ]]
shape: (1, 512), dtype: uint8
[[246 215  87 107 249 190 101  32 187  18 124  90  57 233 245 243 245  54
  229  47 188 147 161 149 149 232  59 217 117 112 243  78  78  39  71  45
  235  53  49  77  49 211  93 136 235  22 150 195 131 172 141 253 220 104
  163 220 110  30  59 182 252 253  70 178 148 152 119 239 167 226 202  58
  179 198  67 117 226  13 204 246 215 163  45 150 158  21 244 214 245 251
  124 155  86 250 183  96 182  90 199  56  31 111 123 123 190  79 247  99
   89 233  61 105  58  13 215 159 198  92 121  39 170 223  79 245  83 143
  175 229 119 127 194 217 207 242  27 251 226  38 204 217 125 175 215 165
  251 197 234  94 221 188 147 247 143 247 124 230 239  34  47 195  36  39
  111 244  43 166 118  15  81 177   7  56 132  50 239 134  78 207 232 188
  194 122 169 215 124 152 187 150  14  45 245  27 198 120 146 108 120 250
  199 178  22  86 175 102   6 237 111 254 214 107 219  37 102 104 255 226
  206 172  75 109 239 189 211  48 105  62 199 238 211 254 255 228 178 189
  116  86 135 224   6 253  98  54 252 168  62  23 163 177 255  58  84 173
  156  84  95 205 140  33 176 150 210 231 221  32  43 201  73 126   4 127
  190 123 115 154 223  79 229 123 241 154  94 250   8 236  76 175 253 247
  240 191 120 174 116 229  37 117 222 214 232 175 255 176 154 207 135 183
  158 136 189  84 155  20  64  76 201  28 109  79 141 188  21 222  71 197
  228 155  94  47 137 250  91 195 201 235 249 255 176 245 112 228 207 229
  111 232 157   6 216 228  55 153 202 249 164  76  65 184 191 188 175  83
  231 174 158  45 128  61 246 191 210 189 120 110 198 126  98 227  94 127
  104 214  77 237  91 235 249  11 246 247  30 152  19 118 142 223   9 245
  196 249 255   0 113   2 115 149 196  59 157 117 252 190 120  93 213  77
  222 215  43 223 222 106 138 251  68 213 163  57  54 252 177 250 172  27
   92 115 104 231  54 240 231  74  60 247  23 242 238 176 136 188  23 165
  118  10 197 183  89 199 220  95 231  61 214  49  19  85  93  41 199  21
  254  28 205 181 118 153 170 155 187  60  90 148 189 218 187 172  95 182
  250 255 147 137 157 225 127 127  42  55 191 114  45 238 228 222  53  94
   42 181  38 254 177 232 150  99]]

Chainer版と同じbinary vectorが出力されていました。

次回はこれを使ってイラストの機械学習をします。

DjangoでTwitter認証を行う

前回の記事のつづきです。 kivantium.hateblo.jp

今回はTwitterでのOAuth認証を実装して、タイムラインを読み込むところまで進めます。

ライブラリのインストール

認証に必要なsocial-auth-app-djangoというライブラリをインストールします。前回作ったプロジェクトのルートディレクトリで以下を実行します。

pipenv shell
pipenv install social-auth-app-django

Twitterアプリケーションの作成

Twitter Developerのページからアプリを作成して、API keyとAPI secret keyを入手します。他に解説記事がたくさんあると思うので詳しくは書きません。

コールバックURLの追加

Twitterアプリケーションに適切なコールバックURLを設定しないと403エラーが出ます。今回のディレクトリ構成だと、ローカルで実行する場合は

http://localhost:8000/complete/twitter/

Herokuで実行する場合は

https://<appname>.herokuapp.com/complete/twitter/

を追加する必要があります。複数追加できるので両方追加しておくといいと思います。

設定の変更

playground/settings.pyを以下のように変更します。

 (略)

 INSTALLED_APPS = [
     'django.contrib.sessions',
     'django.contrib.messages',
     'django.contrib.staticfiles',
+    'hello',
+    'social_django',
 ]
 
 (略)

 TEMPLATES = [
     {
         'BACKEND': 'django.template.backends.django.DjangoTemplates',
         'DIRS': [],
         'APP_DIRS': True,
         'OPTIONS': {
             'context_processors': [
                 'django.template.context_processors.debug',
                 'django.template.context_processors.request',
                 'django.contrib.auth.context_processors.auth',
                 'django.contrib.messages.context_processors.messages',
+                'social_django.context_processors.backends',
+                'social_django.context_processors.login_redirect',
             ],
         },
     },
 ]
 
 (略)
 
 AUTH_PASSWORD_VALIDATORS = [
 
 LANGUAGE_CODE = 'en-us'
 
-TIME_ZONE = 'UTC'
+TIME_ZONE = 'Asia/Tokyo'
 
 USE_I18N = True
 
 (略)
 
 STATIC_URL = '/static/'
 
 STATICFILES_STORAGE = 'whitenoise.storage.CompressedManifestStaticFilesStorage'
 
+AUTHENTICATION_BACKENDS = [
+    'social_core.backends.twitter.TwitterOAuth',
+    'django.contrib.auth.backends.ModelBackend',
+]
+
+SOCIAL_AUTH_LOGIN_REDIRECT_URL = '/'
+
 # Use different settings in local environment and Heroku
 # https://qiita.com/miler0528/items/1926e93ed97979f8e9fa

 (略)
 
 if not DEBUG:
     import django_heroku
     django_heroku.settings(locals())
+
+    SOCIAL_AUTH_TWITTER_KEY = os.environ['TWITTER_CONSUMER_KEY']
+    SOCIAL_AUTH_TWITTER_SECRET = os.environ['TWITTER_CONSUMER_SECRET']

playground/local_settings.pyの末尾に以下の内容を追記します。

SOCIAL_AUTH_TWITTER_KEY = 'XXXXXXXXXXXXXXXXXXXXXXXXX'
SOCIAL_AUTH_TWITTER_SECRET = 'XXXXXXXXXXXXXXXXXXXXXXXXX'

TWITTER_CONSUMER_KEYTWITTER_CONSUMER_SECRETには先ほど作成したAPI keyおよびAPI secret keyを入れます。

設定を変更したのでマイグレートを行います。

python manage.py migrate

認証システムの実装

hello/urls.pyを以下のように変更してログアウト画面を追加します。

-from django.urls import path
+from django.urls import path, include
+import django.contrib.auth.views
 
 from . import views
 
 urlpatterns = [
     path('', views.index, name='index'),
+    path('logout/',
+        django.contrib.auth.views.LogoutView.as_view(template_name = 'hello/logout.html'),
+        name='logout'),
 ]

playground/urls.pyに一行追記します。

 urlpatterns = [
     path('', include('hello.urls')),
+    path('', include('social_django.urls', namespace='social')),
 ]

hello/views.pyを次のように変更します。こうすることで、ログインしているかどうかでトップページにアクセスしたときの表示内容を変えることができます。

-from django.http import HttpResponse
-
+from django.shortcuts import render
+from social_django.models import UserSocialAuth
 
 def index(request):
-    return HttpResponse("Hello, world!")
+    if request.user.is_authenticated:
+        user = UserSocialAuth.objects.get(user_id=request.user.id)
+        return render(request,'hello/index.html', {'user': user})
+    else:
+        return render(request,'hello/index.html')

テンプレートの追加

テンプレートを入れるフォルダを作成します。

mkdir -p hello/templates/hello

hello/templates/hello/index.htmlを以下の内容で作成します。

<html>
  <head>
    <title>index</title>
  </head>
  <body>
  {% if request.user.is_authenticated %}
    <p>あなたはログインしています (screen_name: {{ user.access_token.screen_name }})</p>
    <a href="/logout"><button type="button">ログアウト</button></a>
  {% else %}
    <p>あなたはログインしていません</p>
    <button type="button" onclick="location.href='{% url 'social:begin' 'twitter' %}'">Twitterでログイン</button>
  {% endif %}
  </body>
</html>

hello/templates/hello/logout.htmlを以下の内容で作成します。

<html>
  <head>
    <title>ログアウト</title>
  </head>
  <body>
    <p>ログアウトしました</p>
    <p><a href="/"><button type="button">トップへ</button></a></p>
  </body>
</html>

以上の作業の結果、以下のようなディレクトリ構成になります。

playground/
├── Pipfile
├── Pipfile.lock
├── Procfile
├── db.sqlite3
├── hello
│   ├── __init__.py
│   ├── admin.py
│   ├── apps.py
│   ├── migrations
│   │   └── __init__.py
│   ├── models.py
│   ├── templates
│   │   └── hello
│   │       ├── index.html
│   │       └── logout.html
│   ├── tests.py
│   ├── urls.py
│   └── views.py
├── manage.py
├── playground
│   ├── __init__.py
│   ├── asgi.py
│   ├── local_settings.py
│   ├── settings.py
│   ├── urls.py
│   └── wsgi.py
└── runtime.txt

ローカルでのテスト

python manage.py runserver

ブラウザでhttp://localhost:8000/にアクセスすると、ログインしていない場合の画面が表示されます。

f:id:kivantium:20200412154331p:plain:w250

Twitterでログイン」ボタンを押すとこんな感じの認証画面に遷移します。 f:id:kivantium:20200412153257p:plain:w600

「連携アプリを認証」ボタンを押すと、ログインしてトップ画面に戻ります。ログインしている場合はscreen_nameが表示されます。

f:id:kivantium:20200412154738p:plain:w400

ログアウトボタンを押すとログアウトされます。

f:id:kivantium:20200412154920p:plain:w200

デプロイ

ローカル環境で正しく動いたらHerokuにデプロイします。

heroku config:set TWITTER_CONSUMER_KEY='XXXXXXXXXXXXXXXXXXXXXXXXX'
heroku config:set TWITTER_CONSUMER_SECRET='XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
git add .
git commit -m 'Add login with Twitter'
git push heroku master
heroku run python manage.py migrate

https://<appname>.herokuapp.com/にアクセスすると、ローカルと同じようなアプリが動くはずです。

タイムラインの読み込み

これでTwitter APIを使う準備が整ったので、試しにタイムラインの読み込みを行ってみます。 Twitter APIを扱うライブラリには以前の記事でも取り上げたtweepyを使うことにします。

tweepyのインストール

pipenv install tweepy

hello/views.pyの変更内容

 from django.shortcuts import render
 from social_django.models import UserSocialAuth
+from django.conf import settings
+
+import tweepy
 
 def index(request):
     if request.user.is_authenticated:
         user = UserSocialAuth.objects.get(user_id=request.user.id)
-        return render(request,'hello/index.html', {'user': user})
+        consumer_key = settings.SOCIAL_AUTH_TWITTER_KEY
+        consumer_secret = settings.SOCIAL_AUTH_TWITTER_SECRET
+        access_token = user.extra_data['access_token']['oauth_token']
+        access_secret = user.extra_data['access_token']['oauth_token_secret']
+        auth = tweepy.OAuthHandler(consumer_key, consumer_secret)
+        auth.set_access_token(access_token, access_secret)
+        api = tweepy.API(auth)
+        timeline = api.home_timeline()
+        return render(request,'hello/index.html', {'user': user, 'timeline': timeline})
     else:
         return render(request,'hello/index.html')

/hello/templates/hello/index.htmlの変更内容

   <body>
   {% if request.user.is_authenticated %}
     <p>あなたはログインしています (screen_name: {{ user.access_token.screen_name }})</p>
+    <ul>
+    {% for tweet in timeline %}
+    <li>@{{ tweet.user.screen_name }}: {{ tweet.text }}</li>
+    {% endfor %}
+    </ul>
     <a href="/logout"><button type="button">ログアウト</button></a>
   {% else %}
     <p>あなたはログインしていません</p>

変更をコミットしてgit push heroku masterでデプロイすると、タイムラインから最新20件のツイートを読み込んで表示することができるようになります。

f:id:kivantium:20200412153845p:plain

次回はいつになるか分かりませんが、Twitterを使った何かしらのアプリを作ってみようと思います。

参考文献

DjangoアプリをHerokuにデプロイする

インターンを始めたら労働のつらさを思い出しました。

というわけで、Webサービスの作り方を勉強していきたいと思います。

WebサービスといえばRuby on Railsというイメージがあったのですが、最近は人気が落ちているという話もよく聞きます。Rubyの経験が全くない自分が今から勉強する必要はないのではないかと思って調べてみたら、Rails・Django・Laravelを3大Webフレームワークと呼んでいる記事を発見しました。このうちのDjangoPythonで書かれているため今までの経験を活かすことができ、機械学習ライブラリなどと組み合わせるのもラクそうだったのでこれを使ってみようと思いました。 f:id:kivantium:20200412000515p:plain:w600

(グラフはStack Overflow Trendsより)

Ruby on Railsチュートリアルの第1章ではHello, world!をHerokuにデプロイしていたため、Djangoでも同じことをやってみようとしたら意外と情報が錯綜していて時間がかかったので自分がやった方法を記録しておきます。DjangoもHerokuも今週始めたばかりです。あまり信用しないでください。

実行環境

環境構築

Deploying Python and Django Apps on Heroku によれば、ルートディレクトリにrequirements.txt, setup.py, Pipfile のいずれかがある場合にPythonによるアプリだと認識されるとのことです。一番モダンそうなpipenvによるパッケージ管理を使うことにしました。

適当なディレクトリで以下のコマンドを実行します。playgroundは好きな名前にしてください。

$ django-admin startproject playground
$ cd playground
$ pipenv --python 3.6
$ pipenv shell
$ pipenv install django-heroku gunicorn

ホームディレクトリで上のコマンドを実行したとすると、この時点での~/playground/ディレクトリの構成は以下のようになっているはずです。

~/playground
├── Pipfile
├── Pipfile.lock
├── manage.py
└── playground
    ├── __init__.py
    ├── asgi.py
    ├── settings.py
    ├── urls.py
    └── wsgi.py

以下、manage.pyがあるディレクトリを「ルートディレクトリ」と呼び、ファイルパスは全てそこからの相対パスで表します。playgroundではない名前をつけた場合は適当に読み替えてください。

ファイルを書き換える

まず、settings.pyと同じフォルダに以下の内容のlocal_settings.pyを作成します。こうすることで開発環境と本番環境でデバッグの有無などを切り替えることができるようになります。SECRET_KEYの行は環境ごとに違うので、playground/settings.pyからコピペしてください。(参考: pipenvでDjangoアプリをherokuにでデプロイしたので手順をメモ

DEBUG = True

SECRET_KEY = '_or(2t#$a9fm9i0_-620e(kv__qz6%%wzwx$j0*+ayo1xd&2d+'

次に、playground/settings.py を書き換えます。編集前後のdiffを示します。赤が削除した行で、青が追加した行です。

+import os
 (略)

 BASE_DIR = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
 # Quick-start development settings - unsuitable for production
 # See https://docs.djangoproject.com/en/3.0/howto/deployment/checklist/
 
-# SECURITY WARNING: keep the secret key used in production secret!
-SECRET_KEY = '_or(2t#$a9fm9i0_-620e(kv__qz6%%wzwx$j0*+ayo1xd&2d+'
 
 # SECURITY WARNING: don't run with debug turned on in production!
-DEBUG = True
+DEBUG = False
 
 ALLOWED_HOSTS = []
 
 (略)
 
 INSTALLED_APPS = []
 
 MIDDLEWARE = [
     'django.middleware.security.SecurityMiddleware',
+     'whitenoise.middleware.WhiteNoiseMiddleware',
     'django.contrib.sessions.middleware.SessionMiddleware',
     'django.middleware.common.CommonMiddleware',
     'django.middleware.csrf.CsrfViewMiddleware',

 (略)
 
USE_TZ = True
 # https://docs.djangoproject.com/en/3.0/howto/static-files/
 
 STATIC_URL = '/static/'
+# 2020/04/26 追記
+# http://whitenoise.evans.io/en/stable/django.html#make-sure-staticfiles-is-configured-correctly
+STATIC_ROOT = os.path.join(BASE_DIR, 'staticfiles')
+
+# Simplified static file serving.
+# https://warehouse.python.org/project/whitenoise/
+
+STATICFILES_STORAGE = 'whitenoise.storage.CompressedManifestStaticFilesStorage'
+
+# Use different settings in local environment and Heroku
+# https://qiita.com/miler0528/items/1926e93ed97979f8e9fa
+try:
+    from .local_settings import *
+except ImportError:
+    pass
+
+if not DEBUG:
+    import django_heroku
+    django_heroku.settings(locals())
  • SECRET_KEYは公開してはいけない情報なのでコードからは削除します。実行するたびに値が違うはずなので適宜読み替えてください。
  • DEBUGは本番環境ではFalseにします。
  • MIDDLEWARESTATICFILES_STORAGE関係の設定はWhiteNoiseを使って静的ファイルを効率的に扱えるようにするためのものです。(参考: Django and Static Assets
  • MIDDLEWAREへの追記はSecurityMiddlewareの直後にする必要があるので注意してください。(参考: enable-whitenoise
  • コードの末尾に追加したのは開発環境と本番環境の設定を切り替えるための処理です。

最後にデプロイに必要なProcfileをルートディレクトリに用意します。playground.wsgiの部分はプロジェクト名に合わせて適宜変更してください。(参考: Configuring Django Apps for Heroku

web: gunicorn playground.wsgi

Pythonのバージョンを指定する場合はルートディレクトリにruntime.txtを用意します。用意しない場合はデフォルトのバージョンが使われます。執筆時点でのデフォルトは3.6.10です(参考: specifying-a-python-version

python-3.6.10

この時点でのディレクトリ構成は以下のようになっているはずです。

playground
├── Pipfile
├── Pipfile.lock
├── Procfile
├── manage.py
├── playground
│   ├── __init__.py
│   ├── asgi.py
│   ├── local_settings.py
│   ├── settings.py
│   ├── urls.py
│   └── wsgi.py
└── runtime.txt

一度ローカル環境で動作を確認してみます。ルートディレクトリで以下を実行します。

python manage.py runserver

ブラウザでhttp://127.0.0.1:8000/にアクセスすると、いい感じのロケットの画像が表示されるはずです。 f:id:kivantium:20200411235809p:plain:w600

Hello, world!

Hello, worldを実行するためのアプリを作成します。

python manage.py startapp hello

hello/views.pyを以下のように編集します。

from django.http import HttpResponse


def index(request):
    return HttpResponse("Hello, world!")

以下の内容のhello/urls.pyを作成します。

from django.urls import path

from . import views

urlpatterns = [
    path('', views.index, name='index'),
]

さらに、playground/urls.pyを以下のように編集します。

from django.urls import include, path

urlpatterns = [
    path('', include('hello.urls')),
]

python manage.py runserverを実行してブラウザでhttp://127.0.0.1:8000/にアクセスすると、Hello, world!が表示されるはずです。

デプロイの準備

Herokuのインストール

The Heroku CLIに従ってインストールします。Ubuntuではsnapでインストールできます。

sudo snap install --classic heroku

初めての場合はユーザー登録をする必要があります。公開鍵を登録しておくと便利です。

heroku login --interactive
heroku keys:add

公開鍵を登録したら、以下を実行します。appnameには好きな名前を設定してください。

heroku login
heroku create <appname>

settings.pyから削除したSECRET_KEY環境変数として設定しておきます。 以下のコマンドを実行してください。(SECRET_KEYの中身は環境に応じて変えてください)

heroku config:set SECRET_KEY='_or(2t#$a9fm9i0_-620e(kv__qz6%%wzwx$j0*+ayo1xd&2d+'

デプロイ

HerokuにはGit経由でデプロイします。ルートディレクトリに以下のような.gitignoreを用意します。

__pycache__
*.pyc
db.sqlite3
staticfiles
local_settings.py

ルートディレクトリで以下を実行します。

git init
git add .
git commit -m 'first commit'
git push heroku master
heroku ps:scale web=1

https://<appname>.herokuapp.com/にアクセスするとHello, world!が表示されたページが確認できると思います。

(pushできなかった場合: Herokuにpush時にdoes not appear to be a git repository出た時の対処 - Qiita

ここまでの作業をまとめたリポジトリはこのコミットなので参考にしてください。 github.com

次回Twitterによるログイン機能を実装します。

参考文献

NWChemをコンパイルして動かす

なにも分からない。

NWChemは計算化学ソフトで、GSoCのOpenChemistryに参加していたので動かしてみたかった。

とりあえず動いた方法

環境はUbuntu 18.04。

ドキュメントのCompiling NWChemにあるNWChem 6.6 on Ubuntu 14.04 (Trusty Tahr)の項目を見ればいいと思っていたが微妙に動かなかった。そもそも現在の最新リリースバージョンは6.8.1なのでドキュメントが更新されてなさそう。以下のように改変したらとりあえず動くようになった。

必要なパッケージのインストール

python-dev gfortran libopenblas-dev libopenmpi-dev openmpi-bin tcsh make 

最新リリースをクローンする

git clone -b hotfix/release-6-8 https://github.com/nwchemgit/nwchem.git nwchem-6.8.1
cd nwchem-6.8.1/src

環境変数の設定

ドキュメントでは

export USE_MPI=y  
export NWCHEM_TARGET=LINUX64  
export USE_PYTHONCONFIG=y  
export PYTHONVERSION=2.7  
export PYTHONHOME=/usr 
export BLASOPT="-lopenblas -lpthread -lrt"  
export BLAS_SIZE=4  
export USE_64TO32=y

とするように指示されているが、BLASOPTを指定するとmake時にLAPACK_LIBを指定しろというエラーが発生する。しかし、その正しい値についてはドキュメントに言及がなく分からなかった。そこで

unset BLASOPT

を実行したところとりあえずエラーが消えたので良しとした。(絶対良くないんだよなぁ)

コンパイル

make nwchem_config NWCHEM_MODULES="all python"
make 64_to_32
make

定石どおり make install しようとしたら動かなかったのでコンパイル後にできたbin/以下をPATHに加えてインストール終了とした。

動作チェック

Getting Started にあるコードが動くことを確認した。(aptでインストールしたバージョンではこれが動かなかったのがそもそもコンパイルしないといけなくなった原因なのでNWChemのメンテナンスが全般的に間に合ってないっぽい)

ついでに閉殻Hartree-Fock法によるHeH+のエネルギー計算 - kivantium活動日記と同じ設定で走らせて結果がどうなるか見てみた。

title "HeH+ STO-3G SCF energy calculation"
charge 1
geometry  
  H 0 0 0
  He 0 0 1.4632
end
basis
  H  library sto-3g
  He library sto-3g
end
task scf energy

結果は

       Final RHF  results 
       ------------------ 

         Total SCF energy =     -2.825194209471
      One-electron energy =     -4.549711559836
      Two-electron energy =      1.001202357200
 Nuclear repulsion energy =      0.723314993165

        Time for solution =      0.0s

となり、前回の結果のTOTAL ENERGY = -2.860662163500と近い値が出ることが確認できた。

閉殻Hartree-Fock法によるHeH+のエネルギー計算

量子化学計算を勉強するために新しい量子化学―電子構造の理論入門〈上〉を読んでいます。
サンプルコードがFortranだったので勉強がてらC言語に移植しました。

新しい量子化学―電子構造の理論入門〈上〉

新しい量子化学―電子構造の理論入門〈上〉

理論

本で200ページくらいかけて説明されている内容をブログに書くのは大変なので省略します。

ちなみに、以下のツイートの「100ページかけて説明してる別の本」はこの本を指します。
難しい内容を易しく簡潔に説明しようとするのは大変で、難しい内容は難しいまま長い時間をかけて説明する必要があるのだなぁという学びがありました。


実装

本のFortran実装を素直にC言語で書き直しました。
本のコードを打ち込むのが面倒だったので、検索して見つけたページにあったコードと見比べながら書きました。

最終的に得られたエネルギーはFortran実装の結果と小数第5位まで一致しました。
この差がどこで生じているのか謎ですが、実装は正しいんじゃないでしょうか……

#include <math.h>
#include <stdio.h>

// プロトタイプ宣言
void HFCALC(int, int, double, double, double, double, double);
void INTGRL(int, int, double, double, double, double, double);
void COLECT(int, int, double, double, double, double, double);
void SCF(int, int, double, double, double, double, double);
double calcS(double, double, double);
double calcT(double, double, double);
double F0(double);
double calcV(double, double, double, double, double);
double TWOE(double, double, double, double, double, double, double);
void FORMG(void);
void DIAG(double [][2], double [][2], double [][2]);
void MULT(double [][2], double [][2], double [][2], int, int);
void MATOUT(double [][2], int, int, int, int, char*);

double S12, T11, T12, T22, V11A, V12A, V22A, V11B, V12B, V22B,
       V1111, V2111, V2121, V2211, V2221, V2222;

// Slater型関数に対する最小2乗近似の短縮係数と軌道指数
double COEF[3][3] = {{1.0, 0.678914, 0.444635}, 
                     {0.0, 0.430129, 0.535328},
                     {0.0, 0.0, 0.154329}};
double EXPON[3][3] = {{0.27095, 0.151623, 0.109818},
                      {0.0, 0.851819, 0.405771},
                      {0.0, 0.0, 2.22766}};
double D1[3], A1[3], D2[3], A2[3];
double PI = 3.1415926535898;

double S[2][2], X[2][2], XT[2][2], H[2][2], F[2][2], G[2][2], C[2][2],
       FPRIME[2][2], CPRIME[2][2], P[2][2], OLDP[2][2], TT[2][2][2][2], E[2][2];

int main(void) {
  int IOP = 2;
  // 基底関数STO-NG
  int N = 3;
  // 基底状態の平衡結合長(a.u.)
  double R = 1.4632;
  // Heに対するSTO-3G軌道指数値
  // He原子に対する最良の軌道指数は1.6875
  // HeH+では電子雲が縮んでいることを反映して1.24倍
  double ZETA1 = 2.0925;
  // Hに対するSTO-3G軌道指数値
  // H原子に対する最良の軌道指数は1.0
  // HeH+では電子雲が縮んでいることを反映して1.24倍
  double ZETA2 = 1.24;
  double ZA = 2.0;
  double ZB = 1.0;
  HFCALC(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}

/*
 * 1s基底関数に対応するSTO-NG最小基底を用いて2電子系2原子分子に対する
 * Hartree-Fock計算を行う関数
 *
 * IOP=0 何も表示しない
 * IOP=1 収束の結果のみ表示
 * IOP=2 全繰り返しの結果を表示
 * N STO-NG を使う (N=1, 2, 3)
 * R 結合長 (a.u.)
 * ZETA1 原子Aの基底関数のSlater軌道指数値
 * ZETA2 原子Bの基底関数のSlater軌道指数値
 * ZA 原子Aの原子番号
 * ZB 原子Bの原子番号
 */
void HFCALC(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  if (IOP != 0) {
    printf("   STO-%dG for atomic numbers %5.2lf and %5.2lf\n\n\n", N, ZA, ZB);
  }
  INTGRL(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
  COLECT(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
  SCF(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}

/* SCF計算に必要な積分を計算する関数*/
void INTGRL(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  double R2 = R * R;
  // 軌道指数のスケーリングと短縮係数の規格化
  for (int i = 0; i < N; ++i) {
    A1[i] = EXPON[i][N-1] * pow(ZETA1, 2.0);
    D1[i] = COEF[i][N-1] * pow(2.0 * A1[i] / PI, 0.75);
    A2[i] = EXPON[i][N-1] * pow(ZETA2, 2.0);
    D2[i] = COEF[i][N-1] * pow(2.0 * A2[i] / PI, 0.75);
  }
  S12 = 0.0;
  T11 = 0.0;
  T12 = 0.0;
  T22 = 0.0;
  V11A = 0.0;
  V12A = 0.0;
  V22A = 0.0;
  V11B = 0.0;
  V12B = 0.0;
  V22B = 0.0;
  V1111 = 0.0;
  V2111 = 0.0;
  V2121 = 0.0;
  V2211 = 0.0;
  V2221 = 0.0;
  V2222 = 0.0;

  // 1電子積分の計算
  // V12A = 原子核Aへの各引力の非対角項
  // RAP2 = 中心Aと中心Pの距離の2乗
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
      double RAP = A2[j] * R / (A1[i] + A2[j]);
      double RAP2 = pow(RAP, 2.0);
      double RBP2 = pow(R - RAP, 2.0);
      S12 += calcS(A1[i], A2[j], R2) * D1[i] * D2[j];
      T11 += calcT(A1[i], A1[j], 0.0) * D1[i] * D1[j];
      T12 += calcT(A1[i], A2[j], R2) * D1[i] * D2[j];
      T22 += calcT(A2[i], A2[j], 0.0) * D2[i] * D2[j];
      V11A += calcV(A1[i], A1[j], 0.0, 0.0, ZA) * D1[i] * D1[j];
      V12A += calcV(A1[i], A2[j], R2, RAP2, ZA) * D1[i] * D2[j];
      V22A += calcV(A2[i], A2[j], 0.0, R2, ZA) * D2[i] * D2[j];
      V11B += calcV(A1[i], A1[j], 0.0, R2, ZB) * D1[i] * D1[j];
      V12B += calcV(A1[i], A2[j], R2, RBP2, ZB) * D1[i] * D2[j];
      V22B += calcV(A2[i], A2[j], 0.0, 0.0, ZB) * D2[i] * D2[j];
    }
  }
  // 2電子積分の計算
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
      for (int k = 0; k < N; ++k) {
        for (int l = 0; l < N; ++l) {
          double RAP = A2[i] * R / (A2[i] + A1[j]);
          //double RBP = R - RAP;
          double RAQ = A2[k] * R / (A2[k] + A1[l]);
          double RBQ = R - RAQ;
          double RPQ = RAP - RAQ;
          double RAP2 = RAP * RAP;
          //double RBP2 = RBP * RBP;
          //double RAQ2 = RAQ * RAQ;
          double RBQ2 = RBQ * RBQ;
          double RPQ2 = RPQ * RPQ;
          V1111 += TWOE(A1[i], A1[j], A1[k], A1[l], 0.0, 0.0, 0.0) 
                   * D1[i] * D1[j] * D1[k] * D1[l];
          V2111 += TWOE(A2[i], A1[j], A1[k], A1[l], R2, 0.0, RAP2)
                   * D2[i] * D1[j] * D1[k] * D1[l];
          V2121 += TWOE(A2[i], A1[j], A2[k], A1[l], R2, R2, RPQ2)
                   * D2[i] * D1[j] * D2[k] * D1[l];
          V2211 += TWOE(A2[i], A2[j], A1[k], A1[l], 0.0, 0.0, R2)
                   * D2[i] * D2[j] * D1[k] * D1[l];
          V2221 += TWOE(A2[i], A2[j], A2[k], A1[l], 0.0, R2, RBQ2)
                   * D2[i] * D2[j] * D2[k] * D1[l];
          V2222 += TWOE(A2[i], A2[j], A2[k], A2[l], 0.0, 0.0, 0.0) * D2[i] *
                   D2[j] * D2[k] * D2[l];
        }
      }
    }
  }
  if (IOP != 0) {
    printf("   R          ZETA1      ZETA2      S12       T11\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", R, ZETA1, ZETA2, S12, T11);
    printf("   T12        T22        V11A       V12A        V22A\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", T12, T22, V11A, V12A, V22A);
    printf("   V11B       V12B       V22B       V1111        V2111\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", V11B, V12B, V22B, V1111,
           V2111);
    printf("   V2121      V2211       V2221     V2222\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf\n", V2121, V2211, V2221, V2222);
  }
}

// 重なり積分
double calcS(double A, double B, double RAB2) {
  return pow(PI / (A + B), 1.5) * exp(-A * B * RAB2 / (A + B));
}

// 運動エネルギー積分
double calcT(double A, double B, double RAB2) {
  return A * B / (A + B) * (3.0 - 2.0 * A * B * RAB2 / (A + B)) *
         pow(PI / (A + B), 1.5) * exp(-A*B*RAB2/(A+B));
}

// F0関数
double F0(double ARG) {
  double ret;
  if (ARG < 1.0e-6) {
    ret = 1.0 - ARG / 3.0;
  } else {
    ret = sqrt(PI / ARG) * erf(sqrt(ARG)) / 2.0;
  }
  return ret;
}

// 核引力積分
double calcV(double A, double B, double RAB2, double RCP2, double ZC) {
  double v = 2.0 * PI / (A + B) * F0((A + B) * RCP2) * exp(-A * B * RAB2 / (A + B));
  return -v * ZC;
}

// 2電子積分
double TWOE(double A, double B, double C, double D, double RAB2, double RCD2,
            double RPQ2) {
  return 2.0 * pow(PI, 2.5) / ((A + B) * (C + D) * sqrt(A + B + C + D))
         * F0((A + B) * (C + D) * RPQ2 / (A + B + C + D))
         * exp(-A * B * RAB2 / (A + B) - C * D * RCD2 / (C + D));
}

// 積分の結果から行列を求める
void COLECT(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  // 核ハミルトニアン cf. 式(3.153)
  H[0][0] = T11 + V11A + V11B;
  H[0][1] = T12 + V12A + V12B;
  H[1][0] = H[0][1];
  H[1][1] = T22 + V22A + V22B;

  // 重なり行列 cf. 式(3.136)
  S[0][0] = 1.0;
  S[0][1] = S12;
  S[1][0] = S[0][1];
  S[1][1] = 1.0;

  // 正準直交化 cf. 式(3.169), p.191
  X[0][0] = 1.0 / sqrt(2.0 * (1.0 + S12));
  X[1][0] = X[0][0];
  X[0][1] = 1.0 / sqrt(2.0 * (1.0 - S12));
  X[1][1] = -X[0][1];

  // 変換行列の転置
  XT[0][0] = X[0][0];
  XT[0][1] = X[1][0];
  XT[1][0] = X[0][1];
  XT[1][1] = X[1][1];

  // 2電子積分の行列
  TT[0][0][0][0] = V1111;
  TT[1][0][0][0] = V2111;
  TT[0][1][0][0] = V2111;
  TT[0][0][1][0] = V2111;
  TT[0][0][0][1] = V2111;
  TT[1][0][1][0] = V2121;
  TT[0][1][1][0] = V2121;
  TT[1][0][0][1] = V2121;
  TT[0][1][0][1] = V2121;
  TT[1][1][0][0] = V2211;
  TT[0][0][1][1] = V2211;
  TT[1][1][1][0] = V2221;
  TT[1][1][0][1] = V2221;
  TT[1][0][1][1] = V2221;
  TT[0][1][1][1] = V2221;
  TT[1][1][1][1] = V2222;

  if (IOP != 0) {
    MATOUT(S, 2, 2, 2, 2, "S   ");
    MATOUT(X, 2, 2, 2, 2, "X   ");
    MATOUT(H, 2, 2, 2, 2, "H   ");
    printf("\n");
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        for (int k = 0; k < 2; ++k) {
          for (int l = 0; l < 2; ++l) {
            printf("   (%2d%2d%2d%2d )%10.6lf\n", i+1, j+1, k+1, l+1, TT[i][j][k][l]);
          }
        }
      }
    }
  }
}

// SCFを実行する
void SCF(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
         double ZB) {
  // 密度行列の収束規準
  double CRIT = 1.0e-4;
  // 繰り返しの最大回数
  int MAXIT = 25;
  // 核ハミルトニアンを最初の推定に用いる(つまりPを零行列で初期化)
  for (int i = 0; i < 2; ++i) {
    for (int j = 0; j < 2; ++j) {
      P[i][j] = 0.0;
    }
  }
  if (IOP == 2) {
    MATOUT(P, 2, 2, 2, 2, "P    ");
  }
  int ITER = 0;
  double EN;
  for (;;) {
    ITER++;
    if (IOP == 2) {
      printf("\n    START OF ITERATION NUMBER = %2d\n", ITER);
    }
    // Fock行列の2電子部分を計算
    FORMG();

    if (IOP == 2) {
      MATOUT(G, 2, 2, 2, 2, "G    ");
    }

    // 核ハミルトニアンを足してFock行列を得る
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        F[i][j] = H[i][j] + G[i][j];
      }
    }

    // 電子エネルギーを計算する
    EN = 0.0;
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        EN += 0.5 * P[i][j] * (H[i][j] + F[i][j]);
      }
    }

    if (IOP == 2) {
      MATOUT(F, 2, 2, 2, 2, "F   ");
      printf("\n\n\n    ELECTRONIC ENERGY = %20.12lf\n", EN);
    }

    // Gを用いてFock行列を変換
    MULT(F, X, G, 2, 2);
    MULT(XT, G, FPRIME, 2, 2);
    // Fock行列の対角化
    DIAG(FPRIME, CPRIME, E);
    // 固有ベクトルから行列Cを得る
    MULT(X, CPRIME, C, 2, 2);
    // 新しい密度行列を計算する
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        OLDP[i][j] = P[i][j];
        P[i][j] = 0.0;
        for (int k = 0; k < 1; ++k) { // ここはk<1で本当にいいのか不明
          P[i][j] += 2.0 * C[i][k] * C[j][k];
        }
      }
    }
    if (IOP == 2) {
      MATOUT(FPRIME, 2, 2, 2, 2, "F'   ");
      MATOUT(CPRIME, 2, 2, 2, 2, "C'   ");
      MATOUT(E, 2, 2, 2, 2, "E   ");
      MATOUT(C, 2, 2, 2, 2, "C   ");
      MATOUT(P, 2, 2, 2, 2, "P   ");
    }
    // deltaを計算する
    double DELTA = 0.0;
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        DELTA += pow(P[i][j] - OLDP[i][j], 2.0);
      }
    }
    DELTA = sqrt(DELTA / 4.0);
    if (IOP != 0) {
      printf("\n    DELTA(CONVERGENCE OF DENSITY MATRIX) = %10.6lf\n", DELTA);
    }

    // 収束判定
    if (DELTA < CRIT) break;
    if (ITER < MAXIT) continue;
    // 収束しなかった場合
    printf("    NO CONVERGENCE IN SCF");
    return;
  }

  double ENT = EN + ZA * ZB / R;
  if (IOP != 0) {
    printf("\n\n    CALCULATION CONVERGED"
           "\n\n    ELECTORONIC ENERGY = %20.12lf"
           "\n\n    TOTAL ENERGY =       %20.12f", EN, ENT);
  }
  if (IOP == 1) {
    MATOUT(G, 2, 2, 2, 2, "G   ");
    MATOUT(F, 2, 2, 2, 2, "F   ");
    MATOUT(E, 2, 2, 2, 2, "E   ");
    MATOUT(C, 2, 2, 2, 2, "C   ");
    MATOUT(P, 2, 2, 2, 2, "P   ");
  }
  MULT(P, S, OLDP, 2, 2);
  if (IOP != 0) {
    MATOUT(OLDP, 2, 2, 2, 2, "PS  ");
  }
}

void FORMG(void) {
  for (int i = 0; i < 2; ++i) {
    for (int j = 0; j < 2; ++j) {
      G[i][j] = 0.0;
      for (int k = 0; k < 2; ++k) {
        for (int l = 0; l < 2; ++l) {
          G[i][j] += P[k][l] * (TT[i][j][k][l] - 0.5 * TT[i][l][k][j]);
        }
      }
    }
  }
}

// Fを対角化し、Cに固有ベクトル、Eに固有値を入れる
void DIAG(double F[][2], double C[][2], double E[][2]) {
  double theta;
  if (fabs(F[0][0] - F[1][1]) > 1.0e-20) {
    theta = 0.5 * atan(2.0 * F[0][1] / (F[0][0] - F[1][1]));
  } else {
    theta = PI / 4.0;
  }
  C[0][0] = cos(theta);
  C[1][0] = sin(theta);
  C[0][1] = sin(theta);
  C[1][1] = -cos(theta);
  E[0][0] = F[0][0] * pow(cos(theta), 2.0) + F[1][1] * pow(sin(theta), 2.0)
            + F[0][1] * sin(2.0 * theta);
  E[1][1] = F[1][1] * pow(cos(theta), 2.0) + F[0][0] * pow(sin(theta), 2.0)
            - F[0][1] * sin(2.0 * theta);
  E[1][0] = 0.0;
  E[0][1] = 0.0;
  // 固有値と固有ベクトルを整列する
  if (E[1][1] > E[0][0]) return;
  double temp = E[1][1];
  E[1][1] = E[0][0];
  E[0][0] = temp;
  temp = C[0][1];
  C[0][1] = C[0][0];
  C[0][0] = temp;
  temp = C[1][1];
  C[1][1] = C[1][0];
  C[1][0] = temp;
}

// AとBの積をCに入れる
void MULT(double A[][2], double B[][2], double C[][2], int IM, int M) {
  for (int i = 0; i < M; ++i) {
    for (int j = 0; j < M; ++j) {
      C[i][j] = 0.0;
      for (int k = 0; k < M; ++k) {
        C[i][j] += A[i][k] * B[k][j];
      }
    }
  }
}

void MATOUT(double A[][2], int IM, int IN, int M, int N, char *LABEL) {
  int IHIGH = 0;
  for (;;) {
    int LOW = IHIGH + 1;
    IHIGH = IHIGH + 5;
    IHIGH = (IHIGH < N) ? IHIGH : N;
    printf("\n\n\n    THE %s ARRAY\n", LABEL);
    printf("               ");
    for (int i = LOW; i <= IHIGH; ++i) {
      printf("          %3d      ", i);
    }
    printf("\n");
    for (int i = 1; i <= M; ++i) {
      printf("%10d     ", i);
      for (int j = LOW; j <= IHIGH; ++j) {
        printf(" %18.10lf", A[i - 1][j - 1]);
      }
      printf("\n");
    }
    if (N <= IHIGH) break;
  }
}

最終的な出力は

    CALCULATION CONVERGED

    ELECTORONIC ENERGY =      -4.227529304014

    TOTAL ENERGY =            -2.860662163500


    THE PS   ARRAY
                           1                  2      
         1            1.5296370381       1.1199277981
         2            0.6424383868       0.4703629619

でした。参照したFortranコードの出力は

    CALCULATION CONVERGED

    ELECTRONIC ENERGY =  -0.422752913203D+01

    TOTAL ENERGY =       -0.286066199152D+01



    THE PS   ARRAY
                           1                  2
         1        0.1529637121D+01   0.1119927796D+01
         2        0.6424383099D+00   0.4703628793D+00

でした。

FortranからCに移植するときのtips

詰まったことを挙げておきます。

配列の添字

Fortranでは添字が1から始まるので配列にアクセスするときに添字をうまく調整する必要があります。

多次元配列の並び順

このページに詳しく書いてありますが、C言語の多次元配列はメモリ上でa[0][0], a[0][1], a[1][0], a[1][1], a[2][0], a[2][1]のような順番で並びますが、Fortranではa(1,1), a(2,1), a(3,1), a(1,2), a(2,2), a(3,2)のような順番に並びます。多次元配列の初期化を書き直すときに気をつける必要があります。

整合配列

Fortranではサブルーチンに配列を渡した配列のサイズをサブルーチン側で定義することができます。(整合配列という)

      dimension a(6)
      a(1) = 1
      a(2) = 2
      a(3) = 3
      a(4) = 4
      a(5) = 5
      a(6) = 6
      call sub(a,3,2)
      end

      subroutine sub(a,n,m)
      dimension a(n,m) <---- a(n,*) でも同じ結果を得る
      do i=1,n
         write(6,*) (a(i,j),j=1,m)
      end do
      return
      end

コードはFortran プログラミングの基礎知識から引用しました。

C99で導入された可変長配列の機能を使うと同じことができるそうです。cf. 3.3 整合配列、可変長配列
今回は配列のサイズが決め打ちできたのでvoid MULT(double A[][2], double B[][2], double C[][2], int IM, int M) { のように関数を定義しました。

C言語のabs

特にFortranとは関係がないのですが、C言語のabsはint型に対して使う関数で、double型に対してはfabsという別の関数が用意されています。(これはC言語で関数オーバーロードができないことが原因です。C11で導入されたGenericという機能を使うとそれっぽいことが実現できるらしいです。cf. Generic - C言語入門

広告コーナー