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表記の6つが返る.
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')