Call pmd from python using nappy

python上でnappyを使ってFortranプログラムであるpmdのコア部分（のライブラリ）を呼び出すことができる．ただし，そのためには，nap/nappy/pmd/ 内で make pmd_wrapper を実行して，pmd_wrapperからpmdを呼び出す準備が必要．

nappy.pmd.PMD を用いる方法

次のようなコードで，そのディレクトリにある in.pmd を読み込んで，1000 stepのMDを実行する（しかしなぜかjupyter notebook上では途中でメッセージもなく止まってしまう問題が解決していない）．

nsys = nappy.io.read('pmdini')
pmd = nappy.pmd.PMD(nsys)
pmd.load_inpmd()
print(pmd.params)
pmd.run(nstp=1000, dt=1.0, ifdmp=0, iprint=0)
print(' Potential energy = ', pmd.get_potential_energy())
print(' Kinetic energy = ', pmd.get_kinetic_energy())
print(' Stress = ', pmd.get_stress())
nsys = pmd.get_system()

上記の出力は次のようになり，stressはVoigt表記の６つが返る．

 Potential energy =  -3450.1590261927245
 Kinetic energy =  0.001166161049854613
 Stress =  [-0.00056849 -0.00121267 -0.00180787  0.0001973   0.00035736  0.00028713]

ASEを用いる方法

nappyには nappy.interface.ase.pmdrun 内に PMD_ASE_Calculator というASE用のcalculatorがあるので，それを使ってASEでMDなどを実行することが可能．

from nappy.interface.ase.pmdrun import PMD_ASE_Calculator as PMDCalc
from nappy.elements import elements
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet
from ase import units
 
atoms = nsys.to_ase_atoms()  # NAPSystemオブジェクトからASE atomsオブジェクトに変換
atoms.calc = PMDCalc(force_type=['Morse', 'Coulomb', 'angular'],
                     cutoff_radius=6.0)  # calculatorを設定
MaxwellBoltzmannDistribution(atoms, temperature_K=300)
dyn = VelocityVerlet(atoms, 1*units.fs)
 
for i in range(10):
    dyn.run(10)
    epot = atoms.get_potential_energy() /len(atoms)
    ekin = atoms.get_potential_energy() /len(atoms)
    print(f' {i:5d}: Epot = {epot:6.3f} eV/atom, Ekin = {ekin:6.3f} eV/atom')