LAMMPSの初期座標ファイル（dataファイル）の作成【バルク水のMD計算】
前回LAMMPSのインプットファイルを作成したので、今回は初期座標ファイル（dataファイルの作成）です。
64の水分子からなるバルク水モデルのMD計算の初期座標を作成します。
192 atoms
128 bonds
64 angles
0 dihedrals
0 impropers
2 atom types
2 bond types
1 angle types
0 dihedral types
0 improper types
0.000000000000000 12.417000000000000 xlo xhi
0.000000000000000 12.417000000000000 ylo yhi
0.000000000000000 12.417000000000000 zlo zhi
Masses
1 15.99940 # OW_spc
2 1.00800 # HW_spc
Pair Coeffs
1 0.155504063097514 3.165570000000000 # OW_spc
2 0.000000000000000 0.000000000000000 # HW_spc
Bond Coeffs
1 0.000000000000000 1.000000000000000
2 0.000000000000000 1.000000000000000
Angle Coeffs
1 0.000000000000000 109.470000000000000
Atoms
1 1 1 -0.848 2.300 6.280 1.130
2 1 2 0.424 1.370 6.260 1.500
3 1 2 0.424 2.310 5.890 0.210
4 2 1 -0.848 2.250 2.750 9.960
5 2 2 0.424 2.600 2.580 10.880
6 2 2 0.424 1.370 2.300 9.840
7 3 1 -0.848 0.190 3.680 6.470
8 3 2 0.424 -0.630 4.110 6.860
9 3 2 0.424 -0.090 2.950 5.840
10 4 1 -0.848 5.690 0.330 11.650
11 4 2 0.424 4.760 0.260 11.280
12 4 2 0.424 5.800 1.220 12.090
13 5 1 -0.848 11.350 7.030 7.170
14 5 2 0.424 11.920 7.810 6.920
15 5 2 0.424 10.750 7.290 7.930
16 6 1 -0.848 7.680 11.440 10.230
17 6 2 0.424 6.900 11.610 10.830
18 6 2 0.424 8.020 12.310 9.870
19 7 1 -0.848 6.850 10.120 6.650
20 7 2 0.424 7.540 9.960 7.350
21 7 2 0.424 6.120 10.690 7.030
22 8 1 -0.848 2.400 10.910 8.860
23 8 2 0.424 2.540 10.070 9.380
24 8 2 0.424 1.850 11.550 9.410
25 9 1 -0.848 1.220 6.430 5.630
26 9 2 0.424 0.770 5.550 5.800
27 9 2 0.424 1.210 6.970 6.470
28 10 1 -0.848 6.130 1.230 7.260
29 10 2 0.424 5.640 0.360 7.350
30 10 2 0.424 5.900 1.660 6.390
31 11 1 -0.848 8.090 0.040 5.020
32 11 2 0.424 8.490 0.950 4.930
33 11 2 0.424 7.090 0.120 5.080
34 12 1 -0.848 1.970 9.760 0.220
35 12 2 0.424 2.860 9.310 0.080
36 12 2 0.424 1.240 9.110 0.030
37 13 1 -0.848 9.910 4.100 0.000
38 13 2 0.424 9.140 4.440 0.540
39 13 2 0.424 9.570 3.590 -0.790
40 14 1 -0.848 10.410 7.010 4.290
41 14 2 0.424 10.670 6.970 5.250
42 14 2 0.424 9.560 6.500 4.150
43 15 1 -0.848 5.500 1.960 8.850
44 15 2 0.424 5.450 1.910 9.850
45 15 2 0.424 5.520 2.920 8.560
46 16 1 -0.848 3.210 9.430 2.420
47 16 2 0.424 4.030 9.820 2.000
48 16 2 0.424 2.940 8.610 1.930
49 17 1 -0.848 3.820 7.000 4.800
50 17 2 0.424 4.270 6.100 4.770
51 17 2 0.424 2.880 6.890 5.130
52 18 1 -0.848 8.030 1.600 9.240
53 18 2 0.424 8.930 1.740 8.820
54 18 2 0.424 7.320 1.620 8.530
55 19 1 -0.848 9.220 5.030 8.990
56 19 2 0.424 8.970 4.940 8.030
57 19 2 0.424 9.700 4.210 9.300
58 20 1 -0.848 5.390 0.640 5.120
59 20 2 0.424 4.580 0.650 5.700
60 20 2 0.424 5.420 1.470 4.570
61 21 1 -0.848 2.970 0.350 1.710
62 21 2 0.424 3.460 1.190 1.500
63 21 2 0.424 3.590 -0.300 2.160
64 22 1 -0.848 9.350 2.360 4.800
65 22 2 0.424 8.870 2.770 4.020
66 22 2 0.424 10.340 2.340 4.610
67 23 1 -0.848 10.760 6.830 2.220
68 23 2 0.424 9.960 6.220 2.250
69 23 2 0.424 11.570 6.300 1.980
70 24 1 -0.848 4.590 11.520 7.410
71 24 2 0.424 3.880 11.250 8.060
72 24 2 0.424 4.330 11.240 6.480
73 25 1 -0.848 8.660 4.540 6.420
74 25 2 0.424 8.340 5.260 5.800
75 25 2 0.424 8.900 3.730 5.890
76 26 1 -0.848 10.170 0.390 7.530
77 26 2 0.424 9.450 0.440 6.840
78 26 2 0.424 9.930 -0.300 8.220
79 27 1 -0.848 5.270 2.560 3.280
80 27 2 0.424 5.540 1.970 2.530
81 27 2 0.424 5.270 3.510 2.970
82 28 1 -0.848 10.640 1.050 0.980
83 28 2 0.424 9.840 0.820 1.530
84 28 2 0.424 11.470 0.790 1.470
85 29 1 -0.848 7.050 10.500 3.680
86 29 2 0.424 6.910 10.570 4.670
87 29 2 0.424 7.890 9.990 3.500
88 30 1 -0.848 4.100 8.130 0.090
89 30 2 0.424 4.960 8.250 0.590
90 30 2 0.424 3.680 7.260 0.360
91 31 1 -0.848 0.320 3.860 0.200
92 31 2 0.424 0.530 4.600 0.840
93 31 2 0.424 -0.570 4.030 -0.230
94 32 1 -0.848 3.670 11.000 5.010
95 32 2 0.424 3.600 11.830 4.450
96 32 2 0.424 3.710 10.200 4.410
97 33 1 -0.848 5.660 5.370 8.650
98 33 2 0.424 5.780 6.030 7.910
99 33 2 0.424 6.120 5.710 9.480
100 34 1 -0.848 0.300 2.030 11.420
101 34 2 0.424 0.770 2.160 12.290
102 34 2 0.424 -0.360 1.290 11.510
103 35 1 -0.848 5.730 8.700 10.290
104 35 2 0.424 6.170 9.590 10.200
105 35 2 0.424 5.100 8.700 11.060
106 36 1 -0.848 1.180 8.620 10.450
107 36 2 0.424 0.430 8.620 9.790
108 36 2 0.424 1.550 7.700 10.540
109 37 1 -0.848 3.070 2.130 12.310
110 37 2 0.424 2.840 2.500 13.210
111 37 2 0.424 2.770 1.180 12.250
112 38 1 -0.848 7.320 6.340 10.640
113 38 2 0.424 7.910 6.080 9.880
114 38 2 0.424 7.040 7.300 10.530
115 39 1 -0.848 3.070 0.630 6.180
116 39 2 0.424 2.960 1.570 6.510
117 39 2 0.424 3.020 0.000 6.950
118 40 1 -0.848 10.230 7.660 9.660
119 40 2 0.424 10.380 7.870 10.620
120 40 2 0.424 9.930 6.710 9.570
121 41 1 -0.848 3.190 8.100 9.490
122 41 2 0.424 4.120 8.460 9.540
123 41 2 0.424 3.130 7.250 10.010
124 42 1 -0.848 3.390 5.090 10.060
125 42 2 0.424 2.870 4.260 9.890
126 42 2 0.424 4.160 5.140 9.420
127 43 1 -0.848 10.020 7.960 12.380
128 43 2 0.424 10.430 7.640 13.240
129 43 2 0.424 9.060 7.690 12.350
130 44 1 -0.848 0.400 5.440 11.140
131 44 2 0.424 1.250 5.110 10.730
132 44 2 0.424 0.530 5.590 12.120
133 45 1 -0.848 8.150 5.720 3.250
134 45 2 0.424 8.220 4.830 2.790
135 45 2 0.424 7.210 6.060 3.170
136 46 1 -0.848 6.710 4.640 0.270
137 46 2 0.424 6.370 3.750 -0.030
138 46 2 0.424 6.970 5.180 -0.530
139 47 1 -0.848 4.730 5.000 1.910
140 47 2 0.424 5.340 5.800 1.950
141 47 2 0.424 3.780 5.310 1.980
142 48 1 -0.848 0.600 8.550 3.090
143 48 2 0.424 -0.260 8.240 3.510
144 48 2 0.424 0.560 8.410 2.100
145 49 1 -0.848 0.830 0.160 10.220
146 49 2 0.424 0.780 0.150 11.220
147 49 2 0.424 0.000 -0.250 9.840
148 50 1 -0.848 0.380 9.270 6.720
149 50 2 0.424 0.980 8.460 6.740
150 50 2 0.424 0.780 9.960 6.120
151 51 1 -0.848 2.630 3.260 7.200
152 51 2 0.424 1.840 3.770 6.860
153 51 2 0.424 2.540 3.110 8.180
154 52 1 -0.848 8.220 10.020 1.300
155 52 2 0.424 8.620 10.010 0.380
156 52 2 0.424 8.320 10.940 1.700
157 53 1 -0.848 9.160 9.100 2.910
158 53 2 0.424 9.790 9.480 2.230
159 53 2 0.424 9.560 8.270 3.300
160 54 1 -0.848 10.390 10.980 6.960
161 54 2 0.424 9.690 10.510 7.500
162 54 2 0.424 10.980 10.300 6.530
163 55 1 -0.848 0.720 1.660 3.180
164 55 2 0.424 0.550 2.490 2.640
165 55 2 0.424 1.620 1.290 2.960
166 56 1 -0.848 6.130 8.420 1.890
167 56 2 0.424 6.690 9.230 1.720
168 56 2 0.424 6.720 7.620 1.920
169 57 1 -0.848 5.940 7.450 6.520
170 57 2 0.424 6.440 8.300 6.330
171 57 2 0.424 5.060 7.470 6.040
172 58 1 -0.848 8.590 1.320 0.160
173 58 2 0.424 8.130 1.470 1.040
174 58 2 0.424 9.030 2.170 -0.140
175 59 1 -0.848 8.590 9.560 8.610
176 59 2 0.424 9.130 8.870 9.090
177 59 2 0.424 8.270 10.250 9.270
178 60 1 -0.848 10.830 9.840 0.870
179 60 2 0.424 10.600 10.370 0.050
180 60 2 0.424 11.640 9.280 0.680
181 61 1 -0.848 0.790 12.400 6.530
182 61 2 0.424 0.780 11.930 7.410
183 61 2 0.424 1.610 12.120 6.020
184 62 1 -0.848 4.280 4.240 5.200
185 62 2 0.424 4.580 3.520 4.580
186 62 2 0.424 3.890 3.840 6.030
187 63 1 -0.848 10.010 0.340 11.540
188 63 2 0.424 9.380 -0.380 11.230
189 63 2 0.424 10.940 -0.020 11.540
190 64 1 -0.848 7.700 0.880 3.010
191 64 2 0.424 7.240 0.010 3.180
192 64 2 0.424 8.610 0.850 3.420
Bonds
1 1 1 2
2 2 1 3
3 1 4 5
4 2 4 6
5 1 7 8
6 2 7 9
7 1 10 11
8 2 10 12
9 1 13 14
10 2 13 15
11 1 16 17
12 2 16 18
13 1 19 20
14 2 19 21
15 1 22 23
16 2 22 24
17 1 25 26
18 2 25 27
19 1 28 29
20 2 28 30
21 1 31 32
22 2 31 33
23 1 34 35
24 2 34 36
25 1 37 38
26 2 37 39
27 1 40 41
28 2 40 42
29 1 43 44
30 2 43 45
31 1 46 47
32 2 46 48
33 1 49 50
34 2 49 51
35 1 52 53
36 2 52 54
37 1 55 56
38 2 55 57
39 1 58 59
40 2 58 60
41 1 61 62
42 2 61 63
43 1 64 65
44 2 64 66
45 1 67 68
46 2 67 69
47 1 70 71
48 2 70 72
49 1 73 74
50 2 73 75
51 1 76 77
52 2 76 78
53 1 79 80
54 2 79 81
55 1 82 83
56 2 82 84
57 1 85 86
58 2 85 87
59 1 88 89
60 2 88 90
61 1 91 92
62 2 91 93
63 1 94 95
64 2 94 96
65 1 97 98
66 2 97 99
67 1 100 101
68 2 100 102
69 1 103 104
70 2 103 105
71 1 106 107
72 2 106 108
73 1 109 110
74 2 109 111
75 1 112 113
76 2 112 114
77 1 115 116
78 2 115 117
79 1 118 119
80 2 118 120
81 1 121 122
82 2 121 123
83 1 124 125
84 2 124 126
85 1 127 128
86 2 127 129
87 1 130 131
88 2 130 132
89 1 133 134
90 2 133 135
91 1 136 137
92 2 136 138
93 1 139 140
94 2 139 141
95 1 142 143
96 2 142 144
97 1 145 146
98 2 145 147
99 1 148 149
100 2 148 150
101 1 151 152
102 2 151 153
103 1 154 155
104 2 154 156
105 1 157 158
106 2 157 159
107 1 160 161
108 2 160 162
109 1 163 164
110 2 163 165
111 1 166 167
112 2 166 168
113 1 169 170
114 2 169 171
115 1 172 173
116 2 172 174
117 1 175 176
118 2 175 177
119 1 178 179
120 2 178 180
121 1 181 182
122 2 181 183
123 1 184 185
124 2 184 186
125 1 187 188
126 2 187 189
127 1 190 191
128 2 190 192
Angles
1 1 2 1 3
2 1 5 4 6
3 1 8 7 9
4 1 11 10 12
5 1 14 13 15
6 1 17 16 18
7 1 20 19 21
8 1 23 22 24
9 1 26 25 27
10 1 29 28 30
11 1 32 31 33
12 1 35 34 36
13 1 38 37 39
14 1 41 40 42
15 1 44 43 45
16 1 47 46 48
17 1 50 49 51
18 1 53 52 54
19 1 56 55 57
20 1 59 58 60
21 1 62 61 63
22 1 65 64 66
23 1 68 67 69
24 1 71 70 72
25 1 74 73 75
26 1 77 76 78
27 1 80 79 81
28 1 83 82 84
29 1 86 85 87
30 1 89 88 90
31 1 92 91 93
32 1 95 94 96
33 1 98 97 99
34 1 101 100 102
35 1 104 103 105
36 1 107 106 108
37 1 110 109 111
38 1 113 112 114
39 1 116 115 117
40 1 119 118 120
41 1 122 121 123
42 1 125 124 126
43 1 128 127 129
44 1 131 130 132
45 1 134 133 135
46 1 137 136 138
47 1 140 139 141
48 1 143 142 144
49 1 146 145 147
50 1 149 148 150
51 1 152 151 153
52 1 155 154 156
53 1 158 157 159
54 1 161 160 162
55 1 164 163 165
56 1 167 166 168
57 1 170 169 171
58 1 173 172 174
59 1 176 175 177
60 1 179 178 180
61 1 182 181 183
62 1 185 184 186
63 1 188 187 189
64 1 191 190 192
Header
最初の行は1行空行にする。
atomsは原子数、bondsは結合の数、anglesは結合角の数、dihedralsは二面角の数、impropersはimproperな二面角の数。
atom typesは原子種の数、bond typesは結合種の数、angle typesは結合角の種類の数、dihedral typesは二面角の種類の数、improper typesはimproperな二面角の種類の数を示している。
0.000000000000000 12.417000000000000 xlo xhi
0.000000000000000 12.417000000000000 ylo yhi
0.000000000000000 12.417000000000000 zlo zhi
はセル境界の位置。0≦x,y,z≦12.417が計算領域となる。
Masses
原子質量
Pair Coeffs
pair_styleで指定した関数系における係数を設定する。今回はpair_style lj/cut/coul/longだったのでLJポテンシャルの係数を意味している。
今回はSPC/E力場を使うので、σOO = 3.166, εOO = 0.1553, σHH = 0.000, εHH = 0.000となっている。
pair_coeffはインプットファイル内にも記述することができるがその場合は、pair_coeff 1 1 3.166 0.1553と記述する。dataファイル内に記述する場合は同一原子対のみのポテンシャルとなり、異種原子間のポテンシャルは総和平均（相乗平均）を取って計算することになる（計算方法はインプットファイル内mixで指定）。
Bond Coeffs
bond_styleで指定した関数系における結合長の力場係数を設定する。SPC/E力場はOH間の距離が1Åの剛体モデルなのでそのように設定する。
Angle Coeffs
angle_styleで指定した関数系における結合角の力場係数を設定する。SPC/E力場はHOHの角度が109.47°の剛体モデルなのでそのように設定する。
Atoms
原子座標。順番に原子番号（particle number）、分子番号、原子種番号、電荷、x座標、y座標、z座標が記されている。
Bonds
結合の設定。AとBが結合しているとしたら、順番に結合番号、結合の種類番号、Aの原子番号、Bの原子番号である。
Angles
結合角の設定。ABCの角度を記述する場合は、順番に結合角番号、結合角の種類番号、Aの原子番号、Bの原子番号、Cの原子番号である。
LAMMPSのインプットファイルの構造【バルクの水のMD計算】
更新サボっていて申し訳ありませんでした。久々の更新になります。今回は前回ダウンロードしたLAMMPSを使ってバルク水の分子動力学計算を行ってみます。
これはLAMMPSのインプットファイルの構造を学ぶだけなので、インプットファイルをどうやって作成するかというところは今のところは無視して構いません。
バルク水のMD計算のインプットファイル
units real atom_style full boundary p p p pair_style lj/cut/coul/long 10. 10. kspace_style pppm 1e-5 bond_style harmonic angle_style harmonic dihedral_style charmm read_data bulkwater.data neighbor 2.0 bin neigh_modify every 1 delay 0 check yes # 出力 dump 1 all custom 100 bulkwater.dump id type xs ys zs ix iy iz dump 2 all xtc 100 bulkwater.xtc # 速度の割り当て velocity all create 300.0 12345 run 0 velocity all scale 300.0 # 温度制御 fix 1 all shake 1e-4 1000 0 b 1 b 2 a 1 fix 2 all nvt temp 300.0 300.0 100. tchain 3 thermo_style multi thermo 100 # MDメイン計算 timestep 1.0 run 500000 write_restart bulkwater.restart
各コマンドの意味について以下で解説していきます。
units
シミュレーションに使用する単位の設定。例えばrealなら質量はg/mol、距離はÅを使用する。他の単位については詳しくは以下のリンクを参照してください。
real以外によく使用するコマンドはmetalとljがあるが、MD計算に使用するのは基本的にrealとmetalです。ソフトマターはrealで金属はmetalという考え方でいいと思います。
ここでは水なのでrealを使用しますが、いずれはmetalの場合も例示したいと思います。
http://lammps.sandia.gov/doc/units.html
atom_style
シミュレーションに使用する原子の型を決める。これも様々な設定方法がある（詳細は下記リンク参照）が、基本的に使用するのはunits realの場合はfull、units metalの場合はatomicで良いと思います。
大は小を兼ねるなので、あまり制限されたものを使用するとあとで大変なので大きい型（full）を使うのがよいはず。
http://lammps.sandia.gov/doc/atom_style.html
boundary
境界条件の決定。boundary x y zで各方向の境界条件を定める。pは周期境界であり基本はこれを使用することが多いです。他にはfやmなどがあり、前者は壁に到達すると原子を消去し、後者は原子が遠く離れてもそこを境界に再設定（Shrink-wrapping）する手法。
http://lammps.sandia.gov/doc/boundary.html
pair_style
原子間のポテンシャルを記述するのだが、ポテンシャルはMD計算の要なので、ここはかなり多数のコマンドが用意されているのでこれも詳しくは以下のマニュアルを参照してください。
http://lammps.sandia.gov/doc/pair_style.html
今回はlj/cut/coul/longを使用します。この意味は短距離相互作用にはLennard-Jonesポテンシャルを、長距離相互作用のクーロン力にはEwald法、もしくはPPPM法を使用するという意味です。
数字はカットオフ値でlj/cut/coul/long 10. 10.の場合、LJポテンシャルを10Åまで計算、クーロン力は10Åまで直接計算した後、それ以外のところでは上記手法などで計算することを意味しています。
原子のMDの場合はLJポテンシャルは10Å程度、クーロン力はEwald(PPPM)法を用いるのが最も多いです。
kspace_style
前項で述べた長距離相互作用の計算はｋ空間で行うので、ここはｋ空間の設定をします。pppmはparticle-particle-particle-mesh Ewald法を意味しており、原子数が大きな系では計算時間がEwald法よりも速くなるという利点があります。Ewaldは通常のEwald法を使用します。その次の数字は許容される計算誤差で小さい方が精度が高く計算できます。
出先でパソコンの充電がやばいのでまたあとで追記します。
充電復活しましたので続きです。
bond_style
原子内のボンド間のポテンシャルを設定する。これも様々なコマンドが存在するので詳しくは以下のリンクを参照してください。
http://lammps.sandia.gov/doc/bond_style.html
使用頻度が高いのはharmonicでE=K(r-r0)で表される調和振動子ポテンシャルです。
その他に比較的よく使われるのはMorseポテンシャルmorseでE=D[1-exp{-α(r-r0)}]と表されます。
angle_style
原子内の結合角ポテンシャルの設定。こちらもbond_styleと同様にharmonicを使用することが多い。表式はE=K(θ-θ0)となる。
他にも様々な記述方法があるのでリンク参照。
http://lammps.sandia.gov/doc/angle_style.html
dihedral_style
二面角ポテンシャルの設定。二面角の場合はcharmmを使用することが多い。式はE=K{1-cos(nφ-d)}と表される。
harmonicも存在するが、あまり使われない。（表式としてはcharmmに似ている）
http://lammps.sandia.gov/doc/dihedral_style.html
read_data
Geometricalファイルの読み込み。Geomファイルの作り方は後日。
neighbor
粒子からneighbor listを構成するまでの距離と方法の指定を行う。
これをいじったらどう変化するのかよく分からん。色んな例を見ているとどれもneighbor 2.0 binになっているます。
http://lammps.sandia.gov/doc/neighbor.html
neigh_modify
neighbor listをどのタイミングで更新するかの設定。
例えばevery 1 delay 0 check yesなら0ステップ後に1ステップ毎に更新。（ただし実際はneighborで設定した距離の半分を移動したら更新する。）
http://lammps.sandia.gov/doc/neigh_modify.html
dump
dumpはアウトプットのフォーマットを指定します。Syntaxはdump ID group-ID style N file argsであり簡単な例はdump 2 all xtc 100 bulkwater.xtcです。
allは全原子に対してのアウトプットを意味しており、xtcはアウトプットのファイルフォーマット、100はアウトプット間隔、bulkwater.xtcはアウトプットファイル名です。
ファイル名の後に、アウトプットする属性を指定することもできる。上記例ではid type xs ys zs ix iy izとなっているが、idは原子番号、typeは原子種番号、xs ys zsは原子位置、ix iy izは周期境界条件のどこの仮想ボックスに位置しているかを出力します。
それ以外にも例えばvx vy vzは原子の速度、fx fy fzは原子に働く力、qは原子電荷を出力します。それ以外のパラメータに関しては下記リンク参照。
http://lammps.sandia.gov/doc/dump.html
やったことはないが画像のアウトプットもできるらしい。構文はdump imageです。
http://lammps.sandia.gov/doc/dump_image.html
velocity
速度の割り当て。グループallの原子に対して300.0Kの温度を乱数12345で割り当てることを意味しています。
ここでshake法（後述）を使用する際は、速度を上手く割り当てられないことに注意してください。それを回避するには、
velocity all create 300.0 12345 run 0 velocity all scale 300.0
と、一度速度を割り当ててから計算を0ステップ実行して、その後300Kに再スケールを行うという方法があります。
http://lammps.sandia.gov/doc/velocity.html
fix
fixコマンドはかなり多くのことができるので個別に記事を書きたいと思います。ここではSHAKE法の設定と、温度制御を行っています。
SHAKE法は水素原子を含む結合を固定する方法で、今回は水の計算にSPC/E力場を用いているので、shake法を使用します。fix 1 all shake 1e-4 1000 0 b 1 b 2 a 1はall分子に対して1e-4の精度でSHAKE法を実行することを意味している。1000はSHAKE方の反復計算回数。0はSHAKE法を使用している結合と角度の情報の出力を何ステップごとに行うかということです。（0は出力しない。）
bは結合長の制限、aは角度の制限です。ここでは結合1と結合2を固定、角度1を固定しています。
http://lammps.sandia.gov/doc/fix_shake.html
また、fix 2 all nvt temp 300.0 300.0 100. tchain 3は熱浴の設定です。今回はNVTアンサンブルを使用していますが、NPTアンサンブルnptやNPHアンサンブルnphも使用できます。
tempは温度制御の宣言。NPTの場合は圧力も制御するのでisoで指定します。
300.0 300.0は初期温度と最終温度で100.は温度緩和の時間設定です。NVTアンサンブルの場合LAMMPSはNose-Hoover法を使用しているので温度緩和は100にすると良いとのこと。
tchainは熱浴における振動制御のパラメータでNose-Hoover法の場合デフォルトは3。
http://lammps.sandia.gov/doc/fix_nh.html
thermo_style
出力する熱力学量の設定。multiにするとetotal ke temp pe ebond eangle edihed eimp evdwl ecoul elong pressが出力される（全エネルギー、運動エネルギー、温度、全ポテンシャルエネルギー、結合エネルギー、結合角のエネルギー、二面角のエネルギー、improperな二面角のエネルギー、ファンデルワールスエネルギー、クーロンエネルギー、k空間のエネルギー、圧力）。
http://lammps.sandia.gov/doc/thermo_style.html
themo
アウトプット間隔の設定。
timestep
計算の時間刻みの設定。unitsで指定した単位系で設定する。realの場合はfs単位。水素原子が存在するので1fs刻みにする。
run
計算時間の設定。realの場合はfs単位。
write_restart
リスタートファイルの作成。
LAMMPSのインストール
今日からLAMMPSを使っていきます。
LAMMPSとはLarge-scale Atomic/Molecular Massively Parallel Simulatorの略で分子動力学法の無償プログラムです。分子動力学法の無償プログラムは色々ありますが、LAMMPSが一番広範囲に使用できるので良く使われています。しかしながらLAMMPSの日本語のチュートリアルは存在しませんし、日本語で使用方法を検索しようとしてもほぼ必要な情報が出てきません。
そこでこのブログで使用方法をまとめていきます。ただ、私自身もLAMMPSを初めて使うので分からないところが山程あります。誰か教えてくれる人が近くにいれば良いのですが、いないので完全に独学になりますので、ここに書かれていること自体が正しくない可能性もあります。
その点に注意して、参考にしていただければうれしいです。
ではLAMMPSをインストールします。LAMMPSはWindowsでも使用できるので、今回はWindowsでインストールします。LinuxやMacOSでもインストールできるので、お使いのOSに合わせてダウンロードしてください。
Windowsのダウンロードは以下のサイトから
http://rpm.lammps.org/windows.html
最新のものをダウンロードすればよいかと思います。MPIは並列計算を使用する場合で、今回は練習のために通常バージョンをダウンロードしておきます。
exeファイルを実行すればインストールがスタートします。（画像はMPIバージョンのもの）
Completedが出れば完了です。数分で終了します。
MPIバージョンの場合は以下のサイトからMPICHもダウンロードしてインストールしてください。
http://www.mpich.org/
今日はここまで。
Pythonで常微分方程式を解く（ルンゲクッタ法、オイラー法、解析的手法）
今日はPythonで常微分方程式を解いてみます。
空気抵抗がある状態でのボールの自由落下における速度と軌跡を求めましょう。
数値計算
数値計算では4次のルンゲクッタ法を使用して計算します。
ルンゲクッタ法は数値解析において微分方程式の初期値問題に対して比較的良い近似解を与える方法ですが、それに関しては詳しくはWikipediaなどをご覧ください。
import numpy as np import matplotlib.pyplot as plt g = 9.81 #重力加速度(m/s^2) k = 0.24 #空気抵抗係数(kg/m) m = 10.0 #ボール重さ(kg) def f(v): #関数の定義 return (m*g-k*v**2)/m vtrj = [0.0] #速度のtrajectory配列 xtrj = [100.0] #位置のtrajectory配列 array_t = [0.0] #時間発展配列 v = vtrj[0] #初期速度 x = xtrj[0] #初期位置 t = array_t[0] #初期時間 h = 0.1e0 #ルンゲクッタ法の刻み幅 #ルンゲクッタ法のメイン計算 while x >= 0.0e0: k1 = h*f(v) k2 = h*f(v+0.5*k1) k3 = h*f(v+0.5*k2) k4 = h*f(v+k3) v = v + (k1+2*k2+2*k3+k4)/6 x = x - v*h t = t + h vtrj = np.append(vtrj,v) xtrj = np.append(xtrj,x) array_t = np.append(array_t,t) #グラフ描写 plt.plot(array_t,vtrj) plt.plot(array_t,xtrj)
終端速度が20.2m/s、100m地点から落とすと6.4秒後に地面に到達することが求められた。
位置の時間発展にはオイラー法を使用している。（ただしこの場合はルンゲクッタ法でも同じ結果が得られる。）
オイラー法の詳しい解説はWikipediaなどを参照してください。
解析的手法
Sympyを使用して解析的に常微分方程式を解くことも可能である。
import numpy as np import matplotlib.pyplot as plt from sympy import * v = Function('v') t = Symbol('t') g = 9.81 #重力加速度(m/s^2) k = 0.24 #空気抵抗係数(kg/m) m = 10.0 #ボール重さ(kg) eq1 = dsolve(v(t).diff(t,1)-(m*g-k*v(t)**2)/m) print(eq1)
dsolveを使用して微分方程式を解いている。
Eq(t + 1.03045701422812*log(1.0*v(t) - 20.2175666191557) - 1.03045701422812*log(1.0*v(t) + 20.2175666191557), C1)
しかし解析解は上のように直感的に理解しにくいのと、場合によっては求められないときもあるので基本的には数値的解法をするのが良いと思う。
今日はここまで。
Pythonでグラフの描画（matplotlib）
今日はグラフの描写を簡単に勉強します。
たぶん科学技術計算をやるような人はグラフはKaleidaGraphとかIGORとかをお使いになっている方が多いような気がするので、とりあえずさわりの部分だけここには書きます。
import numpy as np import matplotlib.pyplot as plt for i in range (0,101): x = i * 0.1e0 y = np.sin(x) plt.plot(x,y,'ok')
まずimport matplotlib.pyplot as pltでmatplotをインポートします。
plt.plot(x,y,'ok')でプロットします。クォート内はプロット方法の指定で、'o'は丸マーカー、'k'は黒色を意味します。
他にも様々な指定の方法があるので、以下のサイトを参考にしてください。
matplotlib で指定可能なマーカーの名前 – Python でデータサイエンス
matplotlib で指定可能な色の名前と一覧 – Python でデータサイエンス
何も指定しない場合は自動的に実線になり、色も決定されます。
import numpy as np import matplotlib.pyplot as plt xpl = [] ypl = [] for i in range (0,101): x = i * 0.1e0 y = np.sin(x) xpl = np.append(x,xpl) ypl = np.append(y,ypl) plt.plot(xpl,ypl)
画面に表示されない場合は
最後にplt.show(xpl,ypl)を付け加えてみてください。
一応、Sympyを使用して変数のままでもグラフは描ける。
import matplotlib.pyplot as plt from sympy import * x = Symbol('x') y = Symbol('y') y = sin(x) plotting.plot(y,(x,0,10)) #0<x<10で描写
今日はここまで。
Pythonで極限、微分、積分の計算（Sympy）
今日は前回に引き続きSympyを使った代数計算を勉強します。
極限
まずはsin(x)/xの無限大での極限を調べてみます。
from sympy import * x = Symbol('x') f = limit(sin(x)/x,x,oo) print (f)
極限はlimitを使用し、カッコ内に関数,極限に近づける変数,近づける値の順に書いていきます。
ooは無限大を意味しています。
0
等比級数の和の極限も計算してみよう。an=(1/2)nの極限は
from sympy import * i = Symbol('i',integer=True) #整数型宣言 n = Symbol('n',integer=True) p = Rational(1,2) f = summation(p**i,(i,1,n)) #summationは和を求める関数 l = limit(f,n,oo) print (l)
1
微分
f(x)=log(x)の微分を求めてみる。
微分は曲線の傾きに等しいので
from sympy import * x = Symbol('x') f = (log(x+y)-log(x))/y l = limit(f,y,0) print (l)
ともできるが、diff関数を使えば簡単に計算できる。
from sympy import * x = Symbol('x') g = log(x) f = diff(g) print (f)
答えはどちらも
1/x
偏微分はdiffで文字を指定すれば良い。
from sympy import * x = Symbol('x') y = Symbol('y') g = log(x+y) f = diff(g,x) #ここだけ変化した print (f)
1/(x + y)
n階微分は
from sympy import * x = Symbol('x') g = log(x) f = diff(g,x,3) #今回は3階微分 print (f)
2/x**3
積分
積分もできます。積分ってかなり難しいですがPythonはすごいですね。
今回もlog(x)の積分をしてみます。
from sympy import * x = Symbol('x') g = log(x) f = integrate(g) print (f)
x*log(x) - x
大学入試頻出の積分がこんなに簡単にできちゃいました。
定積分は範囲を指定してあげれば良くて
f = integrate(g,(x,1,exp(1)))
とします。
1
積分範囲の極限を取る場合は
from sympy import * x = Symbol('x') g = exp(-x)*sin(x) f = integrate(g,(x,0,oo)) print (f)
とできる。
1/2
今日はここまで。
Pythonで様々な代数計算（Sympy）
これまで配列の計算を行ってきましたがSympyライブラリを使用すれば、Pythonは代数計算にも威力を発揮します。
面倒くさい方程式などを解くときに利用できますので、今回はSympyを使用して簡単な代数計算を行ってみます。
分数型
Sympyを使用すれば、分数型を使用することができます。
from sympy import * #Sympyをimportします a = Rational(1,3) #分数型宣言Rational print (a)
1/3
Numpyなどでは0.333333…となっていたところが1/3として扱うことができます。
分数計算もこの通り
from sympy import * a = Rational(1,3) b = Rational(2,5) c = a + b print (c)
11/15
変数の取扱い
from sympy import * x = Symbol('x') y = Symbol('y') f = x + y + 2*y - 4*x +2*x*y print (f)
2*x*y - 3*x + 3*y
Symbol関数を使用することで変数を変数のまま計算することも出来る。
展開も
from sympy import * x = Symbol('x') y = Symbol('y') f = expand((x + y)**4) print (f)
x**4 + 4*x**3*y + 6*x**2*y**2 + 4*x*y**3 + y**4
f = expand(cos(x + y), trig=True)
とtrig=Trueを加えれば
-sin(x)*sin(y) + cos(x)*cos(y)
三角関数の展開も出来る（加法定理）。
方程式を解く
1元n次方程式
x2-6x-9=0を解いてみる
from sympy import * x = Symbol('x') f = x**2 + -6 * x -9 s = solve(f,x) #f=0として計算する print (s)
[3 + 3*sqrt(2), -3*sqrt(2) + 3]
solveにより方程式を解ける。
n元n次連立方程式
x2+y2=0
xy=1
の連立方程式を解こう。
from sympy import * x = Symbol('x') y = Symbol('y') f1 = x**2 + y**2 f2 = x*y - 1 s = solve([f1,f2],[x,y]) print (s)
[(-(-sqrt(2)/2 - sqrt(2)*I/2)**3, -sqrt(2)/2 - sqrt(2)*I/2), (-(-sqrt(2)/2 + sqrt(2)*I/2)**3, -sqrt(2)/2 + sqrt(2)*I/2), (-(sqrt(2)/2 - sqrt(2)*I/2)**3, sqrt(2)/2 - sqrt(2)*I/2), (-(sqrt(2)/2 + sqrt(2)*I/2)**3, sqrt(2)/2 + sqrt(2)*I/2)]
solve([f1,f2],[x,y])と書けばx,yのペアで解いてくれる。
ある文字について解きたいとき
exp(x+y)=1をyについて解いてみる。
from sympy import * x = Symbol('x') y = Symbol('y') f = exp(x + y) - 1 s = solve(f,y) print (s)
[log(exp(-x))]
今日はここまで。
Pythonでファイルの読み込み（文字列を含むファイル、連番ファイル）
科学技術計算では他のファイルから読み込んで解析するということが多いのでファイルの読み込みについて勉強していきます。
文字列を含むファイル
,Apple,Orange,Pineapple Day1,5,10,3 Day2,8,13,0 Day3,2,16,3
上記のようなtest.csvファイルを読み込んで、平均値を出してみます。
import numpy as np file = open('test.csv','r') fruit = file.readline() #最初の1行を文字列として読み取る fruit = fruit.replace('\n','') #readlineは改行も読み取ってしまうので置き換え fruit = fruit.split(',') #カンマで分割してリストにする num = np.genfromtxt('test.csv',skip_header=1,delimiter=',') #数値を読み取る(skip_headerはn行目を飛ばす） num = np.transpose(num) for i in range (1,4): ave = np.average(num[i]) #average関数で平均を求める print(fruit[i],ave) file.close
Apple 5.0 Orange 13.0 Pineapple 2.0
skip_headerにより邪魔な文字列を無視して解析することができます。
連番ファイル
,Apple,Orange,Pineapple Day1,1,1,91 Day2,1,2,62 Day3,1,0,120
このようなパイナップル人気が異常に高いパラレルワールドでの販売個数をtest2.csvとします。
これを連続して読み込むためには
import numpy as np file = open('test.csv','r') fruit = file.readline() fruit = fruit.replace('\n','') fruit = fruit.split(',') for i in range (1,3): #2つのファイルを読み込む num = np.genfromtxt('test%d.csv'%i,skip_header=1,delimiter=',') #%dと％iがポイント num = np.transpose(num) for i in range (1,4): ave = np.average(num[i]) print(fruit[i],ave) file.close
for内に'test%d.csv'%iと設定してあげることで順番にファイルを読み込んでいきます。
Apple 5.0 Orange 13.0 Pineapple 2.0 Apple 1.0 Orange 1.0 Pineapple 91.0
もし連番ファイルが01,02,…という形式であれば
for i in range (1,3): num = np.genfromtxt('test%02d.csv'%i,skip_header=1,delimiter=',')
としましょう。
今日はここまで。
Pythonにおける行列の固有値問題の解法、行列の対角化
今日は行列の固有値問題を解いて、さらに行列の対角化をしてみます。
固有値問題
まずは固有値問題から。固有値問題とはAx=λxを満たすような固有ベクトルxと固有値λを求めるような問題のことです。
他のプログラムでは面倒くさいですがPythonなら簡単に求まります。
import numpy as np list_a = (0,14,2) list_b = (-1,9,-1) list_c = (-2,4,8) mat_a = np.array([list_a,list_b,list_c]) eigv_a,eigm_a=np.linalg.eig(mat_a) print(eigv_a) print(eigm_a)
[4. 6. 7.] [[ 9.42809042e-01 -9.12870929e-01 -8.94427191e-01] [ 2.35702260e-01 -3.65148372e-01 -4.47213595e-01] [ 2.35702260e-01 -1.82574186e-01 -5.57327535e-16]]
固有値と固有ベクトルの順番は同じです。そして固有ベクトルは1に規格化されています。
すなわち
固有値4のとき固有ベクトル(9.42809042e-01,2.35702260e-01,2.35702260e-01)
人間にわかりやすくすれば(4,1,1)
固有値6のとき固有ベクトル(-9.12870929e-01,-3.65148372e-01,-1.82574186e-01)
人間にわかりやすくすれば(5,2,1)
固有値6のとき固有ベクトル(-8.94427191e-01,-4.47213595e-01,-5.57327535e-16)
人間にわかりやすくすれば(2,1,0)です。
そしてP-1AP=BのPがeigm_aになっています。
eigm_a = np.transpose(eigm_a)
np.transposeで転置行列を求めましょう。
[[ 9.42809042e-01 2.35702260e-01 2.35702260e-01] [-9.12870929e-01 -3.65148372e-01 -1.82574186e-01] [-8.94427191e-01 -4.47213595e-01 -5.57327535e-16]]
最初から綺麗に出したいときは
import numpy as np list_a = (0,14,2) list_b = (-1,9,-1) list_c = (-2,4,8) mat_a = np.array([list_a,list_b,list_c]) eigv_a,eigm_a = np.linalg.eig(mat_a) eigm_a[np.abs(eigm_a) <= 1.0e-10] = 0.0 #端数を削除するため10^-10以下は0とみなす eigm_a = np.transpose(eigm_a) eigm_abs = np.abs(eigm_a) #行列の絶対値を求める eigmin = eigm_abs.min(axis=1) #行列の絶対値の各行の最小値を求める for i in range (0,3): if eigmin[i] == 0.0: eigmin[i] = eigm_abs[i,0] #各行の最小値が0の場合はとりあえず1番目 for i in range (0,3): eigm_a[i,] = eigm_a[i,]/eigmin[i] #規格化する print(eigm_a)
というのを考えました。汚いコードになってしまったので、もっと綺麗なプログラムがあるはずです。みなさんで作ってください。
[[ 4. 1. 1. ] [-5. -2. -1. ] [-1. -0.5 0. ]]
今日はここまで。
Pythonで逆行列、行列式、rankの計算
今日はNumpyの機能を最大限に活かして行列の基本的な計算をいくつか解いてみたいと思います。
逆行列計算
import numpy as np list_a = (0,4) list_b = (1,5) mat_a = np.matrix([list_a,list_b]) mat_b = np.linalg.inv(mat_a) print(mat_b) mat_c = mat_a * mat_b print(mat_c)
[[-1.25 1. ] [ 0.25 0. ]] [[1. 0.] [0. 1.]]
簡単に逆行列計算ができました。np.linalg.inv関数で逆関数をたちどころに計算することができます。
しかしながら逆計算をプログラムで解くのは注意が必要です。例えば以下のような場合
import numpy as np list_a = (3.1415,9.2653) list_b = (5.8979,3.2384) mat_a = np.matrix([list_a,list_b]) mat_b = np.linalg.inv(mat_a) print(mat_b) mat_c = mat_a * mat_b print(mat_c)
[[-0.07281823 0.2083383 ] [ 0.13261939 -0.07063935]] [[ 1.00000000e+00 -1.11022302e-16] [ 0.00000000e+00 1.00000000e+00]]
逆行列との積が単位行列になりません。これは計算途中で丸め込みの誤差が生じるためだと思われます。逆行列計算には過信しないようにしましょう。
行列式の計算
import numpy as np list_a = (-1,4,1) list_b = (1,1,2) list_c = (0,9,1) mat_a = np.matrix([list_a,list_b,list_c]) det_a = np.linalg.det(mat_a) print(det_a)
22.000000000000004
np.linalg.detが行列式を求める関数です。
実際の答えは22なのですが、変な端数が出てきています。原因がよく分からないのですが、どうもアルゴリズムがLU分解をしているせいではないかと思っているのですが詳しい方教えてください。とりあえず概算は求めることができます。
rankの計算
import numpy as np list_a = (-1,4,1) list_b = (1,1,2) list_c = (0,9,1) mat_a = np.matrix([list_a,list_b,list_c]) rank_a = np.linalg.matrix_rank(mat_a) print(rank_a)
3
階級(rank)も求められます。rankとはなにかを分かりやすく説明するのは避けますが、簡単に言えば、正則な小行列で大きさが最大なものの列数です。
今日はここまで。
Pythonの基礎（ループ処理、条件分岐）
あらゆるプログラミングで重要なループ処理と条件分岐を勉強します。
for文
Pythonでのループ処理はforを使用する。
1から100までの和を求めるプログラムを作ってみよう。
import numpy as np sum=0 for i in range (1,101): sum = sum + i print(sum)
5050
ループ内に入れたいものをスペース4つで下げます。下がっている文はループ内で処理される。
注意してほしいのはPythonは最初が0からスタートしているので、100まで足すときには範囲を1から101にしなければなりません。
またfor文末のコロンを忘れないようにしましょう。
while文
ある条件の間に限りループを実行したいときはwhileを使用する。
例えば200から7ずつ引いていって、100を切るまでの各値を計算してみよう。
import numpy as np rest = 200 array_r = np.array(rest) #表示を横に並べたいので配列化する while rest > 100: rest = rest - 7 array_r = np.append(array_r,rest) #計算した結果を配列に追加する print(array_r)
[200 193 186 179 172 165 158 151 144 137 130 123 116 109 102 95]
if文
数値が偶数か奇数かを判定するプログラムをつくってみます。
すこし欲張りしてコンソールから数字を読み取って偶数か奇数かを判定するソフトにしましょう。
import numpy as np num = input('enter a value :') #数字の読み取り num = int(num) #整数型への変換 sur = np.mod(num,2) if sur == 0: print('Even') else: print('Odd')
inputはコンソールの入力を受け取る関数です。inputではまず文字列型としてデータを受け取るので一度整数型に直すためにint関数を使用しています。
if文のところは他のプログラムと大きく変わりません。コロンをつけるのとインデントに注意してください。3つ以上に分岐するときはif elif elseと使います。
ちょっと応用
for文やwhile文とif文が混ざったプログラムを書いてみましょう。
50番目までのフィボナッチ数列で偶数だけを抜き出すプログラムを書いてみます。
import numpy as np fib = [0,1] fib2 = [0] for i in range(1,51): num = fib[i] + fib[i-1] #i番目とi-1番目を足す fib = fib + [num] #fib配列に追加（np.appendと同じ） if num % 2 == 0: #%はnp.modと同じ fib2 = fib2 + [num] #偶数だけをfib2配列に追加 print(fib2)
[0, 2, 8, 34, 144, 610, 2584, 10946, 46368, 196418, 832040, 3524578, 14930352, 63245986, 267914296, 1134903170, 4807526976, 20365011074]
数が多いのでフィボナッチ数列が10000を超えたらループを停止するようにしてみます。
import numpy as np fib = [0,1] fib2 = [0] for i in range(1,51): num = fib[i] + fib[i-1] fib = fib + [num] if num % 2 == 0: fib2 = fib2 + [num] if num > 10000: break print(fib2)
[0, 2, 8, 34, 144, 610, 2584, 10946]
breakによりif文の中の条件を満たした場合、forループの外に抜けることができます。
今日はここまで。
Pythonによる行列の生成および四則演算
今日は行列の生成と四則演算を行う。
行列の生成
import numpy as np list_a = [0,1,2,3,4] list_b = [5,6,7,8,9] list_c = [10,11,12,13,14] matrix_a = np.array([list_a,list_b,list_c]) print(matrix_a)
[[ 0 1 2 3 4] [ 5 6 7 8 9] [10 11 12 13 14]]
3×5行列が生成できた。
txt,csvファイルからの読み取り
以下のようなcsvファイルから、これを読み取って行列にする。
0,1,2,3 4,5,6,7 8,9,10,11
import numpy as np matrix_a = np.genfromtxt('test.csv', delimiter=',', dtype=float) print(matrix_a)
今回は浮動小数点dtype=floatとして読み込んでみる。
[[ 0. 1. 2. 3.] [ 4. 5. 6. 7.] [ 8. 9. 10. 11.]]
単位行列、零行列の生成
零行列は
import numpy as np matrix_a = np.zeros((3,3)) print(matrix_a)
[[0. 0. 0.] [0. 0. 0.] [0. 0. 0.]]
単位行列は
import numpy as np matrix_a = np.eye(3) print(matrix_a)
[[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]]
で生成できる。
行列の四則演算
基本的には配列の計算と同じ
import numpy as np list_a = [0,1,2,3,4] list_b = [5,6,7,8,9] list_c = [10,11,12,13,14] matrix_a = np.array([list_a,list_b,list_c]) matrix_a = matrix_a + 1 print(matrix_a)
[[ 1 2 3 4 5] [ 6 7 8 9 10] [11 12 13 14 15]]
単純に数値を足すと全要素に対して足し算される。他の演算に関しても同様。
行列同士の演算は、行列の演算則に従うが、積に注意。
import numpy as np list_a = [0,1] list_b = [2,3] matrix_a = np.array([list_a,list_b]) matrix_b = np.array([list_b,list_a]) matrix_c = matrix_a * matrix_b print(matrix_c)
と書くと各要素同士を掛け合わせてしまう。
[[0 3] [0 3]]
行列の演算をするときは
import numpy as np list_a = [0,1] list_b = [2,3] matrix_a = np.array([list_a,list_b]) matrix_b = np.array([list_b,list_a]) matrix_c = np.dot(matrix_a,matrix_b) print(matrix_c)
としなければならない。
[[0 1] [4 9]]
もしくは、
import numpy as np list_a = [0,1] list_b = [2,3] matrix_a = np.matrix([list_a,list_b]) matrix_b = np.matrix([list_b,list_a]) matrix_c = matrix_a * matrix_b print(matrix_c)
とnp.matrixで行列を生成する方法でもよい。
[[0 1] [4 9]]
今日はここまで。
Pythonによる配列の生成および四則演算
数学の基本である配列の計算をまとめておきます。
配列の生成
import numpy as np list_a = [0,1,2,3,4,5] array_a = np.array(list_a) print(array_a)
1行目はNumpyのインポート。これで化学計算のプログラミングが簡単になる。
2行目はリストの設定。
3行目は作ったリストは配列だということをpythonさんサイドに認識させます。
実行はdebugボタンを押してください。(今回はSpyderを使っているが、コンソール上で行うときはpython test.pyなどを実行。)
[0 1 2 3 4 5]
配列が完成しました。 型宣言をすることもできます。例えば浮動小数点にするためにfloatを追加すると
import numpy as np list_a = [0,1,2,3,4,5] array_a = np.array(list_a, dtype = float) print(array_a)
[0. 1. 2. 3. 4. 5.]
小数点がついて浮動小数点扱いになったことが分かります。
txtやcsvからの読み取り
まず以下のようなtest.txtファイルをつくります。
0 1 2 3 4 5 6 7 8 9
これを読み取るときは
import numpy as np array_a = np.genfromtxt('test.txt', delimiter = ' ', dtype = int) print(array_a)
genfromtxtはファイルの読み込みを指示するコードです。
delimiterは区切り文字の指定で、test.csvのときはカンマ区切りなので
array_a = np.genfromtxt('test.csv', delimiter = ',', dtype = int)
とします。
配列の四則演算
配列にそのまま1を加算すると全ての要素に1が加算されます。
import numpy as np list_a = [10,20,30,40,50] array_a = np.array(list_a) array_a = array_a + 1 print(array_a)
[11 21 31 41 51]
これはあらゆる四則演算に共通
例えば剰余計算は
import numpy as np list_a = [15,24,58,32,48] array_a = np.array(list_a) array_a = np.mod(array_a,3) print(array_a)
[0 0 1 2 0]
となる。
配列同士の四則演算も同様に
import numpy as np list_a = [15,24,58,32,48] list_b = [1,3,5,2,7] array_a = np.array(list_a) array_b = np.array(list_b) array_c = array_a / array_b print(array_c)
[15. 8. 11.6 16. 6.85714286]
Fortranでは許されないがPythonでは自動的にfloatにしてくれます。
今日はここまで。
「Anaconda」のインストール
まずプログラミング環境を整備するために「Anaconda」をインストールします。
Anacondaはディストリビューションだなんだと説明しているものもありますが、そういうことは特に触れません。導入したら簡単に計算できるようになるもので良いと思います。
ダウンロード方法は
https://www.anaconda.com/download/
のサイトから自分のパソコンのOSにあったものをダウンロードします。
私はWindows7を使用するので、Windowsの64bitをダウンロードしました。
起動すると上のようなインストーラ画面が出るので指示に従って進めていきます。
途中にチェックボックスが出てきます。上はAnacondaのPATHを設定します。下はAnacondaを標準のPythonとして設定します。
よく分からなかったら全部チェックを外しておいていいです。
インストールが完了して起動すると上のような画面が出ます。
本ブログではSpyderをメインに使用するので、Spyderをクリックしましょう。
Spyderの画面はこんな感じ。左側はコードを書くエディタで、右下のがプログラムを走らせるコンソールです。
今日はここまで。
このブログを始めた理由
私は元々Fortranで科学技術計算をしていたのですが、最近の流行りはPythonのようで、ライブラリもPythonの方が多く揃っているようです。
そこでPythonを勉強しようと思ったのですが、科学技術計算に使用するための方法を解説しているようなネットの記事や専門書は、元々Pythonを使っていた人やプログラムに親しい人に向けて書いたものが多いような気がします。
私はFortranを使っていると言いましたが、がっつりプログラムを組んでいるわけではなくて解析に使用しているだけで、プログラミングに関しては初心者です。
無理な計算でなければExcelをポチポチやって解析しているレベルの人なので、そういう本やネットの記事を読んでもよくわかりません。
特にネットのプログラマーの多くの記事は、小難しい横文字を連発して私にはさっぱりです。
そこでプログラミング初心者でも実用上使えるレベルにだけ学べるようなサイトがあればいいと思い作りました。難しい言葉は一切使わないつもりです。私も分からないので。
本当はQiitaでやろうかと思ったのですが、なんかレベルが高そうだし雰囲気も怖かったので、すこし懐が深そうなはてなブログでやります。
私もプログラミング初心者なので私自身の備忘録として基本的には使っていこうと思いますが、他の見てくださった方も参考になれば嬉しいです。
