Addressed inability of python wrapper to handle mulitple chimesFF instances

serial_interface/api/chimescalc_serial_py.py

@@ -6,24 +6,126 @@
 
 		import wrapper_py
 		wrapper_py.chimes_wrapper = wrapper_py.init_chimes_wrapper("/path/to/libchimescalc_dl.so")
-		wrapper_py.set_chimes()
-		wrapper_py.init_chimes()
-		wrapper_py.read_params("some_parameter_file.txt")
+		chimes_ptr = wrapper_py.chimes_open_instance()
+                        where chimes_ptr is a void pointer to the ChimesFF instance; 
+                        this is needed to close and create new instances
+		wrapper_py.set_chimes_instance()
+		wrapper_py.init_chimes_instance()
+
+        To free up the memory for a chimesFF instance call
+
+                wrapper_py.chimes_close_instance(chimes_ptr)
 
     ChIMES Calculator
     Copyright (C) 2020 Rebecca K. Lindsey, Nir Goldman, and Laurence E. Fried
 	Contributing Author: Rebecca K. Lindsey (2020)
 
+
 """
 
 import ctypes
 
-
 chimes_wrapper = None
 
 def init_chimes_wrapper(lib_name):
 	return ctypes.CDLL(lib_name)
 
+def chimes_open_instance():
+	""" Allocates memory for a the chimesFF object and returns pointer to that
+        object"""
+	chimes_wrapper.chimes_open_instance.restype = ctypes.POINTER(ctypes.c_void_p)
+	chimes_ptr = chimes_wrapper.chimes_open_instance()
+	return chimes_ptr
+
+def chimes_close_instance(chimes_ptr):
+	""" Frees memory for the chimesFF object """
+	chimes_wrapper.chimes_close_instance(chimes_ptr)
+
+
+def set_chimes_instance(chimes_ptr, small=False):
+	""" Instantiates the chimesFF object """
+
+	chimes_wrapper.set_chimes_serial_instance(chimes_ptr, small)
+	return
+
+def init_chimes_instance(chimes_ptr, param_file, rank):
+	""" 
+	Initializes the chimesFF object (sets MPI rank)
+	Optionally takes an  MPI rank as input
+	"""
+	in_paramfile = ctypes.c_char_p(param_file.encode())
+	in_rank      = ctypes.c_int(rank)
+	chimes_wrapper.init_chimes_serial_instance(chimes_ptr, in_paramfile, ctypes.byref(in_rank))
+	return
+
+def calculate_chimes_instance(chimes_ptr, natoms,xcrd,ycrd,zcrd,atmtyps,cell_a,cell_b,cell_c,energy,fx,fy,fz,stress):
+	""" 
+	Computes the ChIMES forces, energy, and stress tensor for a given system
+
+	Inputs:	
+	natoms:	Number of atoms in system 
+	xcrd:	System x-coordinates
+	ycrd:	System y-coordinates
+	zcrd:	System z-coordinates
+	atmtyps:System atom types
+	cell_a:	System a lattice vector
+	cell_b:	System b lattice vector
+	cell_c:	System c lattice vector
+	energy:	System energy
+	fx:	X force components for system atoms
+	fy:	Y force components for system atoms
+	fz:	Z force components for system atoms
+	stress:	System stress tensor
+
+	Returns updated fx, fy, fz, stress, and energy
+
+	"""
+
+	encoded = []
+
+	for i in range(natoms):
+		encoded.append(atmtyps[i].encode())
+
+	in_natom   = ctypes.c_int(natoms) 
+	in_xcrd    = (ctypes.c_double * natoms) (*xcrd)
+	in_ycrd    = (ctypes.c_double * natoms) (*ycrd)
+	in_zcrd    = (ctypes.c_double * natoms) (*zcrd)
+	in_atmtyps = (ctypes.c_char_p * natoms) (*encoded)
+	in_cell_a  = (ctypes.c_double * 3) (*cell_a)
+	in_cell_b  = (ctypes.c_double * 3) (*cell_b)
+	in_cell_c  = (ctypes.c_double * 3) (*cell_c)
+	in_energy  = ctypes.c_double(energy)
+	in_fx      = (ctypes.c_double * natoms) (*fx)
+	in_fy      = (ctypes.c_double * natoms) (*fy)
+	in_fz      = (ctypes.c_double * natoms) (*fz)
+	in_stress  = (ctypes.c_double * 9) (*stress)
+
+	chimes_wrapper.calculate_chimes_instance(chimes_ptr,
+					in_natom,   
+					in_xcrd,   
+					in_ycrd,   
+					in_zcrd,   
+					in_atmtyps,   
+					in_cell_a,
+					in_cell_b, 
+					in_cell_c, 
+					ctypes.byref (in_energy), 
+					in_fx,     
+					in_fy,     
+					in_fz,     
+					in_stress)
+
+	return in_fx, in_fy, in_fz, in_stress, in_energy.value
+
+"""
+    Below are depricated functions for the chimesFF wrapper kept for
+    backwards compatability.
+
+    They do not allow for multiple instances of chimesFF to be intialized
+    for a single system call.
+
+"""
+
 def set_chimes(small=False):
 	""" Instantiates the chimesFF object """
 	chimes_wrapper.set_chimes_serial(small)
@@ -99,3 +201,4 @@ def calculate_chimes(natoms,xcrd,ycrd,zcrd,atmtyps,cell_a,cell_b,cell_c,energy,f
 	return in_fx, in_fy, in_fz, in_stress, in_energy.value
 
 
+

serial_interface/examples/python_instance/Makefile

@@ -0,0 +1,49 @@
+C_LOC = $(realpath .)
+
+CXX     = g++ -O3 -std=c++11
+
+# Rudimentary OS detection
+UNAME := $(shell uname)
+ifeq (${UNAME},Linux)
+	CXX += -fPIC
+endif
+
+CHIMESFF_LOC=$(C_LOC)/../../../chimesFF/src
+CHIMESFF_SRC=$(CHIMESFF_LOC)/chimesFF.cpp
+CHIMESFF_HDR=$(CHIMESFF_LOC)/chimesFF.h
+
+chimesFF.o : $(CHIMESFF_SRC)
+	$(CXX) -c $(CHIMESFF_SRC) -I $(CHIMESFF_LOC)
+
+SERIAL_LOC=$(C_LOC)/../../src
+SERIAL_SRC=$(SERIAL_LOC)/serial_chimes_interface.cpp
+SERIAL_HDR=$(SERIAL_LOC)/serial_chimes_interface.h
+
+serial_chimes_interface.o : $(SERIAL_SRC) 
+	$(CXX) -c $(SERIAL_SRC) -I $(SERIAL_LOC) -I $(CHIMESFF_LOC) 
+
+WRAPPER_LOC=$(C_LOC)/../../api
+WRAPPER_SRC=$(WRAPPER_LOC)/chimescalc_serial_C.cpp
+WRAPPER_HDR=$(WRAPPER_LOC)/chimescalc_serial_C.h
+
+chimescalc_serial_C.o : $(WRAPPER_SRC)
+	$(CXX) -c $(WRAPPER_SRC) -I $(WRAPPER_LOC) -I $(SERIAL_LOC) -I $(CHIMESFF_LOC)	
+
+
+chimescalc_serial_C.so:
+	$(CXX) -shared -o libchimescalc-serial_dl.so chimescalc_serial_C.o serial_chimes_interface.o chimesFF.o
+
+clean:
+	rm -f *.o
+
+clean-all:
+
+	make clean
+	rm -f *.so
+
+all:
+	make chimesFF.o
+	make serial_chimes_interface.o
+	make chimescalc_serial_C.o
+	make chimescalc_serial_C.so
+	make clean

serial_interface/examples/python_instance/README

@@ -0,0 +1,23 @@
+This code demonstrates how the serial_chimes_interface can be used to obtain 
+stress tensor, energy, and per-atom forces for a given system using python. 
+See main.py for a usage example.
+
+Note: This script takes as input a standard ChIMES parameter file,
+and a .xyz file with a, b, and c cell vectors  in the comment line.
+This implementation uses ghost atoms/layering thus a number of  layers
+for the system must be specified such that the effective box lengths
+aregreater than two times the largest outer cutoff, or results will not be 
+correct.
+
+Note: This code requires libwrapper-C.so in the same directory, which can be 
+generated via "make all"
+
+Note: Expects to be run with python version 3.X
+
+Compile with:
+	make all
+To test the executable:
+
+python3 main.py <parameter file> <coordinate file> <nlayers>
+
+e.g. python3 main.py  ../../tests/force_fields/published_params.liqC.2b.cubic.txt ../../tests/configurations/liqC.2.5gcc_6000K.OUTCAR_#000.xyzf 2

serial_interface/examples/python_instance/main.py

@@ -0,0 +1,171 @@
+"""
+
+	A simple usage example for computing system force, stress, and energy via 
+	the serial_chimes_calculator python wrapper.
+
+	Expects "libchimescalc_dl.so" in the same directory as this script. This file
+	can be generated by typing "make all" in the directory containing this file
+
+	Depends on chimescalc_serial_py.py from the chimes_serial_interface's api files
+
+	Expects to be run with python version 3.X
+
+	Run with: "python3 <this file> <parameter file> <coordinate file>"  or 
+               python3 Run with: python <this file> <parameter file> <xyz file> <allow_replicates(0/1)> <debug flag (0/1)> <path to wrapper_py.py>
+
+    ChIMES Calculator
+    Copyright (C) 2020 Rebecca K. Lindsey, Nir Goldman, and Laurence E. Fried
+	Contributing Author: Rebecca K. Lindsey (2020)
+
+"""
+
+import os
+import sys
+import math
+
+if (len(sys.argv) != 4) and (len(sys.argv) != 6):
+
+	print( "ERROR: Wrong number of commandline args")
+	print( "       Run with: python <this file> <parameter file> <xyz file> <allow_replicates(0/1)>")
+	print( "       or")	
+	print( "       Run with: python3 <this file> <parameter file> <xyz file> <allow_replicates(0/1)> <debug flag (0/1)> <path to chimescalc_serial_py.py>")
+	exit()
+
+# A small helper function
+
+def str2bool(v):
+  return v.lower() in ("yes", "true", "t", "1")
+
+# Import ChIMES modules
+
+curr_path = os.getcwd()
+
+chimes_module_path = os.path.abspath( curr_path + "/../../api/")
+if len(sys.argv) == 6:
+	chimes_module_path = os.path.abspath(sys.argv[4])
+sys.path.append(chimes_module_path)
+
+import chimescalc_serial_py 
+
+# Read in the parameter and coordinate filename 
+
+small = False
+
+
+
+param_file =          sys.argv[1]  # parameter file
+coord_file =          sys.argv[2]  # coordinate file
+small      = str2bool(sys.argv[3]) # allow replicates?
+
+print("Read args:")
+
+for i in range(len(sys.argv)-1):
+	print (i+1,sys.argv[i+1])
+
+# Initialize the ChIMES calculator curr_path should be /.../usr/WS2/rlindsey/chimes_calculator-fork/serial_interface/tests/
+
+print("CURR PATH:",curr_path)
+
+chimescalc_serial_py.chimes_wrapper = chimescalc_serial_py.init_chimes_wrapper(curr_path + "/libchimescalc_dl.so")
+chimes_ptr = chimescalc_serial_py.chimes_open_instance()
+chimescalc_serial_py.set_chimes_instance(chimes_ptr, small)
+
+rank = 0
+
+chimescalc_serial_py.init_chimes_instance(chimes_ptr, param_file, rank)
+
+# Read the coordinates, set up the force, stress, and energy vars
+
+natoms  = None
+lx      = None
+ly      = None
+lz      = None
+atmtyps = []
+xcrd    = []
+ycrd    = []
+zcrd    = []
+
+ifstream = open(coord_file,'r')
+natoms   = int(ifstream.readline())
+cell_a       = ifstream.readline().split()
+cell_b       = [float(cell_a[3]), float(cell_a[4]), float(cell_a[5])]
+cell_c       = [float(cell_a[6]), float(cell_a[7]), float(cell_a[8])]
+cell_a       = [float(cell_a[0]), float(cell_a[1]), float(cell_a[2])]
+
+energy = 0.0
+stress = [0.0]*9
+fx     = [] 
+fy     = []
+fz     = []
+
+for i in range(natoms):
+
+	atmtyps.append(ifstream.readline().split())
+
+	xcrd.append(float(atmtyps[-1][1]))
+	ycrd.append(float(atmtyps[-1][2]))
+	zcrd.append(float(atmtyps[-1][3]))
+
+	atmtyps[-1] = atmtyps[-1][0]
+
+	fx.append(0.0)
+	fy.append(0.0)
+	fz.append(0.0)
+
+# Do the calculations
+
+fx, fy, fz, stress, energy = chimescalc_serial_py.calculate_chimes_instance(
+                           chimes_ptr, 
+                           natoms, 
+			   xcrd, 
+			   ycrd, 
+			   zcrd, 
+			   atmtyps, 
+			   cell_a,
+			   cell_b,
+			   cell_c,
+			   energy,
+			   fx,
+			   fy,
+			   fz,
+			   stress)
+
+# close calculator instance
+chimescalc_serial_py.chimes_close_instance(chimes_ptr)
+
+print ("Success!")
+print ("Energy (kcal/mol)",energy)
+print("Stress tensors (GPa): ")
+print(stress[0]*6.9479)
+print(stress[4]*6.9479)
+print(stress[8]*6.9479)
+print(stress[1]*6.9479)
+print(stress[2]*6.9479)
+print(stress[5]*6.9479)
+print("Forces (kcal/mol/A): ")
+for i in range(natoms):
+	print(fx[i])
+	print(fy[i])
+	print(fz[i])
+
+debug = 0
+
+
+if len(sys.argv) == 6:
+	debug = int(sys.argv[-1])
+
+if debug == 1:
+
+	ofstream = open("debug.dat",'w')
+
+	ofstream.write("{0:0.6f}\n".format(energy))
+	ofstream.write("{0:0.6f}\n".format(stress[0]*6.9479))
+	ofstream.write("{0:0.6f}\n".format(stress[4]*6.9479))
+	ofstream.write("{0:0.6f}\n".format(stress[8]*6.9479))
+	ofstream.write("{0:0.6f}\n".format(stress[1]*6.9479))
+	ofstream.write("{0:0.6f}\n".format(stress[2]*6.9479))
+	ofstream.write("{0:0.6f}\n".format(stress[5]*6.9479))
+	for i in range(natoms):
+		ofstream.write("{0:0.6e}\n{1:0.6e}\n{2:0.6e}\n".format(fx[i],fy[i],fz[i]))
+
+	ofstream.close()