nogawaparkのブログ

書く練習。

二次元氷の場合の数の計算 1

初めに

二次元氷というものを前から知っていて興味を持っていた。以前よりは能力が上がった部分もあり、取り組んでみようと思う。二次元氷の数理について調査、計算していく。

二次元氷について

二次元氷というものがある。普通氷は三次元のものであり、二次元氷は仮想的な結晶である。性質は

  • 酸素Oが二次元格子点上に配置される

  • 水素Hが線分OO間に配置される

  • 各HはどちらかのOに属する

  • 各Oは二つのHを持つ

というルールがある。下の絵が分かりやすい。(この論文中の図を拝借https://www.researchgate.net/publication/235466723_Dichotomous_collective_proton_dynamics_in_ice)

f:id:nogawapark:20190507090641p:plain
2次元氷の模式図

この図の白丸がOで小さな黒丸がHである。(図中のポテンシャル図のようなものは無視してください)

上の図ではHOHの向きがすべて同じであるが、そうではない組み合わせも構成することができる。

また、今はOO間OH間距離は一定で、∠HOH=90°, 180° という抽象的な世界を考えている。

二次元氷の場合の数

このような二次元氷の場合の数は何通りであるか、それについてはwikipedia (https://en.wikipedia.org/wiki/Ice-type_model)に書かれているように、1967年に Lieb によって厳密に解かれている。解かれているらしいが分子数無限大での値のみがwikiにのっている。(論文は読んでいない。一般の m×n の氷の場合の数って厳密解ってあった?)。場合の数を  W、一辺の長さを n とすると、分子数  n×n の二次元氷の時  \lim _ {n \to \infty}W^ {1/n^ 2} = (4/3)^ {3/2}=1.53960... となる。

境界条件

境界条件は周期境界条件と境界制約なしの二つがあったと記憶している。

プログラムで計算

今回はプログラムで何通りになるか計算を行った。境界条件の制約はなしで行った。 結果は以下のようになった。場合の数が多すぎて、3×3までしか計算できなかった。

n が 3 程度では  W^ {1/n^ 2} の 1.53960... への収束は分からなかった。

n×n の氷の場合の数  W の n 依存性

 n  W  W^ {1/n^ 2}
1 6 6
2 343 4.304
3 90732 3.555

コード

上の計算は以下の python コードで行った。

コーディングの雑感

  • オブジェクトのアトリビュートを別の新規のオブジェクトのアトリビュートとしてはいけない。書き換えが同時に起こるから。
  • numpy array を要素に持つリスト中に指定された array があるかを、単純に in で調べたらエラーになった。理由はわかっていない。
  • 一応 1×2 の氷の場合の数があっていたから OK とした。
  • 並列化で高速化はできそう。
  • データ構造の持ち方に無駄があるかもしれない。
  • チェック機能が重そうである。
'''
2次元氷が与えられた大きさに対して何通りあるか求める
'''
import numpy as np
import itertools

NORTH = np.array([0, -1])
EAST = np.array([1, 0])
SOUTH = np.array([0, 1])
WEST = np.array([-1, 0])
NEWS = [NORTH, EAST, SOUTH, WEST]

class H2O:
    def __init__(self, place, gotoes, h2o_num):
        self.place = place
        self.gotoes = gotoes
        self.h2o_num = h2o_num
    def __repr__(self):
        return f'{self.place}, {self.gotoes}'

class Ice:
    '''h2oを管理するクラス
    '''
    def __init__(self, ice=None):
        if ice:
            self.h2os = [h2o for h2o in ice.h2os]
            self.places = [h2o.place for h2o in ice.h2os]
        else:
            self.h2os = []
            self.places = []
    
    def __repr__(self):
        h2o_nums = [h2o.h2o_num for h2o in self.h2os]
        return f'{h2o_nums}'

    def add(self, h2o):
        self.h2os.append(h2o)
        self.places.append(h2o.place)
    
    def get_h2o(self, place):
        return [h2o for h2o in self.h2os if np.array_equal(h2o.place, place)][0]

    def view_ice(self):
        pass
    
def first_ice():
    h2os = []
    ices = []
    h2o_num = 0
    place = np.array([0, 0])
    for arrow1, arrow2 in list(itertools.combinations(NEWS, 2)):
        h2o = H2O(place, [place + arrow1, place + arrow2], h2o_num)
        h2os.append(h2o)
        ice = Ice()
        ice.add(h2o)
        ices.append(ice)
        h2o_num += 1
    return ices, h2o_num

def arr_in_list(arr, lst):
    for elm in lst:
        if np.array_equal(arr, elm):
            return True
    return False

def check_h2o(ice, h2o):
    '''check if can insert h2o in ice
    '''
    for goto_place in h2o.gotoes:
        if not arr_in_list(goto_place, ice.places):
            continue
        h2o_goto = ice.get_h2o(goto_place)
        if arr_in_list(h2o.place, h2o_goto.gotoes):
            return False
    return True

def get_next_h2os(place, h2o_num):
    h2os = []
    for arrow1, arrow2 in list(itertools.combinations(NEWS, 2)):
        h2o = H2O(place, [place + arrow1, place + arrow2], h2o_num)    
        h2os.append(h2o)
        h2o_num += 1
    return h2os, h2o_num

def num_ice(m, n):
    '''(m,n)の大きさの2次元氷が与えられたときの場合の数を求める。
    境界は自由とする
    '''
    ices, h2o_num = first_ice()
    for i in range(1, m*n):
        place = np.array([i%n, i//n])
        next_h2os, h2o_num = get_next_h2os(place, h2o_num)
        ices_next = []
              
        for ice in ices:
            for h2o in next_h2os:
                ice_next = Ice(ice)
                if check_h2o(ice_next, h2o):
                    ice_next.add(h2o)
                    ices_next.append(ice_next)

        ices = ices_next
    '''
    for ice in ices:
        print(f'------------------')
        for h2o in ice.h2os:
            print(h2o)
    '''
    #print(len(ices))
    return len(ices)

if __name__=='__main__':
    for n in range(1, 4):
        num = num_ice(n, n)
        print(n, num, num**(1/n**2))