Skip to content

Commit

Permalink
fix: a fatal performace issue in PMD_ASE_Calculator.
Browse files Browse the repository at this point in the history
  • Loading branch information
ryokbys committed May 5, 2024
1 parent 5da437b commit 04655b6
Showing 7 changed files with 99 additions and 48 deletions.
87 changes: 61 additions & 26 deletions nappy/interface/ase/pmdrun.py
Original file line number Diff line number Diff line change
@@ -1,8 +1,9 @@
"""
ASE interface to pmd.
"""
from __future__ import print_function
Changes:
<2024.05.04 RK> Changed to call pmd via f2py and dynamic library instead of via subprocess.
"""
import os
import subprocess
import numpy as np
@@ -12,15 +13,30 @@
from nappy.interface.ase.pmdio import get_fmvs
from nappy.napsys import NAPSystem

try:
import nappy.pmd.pmd_wrapper as pw
except:
pass

__author__ = "Ryo KOBAYASHI"
__version__ = "160605"
__version__ = "240504"
__LICENSE__ = "MIT"

CALC_END_MARK = "Job finished "
some_changes = ['positions', 'numbers', 'cell',]

def str2char(string,nlen):
slen = len(string)
if slen > nlen:
raise ValueError('slen > nlen !')
c = np.empty(nlen,dtype='c')
c[:] = ' '
c[0:slen] = [ string[i] for i in range(slen) ]
c[slen:] = [ ' ' for i in range(nlen-slen) ]
return c

class PMD(FileIOCalculator):

class PMD_ASE_Calculator(Calculator):
"""
Class for PMD calculation for ASE.
@@ -76,7 +92,7 @@ class PMD(FileIOCalculator):
'boundary': 'ppp',
}

def __init__(self, restart=None, ignore_bad_restart_file=False,
def __init__(self, restart=None,
label='pmd', atoms=None,
command='pmd > out.pmd',
dimension=(True,True,True),
@@ -104,10 +120,11 @@ def __init__(self, restart=None, ignore_bad_restart_file=False,
"""

Calculator.__init__(self, restart, ignore_bad_restart_file,
label, atoms, **kwargs)
Calculator.__init__(self, restart=restart,
label=label, atoms=atoms,
directory='.', **kwargs)

if label not in ['pmd']:
if label != 'pmd':
raise RuntimeError('label must be pmd.')

if self.parameters['force_type'] is None:
@@ -120,7 +137,6 @@ def __init__(self, restart=None, ignore_bad_restart_file=False,
else:
self.command = command +' > out.'+self.label

self.specorder= specorder
self.dimension = dimension

def set(self, **kwargs):
@@ -131,19 +147,36 @@ def set(self, **kwargs):
def calculate(self, atoms=None, properties=['energy'],
system_changes=some_changes):
Calculator.calculate(self, atoms, properties, system_changes)
self.write_input(self.atoms, properties, system_changes)

olddir = os.getcwd()
try:
os.chdir(self.directory)
errorcode = subprocess.call(self.command, shell=True)
finally:
os.chdir(olddir)

if errorcode:
raise RuntimeError('%s returned an error: %d' %
(self.name, errorcode))
self.read_results()
# # Write in.pmd to be read from pmd executable
# self.write_input(self.atoms, properties, system_changes)

# olddir = os.getcwd()
# try:
# os.chdir(self.directory)
# errorcode = subprocess.call(self.command, shell=True)
# finally:
# os.chdir(olddir)

# if errorcode:
# raise RuntimeError('%s returned an error: %d' %
# (self.name, errorcode))
# self.read_results()

#...The code below is for calling pmd via f2py and dynamic library
if atoms is not None:
self.nsys = nappy.io.from_ase(atoms, get_forces=False)

#...Call pmd via f2py and dynamic library
pmd = nappy.pmd.PMD(self.nsys)
pmd.set_params(force_type = self.parameters['force_type'],
cutoff_radius = self.parameters['cutoff_radius'],)
pmd.run(nstp=0, iprint=0)

self.results['energy'] = pmd.get_potential_energy()
self.results['forces'] = pmd.get_forces().T
self.results['stresses'] = pmd.get_stress()
self.results['stress'] = pmd.get_stress().sum() /3
return None

def relax(self, atoms=None, properties=['energy'],
system_changes=some_changes,
@@ -293,9 +326,9 @@ def get_input_txt(params,fmvs):
'temperature_relax_time','flag_temp_dist','',
'factor_direction','',
'stress_control','pressure_target','stress_target',
'stress_relax_time','flag_compute_stress','',
'stress_relax_time','',
'overlay','overlay_type','',
'zload_type','final_strain','',
'final_strain','',
'boundary']

int_keys=['num_nodes_x','num_nodes_y','num_nodes_z','num_omp_threads',
@@ -310,7 +343,7 @@ def get_input_txt(params,fmvs):
'converge_eps','final_strain']
str_keys=['io_format','force_type','temperature_control',
'stress_control','flag_temp_dist',
'flag_compute_stress','zload_type','boundary',
'boundary',
'overlay_type']

for key in order:
@@ -332,7 +365,9 @@ def get_input_txt(params,fmvs):
elif key == 'temperature_target':
vals = params[key]
for i,v in enumerate(vals):
txt += '{0:25s} {1:2d} {2:6.1f}\n'.format(key,i+1,v)
#...Format of temperature target has been changed...
# txt += '{0:25s} {1:2d} {2:6.1f}\n'.format(key,i+1,v)
txt += f'{key:25s} {v:6.1f}\n'
elif key == 'factor_direction':
# vals = params[key]
vals= fmvs
12 changes: 7 additions & 5 deletions nappy/io.py
Original file line number Diff line number Diff line change
@@ -1413,7 +1413,7 @@ def get_PDB_txt(nsys,**kwargs):

return txt

def from_ase(atoms,specorder=None):
def from_ase(atoms, specorder=None, get_forces=True):
"""
Convert ASE Atoms object to NAPSystem object.
"""
@@ -1445,10 +1445,12 @@ def from_ase(atoms,specorder=None):
vels = np.zeros((natm,3))
else:
vels = np.array(vels)
try:
frcs = atoms.get_forces()
except:
frcs = np.zeros((natm,3))
frcs = np.zeros((natm,3))
if get_forces:
try:
frcs = atoms.get_forces()
except:
pass

#...Create arrays to be installed into nsys.atoms
sids = [ nsys.specorder.index(si)+1 for si in symbols ]
12 changes: 11 additions & 1 deletion nappy/pmd/__init__.py
Original file line number Diff line number Diff line change
@@ -74,7 +74,10 @@ def run(self, nstp=0, dt=1.0, ifdmp=0, dmp=0.99, conveps=1e-5, convnum=3,
ispcs = self.nsys.atoms.sid.values
#...Run pmd by calling fortran-compiled library
#print('calling pw.run')
res = pw.run(rtot.T,vtot.T,naux,hmat,ispcs,initialize)
try:
res = pw.run(rtot.T,vtot.T,naux,hmat,ispcs,initialize)
except:
raise
#print('out from pw.run')
self.result = {}
self.result['rtot'] = res[0]
@@ -259,6 +262,13 @@ def get_stress_tensor(self):
return None
return self.result['stnsr']

def get_stress(self):
if not hasattr(self,'result'):
return None
stnsr = self.result['stnsr']
return np.array([stnsr[0,0], stnsr[1,1], stnsr[2,2],
stnsr[1,2], stnsr[0,2], stnsr[0,1]])

def get_forces(self):
if not hasattr(self,'result'):
return None
8 changes: 5 additions & 3 deletions nappy/pmd/pmd_wrapper.F90
Original file line number Diff line number Diff line change
@@ -134,6 +134,11 @@ subroutine set_pmdvars(nsp0,ns,ls,cspcs,nf,lf,cfrcs,rc0,rbuf0, &
force_list(i) = trim(c128)
end do

rc1nn = 3d0
rbuf = rbuf0
!.....RC should be set before calling init_force
rc = rc0

!.....It is required to call init_force and read some in.params.XXX to define aux array
call init_force(.true.)
!.....Before allocating auxiliary array, set naux (num of auxiliary data)
@@ -154,9 +159,6 @@ subroutine set_pmdvars(nsp0,ns,ls,cspcs,nf,lf,cfrcs,rc0,rbuf0, &
endif
enddo
dt = dt0
rbuf = rbuf0
rc1nn = 3d0
rc = rc0
nx = 1
ny = 1
nz = 1
2 changes: 1 addition & 1 deletion pmd/force_Morse.F90
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
module Morse
!-----------------------------------------------------------------------
! Last modified: <2023-12-18 16:13:18 KOBAYASHI Ryo>
! Last modified: <2024-05-04 22:43:08 KOBAYASHI Ryo>
!-----------------------------------------------------------------------
! Parallel implementation of Morse pontential.
! - For BVS, see Adams & Rao, Phys. Status Solidi A 208, No.8 (2011)
9 changes: 5 additions & 4 deletions pmd/mod_pairlist.F90
Original file line number Diff line number Diff line change
@@ -171,7 +171,7 @@ subroutine mk_lspr_para(l1st)

call mk_lscl_para()

call set_nnmax()
call set_nnmax(l1st)

if( .not.allocated(lspr) ) then
allocate(lspr(0:nnmax,namax))
@@ -954,7 +954,7 @@ subroutine sort_lspr(namax,natm,ra,nnmax,h,lspr)
return
end subroutine sort_lspr
!=======================================================================
subroutine set_nnmax()
subroutine set_nnmax(l1st)
!
! Estimate and set nnmax from lscl. Thus this must be called after
! lscl and lshd are computed.
@@ -965,13 +965,14 @@ subroutine set_nnmax()
use pmdvars,only: vol,myid_md,mpi_md_world,nnmax,nxyz,rc,rbuf, &
iprint,ratio_nnmax_update,ntot
include "mpif.h"
logical,intent(in):: l1st

integer:: ierr,ic,i,inc,nmaxl,nmax,nnmax_estimate,nnmax_prev, &
ncell,ncelltot
integer:: ix,iy,iz
real(8):: volc,rho
real(8),parameter:: pi = 3.14159265358979d0
logical,save:: l1st = .true.
!!$ logical,save:: l1st = .true.

nmaxl = 0
ncell = 0 ! num. of cells that contain at least one atom
@@ -1026,7 +1027,7 @@ subroutine set_nnmax()
call flush(6)
endif
endif
l1st = .false.
!!$ l1st = .false.
return
end subroutine set_nnmax
!=======================================================================
17 changes: 9 additions & 8 deletions pmd/pmd_core.F90
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
!-----------------------------------------------------------------------
! Last-modified: <2024-04-10 16:48:46 KOBAYASHI Ryo>
! Last-modified: <2024-05-04 22:00:54 KOBAYASHI Ryo>
!-----------------------------------------------------------------------
! Core subroutines/functions needed for pmd.
!-----------------------------------------------------------------------
@@ -164,7 +164,7 @@ subroutine pmd_core(hunit,hmat,ntot0,tagtot,rtot,vtot,atot,stot &
endif
!.....perform space decomposition after reading atomic configuration
tmp = mpi_wtime()
call space_decomp(ntot0,tagtot,rtot,vtot,auxtot)
call space_decomp(ntot0,tagtot,rtot,vtot,auxtot,.true.)
call accum_time('space_decomp',mpi_wtime()-tmp)
!.....Some conversions
do i=1,natm
@@ -1092,7 +1092,7 @@ subroutine oneshot(hunit,hmat,ntot0,tagtot,rtot,vtot,atot,stot, &
write(6,'(a,"[ ",3f12.3," ]")') ' b = ',h(1:3,2,0)
write(6,'(a,"[ ",3f12.3," ]")') ' c = ',h(1:3,3,0)
endif
call space_decomp(ntot0,tagtot,rtot,vtot,auxtot)
call space_decomp(ntot0,tagtot,rtot,vtot,auxtot,.true.)

!.....Some conversions
nsp= 0
@@ -1205,7 +1205,7 @@ subroutine oneshot4fitpot(hunit,hmat,ntot0,tagtot,rtot,vtot,atot,stot, &
write(6,'(a,"[ ",3f12.3," ]")') ' c = ',h(1:3,3,0)
endif
call boxmat(h,hi,ht,g,gi,gt,vol,sgm)
call space_decomp(ntot0,tagtot,rtot,vtot,auxtot)
call space_decomp(ntot0,tagtot,rtot,vtot,auxtot,.true.)

!.....Some conversions
nsp= 0
@@ -1379,7 +1379,7 @@ subroutine min_core(hunit,hmat,ntot0,tagtot,rtot,vtot,atot,stot &

!.....perform space decomposition after reading atomic configuration
tmp = mpi_wtime()
call space_decomp(ntot0,tagtot,rtot,vtot,auxtot)
call space_decomp(ntot0,tagtot,rtot,vtot,auxtot,.true.)
call accum_time('space_decomp',mpi_wtime()-tmp)
!.....Some conversions
do i=1,natm
@@ -2839,7 +2839,7 @@ subroutine vfire(num_fire,alp0_fire,alp_fire,falp_fire,dtmax_fire &
! & ,ierr)
end subroutine vfire
!=======================================================================
subroutine space_decomp(ntot0,tagtot,rtot,vtot,auxtot)
subroutine space_decomp(ntot0,tagtot,rtot,vtot,auxtot,l1st)
!
! Decompose the system and scatter atoms to every process.
!
@@ -2851,6 +2851,7 @@ subroutine space_decomp(ntot0,tagtot,rtot,vtot,auxtot)
real(8),intent(inout):: rtot(3,ntot0),tagtot(ntot0)
real(8),intent(in):: vtot(3,ntot0)
real(8),intent(in):: auxtot(naux,ntot0)
logical,intent(in):: l1st
!!$ real(8),intent(in):: hunit,h(3,3,0:1)
!!$ real(8),intent(inout):: rtot(3,ntot0),tagtot(ntot0)
!!$ integer,intent(in):: ntot0,naux
@@ -2862,7 +2863,7 @@ subroutine space_decomp(ntot0,tagtot,rtot,vtot,auxtot)
integer:: myxt,myyt,myzt,nmin
real(8):: sxogt,syogt,szogt
real(8):: t0
logical,save:: l1st = .true.
!!$ logical,save:: l1st = .true.

t0 = mpi_wtime()

@@ -2914,7 +2915,7 @@ subroutine space_decomp(ntot0,tagtot,rtot,vtot,auxtot)
call estimate_nbmax(nalmax,h,nx,ny,nz,vol,rc,rbuf,nbmax,boundary)
namax = namax +nbmax
if( iprint.ne.0 .and. l1st ) then
l1st = .false.
!!$ l1st = .false.
print *,''
print '(a)', ' space_decomp:'
print '(a,2f6.3)',' rcut, rbuf = ',rc,rbuf

0 comments on commit 04655b6

Please sign in to comment.