LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Problem Calling Certain Python Modules from LAMMPS Input File
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Problem Calling Certain Python Modules from LAMMPS Input File

From: Axel Kohlmeyer <akohlmey@...24...>
Date: Fri, 1 Dec 2017 11:46:14 -0500

On Fri, Dec 1, 2017 at 11:25 AM, Brenda Rubenstein <rubenstein.brenda@...24...> wrote:
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 "", line 13, in <module>
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!

​my first suspicion would be that your PYTHONPATH environment variable is not properly configured.
at any rate, below is an absolute minimum stripped down version of your example.

file in.grid:

variable tx equal 1.0
variable ty equal 1.0
variable tz equal 1.0
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
python force_CH4_x invoke


from scipy.interpolate import griddata

def force_CH4_x(x,y,z):
    print("in force_CH4_x")
    return 0.0

​...and those work fine for me with the current development version of LAMMPS:​

]$ ~/compile/lammps-icms/src/lmp_omp -in in.grid 
LAMMPS (23 Oct 2017)
  using 1 OpenMP thread(s) per MPI task
in force_CH4_x
Total wall time: 0:00:00

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

​using a two year old version of LAMMPS is never a good idea, but it is particularly bad for interfacing to python since a lot of bugfixes and improvements have been added since.​



Relevant Parts of Python Script (name of script is


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

Dr. Brenda Rubenstein
Assistant Professor of Chemistry
Brown University

Check out the vibrant tech community on one of the world's most
engaging tech sites,!
lammps-users mailing list

Dr. Axel Kohlmeyer  akohlmey@...24...
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.