Dear LAMMPS Users,

I have been trying to call a python script from lammps (relevant portions shown below) in order to evaluate forces from a very specific potential energy surface. Even though the python script alone functions properly, when called from lammps, I run into errors of the form:

Traceback (most recent call last):

File “interp_CH4_pot_short_Brenda_modified.py”, line 13, in

from scipy.interpolate import griddata

ImportError: cannot import name griddata

ERROR: Could not process Python file (…/python.cpp:176)

Based upon a number of tests (using mock python scripts and by playing with which scipy functions I am calling), it seems as though lammps has trouble whenever I try to use the scipy.interpolate package. It does not want to load the required griddata (or other interpolate) modules. However, it does not incur any problems running non scipy interpolate functions. Is this a systematic problem loading certain modules? Are there any ideas as to why this won’t work for me? Any help or advice would be greatly appreciated!

Details:

I am using the LAMMPS Dec 7, 2015 version. Maybe later versions can do this?

Relevant Parts of Python Script (name of script is interp_CH4_pot_short_Brenda_modified.py):

#!/usr/bin/python

import sys

import os

import time

import math

import numpy as np

import shutil

import scipy

import subprocess

from scipy.interpolate import griddata

potgrad = np.load(‘potgrad.npy’)

gridsize = .1

gridx = np.ceil(19.2/gridsize) + 1

gridy = np.ceil(3.056/gridsize) + 1

gridz = np.ceil(5.15/gridsize) + 1

gridx = int(gridx)

gridy = int(gridy)

gridz = int(gridz)

gridcoords = np.zeros((gridx,gridy,gridz,3))

gridpot = np.zeros((gridx,gridy,gridz))

gridforce = np.zeros((gridx,gridy,gridz,3))

xc = np.linspace(-1.1,18.1,gridx) #.1 A extra on each side

yc = np.linspace(-.1,2.956231451,gridy) #.1 A extra on each side

zc = np.linspace(-2.573458221,2.573458221,gridz) #.1 A extra on each side

spacex = 19.2/(gridx-1)

spacey = 3.056231451/(gridy-1)

spacez = 5.146916442/(gridz-1)

gridcoords = np.array(gridcoords)

forcey = scipy.interpolate.RegularGridInterpolator((xc,yc,zc), -potgrad[1,:,:,:])

def force_CH4_x(x,y,z):

x_p = np.mod(x,2.856231451)

y_p = np.mod(y,4.946916442) - 2.47345822

force_x = forcey((z,x_p,y_p))

return force_x

force_2 = force_CH4_x(1,1,1)

print force_2

Relevant Call in LAMMPS Input File:

variable tz equal z[1]

variable tx equal x[1]

variable ty equal y[1]

variable force_z2 python force_CH4_x

python force_CH4_x input 3 v_tx v_ty v_tz return v_force_z2 format ffff file interp_CH4_pot_short_Brenda_modified.py

Best,

Dr. Brenda Rubenstein

Assistant Professor of Chemistry

Brown University

https://www.brown.edu/research/labs/rubenstein/home