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
前項で述べた長距離相互作用の計算はk空間で行うので、ここはk空間の設定をします。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でやろうかと思ったのですが、なんかレベルが高そうだし雰囲気も怖かったので、すこし懐が深そうなはてなブログでやります。
私もプログラミング初心者なので私自身の備忘録として基本的には使っていこうと思いますが、他の見てくださった方も参考になれば嬉しいです。