弾性係数（elastic constants）を得るには，計算したい弾性係数に対応した微小な変形を与えた場合のエネルギー変化を計算すれば良い．
ここで紹介する方法は cubic な系に限るが，上記を自動化するスクリプトでの例である．
Lattice constant のページで，平衡体積（格子定数）が得られていると仮定している．変形後のエネルギーの方が低くなってしまうと，正値の弾性係数が得られない．
pmdiniファイルにて正しい格子定数を設定する．elasticity.pyスクリプトを実行し，複数pmd実行のための準備を行う$ python /path/to/nap/nappy/elasticity.py prepare pmdinielast_##というディレクトリが作られる．その各ディレクトリの中に微小変形が加わったpmdiniファイルと，いくつかの設定が記載されているconf.elast.yamlファイルが入っている．- 各ディレクトリで
pmdを実行する．$ for d in elast_*; do echo $d ; cd $d; ln -s ../in.p* .; pmd > out.pmd ; cd ..; done - 最後にもう一度
elasticity.pyを実行し，弾性係数を計算する．$ python /path/to/nap/nappy/elasticity.py analyze pmdini strs.pmd C_ij [GPa]: 108.338 74.200 74.115 0.096 0.000 0.019 74.200 108.263 74.061 0.071 -0.001 0.013 74.115 74.061 108.275 0.039 -0.002 0.005 0.096 0.071 0.039 108.211 -0.039 0.055 0.000 -0.001 -0.002 -0.039 108.211 0.057 0.019 0.013 0.005 0.055 0.057 108.211 Spacegroup number = 1 ==> Triclinic C_ij [GPa]: 108.338 74.200 74.115 0.096 0.000 0.019 74.200 108.263 74.061 0.071 -0.001 0.013 74.115 74.061 108.275 0.039 -0.002 0.005 0.096 0.071 0.039 108.211 -0.039 0.055 0.000 -0.001 -0.002 -0.039 108.211 0.057 0.019 0.013 0.005 0.055 0.057 108.211 Bulk modulus (K) = 85.514 GPa shear modulus (G) = 53.145 GPa Poisson's ratio (nu) = 0.243 efinition of elastic moduli, see https://materialsproject.org/wiki/index.php/Elasticity_calculations 123 = c11 +c22 +c33 = 324.876 c231 = c12 +c13 +c23 = 222.376 c456 = c44 +c55 +c66 = 324.633 s123 = s11 +s22 +s33 = 0.062 s231 = s12 +s13 +s23 = -0.025 s456 = s44 +s55 +s66 = 0.028 Kv = (c123 +2*c231)/9 = 85.514 Kr = 1.0 /(s123 +2*s231) = 85.514 Gv = (c123 -c231 +3*c456)/15 = 71.760 Gr = 15.0 /(4.0*s123 -4.0*s231 +3.0*s456) = 34.531 K = (Kv +Kr)/2 G = (Gv +Gr)/2 nu = (3*K -2*G)/(6*K +2*G) elasticity.py analyze done