diff --git a/external_codes/mpi_wrapper/mpi3/communicator.hpp b/external_codes/mpi_wrapper/mpi3/communicator.hpp
index c5b009b21..df1bc73b7 100644
--- a/external_codes/mpi_wrapper/mpi3/communicator.hpp
+++ b/external_codes/mpi_wrapper/mpi3/communicator.hpp
@@ -17,6 +17,7 @@
 #include "../mpi3/status.hpp"
 #include "../mpi3/type.hpp"
 #include "../mpi3/error.hpp"
+#include "../mpi3/group.hpp"
 
 #include "../mpi3/detail/basic_communicator.hpp"
 #include "../mpi3/detail/buffer.hpp"
@@ -146,8 +147,6 @@ class group;
 
 using boost::optional;
 
-template<class T = void>
-struct window;
 
 struct error_handler;
 struct keyval;
@@ -185,6 +184,8 @@ struct shared_communicator; // intracommunicator
 using boost::any;
 using boost::any_cast;
 
+template<class T> struct window;
+
 class communicator : protected detail::basic_communicator{
 	friend struct detail::package;
 protected:
diff --git a/external_codes/mpi_wrapper/mpi3/dynamic_window.hpp b/external_codes/mpi_wrapper/mpi3/dynamic_window.hpp
index 43ddb184f..53d61b9dd 100644
--- a/external_codes/mpi_wrapper/mpi3/dynamic_window.hpp
+++ b/external_codes/mpi_wrapper/mpi3/dynamic_window.hpp
@@ -16,7 +16,7 @@ struct dynamic_window : window<T>{
 	dynamic_window() : window<T>{}{}
 	public:
 	dynamic_window(communicator& comm){
-		int s = MPI_Win_create_dynamic(MPI_INFO_NULL, comm.impl_, &(this->impl_));
+		int s = MPI_Win_create_dynamic(MPI_INFO_NULL, comm.get(), &(this->impl_));
 		if(s != MPI_SUCCESS) throw std::runtime_error("cannot create dynamic window");
 	}
 	template<class TT = char>
diff --git a/external_codes/mpi_wrapper/mpi3/shm/allocator.hpp b/external_codes/mpi_wrapper/mpi3/shm/allocator.hpp
index 040b1b8ec..a5d50fe1e 100644
--- a/external_codes/mpi_wrapper/mpi3/shm/allocator.hpp
+++ b/external_codes/mpi_wrapper/mpi3/shm/allocator.hpp
@@ -4,22 +4,57 @@
 #ifndef MPI3_SHM_ALLOCATOR_HPP
 #define MPI3_SHM_ALLOCATOR_HPP
 
-#include "../../mpi3/shared_window.hpp"
+//#include "../../mpi3/shared_window.hpp"
+#include "../shm/pointer.hpp"
 
 namespace boost{
 namespace mpi3{
 namespace shm{
 
-template<class... As>
-using allocator = mpi3::intranode::allocator<As...>;
+//template<class... As>
+//using allocator = mpi3::intranode::allocator<As...>;
 
-template<class T>
-using pointer = mpi3::intranode::array_ptr<T>;
+//template<class T>
+//using pointer = mpi3::intranode::array_ptr<T>;
 
-template<class Ptr>
-struct pointer_traits : std::pointer_traits<Ptr>{
-	static auto to_address(Ptr const& p){
-		return std::addressof(*p);
+template<class T = void>
+struct allocator{
+	template<class U> struct rebind{typedef allocator<U> other;};
+	using value_type = T;
+	using pointer = shm::pointer<T>;
+	using const_pointer = shm::pointer<T const>;
+	using size_type = mpi3::size_t;
+	using difference_type = mpi3::ptrdiff_t;
+private:
+	mpi3::shared_communicator* comm_;
+	allocator() = delete;
+public:
+	allocator(mpi3::shared_communicator* comm) : comm_{comm}{}
+	template<class U> allocator(allocator<U> const& o) : comm_{o.comm_}{}
+	~allocator() = default;
+
+	pointer allocate(size_type n, const void* /*hint*/ = 0) const{
+		pointer ret = 0;
+	// this code using reset is to workaround a bug in outdated libc++ https://bugs.llvm.org/show_bug.cgi?id=42021
+		ret.w_.reset(new shared_window<T>{
+			comm_->make_shared_window<T>(comm_->root()?n:0)
+		});
+	//	ret.w_ = std::make_shared<shared_window<T>>(
+	//		comm_->make_shared_window<T>(comm_->root()?n:0)
+	//	);
+		return ret;
+	}
+	void deallocate(pointer p, size_type){p.w_.reset();}
+	allocator& operator=(allocator const& other) = default;
+	bool operator==(allocator const& other) const{return comm_==other.comm_;}
+	bool operator!=(allocator const& other) const{return not(other == *this);}
+	template<class P = pointer, class... As>
+	void construct(P p, As&&... as){
+		if(comm_->root()) ::new((void*)p) typename std::pointer_traits<P>::element_type(std::forward<As>(as)...);
+	}
+	template<class P = pointer> void destroy(P p){
+		using V = typename std::pointer_traits<P>::element_type;
+		p->~V();
 	}
 };
 
@@ -33,6 +68,7 @@ struct pointer_traits : std::pointer_traits<Ptr>{
 
 namespace mpi3 = boost::mpi3;
 
+#if 0
 template<class A>
 struct monotonic{
 	A& a_;
@@ -57,22 +93,13 @@ struct monotonic{
 	}
 	void deallocate(monotonic::pointer, monotonic::size_type){}
 };
+#endif
 
 int mpi3::main(int, char*[], mpi3::communicator world){
 
-	std::allocator<double> stda;
-	std::allocator<double>::rebind<int>::other dd;
-	monotonic<std::allocator<double>> a(1024, stda);
-	std::vector<double, std::allocator<double>> v1(600, stda); 
-
-	std::vector<double, monotonic<std::allocator<double> > > v2(600, a); 
-
-#if 1
-//	std::vector<double, monotonic<std::allocator<double>>> v2(600, a); 
-
 	mpi3::shared_communicator node = world.split_shared();
 
-	mpi3::shm::allocator<double> A1(node);
+	mpi3::shm::allocator<double> A1(std::addressof(node));
 	mpi3::shm::pointer<double> data1 = A1.allocate(80);
 
 	using ptr = decltype(data1);
@@ -90,7 +117,7 @@ int mpi3::main(int, char*[], mpi3::communicator world){
 	assert(data1[3] == 3.4);
 
 	A1.deallocate(data1, 80);
-#endif
+
 	return 0;
 }
 
diff --git a/external_codes/mpi_wrapper/mpi3/shm/memory.hpp b/external_codes/mpi_wrapper/mpi3/shm/memory.hpp
index 60514e4c4..a3b02226e 100644
--- a/external_codes/mpi_wrapper/mpi3/shm/memory.hpp
+++ b/external_codes/mpi_wrapper/mpi3/shm/memory.hpp
@@ -81,6 +81,23 @@ struct mapped_region{
 
 }}}
 
+
+namespace boost{
+namespace mpi3{
+namespace shm{
+
+template<class Alloc, typename T, typename Size, typename TT>
+	pointer<T> uninitialized_fill_n(Alloc const& a, pointer<T> f, Size n, TT const& val){
+	assert(0);
+	using std::uninitialized_fill_n;
+	if(mpi3::group(*f.wSP_).root()) uninitialized_fill_n(to_address(f), n, val);
+	first.wSP_->fence();
+	first.wSP_->fence();
+	return first + n;
+}
+
+}}}
+
 #ifdef _TEST_BOOST_MPI3_SHM_MEMORY
 
 #include "../../mpi3/main.hpp"
diff --git a/external_codes/mpi_wrapper/mpi3/shm/mutex.hpp b/external_codes/mpi_wrapper/mpi3/shm/mutex.hpp
index ecb8cd7a5..4917cf627 100644
--- a/external_codes/mpi_wrapper/mpi3/shm/mutex.hpp
+++ b/external_codes/mpi_wrapper/mpi3/shm/mutex.hpp
@@ -15,14 +15,14 @@ class mutex{
 	mpi3::shm::allocator<std::atomic_flag> alloc_;//(node);
 	mpi3::shm::pointer<std::atomic_flag> f_;
 	public:
-	mutex(mpi3::shared_communicator& scomm) : scomm_(scomm), alloc_(scomm_), f_(alloc_.allocate(1)){
+	mutex(mpi3::shared_communicator& scomm) : scomm_(scomm), alloc_(std::addressof(scomm_)), f_(alloc_.allocate(1)){
 		if(scomm_.root()) alloc_.construct(&*f_, false);
 		scomm_.barrier();
 	}
 	mutex(mutex const&) = delete;
 	mutex(mutex&&) = delete;
-	void lock(){while(f_->test_and_set(std::memory_order_acquire));}
-	void unlock(){f_->clear(std::memory_order_release);}
+	void lock(){while((*f_).test_and_set(std::memory_order_acquire));}
+	void unlock(){(*f_).clear(std::memory_order_release);}
 	~mutex(){
 		if(scomm_.root()) alloc_.destroy(&*f_);
 		scomm_.barrier();
diff --git a/external_codes/mpi_wrapper/mpi3/window.hpp b/external_codes/mpi_wrapper/mpi3/window.hpp
index be1a9171c..a2c62057e 100644
--- a/external_codes/mpi_wrapper/mpi3/window.hpp
+++ b/external_codes/mpi_wrapper/mpi3/window.hpp
@@ -1,12 +1,15 @@
 #if COMPILATION_INSTRUCTIONS
-(echo "#include\""$0"\"" > $0x.cpp) && mpic++ -O3 -std=c++14 -Wfatal-errors -D_TEST_BOOST_MPI3_WINDOW $0x.cpp -o $0x.x && time mpirun -np 4 $0x.x $@ && rm -f $0x.x $0x.cpp; exit
+(echo "#include\""$0"\"" > $0x.cpp) && mpic++ -O3 -std=c++17 `#-Wfatal-errors` -Wall -Wextra -Wpedantic -D_TEST_BOOST_MPI3_WINDOW $0x.cpp -o $0x.x && time mpirun -np 4 $0x.x $@ && rm -f $0x.x $0x.cpp; exit
 #endif
-//  (C) Copyright Alfredo A. Correa 2018.
+//  (C) Copyright Alfredo A. Correa 2019
 #ifndef BOOST_MPI3_WINDOW_HPP
 #define BOOST_MPI3_WINDOW_HPP
 
-#include "../mpi3/communicator.hpp"
+#include "../mpi3/error.hpp"
+#include "../mpi3/type.hpp"
 #include "../mpi3/detail/datatype.hpp"
+#include "../mpi3/detail/call.hpp"
+#include "../mpi3/communicator.hpp"
 
 #define OMPI_SKIP_MPICXX 1  // https://github.com/open-mpi/ompi/issues/5157
 #include<mpi.h>
@@ -14,51 +17,42 @@
 namespace boost{
 namespace mpi3{
 
-class group;
+struct target{
+	int rank;
+	mpi3::size_t disp;
+};
 
-template<class T = void>
-class panel;
+template<class T = void> class panel;
 
-//template<class T = void> struct window;
+template<class T = void> struct window;
 
-struct basic_window{
+template<>
+struct window<void>{
 protected:
 	MPI_Win impl_ = MPI_WIN_NULL;
-	basic_window() = default;
-	basic_window(MPI_Win w) : impl_{w}{}
 public:
 	MPI_Win& operator&(){return impl_;}
 	MPI_Win const& operator&() const{return impl_;}
-};
-
-template<>
-struct window<void> : basic_window{
-	public:
-	window(){}
-	template<class T>
-	window(T* base, mpi3::size_t size, communicator& comm){
-		int s = MPI_Win_create(
-			(void*)base, size*sizeof(T), sizeof(T), MPI_INFO_NULL, 
-			&comm, &impl_
-		); // TODO probably need alignof
-		if(s != MPI_SUCCESS) throw std::runtime_error{"cannot create window"};
+	void clear(){
+		if(impl_ != MPI_WIN_NULL) MPI3_CALL(MPI_Win_free)(&impl_);
+		assert(impl_ == MPI_WIN_NULL);
 	}
-	window(void* base, mpi3::size_t size, communicator& comm){
-		int s = MPI_Win_create(base, size, 1, MPI_INFO_NULL, &comm, &impl_);
-		if(s != MPI_SUCCESS) throw std::runtime_error{"cannot create window"};
+protected:
+	window() = default;
+public:
+	template<class T, class Size = mpi3::size_t>
+	window(communicator const& c, T* b, Size n = 0){
+		MPI3_CALL(MPI_Win_create)(b, n*sizeof(T), alignof(T), MPI_INFO_NULL, c.get(), &impl_);
+		assert( alignof(T) == sizeof(T) ); // to see in what type it is different
 	}
-	window(communicator& comm) : window((void*)nullptr, 0, comm){}
-	window(window const&) = delete; // cannot be duplicated, see text before sec. 4.5 in Using Adv. MPI
-	window(window&& other) : basic_window{std::exchange(other.impl_, MPI_WIN_NULL)}{//is movable if null is not a correct state?
-//		other.impl_ = MPI_WIN_NULL;
+	window(window const&) = delete;// see text before §4.5 in Using Adv. MPI
+	window(window&& o) noexcept : impl_{std::exchange(o.impl_, MPI_WIN_NULL)}{}
+	window& operator=(window const&) = delete; // see cctor
+	window& operator=(window&& other){// self assignment is undefined
+		clear(); swap(*this, other); return *this;
 	}
-	window& operator=(window const&) = delete; // see comment in cctor
-	window& operator=(window&& other){ // guard self assignment? probably worthless
-		if(impl_ != MPI_WIN_NULL) MPI_Win_free(&impl_);
-		impl_ = std::exchange(other.impl_, MPI_WIN_NULL);
-		return *this;
-	}
-	~window(){if(impl_ != MPI_WIN_NULL) MPI_Win_free(&impl_);}
+	friend void swap(window& a, window& b){std::swap(a.impl_, b.impl_);}
+	~window(){clear();}
 	template<typename It1, typename Size, typename V = typename std::iterator_traits<It1>::value_type>
 	void accumulate_n(It1 first, Size count, int target_rank, int target_disp = 0){
 		using detail::data;
@@ -84,37 +78,40 @@ struct window<void> : basic_window{
 	void flush_local(){return flush_local_all();}
 	void* base() const{
 		void* base; int flag;
-		int s = MPI_Win_get_attr(impl_, MPI_WIN_BASE, &base, &flag);
-		if(s != MPI_SUCCESS) throw std::runtime_error("cannot get base");
+		MPI3_CALL(MPI_Win_get_attr)(impl_, MPI_WIN_BASE, &base, &flag);
 		assert(flag);
 		return base;
 	}
 	mpi3::size_t const& size() const{
 		MPI_Aint* size_p; int flag;
-		int s = MPI_Win_get_attr(impl_, MPI_WIN_SIZE, &size_p, &flag);
-		if(s != MPI_SUCCESS) throw std::runtime_error("cannot get base");
+		MPI3_CALL(MPI_Win_get_attr)(impl_, MPI_WIN_SIZE, &size_p, &flag);
 		assert(flag);
 		return *size_p;
 	}
 	int const& disp_unit() const{
 		int* disp_unit_p; int flag;
-		int s = MPI_Win_get_attr(impl_, MPI_WIN_DISP_UNIT, &disp_unit_p, &flag);
-		if(s != MPI_SUCCESS) throw std::runtime_error("cannot get base");
+		MPI3_CALL(MPI_Win_get_attr)(impl_, MPI_WIN_DISP_UNIT, &disp_unit_p, &flag);
 		assert(flag);
 		return *disp_unit_p;
 	}
 //	get_errhandler(...);
+	group get_group() const{
+		group ret;
+		MPI3_CALL(MPI_Win_get_group)(impl_, &(&ret));
+		return ret;
+	}
 //	group get_group(){use reinterpret_cast?}
 //	... get_info
 //	... get_name
+	// lock arguments are reversed
 	void lock(int rank, int lock_type = MPI_LOCK_EXCLUSIVE, int assert = MPI_MODE_NOCHECK){
-		MPI_Win_lock(lock_type, rank, assert, impl_);
+		MPI3_CALL(MPI_Win_lock)(lock_type, rank, assert, impl_);
 	}
 	void lock_exclusive(int rank, int assert = MPI_MODE_NOCHECK){
-		MPI_Win_lock(MPI_LOCK_EXCLUSIVE, rank, assert, impl_);
+		MPI3_CALL(MPI_Win_lock)(MPI_LOCK_EXCLUSIVE, rank, assert, impl_);
 	}
 	void lock_shared(int rank, int assert = MPI_MODE_NOCHECK){
-		MPI_Win_lock(MPI_LOCK_SHARED, rank, assert, impl_);
+		MPI3_CALL(MPI_Win_lock)(MPI_LOCK_SHARED, rank, assert, impl_);
 	}
 	void lock_all(int assert = MPI_MODE_NOCHECK){MPI_Win_lock_all(assert, impl_);}
 //	void post(group const& g, int assert = MPI_MODE_NOCHECK) const{
@@ -133,39 +130,33 @@ struct window<void> : basic_window{
 	void unlock(int rank) const{MPI_Win_unlock(rank, impl_);}
 	void unlock_all(){MPI_Win_unlock_all(impl_);}
 	void wait() const{MPI_Win_wait(impl_);}
-
 //	void fetch_and_op(T const*  origin, T* target, int target_rank, int target_disp = 0) const{
 //		MPI_Fetch_and_op(origin, target, detail::datatype<T>{}, target_rank, target_disp, , impl_);
 //	}
-
 //	template<class T, class Op, class datatype = detail::datatype<T>, >
 //	void fetch_and_op(T const*  origin, T* target, int target_rank, int target_disp = 0) const{
 //		MPI_Fetch_and_op(origin, target, datatype{}, target_rank, target_disp, , impl_);
 //	}
-
 //	void fetch_exchange(T const*  origin, T* target, int target_rank, int target_disp = 0) const{
 //		MPI_Fetch_and_op(origin, target,detail::datatype<T>{}, target_rank, target_disp, MPI_REPLACE, impl_);
 //	}
-
 //	maybe this goes to a pointer impl
-
 	template<class T>
-	void fetch_sum_value(T const& origin, T& target, int target_rank, int target_disp = 0) const{
-		MPI_Fetch_and_op(&origin, &target, detail::basic_datatype<T>{}, target_rank, target_disp, MPI_SUM, impl_);
+	void fetch_sum_value(T const& origin, T& target, int target_rank, int target_disp=0) const{
+		MPI3_CALL(MPI_Fetch_and_op)(&origin, &target, detail::basic_datatype<T>{}, target_rank, target_disp, MPI_SUM, impl_);
 	}
 	template<class T>
 	void fetch_prod_value(T const& origin, T& target, int target_rank, int target_disp = 0) const{
-		MPI_Fetch_and_op(&origin, &target, detail::basic_datatype<T>{}, target_rank, target_disp, MPI_PROD, impl_);
+		MPI3_CALL(MPI_Fetch_and_op)(&origin, &target, detail::basic_datatype<T>{}, target_rank, target_disp, MPI_PROD, impl_);
 	}
 	template<class T>
 	void fetch_replace_value(T const&  origin, T& target, int target_rank, int target_disp = 0) const{
-		MPI_Fetch_and_op(&origin, &target, detail::basic_datatype<T>{}, target_rank, target_disp, MPI_REPLACE, impl_);
+		MPI3_CALL(MPI_Fetch_and_op)(&origin, &target, detail::basic_datatype<T>{}, target_rank, target_disp, MPI_REPLACE, impl_);
 	}
 	template<class CI1, class CI2, class datatype = detail::basic_datatype<typename std::iterator_traits<CI1>::value_type> >
 	void fetch_replace(CI1 it1, CI2 it2, int target_rank, int target_disp = 0) const{
-		MPI_Fetch_and_op(std::addressof(*it1), std::addressof(*it2), datatype{}, target_rank, target_disp, MPI_REPLACE, impl_); 
+		MPI3_CALL(MPI_Fetch_and_op)(std::addressof(*it1), std::addressof(*it2), datatype{}, target_rank, target_disp, MPI_REPLACE, impl_); 
 	}
-
 	template<class ContiguousIterator>
 	void blocking_put_n(ContiguousIterator it, int count, int target_rank, int target_offset = 0){
 		lock_shared(target_rank, 0);
@@ -175,7 +166,7 @@ struct window<void> : basic_window{
 	template<class ContiguousIterator>
 	void put_n(ContiguousIterator it, std::size_t n, int target_rank, int target_disp = 0) const{
 		using detail::data;
-		MPI_Put(
+		MPI3_CALL(MPI_Put)(
 			data(it), /* void* origin_address = a + i*/ 
 			n, /*int origin_count = 1 */
 			detail::basic_datatype<typename std::iterator_traits<ContiguousIterator>::value_type>{}, 
@@ -198,7 +189,7 @@ struct window<void> : basic_window{
 	template<typename ContiguousIterator, typename Size>
 	ContiguousIterator get_n(ContiguousIterator it, Size n, int target_rank, int target_disp = 0) const{
 		using detail::data;
-		int s = MPI_Get(
+		MPI3_CALL(MPI_Get)(
 			data(it), /* void* origin_address = b + i*/
 			n, /*int origin_count = 1 */
 			detail::basic_datatype<typename std::iterator_traits<ContiguousIterator>::value_type>{}, 
@@ -208,7 +199,6 @@ struct window<void> : basic_window{
 			detail::basic_datatype<typename std::iterator_traits<ContiguousIterator>::value_type>{}, 
 			impl_
 		);
-		if(s != MPI_SUCCESS) throw std::runtime_error("cannot get_n");
 		return it + n;
 	}
 	template<typename ContiguousIterator>
@@ -224,9 +214,12 @@ struct window<void> : basic_window{
 
 template<class T>
 struct window : window<void>{
+protected:
 	window() = default;
-	window(T* base, mpi3::size_t size, communicator& comm) : window<void>(static_cast<void*>(base), size*sizeof(T), comm){}
-	T* base() const{return window<void>::base();}
+public:
+	template<class Size = mpi3::size_t>
+	window(communicator const& c, T* b, Size n = 0) : window<void>{c, b, n}{}
+	T* base() const{return static_cast<T*>(window<void>::base());}
 	mpi3::size_t size() const{return window<void>::size()/sizeof(T);}
 };
 
@@ -238,16 +231,6 @@ class panel{
 //	friend window;
 };
 
-template<class T /*=void*/>
-window<T> communicator::make_window(T* t, mpi3::size_t n){
-	return window<T>(t, n, *this);
-}
-
-template<class T /*=void*/>
-window<T> communicator::make_window(){
-	return make_window<T>((T*)nullptr, 0);
-}
-
 template<class T> struct reference;
 
 template<class T>
@@ -293,30 +276,6 @@ shm_pointer<T> communicator::allocate_shared(MPI_Aint size) const
 }
 #endif 
 
-template<class T>
-void communicator::deallocate_shared(pointer<T>){
-//	MPI_Free_mem(p.base_ptr(rank()));
-}
-
-template<class T>
-void communicator::deallocate(pointer<T>&, MPI_Aint){
-//	p.pimpl_->fence();
-//	MPI_Free_mem(p.local_ptr());
-//	MPI_Win_free(&p.pimpl_->impl_);
-//	delete p.pimpl_;
-//	p.pimpl_ == nullptr;
-}
-
-template<class T>
-window<T> communicator::make_window(mpi3::size_t size){
-	mpi3::info inf;
-	void* ptr;
-	window<T> ret;
-	int s = MPI_Win_allocate(size*sizeof(T), sizeof(T), inf.impl_, this->impl_, &ptr, &ret.impl_);
-	if(s != MPI_SUCCESS) throw std::runtime_error("cannot window_allocate");
-	return ret;
-}
-
 }}
 
 #ifdef _TEST_BOOST_MPI3_WINDOW
@@ -329,21 +288,29 @@ using std::cout;
 
 int mpi3::main(int, char*[], mpi3::communicator world){
 
-	std::vector<double> darr(world.rank()?0:100);
-	mpi3::window<double> w = world.make_window(darr.data(), darr.size());
-	w.fence();
+	std::vector<double> darr(world.rank()?0:20);
+	std::iota(darr.begin(), darr.end(), 0);
+
+	mpi3::window<double> win{world, darr.data(), darr.size()};
+	if(world.rank() == 0){
+		std::cout << win.size() << std::endl;
+		assert( win.size() == 20 );
+		assert( win.base()[13] == 13 );
+	}else{
+		assert(win.size() == 0);
+		assert(win.base() == nullptr );
+	}
+	win.fence();
 	if(world.rank() == 0){
 		std::vector<double> a = {5., 6.};
-		w.put(a.begin(), a.end(), 0);
+		win.put(a.begin(), a.end(), 0);
 	}
-	world.barrier();
-	w.fence();
+	mpi3::communicator{win.get_group()}.barrier();
+	win.fence();
 	std::vector<double> b(2);
-	w.get(b.begin(), b.end(), 0);
-	w.fence();
-	assert( b[0] == 5.);
-	world.barrier();
-
+	win.get(b.begin(), b.end(), 0);
+	win.fence();
+	assert( b[0] == 5. and b[1] == 6. );
 	return 0;
 }
 
diff --git a/nexus/bin/nxs-test b/nexus/bin/nxs-test
index 7c0e3e49c..524808ef9 100755
--- a/nexus/bin/nxs-test
+++ b/nexus/bin/nxs-test
@@ -35,11 +35,6 @@ for d in dirs_check:
 # enforce association between nxs-test executable and Nexus library
 sys.path.insert(0,library_dir)
 
-# global data
-global_data = dict(
-    job_ref_table = False,
-    )
-
 from generic import generic_settings,object_interface,NexusError
 
 generic_settings.raise_error = True
@@ -1908,643 +1903,52 @@ def settings_operation():
 
 
 def machines():
-    # divert logging function
-    nlog_divert()
 
-    nlabel('imports')
-    from random import randint
-    from generic import obj
-    from machines import Machine,Workstation,InteractiveCluster,Supercomputer
-    from machines import Job
-    from machines import get_machine,get_machine_name
+    nunit('import')
 
+    nunit('cpu_count')
 
-    nlabel('unimplemented_base_interface')
-    arg0 = None
-    arg1 = None
-    try:
-        Machine.query_queue(arg0)
-        nfail()
-    except:
-        npass()
-    #end try
-    try:
-        Machine.submit_jobs(arg0)
-        nfail()
-    except:
-        npass()
-    #end try
-    try:
-        Machine.process_job(arg0,arg1)
-        nfail()
-    except:
-        npass()
-    #end try
-    try:
-        Machine.process_job_options(arg0,arg1)
-        nfail()
-    except:
-        npass()
-    #end try
-    try:
-        Machine.write_job(arg0,arg1,file=False)
-        nfail()
-    except:
-        npass()
-    #end try
-    try:
-        Machine.submit_job(arg0,arg1)
-        nfail()
-    except:
-        npass()
-    #end try
+    nunit('options')
 
+    nunit('job_init')
 
-    nlabel('file_setup')
-    nassert(len(Machine.machines)>0)
-    for m in Machine.machines:
-        nassert(isinstance(m,Machine))
-        exists = m.name in Machine.machines
-        nassert(exists)
-        nassert(Machine.exists(m.name))
-        nassert(Machine.is_unique(m))
-        m.validate()
-    #end for
+    nunit('job_time')
 
+    nunit('job_set_id')
 
-    nlabel('add')
-    mtest = Machine.machines.first()
-    nassert(isinstance(mtest,Machine))
-    try:
-        Machine.add(mtest)
-        nfail()
-    except:
-        npass()
-    #end try
-    try:
-        Machine.add('my_machine')
-        nfail()
-    except:
-        npass()
-    #end try
+    nunit('job_get_machine')
 
+    nunit('job_set_environment')
 
-    nlabel('get')
-    m = Machine.get(mtest.name)
-    nassert(isinstance(m,Machine))
-    nassert(id(m)==id(mtest))
-    try:
-        Machine.get(m)
-        nfail()
-    except:
-        npass()
-    #end try
-    try:
-        Machine.get('some_nonexistant_machine')
-        nfail()
-    except:
-        npass()
-    #end try
+    nunit('job_clone')
 
+    nunit('job_serial_clone')
 
-    nlabel('intantiation')
-    # test guards against empty/invalid instantiation
-    try:
-        Machine()
-        nfail()
-    except:
-        npass()
-    #end try
-    try:
-        Machine(123)
-        nfail()
-    except:
-        npass()
-    #end try
-    # test creation of a new machine
-    test_name = 'test_machine'
-    nassert(not Machine.exists(test_name))
-    m = Machine(name=test_name)
-    nassert(isinstance(m,Machine))
-    nassert(Machine.exists(m.name))
-    nassert(Machine.is_unique(m))
-    m.validate()
+    nunit('machine_virtuals')
 
-    # test guards against multiple instantiation
-    try:
-        Machine(name=test_name)
-        nfail()
-    except:
-        npass()
-    #end try
+    nunit('machine_list')
 
-    # remove test machine
-    del Machine.machines.test_machine
-    nassert(not Machine.exists(test_name))
+    nunit('machine_add')
 
+    nunit('machine_get')
 
-    # sort known machines
-    workstations   = obj()
-    supercomputers = obj()
-    for machine in Machine.machines:
-        if isinstance(machine,Workstation):
-            workstations.append(machine)
-        elif isinstance(machine,Supercomputer):
-            supercomputers[machine.name] = machine
-        else:
-            Machine.class_error('unknown machine type encountered: {0}'.format(machine.__class__.__name__),'check_process_job_idempotency')
-        #end if
-    #end for
+    nunit('machine_instantiation')
 
+    nunit('workstation_init')
 
-    nlabel('process_job_idempotency')
-    def check_process_job_idempotency(nw=10,nwj=10,nsj=10,nij=10):
-        allow_warn = Machine.allow_warnings
-        Machine.allow_warnings = False
+    nunit('workstation_scheduling')
 
-        not_idempotent = obj()
+    nunit('supercomputer_init')
 
-        # check workstations
-        nworkstations = nw
-        if nworkstations is None:
-            nworkstations=len(workstations)
-        #end if
-        nworkstations = min(nworkstations,len(workstations))
-        njobs = nwj
-        for nm in range(nworkstations):
-            if nworkstations<len(workstations):
-                machine = workstations.select_random() # select machine at random
-            else:
-                machine = workstations[nm]
-            #end if
-            cores_min     = 1
-            cores_max     = machine.cores
-            processes_min = 1
-            processes_max = machine.cores
-            threads_min   = 1
-            threads_max   = machine.cores
-            job_inputs = []
-            job_inputs_base = []
-            for nj in range(njobs): # vary cores
-                cores   = randint(cores_min,cores_max)
-                threads = randint(threads_min,threads_max)
-                job_inputs_base.append(obj(cores=cores,threads=threads))
-            #end for
-            for nj in range(njobs): # vary processes
-                processes   = randint(processes_min,processes_max)
-                threads = randint(threads_min,threads_max)
-                job_inputs_base.append(obj(processes=processes,threads=threads))
-            #end for
-            job_inputs.extend(job_inputs_base)
-            for job_input in job_inputs_base: # run in serial
-                ji = job_input.copy()
-                ji.serial = True
-                job_inputs.append(ji)
-            #end for
-            # perform idempotency test
-            machine_idempotent = True
-            for job_input in job_inputs:
-                job = Job(machine=machine.name,**job_input)
-                job2 = obj.copy(job)
-                machine.process_job(job2)
-                machine_idempotent &= job==job2
-            #end for
-            if not machine_idempotent:
-                not_idempotent[machine.name] = machine
-            #end if
-        #end for
+    nunit('supercomputer_scheduling')
 
-        # check supercomputers
-        njobs = nsj
-        small_node_ceiling = 20
-        nodes_min   = 1
-        cores_min   = 1
-        threads_min = 1
-        shared_job_inputs = obj(name='some_job',account='some_account')
-        for machine in supercomputers:
-            job_inputs = []
-            job_inputs_base = []
-            threads_max = 2*machine.cores_per_node
-            # sample small number of nodes more heavily
-            nodes_max   = min(small_node_ceiling,machine.nodes)
-            cores_max   = min(small_node_ceiling*machine.cores_per_node,machine.cores)
-            for nj in range(njobs): # nodes alone
-                nodes   = randint(nodes_min,nodes_max)
-                threads = randint(threads_min,threads_max)
-                job_input = obj(nodes=nodes,threads=threads,**shared_job_inputs)
-                job_inputs_base.append(job_input)
-            #end for
-            for nj in range(njobs): # cores alone
-                cores   = randint(cores_min,cores_max)
-                threads = randint(threads_min,threads_max)
-                job_input = obj(cores=cores,threads=threads,**shared_job_inputs)
-                job_inputs_base.append(job_input)
-            #end for
-            for nj in range(njobs): # nodes and cores
-                nodes   = randint(nodes_min,nodes_max)
-                cores   = randint(cores_min,cores_max)
-                threads = randint(threads_min,threads_max)
-                job_input = obj(nodes=nodes,cores=cores,threads=threads,**shared_job_inputs)
-                job_inputs_base.append(job_input)
-            #end for
-            # sample full node set
-            nodes_max = machine.nodes
-            cores_max = machine.cores
-            for nj in range(njobs): # nodes alone
-                nodes   = randint(nodes_min,nodes_max)
-                threads = randint(threads_min,threads_max)
-                job_input = obj(nodes=nodes,threads=threads,**shared_job_inputs)
-                job_inputs_base.append(job_input)
-            #end for
-            for nj in range(njobs): # cores alone
-                cores   = randint(cores_min,cores_max)
-                threads = randint(threads_min,threads_max)
-                job_input = obj(cores=cores,threads=threads,**shared_job_inputs)
-                job_inputs_base.append(job_input)
-            #end for
-            for nj in range(njobs): # nodes and cores
-                nodes   = randint(nodes_min,nodes_max)
-                cores   = randint(cores_min,cores_max)
-                threads = randint(threads_min,threads_max)
-                job_input = obj(nodes=nodes,cores=cores,threads=threads,**shared_job_inputs)
-                job_inputs_base.append(job_input)
-            #end for
-            job_inputs.extend(job_inputs_base)
-            # now add serial jobs
-            for job_input in job_inputs_base:
-                ji = job_input.copy()
-                ji.serial = True
-                job_inputs.append(ji)
-            #end for
-            # now add local, serial jobs
-            for job_input in job_inputs_base:
-                ji = job_input.copy()
-                ji.serial = True
-                ji.local  = True
-                job_inputs.append(ji)
-            #end for
-            # perform idempotency test
-            machine_idempotent = True
-            for job_input in job_inputs:
-                job = Job(machine=machine.name,**job_input)
-                job2 = obj.copy(job)
-                machine.process_job(job2)
-                #job_idempotent = job==job2
-                job_idempotent = object_eq(job,job2)
-                if not job_idempotent:
-                    d,d1,d2 = object_diff(job,job2,full=True)
-                    change = obj(job_before=obj(d1),job_after=obj(d2))
-                    msg = machine.name+'\n'+str(change)
-                    nfail(msg)
-                #end if
-                machine_idempotent &= job_idempotent
-            #end for
-            if not machine_idempotent:
-                not_idempotent[machine.name] = machine
-            #end if
-        #end for
+    nunit('process_job')
 
-        if len(not_idempotent)>0:
-            mlist = ''
-            for name in sorted(not_idempotent.keys()):
-                mlist+= '\n  '+name
-            #end for
-            Machine.class_error('\n\nsome machines failed process_job idempotency test:{0}'.format(mlist))
-        #end if
-        Machine.class_log('done checking idempotency')
-        Machine.allow_warnings = allow_warn
-    #end def check_process_job_idempotency
+    nunit('job_run_command')
 
-    check_process_job_idempotency()
+    nunit('write_job')
 
-
-    # test Job.run_command
-    nlabel('job_run_command')
-    def parse_job_command(command):
-        tokens = command.replace(':',' ').split()
-        launcher = tokens[0]
-        exe = tokens[-1]
-        args = []
-        options = obj()
-        last_option = None
-        for t in tokens[1:-1]:
-            if t.startswith('-'):
-                options[t] = None
-                last_option = t
-            elif last_option is not None:
-                options[last_option] = t
-                last_option = None
-            else:
-                args.append(t)
-            #end if
-        #end for
-        jc = obj(
-            launcher   = launcher,
-            executable = exe,
-            args       = args,
-            options    = options,
-            )
-        return jc
-    #end def parse_job_command
-
-    def job_commands_equal(c1,c2):
-        jc1 = parse_job_command(c1)
-        jc2 = parse_job_command(c2)
-        return object_eq(jc1,jc2)
-    #end def job_command_equal
-
-    job_run_ref = obj({
-        ('amos'           , 'n1'            ) : 'srun test.x',
-        ('amos'           , 'n1_p1'         ) : 'srun test.x',
-        ('amos'           , 'n2'            ) : 'srun test.x',
-        ('amos'           , 'n2_t2'         ) : 'srun test.x',
-        ('amos'           , 'n2_t2_e'       ) : 'srun test.x',
-        ('amos'           , 'n2_t2_p2'      ) : 'srun test.x',
-        ('bluewaters_xe'  , 'n1'            ) : 'aprun -n 32 test.x',
-        ('bluewaters_xe'  , 'n1_p1'         ) : 'aprun -n 1 test.x',
-        ('bluewaters_xe'  , 'n2'            ) : 'aprun -n 64 test.x',
-        ('bluewaters_xe'  , 'n2_t2'         ) : 'aprun -d 2 -n 32 test.x',
-        ('bluewaters_xe'  , 'n2_t2_e'       ) : 'aprun -d 2 -n 32 test.x',
-        ('bluewaters_xe'  , 'n2_t2_p2'      ) : 'aprun -d 2 -n 4 test.x',
-        ('bluewaters_xk'  , 'n1'            ) : 'aprun -n 16 test.x',
-        ('bluewaters_xk'  , 'n1_p1'         ) : 'aprun -n 1 test.x',
-        ('bluewaters_xk'  , 'n2'            ) : 'aprun -n 32 test.x',
-        ('bluewaters_xk'  , 'n2_t2'         ) : 'aprun -d 2 -n 16 test.x',
-        ('bluewaters_xk'  , 'n2_t2_e'       ) : 'aprun -d 2 -n 16 test.x',
-        ('bluewaters_xk'  , 'n2_t2_p2'      ) : 'aprun -d 2 -n 4 test.x',
-        ('cades'          , 'n1'            ) : 'mpirun -np 36 test.x',
-        ('cades'          , 'n1_p1'         ) : 'mpirun -np 1 test.x',
-        ('cades'          , 'n2'            ) : 'mpirun -np 72 test.x',
-        ('cades'          , 'n2_t2'         ) : 'mpirun -np 36 --npersocket 9 test.x',
-        ('cades'          , 'n2_t2_e'       ) : 'mpirun -np 36 --npersocket 9 test.x',
-        ('cades'          , 'n2_t2_p2'      ) : 'mpirun -np 4 --npersocket 1 test.x',
-        ('cetus'          , 'n1'            ) : 'runjob --np 16 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
-        ('cetus'          , 'n1_p1'         ) : 'runjob --np 1 -p 1 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
-        ('cetus'          , 'n2'            ) : 'runjob --np 32 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
-        ('cetus'          , 'n2_t2'         ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
-        ('cetus'          , 'n2_t2_e'       ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 ENV_VAR=1 : test.x',
-        ('cetus'          , 'n2_t2_p2'      ) : 'runjob --np 4 -p 2 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
-        ('chama'          , 'n1'            ) : 'srun test.x',
-        ('chama'          , 'n1_p1'         ) : 'srun test.x',
-        ('chama'          , 'n2'            ) : 'srun test.x',
-        ('chama'          , 'n2_t2'         ) : 'srun test.x',
-        ('chama'          , 'n2_t2_e'       ) : 'srun test.x',
-        ('chama'          , 'n2_t2_p2'      ) : 'srun test.x',
-        ('cooley'         , 'n1'            ) : 'mpirun -np 12 test.x',
-        ('cooley'         , 'n1_p1'         ) : 'mpirun -np 1 test.x',
-        ('cooley'         , 'n2'            ) : 'mpirun -np 24 test.x',
-        ('cooley'         , 'n2_t2'         ) : 'mpirun -np 12 test.x',
-        ('cooley'         , 'n2_t2_e'       ) : 'mpirun -np 12 test.x',
-        ('cooley'         , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
-        ('cori'           , 'n1'            ) : 'srun test.x',
-        ('cori'           , 'n1_p1'         ) : 'srun test.x',
-        ('cori'           , 'n2'            ) : 'srun test.x',
-        ('cori'           , 'n2_t2'         ) : 'srun test.x',
-        ('cori'           , 'n2_t2_e'       ) : 'srun test.x',
-        ('cori'           , 'n2_t2_p2'      ) : 'srun test.x',
-        ('edison'         , 'n1'            ) : 'srun test.x',
-        ('edison'         , 'n1_p1'         ) : 'srun test.x',
-        ('edison'         , 'n2'            ) : 'srun test.x',
-        ('edison'         , 'n2_t2'         ) : 'srun test.x',
-        ('edison'         , 'n2_t2_e'       ) : 'srun test.x',
-        ('edison'         , 'n2_t2_p2'      ) : 'srun test.x',
-        ('eos'            , 'n1'            ) : 'aprun -n 16 test.x',
-        ('eos'            , 'n1_p1'         ) : 'aprun -n 1 test.x',
-        ('eos'            , 'n2'            ) : 'aprun -n 32 test.x',
-        ('eos'            , 'n2_t2'         ) : 'aprun -ss -cc numa_node -d 2 -n 16 test.x',
-        ('eos'            , 'n2_t2_e'       ) : 'aprun -ss -cc numa_node -d 2 -n 16 test.x',
-        ('eos'            , 'n2_t2_p2'      ) : 'aprun -ss -cc numa_node -d 2 -n 4 test.x',
-        ('jaguar'         , 'n1'            ) : 'aprun -n 16 test.x',
-        ('jaguar'         , 'n1_p1'         ) : 'aprun -n 1 test.x',
-        ('jaguar'         , 'n2'            ) : 'aprun -n 32 test.x',
-        ('jaguar'         , 'n2_t2'         ) : 'aprun -d 2 -n 16 test.x',
-        ('jaguar'         , 'n2_t2_e'       ) : 'aprun -d 2 -n 16 test.x',
-        ('jaguar'         , 'n2_t2_p2'      ) : 'aprun -d 2 -n 4 test.x',
-        ('komodo'         , 'n1'            ) : 'mpirun -np 12 test.x',
-        ('komodo'         , 'n1_p1'         ) : 'mpirun -np 1 test.x',
-        ('komodo'         , 'n2'            ) : 'mpirun -np 24 test.x',
-        ('komodo'         , 'n2_t2'         ) : 'mpirun -np 12 test.x',
-        ('komodo'         , 'n2_t2_e'       ) : 'mpirun -np 12 test.x',
-        ('komodo'         , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
-        ('kraken'         , 'n1'            ) : 'aprun -n 12 test.x',
-        ('kraken'         , 'n1_p1'         ) : 'aprun -n 1 test.x',
-        ('kraken'         , 'n2'            ) : 'aprun -n 24 test.x',
-        ('kraken'         , 'n2_t2'         ) : 'aprun -d 2 -n 12 test.x',
-        ('kraken'         , 'n2_t2_e'       ) : 'aprun -d 2 -n 12 test.x',
-        ('kraken'         , 'n2_t2_p2'      ) : 'aprun -d 2 -n 4 test.x',
-        ('lonestar'       , 'n1'            ) : 'ibrun -n 12 -o 0 test.x',
-        ('lonestar'       , 'n1_p1'         ) : 'ibrun -n 1 -o 0 test.x',
-        ('lonestar'       , 'n2'            ) : 'ibrun -n 24 -o 0 test.x',
-        ('lonestar'       , 'n2_t2'         ) : 'ibrun -n 12 -o 0 test.x',
-        ('lonestar'       , 'n2_t2_e'       ) : 'ibrun -n 12 -o 0 test.x',
-        ('lonestar'       , 'n2_t2_p2'      ) : 'ibrun -n 4 -o 0 test.x',
-        ('matisse'        , 'n1'            ) : 'mpirun -np 16 test.x',
-        ('matisse'        , 'n1_p1'         ) : 'mpirun -np 1 test.x',
-        ('matisse'        , 'n2'            ) : 'mpirun -np 32 test.x',
-        ('matisse'        , 'n2_t2'         ) : 'mpirun -np 16 test.x',
-        ('matisse'        , 'n2_t2_e'       ) : 'mpirun -np 16 test.x',
-        ('matisse'        , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
-        ('mira'           , 'n1'            ) : 'runjob --np 16 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
-        ('mira'           , 'n1_p1'         ) : 'runjob --np 1 -p 1 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
-        ('mira'           , 'n2'            ) : 'runjob --np 32 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
-        ('mira'           , 'n2_t2'         ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
-        ('mira'           , 'n2_t2_e'       ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 ENV_VAR=1 : test.x',
-        ('mira'           , 'n2_t2_p2'      ) : 'runjob --np 4 -p 2 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
-        ('oic5'           , 'n1'            ) : 'mpirun -np 32 test.x',
-        ('oic5'           , 'n1_p1'         ) : 'mpirun -np 1 test.x',
-        ('oic5'           , 'n2'            ) : 'mpirun -np 64 test.x',
-        ('oic5'           , 'n2_t2'         ) : 'mpirun -np 32 test.x',
-        ('oic5'           , 'n2_t2_e'       ) : 'mpirun -np 32 test.x',
-        ('oic5'           , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
-        ('redsky'         , 'n1'            ) : 'srun test.x',
-        ('redsky'         , 'n1_p1'         ) : 'srun test.x',
-        ('redsky'         , 'n2'            ) : 'srun test.x',
-        ('redsky'         , 'n2_t2'         ) : 'srun test.x',
-        ('redsky'         , 'n2_t2_e'       ) : 'srun test.x',
-        ('redsky'         , 'n2_t2_p2'      ) : 'srun test.x',
-        ('serrano'        , 'n1'            ) : 'srun test.x',
-        ('serrano'        , 'n1_p1'         ) : 'srun test.x',
-        ('serrano'        , 'n2'            ) : 'srun test.x',
-        ('serrano'        , 'n2_t2'         ) : 'srun test.x',
-        ('serrano'        , 'n2_t2_e'       ) : 'srun test.x',
-        ('serrano'        , 'n2_t2_p2'      ) : 'srun test.x',
-        ('skybridge'      , 'n1'            ) : 'srun test.x',
-        ('skybridge'      , 'n1_p1'         ) : 'srun test.x',
-        ('skybridge'      , 'n2'            ) : 'srun test.x',
-        ('skybridge'      , 'n2_t2'         ) : 'srun test.x',
-        ('skybridge'      , 'n2_t2_e'       ) : 'srun test.x',
-        ('skybridge'      , 'n2_t2_p2'      ) : 'srun test.x',
-        ('solo'           , 'n1'            ) : 'srun test.x',
-        ('solo'           , 'n1_p1'         ) : 'srun test.x',
-        ('solo'           , 'n2'            ) : 'srun test.x',
-        ('solo'           , 'n2_t2'         ) : 'srun test.x',
-        ('solo'           , 'n2_t2_e'       ) : 'srun test.x',
-        ('solo'           , 'n2_t2_p2'      ) : 'srun test.x',
-        ('stampede2'      , 'n1'            ) : 'ibrun -n 68 -o 0 test.x',
-        ('stampede2'      , 'n1_p1'         ) : 'ibrun -n 1 -o 0 test.x',
-        ('stampede2'      , 'n2'            ) : 'ibrun -n 136 -o 0 test.x',
-        ('stampede2'      , 'n2_t2'         ) : 'ibrun -n 68 -o 0 test.x',
-        ('stampede2'      , 'n2_t2_e'       ) : 'ibrun -n 68 -o 0 test.x',
-        ('stampede2'      , 'n2_t2_p2'      ) : 'ibrun -n 4 -o 0 test.x',
-        ('summit'         , 'n1'            ) : 'jsrun -a 21 -r 2 -b rs -c 21 -d packed -n 2 -g 0 test.x',
-        ('summit'         , 'n1_g6'         ) : 'jsrun -a 7 -r 6 -b rs -c 7 -d packed -n 6 -g 1 test.x',
-        ('summit'         , 'n2'            ) : 'jsrun -a 21 -r 2 -b rs -c 21 -d packed -n 4 -g 0 test.x',
-        ('summit'         , 'n2_g6'         ) : 'jsrun -a 7 -r 6 -b rs -c 7 -d packed -n 12 -g 1 test.x',
-        ('summit'         , 'n2_t2'         ) : 'jsrun -a 10 -r 2 -b rs -c 20 -d packed -n 4 -g 0 test.x',
-        ('summit'         , 'n2_t2_e'       ) : 'jsrun -a 10 -r 2 -b rs -c 20 -d packed -n 4 -g 0 test.x',
-        ('summit'         , 'n2_t2_e_g6'    ) : 'jsrun -a 3 -r 6 -b rs -c 6 -d packed -n 12 -g 1 test.x',
-        ('summit'         , 'n2_t2_g6'      ) : 'jsrun -a 3 -r 6 -b rs -c 6 -d packed -n 12 -g 1 test.x',
-        ('supermuc'       , 'n1'            ) : 'mpiexec -n 28 test.x',
-        ('supermuc'       , 'n1_p1'         ) : 'mpiexec -n 1 test.x',
-        ('supermuc'       , 'n2'            ) : 'mpiexec -n 56 test.x',
-        ('supermuc'       , 'n2_t2'         ) : 'mpiexec -n 28 test.x',
-        ('supermuc'       , 'n2_t2_e'       ) : 'mpiexec -n 28 test.x',
-        ('supermuc'       , 'n2_t2_p2'      ) : 'mpiexec -n 4 test.x',
-        ('supermucng'     , 'n1'            ) : 'mpiexec -n 48 test.x',
-        ('supermucng'     , 'n1_p1'         ) : 'mpiexec -n 1 test.x',
-        ('supermucng'     , 'n2'            ) : 'mpiexec -n 96 test.x',
-        ('supermucng'     , 'n2_t2'         ) : 'mpiexec -n 48 test.x',
-        ('supermucng'     , 'n2_t2_e'       ) : 'mpiexec -n 48 test.x',
-        ('supermucng'     , 'n2_t2_p2'      ) : 'mpiexec -n 4 test.x',
-        ('taub'           , 'n1'            ) : 'mpirun -np 12 test.x',
-        ('taub'           , 'n1_p1'         ) : 'mpirun -np 1 test.x',
-        ('taub'           , 'n2'            ) : 'mpirun -np 24 test.x',
-        ('taub'           , 'n2_t2'         ) : 'mpirun -np 12 test.x',
-        ('taub'           , 'n2_t2_e'       ) : 'mpirun -np 12 test.x',
-        ('taub'           , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
-        ('theta'          , 'n1'            ) : 'aprun -e OMP_NUM_THREADS=1 -d 1 -cc depth -j 1 -n 64 -N 64 test.x',
-        ('theta'          , 'n1_p1'         ) : 'aprun -e OMP_NUM_THREADS=1 -d 1 -cc depth -j 1 -n 1 -N 1 test.x',
-        ('theta'          , 'n2'            ) : 'aprun -e OMP_NUM_THREADS=1 -d 1 -cc depth -j 1 -n 128 -N 64 test.x',
-        ('theta'          , 'n2_t2'         ) : 'aprun -e OMP_NUM_THREADS=2 -d 2 -cc depth -j 1 -n 64 -N 32 test.x',
-        ('theta'          , 'n2_t2_e'       ) : 'aprun -e OMP_NUM_THREADS=2 -d 2 -cc depth -j 1 -n 64 -N 32 test.x',
-        ('theta'          , 'n2_t2_p2'      ) : 'aprun -e OMP_NUM_THREADS=2 -d 2 -cc depth -j 1 -n 4 -N 2 test.x',
-        ('titan'          , 'n1'            ) : 'aprun -n 16 test.x',
-        ('titan'          , 'n1_p1'         ) : 'aprun -n 1 test.x',
-        ('titan'          , 'n2'            ) : 'aprun -n 32 test.x',
-        ('titan'          , 'n2_t2'         ) : 'aprun -d 2 -n 16 test.x',
-        ('titan'          , 'n2_t2_e'       ) : 'aprun -d 2 -n 16 test.x',
-        ('titan'          , 'n2_t2_p2'      ) : 'aprun -d 2 -n 4 test.x',
-        ('tomcat3'        , 'n1'            ) : 'mpirun -np 64 test.x',
-        ('tomcat3'        , 'n1_p1'         ) : 'mpirun -np 1 test.x',
-        ('tomcat3'        , 'n2'            ) : 'mpirun -np 128 test.x',
-        ('tomcat3'        , 'n2_t2'         ) : 'mpirun -np 64 test.x',
-        ('tomcat3'        , 'n2_t2_e'       ) : 'mpirun -np 64 test.x',
-        ('tomcat3'        , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
-        ('uno'            , 'n1'            ) : 'srun test.x',
-        ('uno'            , 'n1_p1'         ) : 'srun test.x',
-        ('uno'            , 'n2'            ) : 'srun test.x',
-        ('uno'            , 'n2_t2'         ) : 'srun test.x',
-        ('uno'            , 'n2_t2_e'       ) : 'srun test.x',
-        ('uno'            , 'n2_t2_p2'      ) : 'srun test.x',
-        ('vesta'          , 'n1'            ) : 'runjob --np 16 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
-        ('vesta'          , 'n1_p1'         ) : 'runjob --np 1 -p 1 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
-        ('vesta'          , 'n2'            ) : 'runjob --np 32 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
-        ('vesta'          , 'n2_t2'         ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
-        ('vesta'          , 'n2_t2_e'       ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 ENV_VAR=1 : test.x',
-        ('vesta'          , 'n2_t2_p2'      ) : 'runjob --np 4 -p 2 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
-        })
-
-    job_inputs_orig = obj(
-        n1        = obj(nodes=1),
-        n1_p1     = obj(nodes=1,processes_per_node=1),
-        n2        = obj(nodes=2),
-        n2_t2     = obj(nodes=2,threads=2),
-        n2_t2_p2  = obj(nodes=2,threads=2,processes_per_node=2),
-        n2_t2_e   = obj(nodes=2,threads=2,env=obj(ENV_VAR=1)),
-        )
-    for name in sorted(supercomputers.keys()):
-        m = supercomputers[name]
-        if m.requires_account:
-            acc = 'ABC123'
-        else:
-            acc = None
-        #end if
-        job_inputs = job_inputs_orig
-        if name=='summit': # exceptional treatment for summit nodes
-            job_inputs = job_inputs_orig.copy()
-            jtypes = list(job_inputs.keys())
-            for jtype in jtypes:
-                if 'p' in jtype:
-                    del job_inputs[jtype]
-                else:
-                    jcpu = job_inputs[jtype]
-                    jcpu.gpus = 0
-                    jgpu = jcpu.copy()
-                    jgpu.gpus = 6
-                    job_inputs[jtype+'_g6'] = jgpu
-                #end if
-            #end for
-        #end if
-        for jtype in sorted(job_inputs.keys()):
-            job = Job(app_command = 'test.x',
-                      machine     = name,
-                      account     = acc,
-                      **job_inputs[jtype]
-                      )
-            command = job.run_command()
-            if global_data['job_ref_table']:
-                sname = "'{0}'".format(name)
-                stype = "'{0}'".format(jtype)
-                nlog("        ({0:<16} , {1:<16}) : '{2}',".format(sname,stype,command))
-                continue
-            #end if
-            ref_command = job_run_ref[name,jtype]
-            if not job_commands_equal(command,ref_command):
-                nfail('Job.run_command for machine "{0}" does not match the reference\njob inputs: {1}\nreference command: {2}\nincorrect command: {3}'.format(name,job_inputs[jtype],ref_command,command))
-            #end for
-        #end for
-    #end for
-
-
-    nlabel('Job.split_nodes')
-    for name in sorted(supercomputers.keys()):
-        m = supercomputers[name]
-        if m.app_launcher=='srun': # no slurm support yet
-            continue
-        #end if
-        if name=='summit': # no summit support
-            continue
-        #end if
-        if m.requires_account:
-            acc = 'ABC123'
-        else:
-            acc = None
-        #end if
-        # make a 4 node job
-        job = Job(app_command = 'test.x',
-                  machine     = name,
-                  account     = acc,
-                  nodes       = 4,
-                  threads     = m.cores_per_node,
-                  )
-        # split the job into 1 and 3 nodes
-        job1,job2 = job.split_nodes(1)
-        # get the split run commands
-        rc  = job.run_command()
-        rc1 = job1.run_command()
-        rc2 = job2.run_command()
-        ns  = ' {0} '.format(job.nodes)
-        ns1 = ' {0} '.format(job1.nodes)
-        ns2 = ' {0} '.format(job2.nodes)
-        # verify that node count is in each command
-        nassert(ns  in rc )
-        nassert(ns1 in rc1)
-        nassert(ns2 in rc2)
-        # verify that text on either side of node count 
-        # agrees for original and split commands
-        rcl ,rcr  = rc.split(ns,1)
-        rc1l,rc1r = rc1.split(ns1,1)
-        rc2l,rc2r = rc2.split(ns2,1)
-        rcf  = rcl+' '+rcr
-        rc1f = rc1l+' '+rc1r
-        rc2f = rc2l+' '+rc2r
-        nassert(job_commands_equal(rcf,rc1f))
-        nassert(job_commands_equal(rcf,rc2f))
-    #end for
-
-
-    # restore logging function
-    nlog_restore()
+    nunit_all()
 #end def machines
 
 
diff --git a/nexus/lib/machines.py b/nexus/lib/machines.py
index 39f9ee3d3..69c556df3 100644
--- a/nexus/lib/machines.py
+++ b/nexus/lib/machines.py
@@ -57,7 +57,7 @@ from imp import load_source
 
 
 import re,subprocess
-def  cpu_count():
+def cpu_count():
     """ Number of virtual or physical CPUs on this system, i.e.
     user/real as output by time(1) when called with an optimally scaling
     userspace-only program"""
@@ -90,6 +90,11 @@ class Options(DevBase):
     #end def __init__
 
 
+    def add(self,**kwargs):
+        self.transfer_from(kwargs)
+    #end def add
+
+
     def read(self,options):
         nopts = 0
         intext = False
@@ -119,11 +124,6 @@ class Options(DevBase):
     #end def read
 
 
-    def add(self,**kwargs):
-        self.transfer_from(kwargs)
-    #end def add
-
-
     def write(self):
         s = ''
         for o in self:
@@ -135,7 +135,7 @@ class Options(DevBase):
 
 
 
-job_defaults = obj(
+job_defaults_assign = obj(
     name               = 'jobname',
     type               = None,
     directory          = None,
@@ -159,7 +159,6 @@ job_defaults = obj(
     hyperthreads       = None,
     ppn                = None,
     gpus               = None, # number of gpus per node
-    compiler           = None,
     serial             = False, # run job serially, no mpi
     local              = False, # run job locally, no queue submission
     days               = 0,
@@ -183,7 +182,10 @@ job_defaults = obj(
     default_cpus_per_task = False, # optionally bypass typical nexus processing for supermucng
     ntasks_per_core    = None,
     cpus_per_task      = None,
+    )
+
     # these are not assigned directly
+job_defaults_nonassign = obj(
     fake               = False,
     app                = None, # name of/path to application
     machine            = None,
@@ -192,8 +194,11 @@ job_defaults = obj(
     app_options        = None,
     run_options        = None,
     sub_options        = None,
+    skip_machine       = False,
     )
 
+job_defaults = obj(job_defaults_assign,job_defaults_nonassign)
+
 class Job(NexusCore):
 
     machine = None #default machine if none is specified in settings
@@ -210,12 +215,6 @@ class Job(NexusCore):
     job_count = 0
 
 
-    intel_compilers = set(['intel','icc','icpc','ifort'])
-    gnu_compilers   = set(['gnu','gcc','g++','gfortran'])
-    pgi_compilers   = set(['pgi'])
-    cray_compilers  = set(['cray'])
-
-
     @staticmethod
     def restore_default_settings():
         Job.machine = None
@@ -247,15 +246,16 @@ class Job(NexusCore):
         kw.set_optional(**job_defaults)
 
         # extract keywords not assigned
-        app         = kw.delete('app')
-        machine     = kw.delete('machine')
-        options     = kw.delete('options')
-        app_flags   = kw.delete('app_flags')
-        app_options = kw.delete('app_options')
-        run_options = kw.delete('run_options')
-        sub_options = kw.delete('sub_options')
-        env         = kw.delete('env')
-        fake        = kw.delete('fake')
+        app          = kw.delete('app')
+        machine      = kw.delete('machine')
+        options      = kw.delete('options')
+        app_flags    = kw.delete('app_flags')
+        app_options  = kw.delete('app_options')
+        run_options  = kw.delete('run_options')
+        sub_options  = kw.delete('sub_options')
+        env          = kw.delete('env')
+        fake         = kw.delete('fake')
+        skip_machine = kw.delete('skip_machine')
 
         # assign keywords
         self.set(**kw)
@@ -303,37 +303,30 @@ class Job(NexusCore):
         if self.app_props is None:
             self.app_props = []
         #end if
-        if self.compiler is not None:
-            if self.compiler in self.intel_compilers:
-                self.compiler = 'intel'
-            elif self.compiler in self.gnu_compilers:
-                self.compiler = 'gnu'
-            elif self.compiler in self.pgi_compilers:
-                self.compiler = 'pgi'
-            elif self.compiler in self.cray_compilers:
-                self.compiler = 'cray'
-            #end if
-        #end if
 
         if self.serial:
             self.cores = 1
             self.nodes = None
         #end if
 
-        if machine is not None:
+        if skip_machine:
             self.machine = machine
-        #end if
-        if env is not None:
-            self.set_environment(**env)
-        #end if
-        #check that the machine exists and have it complete the job info
-        self.process()
+        else:
+            if machine is not None:
+                self.machine = machine
+            #end if
+            if env is not None:
+                self.set_environment(**env)
+            #end if
+            #check that the machine exists and have it complete the job info
+            self.process()
 
-        machine = self.get_machine() 
-        self.batch_mode = machine.in_batch_mode()
+            machine = self.get_machine() 
+            self.batch_mode = machine.in_batch_mode()
 
-        if self.bundled_jobs is not None and not machine.batch_capable:
-            self.error('running batched/bundled jobs on {0} is either not possible or not yet implemented, sorry.'.format(machine.name))
+            if self.bundled_jobs is not None and not machine.batch_capable:
+                self.error('running batched/bundled jobs on {0} is either not possible or not yet implemented, sorry.'.format(machine.name))
+            #end if
         #end if
 
         self.normalize_time()
@@ -353,6 +346,7 @@ class Job(NexusCore):
     #end def process
 
 
+    # test needed
     def process_options(self,machine=None):
         if machine is None:
             machine = self.get_machine()
@@ -361,6 +355,7 @@ class Job(NexusCore):
     #end def process_options
 
 
+    # test needed
     def initialize(self,sim):
         self.set_id()
         self.identifier = sim.identifier
@@ -420,6 +415,7 @@ class Job(NexusCore):
     #end def initialize
 
 
+    # test needed
     def renew_app_command(self,sim):
         if not self.user_app_command:
             self.app_command = sim.app_command()
@@ -432,6 +428,7 @@ class Job(NexusCore):
     #end def set_id
 
 
+    # remove?
     def set_processes(self):
         if self.processes is None:
             self.error('processes should have been set before now\ncontact the developers and have them fix this','Developer')
@@ -462,10 +459,6 @@ class Job(NexusCore):
 
     def divert_out_err(self):
         self.identifier += '_divert'
-        #prefix,ext = self.outfile.split('.',1)
-        #self.outfile = prefix+'_divert.'+ext
-        #prefix,ext = self.errfile.split('.',1)
-        #self.errfile = prefix+'_divert.'+ext
     #end def divert_out_err
 
 
@@ -498,28 +491,22 @@ class Job(NexusCore):
     #end if
 
 
+    # remove?
     def determine_end_status(self,status):
         if not nexus_core.generate_only:
             self.successful = False # not really implemented yet
-            #errfile = os.path.join(self.directory,self.errfile)
-            #if os.path.exists(errfile):
-            #    fobj = open(errfile,'r')
-            #    errors = fobj.read()
-            #    fobj.close()
-            #    if errors.strip()=='':
-            #        self.successful=True
-            #    #end if
-            ##end if
         #end if
     #end def determine_end_status
 
 
+    # test needed
     def write(self,file=False):
         machine = self.get_machine()
         return machine.write_job(self,file=file)
     #end def write
 
 
+    # test needed
     def submit(self):
         machine = self.get_machine()
         machine.add_job(self)
@@ -527,6 +514,7 @@ class Job(NexusCore):
     #end def submit
 
 
+    # test needed
     def reenter_queue(self):
         machine = self.get_machine()
         machine.requeue_job(self)
@@ -856,6 +844,7 @@ class Machine(NexusCore):
     
     def restore_default_settings(self):
         self.account         = None
+        self.user            = None
         self.local_directory = None
         self.app_directory   = None
         self.app_directories = None
@@ -1133,6 +1122,7 @@ class Workstation(Machine):
 
 
 
+# test needed
 class InteractiveCluster(Workstation):
 
     def __init__(self,*args,**kwargs):
@@ -1243,7 +1233,7 @@ class Supercomputer(Machine):
                  sub_launcher   = None,
                  queue_querier  = None,
                  job_remover    = None,
-                 name           = None
+                 name           = None,
                  ):
         if name is None:
             if self.name is not None:
@@ -1372,6 +1362,8 @@ class Supercomputer(Machine):
                                  USUSP = 'suspended',
                                  SSUSP = 'suspended',
                                  )
+        elif self.queue_querier=='test_query':
+            None
         else:
             self.error('ability to query queue with '+self.queue_querier+' has not yet been implemented')
         #end if
@@ -1379,11 +1371,13 @@ class Supercomputer(Machine):
     #end def __init__
 
 
+    # test needed
     def interactive_representation(self,cores):
         return InteractiveCluster(self,cores)
     #end def interactive_representation
 
 
+    # test needed
     def requeue_job(self,job):
         if isinstance(job,Job):
             jid = job.internal_id
@@ -1693,6 +1687,11 @@ class Supercomputer(Machine):
                     #end if
                 #end if
             #end for
+        elif self.queue_querier=='test_query': # for testing
+            # pretend that all jobs have finished
+            for pid in self.processes.keys():
+                self.system_queue[pid] = 'complete'
+            #end for
         else:
             self.error('ability to query queue with '+self.queue_querier+' has not yet been implemented')
         #end if
@@ -3281,7 +3280,6 @@ class Tomcat3(Supercomputer):
 #  workstations
 for cores in range(1,128+1):
     Workstation('ws'+str(cores),cores,'mpirun'),
-    Workstation('node'+str(cores),cores,'mpirun'),
 #end for
 #  supercomputers and clusters
 #            nodes sockets cores ram qslots  qlaunch  qsubmit     qstatus    qdelete
diff --git a/nexus/lib/numerics.py b/nexus/lib/numerics.py
index 914faaab2..78014604c 100644
--- a/nexus/lib/numerics.py
+++ b/nexus/lib/numerics.py
@@ -83,9 +83,10 @@
 
 import sys
 import inspect
+import numpy as np
 from numpy import array,ndarray,zeros,linspace,pi,exp,sqrt,polyfit,polyval
 from numpy import sum,abs,arange,empty,sin,cos,dot,atleast_2d,ogrid
-from numpy import ones_like,sign,random,cross,prod
+from numpy import ones_like,sign,cross,prod
 from numpy.linalg import norm
 from generic import obj
 from developer import unavailable,warn,error
@@ -939,7 +940,7 @@ def equilibration_length(x,tail=.5,plot=False,xlim=None,bounces=2,random=True):
         bounce[1] = max(bounce[1],bounce[0])
         #print len(x),crossings,crossings[1]-crossings[0]+1
         if random:
-            eqlen = bounce[0]+random.randint(bounce[1]-bounce[0]+1)
+            eqlen = bounce[0]+np.random.randint(bounce[1]-bounce[0]+1)
         else:
             eqlen = (bounce[0]+bounce[1])//2
         #end if
diff --git a/nexus/lib/qmcpack_method_analyzers.py b/nexus/lib/qmcpack_method_analyzers.py
index d96773743..5f8d02b51 100644
--- a/nexus/lib/qmcpack_method_analyzers.py
+++ b/nexus/lib/qmcpack_method_analyzers.py
@@ -79,7 +79,7 @@ class MethodAnalyzer(QAanalyzer):
             complete &= 'opt' in files
         #end if
         equil = request.equilibration
-        nblocks_exclude = 0
+        nblocks_exclude = -1
         if isinstance(equil,int):
             nblocks_exclude = equil
         elif isinstance(equil,(dict,obj)):
diff --git a/nexus/lib/qmcpack_quantity_analyzers.py b/nexus/lib/qmcpack_quantity_analyzers.py
index 192445609..8e0489be5 100644
--- a/nexus/lib/qmcpack_quantity_analyzers.py
+++ b/nexus/lib/qmcpack_quantity_analyzers.py
@@ -134,6 +134,7 @@ class DatAnalyzer(QuantityAnalyzer):
         if equilibration!=None and nbe==-1:
             self.load_data()
             nbe = equilibration_length(self.data[equilibration])
+            assert nbe>=0, 'Number of equilibration blocks is negative.'
             self.method_info.nblocks_exclude = nbe
         #end if
     #end def __init__
diff --git a/nexus/lib/structure.py b/nexus/lib/structure.py
index 8d37ade34..e51f9c865 100644
--- a/nexus/lib/structure.py
+++ b/nexus/lib/structure.py
@@ -26,6 +26,100 @@
 
 #! /usr/bin/env python
 
+"""
+The :py:mod:`structure` module provides support for atomic structure I/O,
+generation, and manipulation.  
+
+
+List of module contents
+-----------------------
+
+Read cif file functions:
+
+* :py:func:`read_cif_celldata`
+* :py:func:`read_cif_cell`
+* :py:func:`read_cif`
+
+Operations on logical conditions:
+
+* :py:func:`equate`
+* :py:func:`negate`
+
+Create a Monkhorst-Pack k-point mesh function
+
+* :py:func:`kmesh`
+
+Tile matrix malipulation functions
+
+* :py:func:`reduce_tilematrix`
+* :py:func:`tile_magnetization`
+
+Rotate plane function
+
+* :py:func:`rotate_plane`
+
+Trivial filter function
+
+* :py:func:`trivial_filter`
+
+* :py:class:`MaskFilter`
+
+* :py:func:`optimal_tilematrix`
+
+Base class for :py:class:`Structure` class:
+
+* :py:class:`Sobj`
+
+Base class for :py:class:`DefectStructure`, :py:class:`Crystal`, and :py:class:`Jellium` classes:
+
+* :py:class:`Structure`
+
+SeeK-path functions
+
+* :py:func:`\_getseekpath`
+* :py:func:`get_conventional_cell`
+* :py:func:`get_primitive_cell`
+* :py:func:`get_kpath`
+* :py:func:`get_symmetry`
+* :py:func:`get_structure_with_bands`
+* :py:func:`get_band_tiling`
+* :py:func:`get_seekpath_full`
+
+Interpolate structures functions
+
+* :py:func:`interpolate_structures`
+
+Animate structures functions
+
+* :py:func:`structure_animation`
+
+Concrete :py:class:`Structure` classes:
+
+* :py:class:`DefectStructure`
+* :py:class:`Crystal`
+* :py:class:`Jellium`
+
+Structure generation functions:
+
+* :py:func:`generate_cell`
+* :py:func:`generate_structure`
+* :py:func:`generate_atom_structure`
+* :py:func:`generate_dimer_structure`
+* :py:func:`generate_trimer_structure`
+* :py:func:`generate_jellium_structure`
+* :py:func:`generate_crystal_structure`
+* :py:func:`generate_defect_structure`
+
+Read structure functions
+
+* :py:func:`read_structure`
+
+
+
+Module contents
+---------------
+"""
+
 import os
 import numpy as np
 from copy import deepcopy
@@ -1457,21 +1551,149 @@ class Structure(Sobj):
         self.skew(d)
     #end def stretch
 
-        
+
     # test needed
-    def skew(self,skew):
+    def rotate(self,r,rp=None,passive=False,units="radians",check=True):
+        """
+        Arbitrary rotation of the structure.
+        Parameters
+        ----------
+        r  : `array_like, float, shape (3,3)` or `array_like, float, shape (3)` or `str`
+            If a 3x3 matrix, then code executes rotation consistent with this matrix -- 
+            it is assumed that the matrix acts on a column-major vector (eg, v'=Rv)
+            If a three-dimensional array, then the operation of the function depends
+            on the input type of rp in the following ways:
+                1. If rp is a scalar, then rp is assumed to be an angle and a rotation 
+                   of rp is made about the axis defined by r
+                2. If rp is a vector, then rp is assumed to be an axis and a rotation is made 
+                   such that r aligns with rp
+                3. If rp is a str, then the rotation is such that r aligns with the
+                   axis given by the str ('x', 'y', 'z', 'a0', 'a1', or 'a2')
+            If a str then the axis, r, is defined by the input label (e.g. 'x', 'y', 'z', 'a1', 'a2', or 'a3')
+            and the operation of the function depends on the input type of rp in the following
+            ways (same as above):
+                1. If rp is a scalar, then rp is assumed to be an angle and a rotation 
+                   of rp is made about the axis defined by r
+                2. If rp is a vector, then rp is assumed to be an axis and a rotation is made 
+                   such that r aligns with rp
+                3. If rp is a str, then the rotation is such that r aligns with the
+                   axis given by the str ('x', 'y', 'z', 'a0', 'a1', or 'a2')
+        rp : `array_like, float, shape (3), optional` or `str, optional`
+            If a 3-dimensional vector is given, then rp is assumed to be an axis and a rotation is made
+            such that the axis r is aligned with rp.
+            If a str, then rp is assumed to be an angle and a rotation about the axis defined by r 
+            is made by an angle rp
+            If a str is given, then rp is assumed to be an axis defined by the given label
+            (e.g. 'x', 'y', 'z', 'a1', 'a2', or 'a3') and a rotation is made such that the axis r 
+            is aligned with rp.
+        passive : `bool, optional, default False`
+            If `True`, perform a passive rotation
+            If `False`, perform an active rotation
+        units : `str, optional, default "radians"`
+            Units of rp, if rp is given as an angle (scalar)
+        check : `bool, optional, default True`
+            Perform a check to verify rotation matrix is orthogonal
+        """
+        if rp is not None:
+            dirmap = dict(x=[1,0,0],y=[0,1,0],z=[0,0,1])
+            if isinstance(r,str): 
+                if r[0]=='a': # r= 'a0', 'a1', or 'a2'
+                    r = self.axes[int(r[1])]
+                else: # r= 'x', 'y', or 'z'
+                    r = dirmap[r]
+                #end if
+            else:
+                r = array(r,dtype=float)
+                if len(r.shape)>1:
+                    self.error('r must be given as a 1-d vector or string, if rp is not None')
+                #end if
+            #end if
+            if isinstance(rp,(int,float)):
+                if units=="radians" or units=="rad":
+                    theta = float(rp)
+                else:
+                    theta = float(rp)*np.pi/180.0
+                c = np.cos(theta)
+                s = np.sin(theta)
+            else:
+                if isinstance(rp,str):
+                    if rp[0]=='a': # rp= 'a0', 'a1', or 'a2'
+                        rp = self.axes[int(rp[1])]
+                    else: # rp= 'x', 'y', or 'z'
+                        rp = dirmap[rp]
+                    #end if
+                else:
+                    rp = array(rp,dtype=float)
+                #end if
+                # go from r,rp to r,theta
+                c = np.dot(r,rp)/np.linalg.norm(r)/np.linalg.norm(rp)
+                if abs(c-1)<1e-6:
+                    s = 0.0
+                    r = np.array([1,0,0])
+                else:
+                    s = np.dot(np.cross(r,rp),np.cross(r,rp))/np.linalg.norm(r)/np.linalg.norm(rp)/np.linalg.norm(np.cross(r,rp)) 
+                    r = np.cross(r,rp)/np.linalg.norm(np.cross(r,rp))
+            #end if
+            # make R from r,theta
+            R = [[     c+r[0]**2.0*(1.0-c), r[0]*r[1]*(1.0-c)-r[2]*s, r[0]*r[2]*(1.0-c)+r[1]*s],
+                 [r[1]*r[0]*(1.0-c)+r[2]*s,      c+r[1]**2.0*(1.0-c), r[1]*r[2]*(1.0-c)-r[0]*s],
+                 [r[2]*r[0]*(1.0-c)-r[1]*s, r[2]*r[1]*(1.0-c)+r[0]*s,      c+r[2]**2.0*(1.0-c)]]
+        else:
+            R = r
+        #end if
+        R = array(R,dtype=float)
+        if passive:
+            R = R.T
+        if check:
+            if not np.allclose(dot(R,R.T),identity(len(R))):
+                self.error('the function, rotate, must be given an orthogonal matrix')
+            #end if
+        #end if
+        self.matrix_transform(R)
+    #end def rotate
+
+
+    # test needed
+    def matrix_transform(self,A): 
+        """
+        Arbitrary transformation matrix (column-major).
+
+        Parameters
+        ----------
+        skew  : `array_like, float, shape (3,3)`
+            Transform the structure using the matrix skew. It is assumed that
+            skew is in column-major form, i.e., it transforms a vector v as
+            v' = Tv
+        """
+        A = A.T
         axinv  = inv(self.axes)
-        axnew  = dot(self.axes,skew)
+        axnew  = dot(self.axes,A)
         kaxinv = inv(self.kaxes)
-        kaxnew = dot(self.kaxes,inv(skew).T)
+        kaxnew = dot(self.kaxes,inv(A).T)
         self.pos     = dot(dot(self.pos,axinv),axnew)
         self.center  = dot(dot(self.center,axinv),axnew)
         self.kpoints = dot(dot(self.kpoints,kaxinv),kaxnew)
         self.axes  = axnew
         self.kaxes = kaxnew
         if self.folded_structure!=None:
-            self.folded_structure.skew(skew)
+            self.folded_structure.matrix_transform(A.T)
         #end if
+    #end def matrix_transform
+
+
+    # test needed
+    def skew(self,skew):
+        """
+        Arbitrary transformation matrix (row-major).
+
+        Parameters
+        ----------
+        skew  : `array_like, float, shape (3,3)`
+            Transform the structure using the matrix skew. It is assumed that
+            skew is in row-major form, i.e., it transforms a vector v as
+            v' = vT
+        """
+        self.matrix_transform(skew.T)
     #end def skew
         
     
diff --git a/nexus/lib/testing.py b/nexus/lib/testing.py
index cecb96e3d..e36743c24 100644
--- a/nexus/lib/testing.py
+++ b/nexus/lib/testing.py
@@ -277,3 +277,20 @@ def restore_nexus_log():
     logging_storage.clear()
     assert(len(logging_storage)==0)
 #end def restore_nexus_log
+
+
+# declare test failure
+#   useful inside try/except blocks
+def failed(msg='Test failed.'):
+    assert False,msg
+#end def failed
+
+
+class TestFailed(Exception):
+    None
+#end class TestFailed
+
+
+global_data = dict(
+    job_ref_table = False,
+    )
diff --git a/nexus/tests/unit/test_machines.py b/nexus/tests/unit/test_machines.py
new file mode 100644
index 000000000..93461acbf
--- /dev/null
+++ b/nexus/tests/unit/test_machines.py
@@ -0,0 +1,1815 @@
+
+import testing
+from testing import value_eq,object_eq,failed,TestFailed
+from testing import divert_nexus_log,restore_nexus_log
+
+
+all_machines = []
+machines_data = dict()
+
+
+def get_machine_data():
+    from generic import obj
+    from machines import Machine,Workstation,Supercomputer
+    if len(machines_data)==0:
+        workstations   = obj()
+        supercomputers = obj()
+        for machine in Machine.machines:
+            if isinstance(machine,Workstation):
+                workstations.append(machine)
+            elif isinstance(machine,Supercomputer):
+                supercomputers[machine.name] = machine
+            else:
+                failed()
+            #end if
+        #end for
+        machines_data['ws'] = workstations
+        machines_data['sc'] = supercomputers
+    #end if
+    ws = machines_data['ws']
+    sc = machines_data['sc']
+    return ws,sc
+#end def get_machine_data
+
+
+def get_all_machines():
+    from machines import Machine
+    if len(all_machines)==0:
+        for m in Machine.machines:
+            all_machines.append(m)
+        #end for
+    #end if
+    return all_machines
+#end def get_all_machines
+
+
+def get_supercomputers():
+    ws,sc = get_machine_data()
+    return sc
+#end def get_supercomputers
+
+
+
+
+def test_import():
+    from machines import Machine,Workstation,InteractiveCluster,Supercomputer
+    from machines import Job,job
+    from machines import get_machine,get_machine_name
+#end def test_import
+
+
+
+def test_cpu_count():
+    from machines import cpu_count
+    assert(isinstance(cpu_count(),int))
+#end def test_cpu_count
+
+
+
+def test_options():
+    from generic import obj
+    from machines import Options
+
+    # empty init
+    o = Options()
+    assert(len(o)==0)
+
+    # filled init
+    inputs = dict(
+        n = '-n 1',
+        d = '-d 2',
+        exe = '--exe',
+        )
+    oi = Options(**inputs)
+    assert(oi.to_dict()==inputs)
+
+    # add
+    oa = Options()
+    oa.add(**inputs)
+    assert(object_eq(oa,oi))
+
+    # read
+    opts = '-a -b 1 -c=2 --dname="0 1 2" --evar -fval -gval other --hval'
+    ref = obj(
+        a     = '-a',
+        b     = '-b 1',
+        c     = '-c=2',
+        dname = '--dname="0 1 2"',
+        evar  = '--evar',
+        fval  = '-fval',
+        gval  = '-gval other',
+        hval  = '--hval',
+        )
+    o = Options()
+    o.read(opts)
+    assert(object_eq(o.to_obj(),ref))
+
+    # write
+    opts_write = o.write()
+    o2 = Options()
+    o2.read(opts_write)
+    assert(object_eq(o2.to_obj(),ref))
+#end def test_options
+
+
+
+def test_job_init():
+    from machines import Job,job
+    from machines import job_defaults
+    from machines import job_defaults_assign,job_defaults_nonassign
+
+    jda  = set(job_defaults_assign.keys())
+    jdna = set(job_defaults_nonassign.keys())
+    jd   = set(job_defaults.keys())
+    assert('nodes' in jda)
+    assert('app' in jdna)
+    assert(len(jda&jdna)==0)
+    assert(jd==(jda|jdna))
+
+    assert(id(job)==id(Job))
+
+    # empty init should fail w/o implicit or explicit machine
+    try:
+        job()
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+
+    # empty init should succeed if machine is bypassed
+    j = job(skip_machine=True)
+
+    assert(len(jda-set(j.keys()))==0)
+    assert(len(j.app_props)==0)
+    j.app_props = None # set back to default
+    assert(object_eq(j.obj(*jda),job_defaults_assign))
+#end def test_job_init
+
+
+
+def test_job_time():
+    from generic import obj
+    from machines import job,Job
+
+    def check_time(j,d,h,m,s,td,th,tm,ts):
+        t = obj(days=d,hours=h,minutes=m,seconds=s)
+        assert(object_eq(j.get_time(),t))
+        assert(j.total_seconds()==ts)
+        assert(j.total_minutes()==tm)
+        assert(j.total_hours()==th)
+        assert(j.total_days()==td)
+    #end def check_time
+
+    sm = dict(skip_machine=True)
+
+    j1 = job(**sm)
+    j2 = job(hours=4,minutes=30,**sm)
+    j3 = job(days=2,hours=8,**sm)
+    j4 = job(hours=53,minutes=127,**sm)
+
+    check_time(j1, 0, 0,  0, 0, 0,  0,    0,      0)
+    check_time(j2, 0, 4, 30, 0, 0,  4,  270,  16200)
+    check_time(j3, 2, 8,  0, 0, 2, 56, 3360, 201600)
+    check_time(j4, 2, 7,  7, 0, 2, 55, 3307, 198420)
+
+    # pbs walltime
+    assert(j1.pbs_walltime()=='00:00:00')
+    assert(j2.pbs_walltime()=='04:30:00')
+    assert(j3.pbs_walltime()=='2:08:00:00')
+    assert(j4.pbs_walltime()=='2:07:07:00')
+
+    # sbatch walltime
+    assert(j1.sbatch_walltime()=='00:00:00')
+    assert(j2.sbatch_walltime()=='04:30:00')
+    assert(j3.sbatch_walltime()=='56:00:00')
+    assert(j4.sbatch_walltime()=='55:07:00')
+
+    # ll walltime
+    assert(j1.ll_walltime()=='00:00:00')
+    assert(j2.ll_walltime()=='04:30:00')
+    assert(j3.ll_walltime()=='56:00:00')
+    assert(j4.ll_walltime()=='55:07:00')
+
+    # lsf walsftime
+    assert(j1.lsf_walltime()=='00:00')
+    assert(j2.lsf_walltime()=='04:30')
+    assert(j3.lsf_walltime()=='56:00')
+    assert(j4.lsf_walltime()=='55:07')
+
+    t = Job.zero_time()
+    t = j1.max_time(t)
+    t = j2.max_time(t)
+    t = j3.max_time(t)
+    t = j4.max_time(t)
+    assert(object_eq(t,j3.get_time()))
+
+#end def test_job_time
+
+
+
+def test_job_set_id():
+    from machines import job
+
+    j1 = job(skip_machine=True)
+    j2 = job(skip_machine=True)
+
+    j1.set_id()
+    j2.set_id()
+
+    assert(isinstance(j1.internal_id,int))
+    assert(isinstance(j2.internal_id,int))
+    assert(j2.internal_id-j1.internal_id==1)
+#end def test_job_set_id
+
+
+
+def test_job_get_machine():
+    from machines import job
+
+    machines = get_all_machines()
+    for m in machines:
+        j = job(machine=m.name,skip_machine=True)
+        mj = j.get_machine()
+        assert(id(mj)==id(m))
+    #end for
+
+#end def test_job_get_machine
+
+
+
+def test_job_set_environment():
+    from machines import job
+
+    workstations,supercomputers = get_machine_data()
+
+    machines = []
+    machines.append(workstations.first())
+    machines.append(supercomputers.first())
+
+    for m in machines:
+        j = job(machine=m.name,skip_machine=True)
+        j.set_environment(OMP_NUM_THREADS=1)
+        assert(j.env['OMP_NUM_THREADS']=='1')
+    #end for
+
+#end def test_job_set_environment
+
+
+
+def test_job_clone():
+    from machines import job
+
+    j1 = job(skip_machine=True)
+    j1.set_id()
+
+    j2 = j1.clone()
+    assert(id(j2)!=id(j1))
+    assert(j2.internal_id-j1.internal_id==1)
+    del j1.internal_id
+    del j2.internal_id
+    assert(object_eq(j2,j1))
+#end def_test_job_clone
+
+
+
+def test_job_serial_clone():
+    from machines import job
+
+    j1 = job(skip_machine=True)
+    assert(not j1.serial)
+
+    j2 = j1.serial_clone()
+    assert(j2.serial)
+    assert(j2.cores==1)
+    assert(id(j2)!=id(j1))
+    keys = 'serial cores init_info'.split()
+    j1.delete(keys)
+    j2.delete(keys)
+    assert(object_eq(j2,j1))
+#end def test_job_serial_clone
+
+
+
+def test_machine_virtuals():
+    from machines import Machine
+    arg0 = None
+    arg1 = None
+    try:
+        Machine.query_queue(arg0)
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+    try:
+        Machine.submit_jobs(arg0)
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+    try:
+        Machine.process_job(arg0,arg1)
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+    try:
+        Machine.process_job_options(arg0,arg1)
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+    try:
+        Machine.write_job(arg0,arg1,file=False)
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+    try:
+        Machine.submit_job(arg0,arg1)
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+#end def test_machine_virtuals
+
+
+
+def test_machine_list():
+    from machines import Machine
+
+    assert(len(Machine.machines)>0)
+    for m in Machine.machines:
+        assert(isinstance(m,Machine))
+        exists = m.name in Machine.machines
+        assert(exists)
+        assert(Machine.exists(m.name))
+        assert(Machine.is_unique(m))
+        m.validate()
+    #end for
+#end def test_machine_list
+
+
+
+def test_machine_add():
+    from machines import Machine
+    mtest = Machine.machines.first()
+    assert(isinstance(mtest,Machine))
+    try:
+        Machine.add(mtest)
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+    try:
+        Machine.add('my_machine')
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+#end def test_machine_add
+
+
+
+def test_machine_get():
+    from machines import Machine
+    mtest = Machine.machines.first()
+    assert(isinstance(mtest,Machine))
+
+    m = Machine.get(mtest.name)
+    assert(isinstance(m,Machine))
+    assert(id(m)==id(mtest))
+    try:
+        Machine.get(m)
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+    try:
+        Machine.get('some_nonexistant_machine')
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+#end def test_machine_get
+
+
+
+def test_machine_instantiation():
+    from machines import Machine
+    # test guards against empty/invalid instantiation
+    try:
+        Machine()
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+    try:
+        Machine(123)
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+    # test creation of a new machine
+    test_name = 'test_machine'
+    assert(not Machine.exists(test_name))
+    m = Machine(name=test_name)
+    assert(isinstance(m,Machine))
+    assert(Machine.exists(m.name))
+    assert(Machine.is_unique(m))
+    m.validate()
+
+    # test guards against multiple instantiation
+    try:
+        Machine(name=test_name)
+        raise TestFailed
+    except TestFailed:
+        failed()
+    except:
+        None
+    #end try
+
+    # remove test machine
+    del Machine.machines.test_machine
+    assert(not Machine.exists(test_name))
+#end def test_machine_instantiation
+
+
+
+def test_workstation_init():
+    from generic import obj
+    from machines import Workstation
+
+    ws = Workstation('wsi',16,'mpirun')
+
+    refws = obj(
+        account         = None,
+        app_directories = None,
+        app_directory   = None,
+        app_launcher    = 'mpirun',
+        cores           = 16,
+        finished        = set([]),
+        local_directory = None,
+        name            = 'wsi',
+        process_granularity = 1,
+        queue_size      = 16,
+        running         = set([]),
+        user            = None,
+        waiting         = set([]),
+        jobs            = obj(),
+        processes       = obj(),
+        )
+
+    assert(object_eq(ws.to_obj(),refws))
+
+#end def test_workstation_init
+
+
+
+# imitate Job.initialize w/o involving simulation object
+def init_job(j,
+             id  = 'job_ident',
+             dir = './',
+             ):
+    import os
+    identifier = id
+    directory  = dir
+    j.set_id()
+    j.identifier  = identifier
+    j.directory   = directory
+    j.abs_dir     = os.path.abspath(directory)
+    j.subdir      = j.directory
+    j.abs_subdir  = j.abs_dir
+    j.name        = identifier
+    j.simid       = j.internal_id+1000
+    j.outfile     = identifier+'.out'
+    j.errfile     = identifier+'.err'
+    j.app_command = 'echo run'
+    j.process()
+#end def init_job
+
+
+
+def test_workstation_scheduling():
+    import time
+    from machines import Workstation
+    from machines import job,Job
+
+    tpath = testing.setup_unit_test_output_directory('machines','test_workstation_scheduling')
+
+    # create workstation for testing
+    ws = Workstation('wss',16,'mpirun')
+
+
+    # test process_job(), process_job_options(), write_job()
+    j = job(machine=ws.name,cores=4,threads=2,local=True)
+    assert(j.machine==ws.name)
+    assert(j.cores==4)
+    assert(j.threads==2)
+    assert(j.local)
+    assert(j.processes==2)
+    assert(j.run_options.np=='-np 2')
+    assert(j.batch_mode==False)
+    init_job(j,dir=tpath) # imitate interaction w/ simulation object
+    assert(ws.write_job(j)=='export OMP_NUM_THREADS=2\nmpirun -np 2 echo run')
+
+    j = job(machine=ws.name,serial=True)
+    assert(j.machine==ws.name)
+    assert(j.cores==1)
+    assert(j.threads==1)
+    assert(j.serial)
+    assert(j.processes==1)
+    assert(j.run_options.np=='-np 1')
+    assert(j.batch_mode==False)
+    init_job(j,dir=tpath) # imitate interaction w/ simulation object
+    assert(ws.write_job(j)=='export OMP_NUM_THREADS=1\necho run')
+
+
+    # test add_job()
+    assert(j.status==Job.states.none)
+    assert(not j.submitted)
+    assert(len(ws.waiting)==0)
+
+    j.submit() # calls ws.add_job(j)
+
+    assert(j.status==Job.states.waiting)
+    assert(j.submitted)
+    assert(ws.waiting==set([j.internal_id]))
+    assert(set(ws.jobs.keys())==set([j.internal_id]))
+    assert(id(ws.jobs[j.internal_id])==id(j))
+
+
+    # test submit_jobs() and submit_job()
+    divert_nexus_log()
+    assert(j.system_id is None)
+    assert(len(ws.running)==0)
+    assert(len(ws.processes)==0)
+
+    ws.submit_jobs() # will call system echo
+
+    assert(j.status==Job.states.running)
+    assert(isinstance(j.system_id,int))
+    assert(len(ws.waiting)==0)
+    assert(ws.running==set([j.internal_id]))
+    assert(set(ws.processes.keys())==set([j.system_id]))
+    p = ws.processes[j.system_id]
+    assert(p.popen.pid==j.system_id)
+    assert(id(p.job)==id(j))
+    assert(set(ws.jobs.keys())==set([j.internal_id]))
+    restore_nexus_log()
+
+    # allow a moment for all system calls to resolve
+    time.sleep(0.1)
+
+    # test query_queue()
+    assert(not j.finished)
+
+    ws.query_queue() # echo process complete
+
+    assert(j.finished)
+    assert(j.status==Job.states.finished)
+    assert(len(ws.running)==0)
+    assert(len(ws.processes)==0)
+    assert(ws.finished==set([j.internal_id]))
+    assert(set(ws.jobs.keys())==set([j.internal_id]))
+
+#end def test_workstation_scheduling
+
+
+
+def test_supercomputer_init():
+    from generic import obj
+    from machines import Theta
+
+    class ThetaInit(Theta):
+        name = 'theta_init'
+    #end class ThetaInit
+
+    sc = ThetaInit(4392,1,64,192,1000,'aprun','qsub','qstata','qdel')
+
+    refsc = obj(
+        account         = None,
+        app_directories = None,
+        app_directory   = None,
+        app_launcher    = 'aprun',
+        cores           = 281088,
+        cores_per_node  = 64,
+        cores_per_proc  = 64,
+        finished        = set([]),
+        job_remover     = 'qdel',
+        local_directory = None,
+        name            = 'theta_init',
+        nodes           = 4392,
+        procs           = 4392,
+        procs_per_node  = 1,
+        queue_querier   = 'qstata',
+        queue_size      = 1000,
+        ram             = 843264,
+        ram_per_node    = 192,
+        running         = set([]),
+        sub_launcher    = 'qsub',
+        user            = None,
+        waiting         = set([]),
+        jobs            = obj(),
+        processes       = obj(),
+        system_queue    = obj(),
+        )
+
+    assert(object_eq(sc.to_obj(),refsc))
+
+#end def test_supercomputer_init
+
+
+
+def test_supercomputer_scheduling():
+    import os
+    import time
+    from generic import obj
+    from machines import Theta
+    from machines import job,Job
+
+    tpath = testing.setup_unit_test_output_directory('machines','test_supercomputer_scheduling')
+
+    # create supercomputer for testing
+    class ThetaSched(Theta):
+        name = 'theta_sched'
+    #end class ThetaSched
+
+    sc = ThetaSched(4392,1,64,192,1000,'aprun','echo','test_query','qdel')
+
+
+    # test process_job() and process_job_options()
+    j = job(machine=sc.name,nodes=2,threads=8,hours=3,minutes=30,account='ABC123')
+    assert(j.machine==sc.name)
+    assert(j.nodes==2)
+    assert(j.threads==8)
+    assert(j.processes==16)
+    assert(j.processes_per_node==8)
+    assert(j.cores==128)
+    assert(j.hours==3)
+    assert(j.minutes==30)
+    assert(j.account=='ABC123')
+    refro = obj(
+        N               = '-N 8',
+        cc              = '-cc depth',
+        d               = '-d 8',
+        e               = '-e OMP_NUM_THREADS=8',
+        j               = '-j 1',
+        n               = '-n 16',
+        )
+    assert(object_eq(j.run_options.to_obj(),refro))
+    assert(j.batch_mode==True)
+
+
+    # test write_job()
+    init_job(j,id='123',dir=tpath) # imitate interaction w/ simulation object
+    ref_wj = '''#!/bin/bash
+#COBALT -q default
+#COBALT -A ABC123
+#COBALT -n 2
+#COBALT -t 210
+#COBALT -O 123
+#COBALT --attrs mcdram=cache:numa=quad
+
+export OMP_NUM_THREADS=8
+aprun -e OMP_NUM_THREADS=8 -d 8 -cc depth -j 1 -n 16 -N 8 echo run'''
+    wj = sc.write_job(j)
+    assert(wj.strip()==ref_wj)
+
+
+    # test add_job()
+    assert(j.status==Job.states.none)
+    assert(not j.submitted)
+    assert(len(sc.waiting)==0)
+
+    j.submit() # calls sc.add_job(j)
+
+    assert(j.status==Job.states.waiting)
+    assert(j.submitted)
+    assert(sc.waiting==set([j.internal_id]))
+    assert(set(sc.jobs.keys())==set([j.internal_id]))
+    assert(id(sc.jobs[j.internal_id])==id(j))
+
+    
+    # test write_job() to file
+    sc.write_job(j,file=True)
+
+    subfile_path = os.path.join(tpath,j.subfile)
+    assert(os.path.exists(subfile_path))
+    assert(open(subfile_path,'r').read().strip()==ref_wj)
+
+
+    # test sub_command()
+    assert(sc.sub_command(j)=='echo 123.echo.in')
+
+
+    # test submit_jobs() and submit_job()
+    divert_nexus_log()
+    assert(j.system_id is None)
+
+    sc.submit_jobs() # will call system echo
+
+    assert(j.status==Job.states.running)
+    assert(j.system_id==123)
+    assert(len(sc.waiting)==0)
+    assert(sc.running==set([j.internal_id]))
+    assert(set(sc.processes.keys())==set([123]))
+    assert(set(sc.jobs.keys())==set([j.internal_id]))
+    restore_nexus_log()
+
+
+    # allow a moment for all system calls to resolve
+    time.sleep(0.1)
+
+
+    # test query_queue()
+    assert(not j.finished)
+
+    sc.query_queue() # echo process complete
+
+    assert(j.finished)
+    assert(j.status==Job.states.finished)
+    assert(len(sc.running)==0)
+    assert(len(sc.processes)==0)
+    assert(sc.finished==set([j.internal_id]))
+    assert(set(sc.jobs.keys())==set([j.internal_id]))
+
+#end def test_supercomputer_scheduling
+
+
+
+def test_process_job():
+    from random import randint
+    from generic import obj
+    from machines import Machine,Job
+
+    nw  = 5
+    nwj = 5
+    nsj = 5
+    nij = 5
+
+    workstations,supercomputers = get_machine_data()
+
+    allow_warn = Machine.allow_warnings
+    Machine.allow_warnings = False
+
+    not_idempotent = obj()
+
+    # check workstations
+    nworkstations = nw
+    if nworkstations is None:
+        nworkstations=len(workstations)
+    #end if
+    nworkstations = min(nworkstations,len(workstations))
+    njobs = nwj
+    for nm in range(nworkstations):
+        if nworkstations<len(workstations):
+            machine = workstations.select_random() # select machine at random
+        else:
+            machine = workstations[nm]
+        #end if
+        cores_min     = 1
+        cores_max     = machine.cores
+        processes_min = 1
+        processes_max = machine.cores
+        threads_min   = 1
+        threads_max   = machine.cores
+        job_inputs = []
+        job_inputs_base = []
+        for nj in range(njobs): # vary cores
+            cores   = randint(cores_min,cores_max)
+            threads = randint(threads_min,threads_max)
+            job_inputs_base.append(obj(cores=cores,threads=threads))
+        #end for
+        for nj in range(njobs): # vary processes
+            processes   = randint(processes_min,processes_max)
+            threads = randint(threads_min,threads_max)
+            job_inputs_base.append(obj(processes=processes,threads=threads))
+        #end for
+        job_inputs.extend(job_inputs_base)
+        for job_input in job_inputs_base: # run in serial
+            ji = job_input.copy()
+            ji.serial = True
+            job_inputs.append(ji)
+        #end for
+        # perform idempotency test
+        machine_idempotent = True
+        for job_input in job_inputs:
+            job = Job(machine=machine.name,**job_input)
+            job2 = obj.copy(job)
+            machine.process_job(job2)
+            machine_idempotent &= job==job2
+        #end for
+        if not machine_idempotent:
+            not_idempotent[machine.name] = machine
+        #end if
+    #end for
+
+    # check supercomputers
+    njobs = nsj
+    small_node_ceiling = 20
+    nodes_min   = 1
+    cores_min   = 1
+    threads_min = 1
+    shared_job_inputs = obj(name='some_job',account='some_account')
+    for machine in supercomputers:
+        job_inputs = []
+        job_inputs_base = []
+        threads_max = 2*machine.cores_per_node
+        # sample small number of nodes more heavily
+        nodes_max   = min(small_node_ceiling,machine.nodes)
+        cores_max   = min(small_node_ceiling*machine.cores_per_node,machine.cores)
+        for nj in range(njobs): # nodes alone
+            nodes   = randint(nodes_min,nodes_max)
+            threads = randint(threads_min,threads_max)
+            job_input = obj(nodes=nodes,threads=threads,**shared_job_inputs)
+            job_inputs_base.append(job_input)
+        #end for
+        for nj in range(njobs): # cores alone
+            cores   = randint(cores_min,cores_max)
+            threads = randint(threads_min,threads_max)
+            job_input = obj(cores=cores,threads=threads,**shared_job_inputs)
+            job_inputs_base.append(job_input)
+        #end for
+        for nj in range(njobs): # nodes and cores
+            nodes   = randint(nodes_min,nodes_max)
+            cores   = randint(cores_min,cores_max)
+            threads = randint(threads_min,threads_max)
+            job_input = obj(nodes=nodes,cores=cores,threads=threads,**shared_job_inputs)
+            job_inputs_base.append(job_input)
+        #end for
+        # sample full node set
+        nodes_max = machine.nodes
+        cores_max = machine.cores
+        for nj in range(njobs): # nodes alone
+            nodes   = randint(nodes_min,nodes_max)
+            threads = randint(threads_min,threads_max)
+            job_input = obj(nodes=nodes,threads=threads,**shared_job_inputs)
+            job_inputs_base.append(job_input)
+        #end for
+        for nj in range(njobs): # cores alone
+            cores   = randint(cores_min,cores_max)
+            threads = randint(threads_min,threads_max)
+            job_input = obj(cores=cores,threads=threads,**shared_job_inputs)
+            job_inputs_base.append(job_input)
+        #end for
+        for nj in range(njobs): # nodes and cores
+            nodes   = randint(nodes_min,nodes_max)
+            cores   = randint(cores_min,cores_max)
+            threads = randint(threads_min,threads_max)
+            job_input = obj(nodes=nodes,cores=cores,threads=threads,**shared_job_inputs)
+            job_inputs_base.append(job_input)
+        #end for
+        job_inputs.extend(job_inputs_base)
+        # now add serial jobs
+        for job_input in job_inputs_base:
+            ji = job_input.copy()
+            ji.serial = True
+            job_inputs.append(ji)
+        #end for
+        # now add local, serial jobs
+        for job_input in job_inputs_base:
+            ji = job_input.copy()
+            ji.serial = True
+            ji.local  = True
+            job_inputs.append(ji)
+        #end for
+        # perform idempotency test
+        machine_idempotent = True
+        for job_input in job_inputs:
+            job = Job(machine=machine.name,**job_input)
+            job2 = obj.copy(job)
+            machine.process_job(job2)
+            job_idempotent = object_eq(job,job2)
+            if not job_idempotent:
+                d,d1,d2 = object_diff(job,job2,full=True)
+                change = obj(job_before=obj(d1),job_after=obj(d2))
+                msg = machine.name+'\n'+str(change)
+                failed(msg)
+            #end if
+            machine_idempotent &= job_idempotent
+        #end for
+        if not machine_idempotent:
+            not_idempotent[machine.name] = machine
+        #end if
+    #end for
+
+    if len(not_idempotent)>0:
+        mlist = ''
+        for name in sorted(not_idempotent.keys()):
+            mlist+= '\n  '+name
+        #end for
+        msg='\n\nsome machines failed process_job idempotency test:{0}'.format(mlist)
+        failed(msg)
+    #end if
+    Machine.allow_warnings = allow_warn
+
+#end def test_process_job
+
+
+
+def test_job_run_command():
+    from generic import obj
+    from machines import Machine,Job
+
+    workstations,supercomputers = get_machine_data()
+
+    allow_warn = Machine.allow_warnings
+    Machine.allow_warnings = False
+
+    def parse_job_command(command):
+        tokens = command.replace(':',' ').split()
+        launcher = tokens[0]
+        exe = tokens[-1]
+        args = []
+        options = obj()
+        last_option = None
+        for t in tokens[1:-1]:
+            if t.startswith('-'):
+                options[t] = None
+                last_option = t
+            elif last_option is not None:
+                options[last_option] = t
+                last_option = None
+            else:
+                args.append(t)
+            #end if
+        #end for
+        jc = obj(
+            launcher   = launcher,
+            executable = exe,
+            args       = args,
+            options    = options,
+            )
+        return jc
+    #end def parse_job_command
+
+    def job_commands_equal(c1,c2):
+        jc1 = parse_job_command(c1)
+        jc2 = parse_job_command(c2)
+        return object_eq(jc1,jc2)
+    #end def job_command_equal
+
+    job_run_ref = obj({
+        ('amos'           , 'n1'            ) : 'srun test.x',
+        ('amos'           , 'n1_p1'         ) : 'srun test.x',
+        ('amos'           , 'n2'            ) : 'srun test.x',
+        ('amos'           , 'n2_t2'         ) : 'srun test.x',
+        ('amos'           , 'n2_t2_e'       ) : 'srun test.x',
+        ('amos'           , 'n2_t2_p2'      ) : 'srun test.x',
+        ('bluewaters_xe'  , 'n1'            ) : 'aprun -n 32 test.x',
+        ('bluewaters_xe'  , 'n1_p1'         ) : 'aprun -n 1 test.x',
+        ('bluewaters_xe'  , 'n2'            ) : 'aprun -n 64 test.x',
+        ('bluewaters_xe'  , 'n2_t2'         ) : 'aprun -d 2 -n 32 test.x',
+        ('bluewaters_xe'  , 'n2_t2_e'       ) : 'aprun -d 2 -n 32 test.x',
+        ('bluewaters_xe'  , 'n2_t2_p2'      ) : 'aprun -d 2 -n 4 test.x',
+        ('bluewaters_xk'  , 'n1'            ) : 'aprun -n 16 test.x',
+        ('bluewaters_xk'  , 'n1_p1'         ) : 'aprun -n 1 test.x',
+        ('bluewaters_xk'  , 'n2'            ) : 'aprun -n 32 test.x',
+        ('bluewaters_xk'  , 'n2_t2'         ) : 'aprun -d 2 -n 16 test.x',
+        ('bluewaters_xk'  , 'n2_t2_e'       ) : 'aprun -d 2 -n 16 test.x',
+        ('bluewaters_xk'  , 'n2_t2_p2'      ) : 'aprun -d 2 -n 4 test.x',
+        ('cades'          , 'n1'            ) : 'mpirun -np 36 test.x',
+        ('cades'          , 'n1_p1'         ) : 'mpirun -np 1 test.x',
+        ('cades'          , 'n2'            ) : 'mpirun -np 72 test.x',
+        ('cades'          , 'n2_t2'         ) : 'mpirun -np 36 --npersocket 9 test.x',
+        ('cades'          , 'n2_t2_e'       ) : 'mpirun -np 36 --npersocket 9 test.x',
+        ('cades'          , 'n2_t2_p2'      ) : 'mpirun -np 4 --npersocket 1 test.x',
+        ('cetus'          , 'n1'            ) : 'runjob --np 16 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
+        ('cetus'          , 'n1_p1'         ) : 'runjob --np 1 -p 1 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
+        ('cetus'          , 'n2'            ) : 'runjob --np 32 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
+        ('cetus'          , 'n2_t2'         ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
+        ('cetus'          , 'n2_t2_e'       ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 ENV_VAR=1 : test.x',
+        ('cetus'          , 'n2_t2_p2'      ) : 'runjob --np 4 -p 2 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
+        ('chama'          , 'n1'            ) : 'srun test.x',
+        ('chama'          , 'n1_p1'         ) : 'srun test.x',
+        ('chama'          , 'n2'            ) : 'srun test.x',
+        ('chama'          , 'n2_t2'         ) : 'srun test.x',
+        ('chama'          , 'n2_t2_e'       ) : 'srun test.x',
+        ('chama'          , 'n2_t2_p2'      ) : 'srun test.x',
+        ('cooley'         , 'n1'            ) : 'mpirun -np 12 test.x',
+        ('cooley'         , 'n1_p1'         ) : 'mpirun -np 1 test.x',
+        ('cooley'         , 'n2'            ) : 'mpirun -np 24 test.x',
+        ('cooley'         , 'n2_t2'         ) : 'mpirun -np 12 test.x',
+        ('cooley'         , 'n2_t2_e'       ) : 'mpirun -np 12 test.x',
+        ('cooley'         , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
+        ('cori'           , 'n1'            ) : 'srun test.x',
+        ('cori'           , 'n1_p1'         ) : 'srun test.x',
+        ('cori'           , 'n2'            ) : 'srun test.x',
+        ('cori'           , 'n2_t2'         ) : 'srun test.x',
+        ('cori'           , 'n2_t2_e'       ) : 'srun test.x',
+        ('cori'           , 'n2_t2_p2'      ) : 'srun test.x',
+        ('edison'         , 'n1'            ) : 'srun test.x',
+        ('edison'         , 'n1_p1'         ) : 'srun test.x',
+        ('edison'         , 'n2'            ) : 'srun test.x',
+        ('edison'         , 'n2_t2'         ) : 'srun test.x',
+        ('edison'         , 'n2_t2_e'       ) : 'srun test.x',
+        ('edison'         , 'n2_t2_p2'      ) : 'srun test.x',
+        ('eos'            , 'n1'            ) : 'aprun -n 16 test.x',
+        ('eos'            , 'n1_p1'         ) : 'aprun -n 1 test.x',
+        ('eos'            , 'n2'            ) : 'aprun -n 32 test.x',
+        ('eos'            , 'n2_t2'         ) : 'aprun -ss -cc numa_node -d 2 -n 16 test.x',
+        ('eos'            , 'n2_t2_e'       ) : 'aprun -ss -cc numa_node -d 2 -n 16 test.x',
+        ('eos'            , 'n2_t2_p2'      ) : 'aprun -ss -cc numa_node -d 2 -n 4 test.x',
+        ('jaguar'         , 'n1'            ) : 'aprun -n 16 test.x',
+        ('jaguar'         , 'n1_p1'         ) : 'aprun -n 1 test.x',
+        ('jaguar'         , 'n2'            ) : 'aprun -n 32 test.x',
+        ('jaguar'         , 'n2_t2'         ) : 'aprun -d 2 -n 16 test.x',
+        ('jaguar'         , 'n2_t2_e'       ) : 'aprun -d 2 -n 16 test.x',
+        ('jaguar'         , 'n2_t2_p2'      ) : 'aprun -d 2 -n 4 test.x',
+        ('komodo'         , 'n1'            ) : 'mpirun -np 12 test.x',
+        ('komodo'         , 'n1_p1'         ) : 'mpirun -np 1 test.x',
+        ('komodo'         , 'n2'            ) : 'mpirun -np 24 test.x',
+        ('komodo'         , 'n2_t2'         ) : 'mpirun -np 12 test.x',
+        ('komodo'         , 'n2_t2_e'       ) : 'mpirun -np 12 test.x',
+        ('komodo'         , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
+        ('kraken'         , 'n1'            ) : 'aprun -n 12 test.x',
+        ('kraken'         , 'n1_p1'         ) : 'aprun -n 1 test.x',
+        ('kraken'         , 'n2'            ) : 'aprun -n 24 test.x',
+        ('kraken'         , 'n2_t2'         ) : 'aprun -d 2 -n 12 test.x',
+        ('kraken'         , 'n2_t2_e'       ) : 'aprun -d 2 -n 12 test.x',
+        ('kraken'         , 'n2_t2_p2'      ) : 'aprun -d 2 -n 4 test.x',
+        ('lonestar'       , 'n1'            ) : 'ibrun -n 12 -o 0 test.x',
+        ('lonestar'       , 'n1_p1'         ) : 'ibrun -n 1 -o 0 test.x',
+        ('lonestar'       , 'n2'            ) : 'ibrun -n 24 -o 0 test.x',
+        ('lonestar'       , 'n2_t2'         ) : 'ibrun -n 12 -o 0 test.x',
+        ('lonestar'       , 'n2_t2_e'       ) : 'ibrun -n 12 -o 0 test.x',
+        ('lonestar'       , 'n2_t2_p2'      ) : 'ibrun -n 4 -o 0 test.x',
+        ('matisse'        , 'n1'            ) : 'mpirun -np 16 test.x',
+        ('matisse'        , 'n1_p1'         ) : 'mpirun -np 1 test.x',
+        ('matisse'        , 'n2'            ) : 'mpirun -np 32 test.x',
+        ('matisse'        , 'n2_t2'         ) : 'mpirun -np 16 test.x',
+        ('matisse'        , 'n2_t2_e'       ) : 'mpirun -np 16 test.x',
+        ('matisse'        , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
+        ('mira'           , 'n1'            ) : 'runjob --np 16 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
+        ('mira'           , 'n1_p1'         ) : 'runjob --np 1 -p 1 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
+        ('mira'           , 'n2'            ) : 'runjob --np 32 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
+        ('mira'           , 'n2_t2'         ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
+        ('mira'           , 'n2_t2_e'       ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 ENV_VAR=1 : test.x',
+        ('mira'           , 'n2_t2_p2'      ) : 'runjob --np 4 -p 2 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
+        ('oic5'           , 'n1'            ) : 'mpirun -np 32 test.x',
+        ('oic5'           , 'n1_p1'         ) : 'mpirun -np 1 test.x',
+        ('oic5'           , 'n2'            ) : 'mpirun -np 64 test.x',
+        ('oic5'           , 'n2_t2'         ) : 'mpirun -np 32 test.x',
+        ('oic5'           , 'n2_t2_e'       ) : 'mpirun -np 32 test.x',
+        ('oic5'           , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
+        ('redsky'         , 'n1'            ) : 'srun test.x',
+        ('redsky'         , 'n1_p1'         ) : 'srun test.x',
+        ('redsky'         , 'n2'            ) : 'srun test.x',
+        ('redsky'         , 'n2_t2'         ) : 'srun test.x',
+        ('redsky'         , 'n2_t2_e'       ) : 'srun test.x',
+        ('redsky'         , 'n2_t2_p2'      ) : 'srun test.x',
+        ('serrano'        , 'n1'            ) : 'srun test.x',
+        ('serrano'        , 'n1_p1'         ) : 'srun test.x',
+        ('serrano'        , 'n2'            ) : 'srun test.x',
+        ('serrano'        , 'n2_t2'         ) : 'srun test.x',
+        ('serrano'        , 'n2_t2_e'       ) : 'srun test.x',
+        ('serrano'        , 'n2_t2_p2'      ) : 'srun test.x',
+        ('skybridge'      , 'n1'            ) : 'srun test.x',
+        ('skybridge'      , 'n1_p1'         ) : 'srun test.x',
+        ('skybridge'      , 'n2'            ) : 'srun test.x',
+        ('skybridge'      , 'n2_t2'         ) : 'srun test.x',
+        ('skybridge'      , 'n2_t2_e'       ) : 'srun test.x',
+        ('skybridge'      , 'n2_t2_p2'      ) : 'srun test.x',
+        ('solo'           , 'n1'            ) : 'srun test.x',
+        ('solo'           , 'n1_p1'         ) : 'srun test.x',
+        ('solo'           , 'n2'            ) : 'srun test.x',
+        ('solo'           , 'n2_t2'         ) : 'srun test.x',
+        ('solo'           , 'n2_t2_e'       ) : 'srun test.x',
+        ('solo'           , 'n2_t2_p2'      ) : 'srun test.x',
+        ('stampede2'      , 'n1'            ) : 'ibrun -n 68 -o 0 test.x',
+        ('stampede2'      , 'n1_p1'         ) : 'ibrun -n 1 -o 0 test.x',
+        ('stampede2'      , 'n2'            ) : 'ibrun -n 136 -o 0 test.x',
+        ('stampede2'      , 'n2_t2'         ) : 'ibrun -n 68 -o 0 test.x',
+        ('stampede2'      , 'n2_t2_e'       ) : 'ibrun -n 68 -o 0 test.x',
+        ('stampede2'      , 'n2_t2_p2'      ) : 'ibrun -n 4 -o 0 test.x',
+        ('summit'         , 'n1'            ) : 'jsrun -a 21 -r 2 -b rs -c 21 -d packed -n 2 -g 0 test.x',
+        ('summit'         , 'n1_g6'         ) : 'jsrun -a 7 -r 6 -b rs -c 7 -d packed -n 6 -g 1 test.x',
+        ('summit'         , 'n2'            ) : 'jsrun -a 21 -r 2 -b rs -c 21 -d packed -n 4 -g 0 test.x',
+        ('summit'         , 'n2_g6'         ) : 'jsrun -a 7 -r 6 -b rs -c 7 -d packed -n 12 -g 1 test.x',
+        ('summit'         , 'n2_t2'         ) : 'jsrun -a 10 -r 2 -b rs -c 20 -d packed -n 4 -g 0 test.x',
+        ('summit'         , 'n2_t2_e'       ) : 'jsrun -a 10 -r 2 -b rs -c 20 -d packed -n 4 -g 0 test.x',
+        ('summit'         , 'n2_t2_e_g6'    ) : 'jsrun -a 3 -r 6 -b rs -c 6 -d packed -n 12 -g 1 test.x',
+        ('summit'         , 'n2_t2_g6'      ) : 'jsrun -a 3 -r 6 -b rs -c 6 -d packed -n 12 -g 1 test.x',
+        ('supermuc'       , 'n1'            ) : 'mpiexec -n 28 test.x',
+        ('supermuc'       , 'n1_p1'         ) : 'mpiexec -n 1 test.x',
+        ('supermuc'       , 'n2'            ) : 'mpiexec -n 56 test.x',
+        ('supermuc'       , 'n2_t2'         ) : 'mpiexec -n 28 test.x',
+        ('supermuc'       , 'n2_t2_e'       ) : 'mpiexec -n 28 test.x',
+        ('supermuc'       , 'n2_t2_p2'      ) : 'mpiexec -n 4 test.x',
+        ('supermucng'     , 'n1'            ) : 'mpiexec -n 48 test.x',
+        ('supermucng'     , 'n1_p1'         ) : 'mpiexec -n 1 test.x',
+        ('supermucng'     , 'n2'            ) : 'mpiexec -n 96 test.x',
+        ('supermucng'     , 'n2_t2'         ) : 'mpiexec -n 48 test.x',
+        ('supermucng'     , 'n2_t2_e'       ) : 'mpiexec -n 48 test.x',
+        ('supermucng'     , 'n2_t2_p2'      ) : 'mpiexec -n 4 test.x',
+        ('taub'           , 'n1'            ) : 'mpirun -np 12 test.x',
+        ('taub'           , 'n1_p1'         ) : 'mpirun -np 1 test.x',
+        ('taub'           , 'n2'            ) : 'mpirun -np 24 test.x',
+        ('taub'           , 'n2_t2'         ) : 'mpirun -np 12 test.x',
+        ('taub'           , 'n2_t2_e'       ) : 'mpirun -np 12 test.x',
+        ('taub'           , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
+        ('theta'          , 'n1'            ) : 'aprun -e OMP_NUM_THREADS=1 -d 1 -cc depth -j 1 -n 64 -N 64 test.x',
+        ('theta'          , 'n1_p1'         ) : 'aprun -e OMP_NUM_THREADS=1 -d 1 -cc depth -j 1 -n 1 -N 1 test.x',
+        ('theta'          , 'n2'            ) : 'aprun -e OMP_NUM_THREADS=1 -d 1 -cc depth -j 1 -n 128 -N 64 test.x',
+        ('theta'          , 'n2_t2'         ) : 'aprun -e OMP_NUM_THREADS=2 -d 2 -cc depth -j 1 -n 64 -N 32 test.x',
+        ('theta'          , 'n2_t2_e'       ) : 'aprun -e OMP_NUM_THREADS=2 -d 2 -cc depth -j 1 -n 64 -N 32 test.x',
+        ('theta'          , 'n2_t2_p2'      ) : 'aprun -e OMP_NUM_THREADS=2 -d 2 -cc depth -j 1 -n 4 -N 2 test.x',
+        ('titan'          , 'n1'            ) : 'aprun -n 16 test.x',
+        ('titan'          , 'n1_p1'         ) : 'aprun -n 1 test.x',
+        ('titan'          , 'n2'            ) : 'aprun -n 32 test.x',
+        ('titan'          , 'n2_t2'         ) : 'aprun -d 2 -n 16 test.x',
+        ('titan'          , 'n2_t2_e'       ) : 'aprun -d 2 -n 16 test.x',
+        ('titan'          , 'n2_t2_p2'      ) : 'aprun -d 2 -n 4 test.x',
+        ('tomcat3'        , 'n1'            ) : 'mpirun -np 64 test.x',
+        ('tomcat3'        , 'n1_p1'         ) : 'mpirun -np 1 test.x',
+        ('tomcat3'        , 'n2'            ) : 'mpirun -np 128 test.x',
+        ('tomcat3'        , 'n2_t2'         ) : 'mpirun -np 64 test.x',
+        ('tomcat3'        , 'n2_t2_e'       ) : 'mpirun -np 64 test.x',
+        ('tomcat3'        , 'n2_t2_p2'      ) : 'mpirun -np 4 test.x',
+        ('uno'            , 'n1'            ) : 'srun test.x',
+        ('uno'            , 'n1_p1'         ) : 'srun test.x',
+        ('uno'            , 'n2'            ) : 'srun test.x',
+        ('uno'            , 'n2_t2'         ) : 'srun test.x',
+        ('uno'            , 'n2_t2_e'       ) : 'srun test.x',
+        ('uno'            , 'n2_t2_p2'      ) : 'srun test.x',
+        ('vesta'          , 'n1'            ) : 'runjob --np 16 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
+        ('vesta'          , 'n1_p1'         ) : 'runjob --np 1 -p 1 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
+        ('vesta'          , 'n2'            ) : 'runjob --np 32 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 : test.x',
+        ('vesta'          , 'n2_t2'         ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
+        ('vesta'          , 'n2_t2_e'       ) : 'runjob --np 16 -p 8 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 ENV_VAR=1 : test.x',
+        ('vesta'          , 'n2_t2_p2'      ) : 'runjob --np 4 -p 2 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=2 : test.x',
+        })
+
+    if testing.global_data['job_ref_table']:
+        print('\n\n')
+    #end if
+    job_inputs_orig = obj(
+        n1        = obj(nodes=1),
+        n1_p1     = obj(nodes=1,processes_per_node=1),
+        n2        = obj(nodes=2),
+        n2_t2     = obj(nodes=2,threads=2),
+        n2_t2_p2  = obj(nodes=2,threads=2,processes_per_node=2),
+        n2_t2_e   = obj(nodes=2,threads=2,env=obj(ENV_VAR=1)),
+        )
+    for name in sorted(supercomputers.keys()):
+        m = supercomputers[name]
+        if m.requires_account:
+            acc = 'ABC123'
+        else:
+            acc = None
+        #end if
+        job_inputs = job_inputs_orig
+        if name=='summit': # exceptional treatment for summit nodes
+            job_inputs = job_inputs_orig.copy()
+            jtypes = list(job_inputs.keys())
+            for jtype in jtypes:
+                if 'p' in jtype:
+                    del job_inputs[jtype]
+                else:
+                    jcpu = job_inputs[jtype]
+                    jcpu.gpus = 0
+                    jgpu = jcpu.copy()
+                    jgpu.gpus = 6
+                    job_inputs[jtype+'_g6'] = jgpu
+                #end if
+            #end for
+        #end if
+        for jtype in sorted(job_inputs.keys()):
+            job = Job(app_command = 'test.x',
+                      machine     = name,
+                      account     = acc,
+                      **job_inputs[jtype]
+                      )
+            command = job.run_command()
+            if testing.global_data['job_ref_table']:
+                sname = "'{0}'".format(name)
+                stype = "'{0}'".format(jtype)
+                print("        ({0:<16} , {1:<16}) : '{2}',".format(sname,stype,command))
+                continue
+            #end if
+            ref_command = job_run_ref[name,jtype]
+            if not job_commands_equal(command,ref_command):
+                failed('Job.run_command for machine "{0}" does not match the reference\njob inputs: {1}\nreference command: {2}\nincorrect command: {3}'.format(name,job_inputs[jtype],ref_command,command))
+            #end for
+        #end for
+    #end for
+
+
+    # test split_nodes
+    for name in sorted(supercomputers.keys()):
+        m = supercomputers[name]
+        if m.app_launcher=='srun': # no slurm support yet
+            continue
+        #end if
+        if name=='summit': # no summit support
+            continue
+        #end if
+        if m.requires_account:
+            acc = 'ABC123'
+        else:
+            acc = None
+        #end if
+        # make a 4 node job
+        job = Job(app_command = 'test.x',
+                  machine     = name,
+                  account     = acc,
+                  nodes       = 4,
+                  threads     = m.cores_per_node,
+                  )
+        # split the job into 1 and 3 nodes
+        job1,job2 = job.split_nodes(1)
+        # get the split run commands
+        rc  = job.run_command()
+        rc1 = job1.run_command()
+        rc2 = job2.run_command()
+        ns  = ' {0} '.format(job.nodes)
+        ns1 = ' {0} '.format(job1.nodes)
+        ns2 = ' {0} '.format(job2.nodes)
+        # verify that node count is in each command
+        assert(ns  in rc )
+        assert(ns1 in rc1)
+        assert(ns2 in rc2)
+        # verify that text on either side of node count 
+        # agrees for original and split commands
+        rcl ,rcr  = rc.split(ns,1)
+        rc1l,rc1r = rc1.split(ns1,1)
+        rc2l,rc2r = rc2.split(ns2,1)
+        rcf  = rcl+' '+rcr
+        rc1f = rc1l+' '+rc1r
+        rc2f = rc2l+' '+rc2r
+        assert(job_commands_equal(rcf,rc1f))
+        assert(job_commands_equal(rcf,rc2f))
+    #end for
+
+    Machine.allow_warnings = allow_warn
+#end def test_job_run_command
+
+
+
+def test_write_job():
+    from generic import obj
+    from machines import Machine,job
+
+    supercomputers = get_supercomputers()
+
+    allow_warn = Machine.allow_warnings
+    Machine.allow_warnings = False
+
+    job_write_ref = dict(
+        amos = '''#!/bin/bash -x
+#SBATCH --export=ALL
+#SBATCH -J None
+#SBATCH -p debug
+#SBATCH -o test.out
+#SBATCH -e test.err
+#SBATCH --nodes 2
+#SBATCH --ntasks-per-node=16
+#SBATCH --cpus-per-task=1
+#SBATCH -t 01:00:00
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+srun test.x''',
+        bluewaters_xe = '''#!/bin/bash
+#PBS -N jobname
+#PBS -l walltime=06:30:00
+#PBS -l nodes=2:ppn=32:xe
+#PBS -o test.out
+#PBS -e test.err
+#PBS -V
+
+echo $PBS_O_WORKDIR
+cd $PBS_O_WORKDIR
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+aprun -n 64 test.x''',
+        bluewaters_xk = '''#!/bin/bash
+#PBS -N jobname
+#PBS -l walltime=06:30:00
+#PBS -l nodes=2:ppn=16:xk
+#PBS -o test.out
+#PBS -e test.err
+#PBS -V
+
+echo $PBS_O_WORKDIR
+cd $PBS_O_WORKDIR
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+aprun -n 32 test.x''',
+        cades = '''#!/bin/bash
+#PBS -A ABC123
+#PBS -W group_list=cades-ABC123
+#PBS -q skylake
+#PBS -N jobname
+#PBS -o test.out
+#PBS -e test.err
+#PBS -l qos=std
+#PBS -l walltime=06:30:00
+#PBS -l nodes=2:ppn=36
+
+echo $PBS_O_WORKDIR
+cd $PBS_O_WORKDIR
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+mpirun -np 72 test.x''',
+        cetus = '''#!/bin/bash
+#COBALT -q default
+#COBALT -A ABC123
+#COBALT -n 2
+#COBALT -t 390
+#COBALT -O None
+
+LOCARGS="--block $COBALT_PARTNAME ${COBALT_CORNER:+--corner} $COBALT_CORNER ${COBALT_SHAPE:+--shape} $COBALT_SHAPE"
+echo "Cobalt location args: $LOCARGS" >&2
+
+
+runjob --np 32 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 ENV_VAR=1 : test.x''',
+        chama = '''#!/bin/bash
+#SBATCH --job-name jobname
+#SBATCH --account=ABC123
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=32
+#SBATCH --cpus-per-task=1
+#SBATCH -t 06:30:00
+#SBATCH -o test.out
+#SBATCH -e test.err
+
+module purge
+module add intel/intel-16.0.1.150
+module add libraries/intel-mkl-16.0.1.150
+module add mvapich2-intel-psm/1.7
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+srun test.x''',
+        cooley = '''#!/bin/bash
+#COBALT -q default
+#COBALT -A ABC123
+#COBALT -n 2
+#COBALT -t 390
+#COBALT -O None
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+mpirun -np 24 test.x''',
+        cori = '''#!/bin/bash
+#SBATCH -p regular
+#SBATCH -C knl
+#SBATCH -J jobname
+#SBATCH -t 06:30:00
+#SBATCH -N 2
+#SBATCH --tasks-per-node=68
+#SBATCH --cpus-per-task=4
+#SBATCH -o test.out
+#SBATCH -e test.err
+#SBATCH --export=ALL
+
+echo $SLURM_SUBMIT_DIR
+cd $SLURM_SUBMIT_DIR
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+srun test.x''',
+        edison = '''#!/bin/bash
+#SBATCH -p regular
+#SBATCH -J jobname
+#SBATCH -t 06:30:00
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=24
+#SBATCH --cpus-per-task=1
+#SBATCH -o test.out
+#SBATCH -e test.err
+#SBATCH --export=ALL
+
+echo $SLURM_SUBMIT_DIR
+cd $SLURM_SUBMIT_DIR
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+srun test.x''',
+        eos = '''#!/bin/bash
+#PBS -A ABC123
+#PBS -q batch
+#PBS -N jobname
+#PBS -o test.out
+#PBS -e test.err
+#PBS -l walltime=06:30:00
+#PBS -l nodes=2
+#PBS -l gres=atlas1
+#PBS -V
+
+echo $PBS_O_WORKDIR
+cd $PBS_O_WORKDIR
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+aprun -n 32 test.x''',
+        jaguar = '''#!/bin/bash
+#PBS -A ABC123
+#PBS -q batch
+#PBS -N jobname
+#PBS -l walltime=06:30:00
+#PBS -l size=32
+#PBS -l gres=widow2%widow3
+#PBS -o test.out
+#PBS -e test.err
+#PBS -V
+
+echo $PBS_O_WORKDIR
+cd $PBS_O_WORKDIR
+export MPICH_PTL_SEND_CREDITS=-1
+export MPICH_MAX_SHORT_MSG_SIZE=1024
+export MPICH_PTL_UNEX_EVENTS=800000
+export MPICH_UNEX_BUFFER_SIZE=16M
+export MPI_MSGS_PER_PROC=32768
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+aprun -n 32 test.x''',
+        komodo = '''#!/bin/bash -x
+#SBATCH --export=ALL
+#SBATCH -J None
+#SBATCH -p defq
+#SBATCH -o test.out
+#SBATCH -e test.err
+#SBATCH --nodes 2
+#SBATCH --ntasks-per-node=12
+#SBATCH --cpus-per-task=1
+#SBATCH -t 06:30:00
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+mpirun -np 24 test.x''',
+        kraken = '''#!/bin/bash
+#PBS -A ABC123
+#PBS -N jobname
+#PBS -l walltime=06:30:00
+#PBS -l size=24
+#PBS -o test.out
+#PBS -e test.err
+#PBS -V
+
+cd $PBS_O_WORKDIR
+export MPICH_PTL_SEND_CREDITS=-1
+export MPICH_MAX_SHORT_MSG_SIZE=1024
+export MPICH_PTL_UNEX_EVENTS=800000
+export MPICH_UNEX_BUFFER_SIZE=16M
+export MPI_MSGS_PER_PROC=32768
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+aprun -n 24 test.x''',
+        lonestar = '''#!/bin/bash
+#$ -q batch
+#$ -N jobname
+#$ -o test.out
+#$ -e test.err
+#$ -l h_rt=06:30:00
+#$ -pe 12way 24
+#$ -cwd
+#$ -V
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+ibrun -n 24 -o 0 test.x''',
+        matisse = '''#!/bin/bash -x
+#SBATCH --export=ALL
+#SBATCH -J None
+#SBATCH -p defq
+#SBATCH -o test.out
+#SBATCH -e test.err
+#SBATCH --nodes 2
+#SBATCH --ntasks-per-node=16
+#SBATCH --cpus-per-task=1
+#SBATCH -t 06:30:00
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+mpirun -np 32 test.x''',
+        mira = '''#!/bin/bash
+#COBALT -q default
+#COBALT -A ABC123
+#COBALT -n 2
+#COBALT -t 390
+#COBALT -O None
+
+LOCARGS="--block $COBALT_PARTNAME ${COBALT_CORNER:+--corner} $COBALT_CORNER ${COBALT_SHAPE:+--shape} $COBALT_SHAPE"
+echo "Cobalt location args: $LOCARGS" >&2
+
+
+runjob --np 32 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 ENV_VAR=1 : test.x''',
+        oic5 = '''#!/bin/bash
+#PBS -q mstqmc13q
+#PBS -N jobname
+#PBS -l walltime=06:30:00
+#PBS -l nodes=2:ppn=32
+#PBS -W x="NACCESSPOLICY:SINGLEJOB"
+#PBS -o test.out
+#PBS -e test.err
+#PBS -V
+
+echo $PBS_O_WORKDIR
+cd $PBS_O_WORKDIR
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+mpirun -np 64 test.x''',
+        redsky = '''#!/bin/bash
+#SBATCH --job-name jobname
+#SBATCH --account=ABC123
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=16
+#SBATCH --cpus-per-task=1
+#SBATCH -t 06:30:00
+#SBATCH -o test.out
+#SBATCH -e test.err
+
+module purge
+module add intel/intel-16.0.1.150
+module add libraries/intel-mkl-16.0.1.150
+module add mvapich2-intel-psm/1.7
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+srun test.x''',
+        serrano = '''#!/bin/bash
+#SBATCH --job-name jobname
+#SBATCH --account=ABC123
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=36
+#SBATCH --cpus-per-task=1
+#SBATCH -t 06:30:00
+#SBATCH -o test.out
+#SBATCH -e test.err
+
+module purge
+module add intel/16.0.3
+module add mkl/16.0
+module add mvapich2-intel-psm2/2.2rc1
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+srun test.x''',
+        skybridge = '''#!/bin/bash
+#SBATCH --job-name jobname
+#SBATCH --account=ABC123
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=32
+#SBATCH --cpus-per-task=1
+#SBATCH -t 06:30:00
+#SBATCH -o test.out
+#SBATCH -e test.err
+
+module purge
+module add intel/intel-16.0.1.150
+module add libraries/intel-mkl-16.0.1.150
+module add mvapich2-intel-psm/1.7
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+srun test.x''',
+        solo = '''#!/bin/bash
+#SBATCH --job-name jobname
+#SBATCH --account=ABC123
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=36
+#SBATCH --cpus-per-task=1
+#SBATCH -t 06:30:00
+#SBATCH -o test.out
+#SBATCH -e test.err
+
+module purge
+module add intel/16.0.3
+module add mkl/16.0
+module add mvapich2-intel-psm2/2.2rc1
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+srun test.x''',
+        stampede2 = '''#!/bin/bash
+#SBATCH --job-name jobname
+#SBATCH --account=ABC123
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=68
+#SBATCH --cpus-per-task=1
+#SBATCH -t 06:30:00
+#SBATCH -o test.out
+#SBATCH -e test.err
+#SBATCH -p normal
+
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+ibrun -n 136 -o 0 test.x''',
+        summit = '''#!/bin/bash
+#BSUB -P ABC123
+#BSUB -J jobname
+#BSUB -o test.out
+#BSUB -e test.err
+#BSUB -W 06:30
+#BSUB -nnodes 2
+#BSUB -alloc_flags "smt1"
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+jsrun -a 7 -r 6 -b rs -c 7 -d packed -n 12 -g 1 test.x''',
+        supermuc = '''#!/bin/bash
+#@ job_name         = jobname
+#@ job_type         = MPICH
+#@ class            = general
+#@ node             = 2
+#@ island_count     = 1,1
+#@ total_tasks      = 56
+#@ wall_clock_limit = 06:30:00
+#@ network.MPI      = sn_all,not_shared,us
+#@ initialdir       = /path/on/supermuc
+#@ output           = test.out
+#@ error            = test.err
+#@ energy_policy_tag = my_energy_tag
+#@ minimize_time_to_solution = yes
+#@ notification     = never
+#@ queue
+. /etc/profile
+. /etc/profile.d/modules.sh
+module unload mpi.ibm
+module load mpi.intel
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+mpiexec -n 56 test.x''',
+        supermucng = '''#!/bin/bash
+#SBATCH --account=ABC123
+#SBATCH --partition=general
+#SBATCH -J jobname
+#SBATCH --time=06:30:00
+#SBATCH -o ./test.out
+#SBATCH -e ./test.err
+#SBATCH --nodes=2
+#SBATCH --ntasks-per-node=48
+#SBATCH --ntasks-per-core=1
+#SBATCH --cpus-per-task=1
+#SBATCH -D ./
+#SBATCH --no-requeue
+#SBATCH --export=NONE
+#SBATCH --get-user-env
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+mpiexec -n 96 test.x''',
+        taub = '''#PBS -q cse
+#PBS -N jobname
+#PBS -l nodes=2:ppn=12
+#PBS -l walltime=06:30:00
+#PBS -e test.err
+#PBS -o test.out
+#PBS -V
+
+cd ${PBS_O_WORKDIR}
+
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+mpirun -np 24 test.x''',
+        theta = '''#!/bin/bash
+#COBALT -q default
+#COBALT -A ABC123
+#COBALT -n 2
+#COBALT -t 390
+#COBALT -O None
+#COBALT --attrs mcdram=cache:numa=quad
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+aprun -e OMP_NUM_THREADS=1 -d 1 -cc depth -j 1 -n 128 -N 64 test.x''',
+        titan = '''#!/bin/bash
+#PBS -A ABC123
+#PBS -q batch
+#PBS -N jobname
+#PBS -o test.out
+#PBS -e test.err
+#PBS -l walltime=06:30:00
+#PBS -l nodes=2
+#PBS -l gres=atlas1
+#PBS -V
+
+echo $PBS_O_WORKDIR
+cd $PBS_O_WORKDIR
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+aprun -n 32 test.x''',
+        tomcat3 = '''#!/bin/bash -l
+#SBATCH -J jobname
+#SBATCH -N 2
+#SBATCH -t 06:30:00
+#SBATCH -p tomcat
+#. /home/rcohen/.bashrc
+unalias cd; source /mnt/beegfs/intel/parallel_studio_xe_2019.3.062/bin/psxevars.sh
+ulimit -a
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+mpirun -np 128 test.x >test.out 2>test.err''',
+        uno = '''#!/bin/bash
+#SBATCH --job-name jobname
+#SBATCH --account=ABC123
+#SBATCH -N 2
+#SBATCH --ntasks-per-node=128
+#SBATCH --cpus-per-task=1
+#SBATCH -t 06:30:00
+#SBATCH -o test.out
+#SBATCH -e test.err
+#SBATCH -p quad
+
+module purge
+module add intel/intel-16.0.1.150
+module add libraries/intel-mkl-16.0.1.150
+module add mvapich2-intel-psm/1.7
+
+export OMP_NUM_THREADS=1
+export ENV_VAR=1
+srun test.x''',
+        vesta = '''#!/bin/bash
+#COBALT -q default
+#COBALT -A ABC123
+#COBALT -n 2
+#COBALT -t 390
+#COBALT -O None
+
+LOCARGS="--block $COBALT_PARTNAME ${COBALT_CORNER:+--corner} $COBALT_CORNER ${COBALT_SHAPE:+--shape} $COBALT_SHAPE"
+echo "Cobalt location args: $LOCARGS" >&2
+
+
+runjob --np 32 -p 16 $LOCARGS --verbose=INFO --envs OMP_NUM_THREADS=1 ENV_VAR=1 : test.x''',
+        )
+
+    if testing.global_data['job_ref_table']:
+        print('\n\n')
+    #end if
+    job_inputs = obj(
+        nodes       = 2,
+        hours       = 6,
+        minutes     = 30,
+        env         = obj(ENV_VAR=1),
+        identifier  = 'test',
+        outfile     = 'test.out',
+        errfile     = 'test.err',
+        app_command = 'test.x',
+        )
+    for name in sorted(supercomputers.keys()):
+        m = supercomputers[name]
+        if m.requires_account:
+            acc = 'ABC123'
+        else:
+            acc = None
+        #end if
+        ji = job_inputs.copy()
+        if name=='summit': # exceptional treatment for summit nodes
+            ji.gpus = 6
+        #end if
+        j = job(
+            machine = name,
+            account = acc,
+            **ji
+            )
+        j.abs_dir = '/path/on/'+name
+        wj = m.write_job(j)
+        if testing.global_data['job_ref_table']:
+            print("        {} = '''{}''',".format(name,wj.strip()))
+            continue
+        #end if
+        ref_wj = job_write_ref[name]
+        assert(wj.strip()==ref_wj)
+    #end for
+
+    Machine.allow_warnings = allow_warn
+#end def test_write_job
diff --git a/nexus/tests/unit/test_numerics.py b/nexus/tests/unit/test_numerics.py
index bb951f570..9e1188fbf 100644
--- a/nexus/tests/unit/test_numerics.py
+++ b/nexus/tests/unit/test_numerics.py
@@ -385,6 +385,8 @@ def test_equilibration_length():
 
     eq = equilibration_length(rs,random=False)
     assert(eq==52)
+
+    eq = equilibration_length(rs,random=True)
 #end def test_equilibration_length
 
 
diff --git a/nexus/tests/unit/test_structure.py b/nexus/tests/unit/test_structure.py
index f30ae0b8c..a603eed96 100644
--- a/nexus/tests/unit/test_structure.py
+++ b/nexus/tests/unit/test_structure.py
@@ -306,6 +306,70 @@ def test_change_units():
 
 
 
+def test_rotate():
+    import numpy as np
+    ref = get_reference_structures()
+    s0 = ref.CuO_prim.copy()
+    s1 = ref.CuO_prim.copy()
+
+    # Test the various parameter choices in the case that rp is given
+    # Perform rotation taking x-axis to x-axis (original positions should be found)
+    s1.rotate('x','x')
+    assert(value_eq(s1.pos,s0.pos))
+    assert(value_eq(s1.axes,s0.axes))
+    assert(value_eq(s1.kaxes,s0.kaxes))
+    # Perform active rotation taking x-coords to y-coords
+    s1.rotate([1,0,0],[0,1,0])
+    assert(value_eq(s1.pos,np.array([-s0.pos[:,1],s0.pos[:,0],s0.pos[:,2]]).T))
+    assert(value_eq(s1.axes,np.array([-s0.axes[:,1],s0.axes[:,0],s0.axes[:,2]]).T))
+    assert(value_eq(s1.kaxes,np.array([-s0.kaxes[:,1],s0.kaxes[:,0],s0.kaxes[:,2]]).T))
+    # Perform passive rotation taking x-axis to y-axis (original positions should be found)
+    s1.rotate('x','y',passive=True)
+    assert(value_eq(s1.pos,s0.pos))
+    assert(value_eq(s1.axes,s0.axes))
+    assert(value_eq(s1.kaxes,s0.kaxes))
+    # Perform active rotation about z-axis by an angle pi/2
+    s1.rotate('z',np.pi/2.0)
+    assert(value_eq(s1.pos,np.array([-s0.pos[:,1],s0.pos[:,0],s0.pos[:,2]]).T))
+    assert(value_eq(s1.axes,np.array([-s0.axes[:,1],s0.axes[:,0],s0.axes[:,2]]).T))
+    assert(value_eq(s1.kaxes,np.array([-s0.kaxes[:,1],s0.kaxes[:,0],s0.kaxes[:,2]]).T))
+    # Perform active rotation taking y-coords to x-coords (original positions should be found)
+    s1.rotate('y',[1,0,0])
+    assert(value_eq(s1.pos[-1],s0.pos[-1]))
+    assert(value_eq(s1.axes[-1],s0.axes[-1]))
+    assert(value_eq(s1.kaxes[-1],s0.kaxes[-1]))
+    # Perform active rotation taking a0-coords to a2-coords
+    s1.rotate(s0.axes[0],s0.axes[2])
+    assert(value_eq(s1.pos[-1],np.array([-2.15536928,3.46035669,0.86507139])))
+    assert(value_eq(s1.axes[-1],np.array([-3.91292278,3.02549423,-1.35344154])))
+    assert(value_eq(s1.kaxes[-1],np.array([-0.90768302,0.83458438,-0.15254555])))
+    # Perform active rotation taking a2-coords to a0-coords (original positions should be found)
+    s1.rotate('a2','a0')
+    assert(value_eq(s1.pos,s0.pos))
+    assert(value_eq(s1.axes,s0.axes))
+    assert(value_eq(s1.kaxes,s0.kaxes))
+
+    # Test the case where rp is not given
+    # Perform active rotation taking a2-coords to a0-coords
+    R = [[0.2570157723433977, 0.6326366344635742,-0.7305571719594085], 
+         [0.4370696746690278, 0.5981289557203555, 0.6717230469572912], 
+         [0.8619240060767753,-0.4919478031900122,-0.12277771249328594]]
+    s1.rotate(R)
+    assert(value_eq(s1.pos[-1],np.array([-2.15536928,3.46035669,0.86507139])))
+    assert(value_eq(s1.axes[-1],np.array([-3.91292278,3.02549423,-1.35344154])))
+    assert(value_eq(s1.kaxes[-1],np.array([-0.90768302,0.83458438,-0.15254555])))
+
+    # A final test which places the structure back into its original form
+    # Perform passive rotation taking a2-coords to a0-coords (original positions should be found)
+    s1.rotate([-0.5871158698555267, -0.8034668669004766, -0.09867091342903483],1.7050154439645373,passive=True)
+    assert(value_eq(s1.pos,s0.pos))
+    assert(value_eq(s1.axes,s0.axes))
+    assert(value_eq(s1.kaxes,s0.kaxes))
+
+#end def test_change_units
+
+
+
 def test_diagonal_tiling():
     ref = get_reference_structures()
     diag_tilings = [
diff --git a/src/QMCWaveFunctions/CMakeLists.txt b/src/QMCWaveFunctions/CMakeLists.txt
index b1167d82f..09466bdf8 100644
--- a/src/QMCWaveFunctions/CMakeLists.txt
+++ b/src/QMCWaveFunctions/CMakeLists.txt
@@ -28,7 +28,6 @@ SET(WFBASE_SRCS
   SPOSetInfo.cpp
   SPOSetInputInfo.cpp
   SPOSet.cpp
-  OptimizableSPOSet.cpp
   CompositeSPOSet.cpp
   HarmonicOscillator/SHOSet.cpp
   HarmonicOscillator/SHOSetBuilder.cpp
@@ -56,7 +55,6 @@ SET(JASTROW_SRCS
   Jastrow/RPAJastrow.cpp
   LatticeGaussianProduct.cpp
   LatticeGaussianProductBuilder.cpp
-  OptimizableSPOBuilder.cpp
   Fermion/SPOSetProxy.cpp
   Fermion/SPOSetProxyForMSD.cpp
   )
diff --git a/src/QMCWaveFunctions/Fermion/DiracDeterminantBase.h b/src/QMCWaveFunctions/Fermion/DiracDeterminantBase.h
index 52d7a611a..c999dee48 100644
--- a/src/QMCWaveFunctions/Fermion/DiracDeterminantBase.h
+++ b/src/QMCWaveFunctions/Fermion/DiracDeterminantBase.h
@@ -43,9 +43,9 @@ public:
         SPOVTimer(*TimerManager.createTimer("DiracDeterminantBase::spoval", timer_level_fine)),
         SPOVGLTimer(*TimerManager.createTimer("DiracDeterminantBase::spovgl", timer_level_fine))
   {
-    Optimizable = Phi->Optimizable;
+    Optimizable  = Phi->isOptimizable();
     is_fermionic = true;
-    ClassName   = "DiracDeterminantBase";
+    ClassName    = "DiracDeterminantBase";
     registerTimers();
   }
 
@@ -76,26 +76,16 @@ public:
   virtual inline void checkInVariables(opt_variables_type& active) override
   {
     Phi->checkInVariables(active);
-    Phi->checkInVariables(myVars);
   }
 
   virtual inline void checkOutVariables(const opt_variables_type& active) override
   {
     Phi->checkOutVariables(active);
-    myVars.clear();
-    myVars.insertFrom(Phi->myVars);
-    myVars.getIndex(active);
   }
 
   virtual void resetParameters(const opt_variables_type& active) override
   {
     Phi->resetParameters(active);
-    for (int i = 0; i < myVars.size(); ++i)
-    {
-      int ii = myVars.Index[i];
-      if (ii >= 0)
-        myVars[i] = active[ii];
-    }
   }
 
   // To be removed with AoS
diff --git a/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.cpp b/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.cpp
index 31a465296..fa993cb49 100644
--- a/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.cpp
+++ b/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.cpp
@@ -392,7 +392,7 @@ MultiDiracDeterminant::MultiDiracDeterminant(SPOSetPtr const& spos, int first)
       buildTableGradTimer(*TimerManager.createTimer("MultiDiracDeterminant::buildTableGrad")),
       ExtraStuffTimer(*TimerManager.createTimer("MultiDiracDeterminant::ExtraStuff"))
 {
-  (spos->Optimizable == true) ? Optimizable = true : Optimizable = false;
+  (spos->isOptimizable() == true) ? Optimizable = true : Optimizable = false;
   ClassName                                                      = "MultiDiracDeterminant";
 
   IsCloned = false;
diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp
index fc6b2b179..566fbd738 100644
--- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp
+++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp
@@ -589,7 +589,7 @@ bool SlaterDetBuilder::putDeterminant(xmlNodePtr cur, int spin_group)
   targetPsi.setndelay(delay_rank);
 #endif
   slaterdet_0->add(adet, spin_group);
-  if (psi->Optimizable)
+  if (psi->isOptimizable())
     slaterdet_0->Optimizable = true;
 
   app_log() << std::endl;
diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetOpt.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetOpt.cpp
index afe862f00..79674946a 100644
--- a/src/QMCWaveFunctions/Fermion/SlaterDetOpt.cpp
+++ b/src/QMCWaveFunctions/Fermion/SlaterDetOpt.cpp
@@ -46,8 +46,8 @@ SlaterDetOpt::SlaterDetOpt(ParticleSet& ptcl, SPOSet* spo_ptr, const int up_or_d
 
   resize(m_nel, m_nmo);
 
-  m_nlc = Phi->OrbitalSetSize;
-  m_nb  = Phi->BasisSetSize;
+  m_nlc = Phi->getOrbitalSetSize();
+  m_nb  = Phi->getBasisSetSize();
 
   // make sure we didn't start with a bad m_nlc
   check_index_sanity();
diff --git a/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.cpp b/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.cpp
index 55717e4cd..d1c7670c1 100644
--- a/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.cpp
+++ b/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.cpp
@@ -25,7 +25,7 @@ void fillRadFunWithPhiBar(int curOrb_,
 {
   Psi1      = Phi;
   Psi2      = Eta;
-  int norb  = Psi1->OrbitalSetSize;
+  int norb  = Psi1->getOrbitalSetSize();
   curOrb    = curOrb_;
   curCenter = curCenter_;
 
@@ -56,7 +56,7 @@ void fillRadFunWithPhi(int curOrb_,
 {
   Psi1      = Phi;
   Psi2      = Eta;
-  int norb  = Psi1->OrbitalSetSize;
+  int norb  = Psi1->getOrbitalSetSize();
   curOrb    = curOrb_;
   curCenter = curCenter_;
   Z         = Zion;
@@ -177,7 +177,7 @@ RealType execute(int curOrb_,
   bool writeout = (thisFile != "NULL");
   Psi1          = Phi;
   Psi2          = Eta;
-  int norb      = Psi1->OrbitalSetSize;
+  int norb      = Psi1->getOrbitalSetSize();
   ValueVector_t X(5);
   curOrb    = curOrb_;
   curCenter = curCenter_;
diff --git a/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.h b/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.h
index f6bd31f80..cd63aa716 100644
--- a/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.h
+++ b/src/QMCWaveFunctions/MolecularOrbitals/CuspCorr.h
@@ -104,7 +104,7 @@ public:
   {
     Psi1     = Phi;
     Psi2     = Eta;
-    int norb = Psi1->OrbitalSetSize;
+    int norb = Psi1->getOrbitalSetSize();
 
     val1.resize(norb);
     grad1.resize(norb);
diff --git a/src/QMCWaveFunctions/MolecularOrbitals/LCOrbitalSetOpt.h b/src/QMCWaveFunctions/MolecularOrbitals/LCOrbitalSetOpt.h
index 01abbd756..bc20e4b45 100644
--- a/src/QMCWaveFunctions/MolecularOrbitals/LCOrbitalSetOpt.h
+++ b/src/QMCWaveFunctions/MolecularOrbitals/LCOrbitalSetOpt.h
@@ -274,7 +274,7 @@ public:
       this->setBasisSet(bs);
 
     // initialize number of molecular orbitals as zero
-    this->OrbitalSetSize = 0;
+    setOrbitalSetSize(0);
   }
 
   ///////////////////////////////////////////////////////////////////////////////////////////////////
@@ -291,7 +291,7 @@ public:
       this->setSPOSet(spo);
 
     // initialize number of molecular orbitals as zero
-    this->OrbitalSetSize = 0;
+    setOrbitalSetSize(0);
   }
 
   ///////////////////////////////////////////////////////////////////////////////////////////////////
@@ -358,7 +358,7 @@ public:
     m_spo_set = spo;
 
     // extract the number of single particle orbitals in the basis set
-    this->BasisSetSize = m_spo_set->OrbitalSetSize;
+    this->BasisSetSize = m_spo_set->getOrbitalSetSize();
   }
 
   ///////////////////////////////////////////////////////////////////////////////////////////////////
@@ -380,7 +380,7 @@ public:
       retval = new LCOrbitalSetOpt(m_basis_set->makeClone(), m_report_level);
 
     retval->C = C;
-    retval->setOrbitalSetSize(this->OrbitalSetSize);
+    retval->setOrbitalSetSize(OrbitalSetSize);
     retval->init_LCOrbitalSetOpt(0.0);
 
     retval->m_B      = m_B;
@@ -429,7 +429,7 @@ public:
   void setOrbitalSetSize(int norbs)
   {
     // record the number of linear combinations (i.e. molecular orbitals)
-    this->OrbitalSetSize = norbs;
+    OrbitalSetSize = norbs;
     app_log() << "LCOrbitalSetOpt finished setOrbitalSetSize with norbs = " << norbs << std::endl;
   }
 
@@ -698,7 +698,7 @@ public:
     //app_log() << "this->OrbitalSetSize = " << this->OrbitalSetSize << std::endl;
     // check sanity
     this->check_input_dim("logdet # of columns", "LCOrbitalSetOpt::evaluate_notranspose", logdet.cols(),
-                          this->OrbitalSetSize);
+                          OrbitalSetSize);
     if (logdet.cols() != dlogdet.cols() || logdet.cols() != d2logdet.cols())
       throw std::runtime_error("logdet, dlogdet, and d2logdet should have the same number of columns in "
                                "LCOrbitalSetOpt::evaluate_notranspose");
@@ -724,7 +724,7 @@ public:
   inline void evaluate(const ParticleSet& P, int iat, ValueVector_t& psi)
   {
     // check input vector dimension
-    this->check_input_dim("psi", "LCOrbitalSetOpt::evaluate", psi.size(), this->OrbitalSetSize);
+    this->check_input_dim("psi", "LCOrbitalSetOpt::evaluate", psi.size(), OrbitalSetSize);
 
     // resize temporary arrays if necessary
     if (m_temp_g.size() != BasisSetSize)
@@ -757,9 +757,9 @@ public:
   inline void evaluate(const ParticleSet& P, int iat, ValueVector_t& psi, GradVector_t& dpsi, ValueVector_t& d2psi)
   {
     // check sanity
-    this->check_input_dim("d2psi", "LCOrbitalSetOpt::evaluate", d2psi.size(), this->OrbitalSetSize);
-    this->check_input_dim("dpsi", "LCOrbitalSetOpt::evaluate", dpsi.size(), this->OrbitalSetSize);
-    this->check_input_dim("psi", "LCOrbitalSetOpt::evaluate", psi.size(), this->OrbitalSetSize);
+    this->check_input_dim("d2psi", "LCOrbitalSetOpt::evaluate", d2psi.size(), OrbitalSetSize);
+    this->check_input_dim("dpsi", "LCOrbitalSetOpt::evaluate", dpsi.size(), OrbitalSetSize);
+    this->check_input_dim("psi", "LCOrbitalSetOpt::evaluate", psi.size(), OrbitalSetSize);
     if (psi.size() != dpsi.size() || psi.size() != d2psi.size())
       throw std::runtime_error("psi, dpsi, and d2psi vectors must be the same length in LCOrbitalSetOpt::evaluate");
 
diff --git a/src/QMCWaveFunctions/MolecularOrbitals/LCOrbitalSetWithCorrection.cpp b/src/QMCWaveFunctions/MolecularOrbitals/LCOrbitalSetWithCorrection.cpp
index 4396af6fe..ea9d66632 100644
--- a/src/QMCWaveFunctions/MolecularOrbitals/LCOrbitalSetWithCorrection.cpp
+++ b/src/QMCWaveFunctions/MolecularOrbitals/LCOrbitalSetWithCorrection.cpp
@@ -192,7 +192,6 @@ LCOrbitalSet<BS, false>* LCOrbitalSetWithCorrection<BS, false>::clone2LCOrbitalS
   newSPO->C = C;
   newSPO->Occ.resize(Occ.size());
   newSPO->Occ       = Occ;
-  newSPO->className = "LCOrbitalSet";
   return newSPO;
 }
 
diff --git a/src/QMCWaveFunctions/OptimizableSPOBuilder.cpp b/src/QMCWaveFunctions/OptimizableSPOBuilder.cpp
deleted file mode 100644
index 7f6d57107..000000000
--- a/src/QMCWaveFunctions/OptimizableSPOBuilder.cpp
+++ /dev/null
@@ -1,31 +0,0 @@
-//////////////////////////////////////////////////////////////////////////////////////
-// This file is distributed under the University of Illinois/NCSA Open Source License.
-// See LICENSE file in top directory for details.
-//
-// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
-//
-// File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
-//                    Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
-//                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
-//                    Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
-//
-// File created by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
-//////////////////////////////////////////////////////////////////////////////////////
-
-
-#include "QMCWaveFunctions/OptimizableSPOBuilder.h"
-#include "OhmmsData/AttributeSet.h"
-
-namespace qmcplusplus
-{
-OptimizableSPOBuilder::OptimizableSPOBuilder(ParticleSet& p, PtclPoolType& psets, Communicate* comm, xmlNodePtr cur)
-    : SPOSetBuilder(comm), targetPtcl(&p)
-{}
-
-SPOSet* OptimizableSPOBuilder::createSPOSetFromXML(xmlNodePtr cur)
-{
-  OptimizableSPOSet* spo = new OptimizableSPOSet();
-  return spo;
-}
-
-} // namespace qmcplusplus
diff --git a/src/QMCWaveFunctions/OptimizableSPOBuilder.h b/src/QMCWaveFunctions/OptimizableSPOBuilder.h
deleted file mode 100644
index 6a3971fb1..000000000
--- a/src/QMCWaveFunctions/OptimizableSPOBuilder.h
+++ /dev/null
@@ -1,43 +0,0 @@
-//////////////////////////////////////////////////////////////////////////////////////
-// This file is distributed under the University of Illinois/NCSA Open Source License.
-// See LICENSE file in top directory for details.
-//
-// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
-//
-// File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
-//                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
-//                    Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
-//                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
-//
-// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
-//////////////////////////////////////////////////////////////////////////////////////
-
-
-#ifndef OPTIMIZABLE_SPO_BUILDER_H
-#define OPTIMIZABLE_SPO_BUILDER_H
-
-#include "QMCWaveFunctions/SPOSetBuilder.h"
-#include "QMCWaveFunctions/OptimizableSPOSet.h"
-
-
-namespace qmcplusplus
-{
-class OptimizableSPOBuilder : public SPOSetBuilder
-{
-protected:
-  typedef std::map<std::string, ParticleSet*> PtclPoolType;
-  typedef std::map<std::string, SPOSet*> SPOPoolType;
-  ParticleSet* targetPtcl;
-
-public:
-  OptimizableSPOBuilder(ParticleSet& p, PtclPoolType& psets, Communicate* comm, xmlNodePtr cur);
-
-  /** initialize the Antisymmetric wave function for electrons
-   *@param cur the current xml node
-   */
-  SPOSet* createSPOSetFromXML(xmlNodePtr cur);
-  //    SPOSet* createSPOSetFromXML(xmlNodePtr cur, SPOPool_t& spo_pool);
-};
-} // namespace qmcplusplus
-
-#endif
diff --git a/src/QMCWaveFunctions/OptimizableSPOSet.cpp b/src/QMCWaveFunctions/OptimizableSPOSet.cpp
deleted file mode 100644
index 99f2e3e81..000000000
--- a/src/QMCWaveFunctions/OptimizableSPOSet.cpp
+++ /dev/null
@@ -1,618 +0,0 @@
-//////////////////////////////////////////////////////////////////////////////////////
-// This file is distributed under the University of Illinois/NCSA Open Source License.
-// See LICENSE file in top directory for details.
-//
-// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
-//
-// File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
-//                    Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
-//                    Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
-//                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
-//                    Ye Luo, yeluo@anl.gov, Argonne National Laboratory
-//                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
-//
-// File created by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
-//////////////////////////////////////////////////////////////////////////////////////
-
-
-#include "Message/Communicate.h"
-#include "QMCWaveFunctions/OptimizableSPOSet.h"
-#include "Numerics/OhmmsBlas.h"
-#include "OhmmsData/AttributeSet.h"
-
-namespace qmcplusplus
-{
-void OptimizableSPOSet::addParameter(std::string id, int iorb, int basis) {}
-
-bool OptimizableSPOSet::put(xmlNodePtr node, SPOPool_t& spo_pool)
-{
-  std::string gsName, basisName;
-  bool same_k       = false;
-  bool mapped_k     = false;
-  bool same_orbital = false;
-  bool lapl         = false;
-  derivScale        = 1.0;
-  OhmmsAttributeSet attrib;
-  attrib.add(gsName, "gs_sposet");
-  attrib.add(basisName, "basis_sposet");
-  attrib.add(same_k, "same_k");
-  attrib.add(mapped_k, "mapped_k");
-  attrib.add(same_orbital, "same_orbital");
-  attrib.add(lapl, "lapl");
-  attrib.add(N, "size");
-  attrib.add(derivScale, "scale");
-  attrib.add(thr, "thr");
-  attrib.put(node);
-  if (N == 0)
-  {
-    app_error() << "You must specify \"size\" attribute for linearopt sposet.\n";
-    abort();
-  }
-  /////////////////////////////////////
-  // First, find ground-state SPOSet //
-  /////////////////////////////////////
-  if (gsName == "")
-  {
-    app_error() << "You must supply \"gs_sposet\".  Aborting.\n";
-    abort();
-  }
-  SPOPool_t::iterator iter = spo_pool.find(gsName);
-  if (iter == spo_pool.end())
-  {
-    app_error() << "No sposet named \"" << gsName << "\" found.  Abort.\n";
-    abort();
-  }
-  else
-  {
-    app_log() << "  Found ground-state SPOSet \"" << gsName << "\".\n";
-    GSOrbitals = iter->second;
-  }
-  //////////////////////////////////////
-  // Now, find basis SPOSet from pool //
-  //////////////////////////////////////
-  iter = spo_pool.find(basisName);
-  if (iter != spo_pool.end())
-  {
-    BasisOrbitals = iter->second;
-    app_log() << "  Found basis SPOSet \"" << basisName << "\".\n";
-  }
-  if (BasisOrbitals == GSOrbitals)
-    BasisOrbitals = 0;
-  if (BasisOrbitals)
-    M = BasisOrbitals->getOrbitalSetSize();
-  else
-    M = GSOrbitals->getOrbitalSetSize() - N;
-  resize(N, M);
-  app_log() << "  linearopt sposet has " << N << " ground-state orbitals and " << M << " basis orbitals.\n";
-  //Map K points to orbitals.
-  ValueMatrix_t allowedOrbs;
-  allowedOrbs.resize(N, M);
-  //    if (mapped_k)
-  //    {
-  //      std::cerr <<"Not Done"<< std::endl;
-  //    }
-  //    else
-  if (same_k)
-  {
-    int off         = BasisOrbitals ? 0 : N;
-    SPOSet* basOrbs = BasisOrbitals ? BasisOrbitals : GSOrbitals;
-    for (int igs = 0; igs < N; igs++)
-    {
-      PosType k_gs = GSOrbitals->get_k(igs);
-      for (int ib = 0; ib < M; ib++)
-      {
-        PosType k_b = basOrbs->get_k(ib + off);
-        if (dot(k_gs - k_b, k_gs - k_b) < 1.0e-6)
-        {
-          //             app_log()<<"Allowing ("<<igs<<","<<ib<<") :";
-          //             for(int ix=0;ix<3;ix++) app_log()<<k_gs[ix]<<"  ";
-          //             app_log()<<" : ";
-          //             for(int ix=0;ix<3;ix++) app_log()<<k_b[ix]<<"  ";
-          //             app_log()<< std::endl;
-          allowedOrbs(igs, ib) = 1;
-        }
-        else
-        {
-          //             app_log()<<"Not allowing ("<<igs<<","<<ib<<") :";
-          //             for(int ix=0;ix<3;ix++) app_log()<<k_gs[ix]<<"  ";
-          //             app_log()<<" : ";
-          //             for(int ix=0;ix<3;ix++) app_log()<<k_b[ix]<<"  ";
-          //             app_log()<< std::endl;
-          allowedOrbs(igs, ib) = 0;
-        }
-      }
-    }
-  }
-  else if (mapped_k)
-  {
-    xmlNodePtr xmlCoefs = node->xmlChildrenNode;
-    while (xmlCoefs != NULL)
-    {
-      std::string cname((const char*)xmlCoefs->name);
-      if (cname == "orbitalmap")
-      {
-        std::string type("0");
-        OhmmsAttributeSet cAttrib;
-        cAttrib.add(type, "type");
-        cAttrib.put(xmlCoefs);
-        if (type != "Array")
-        {
-          // app_error() << "Unknown correlation type " + type + " in OptimizableSPOSet." + "Resetting to \"Array\"\n";
-          xmlNewProp(xmlCoefs, (const xmlChar*)"type", (const xmlChar*)"Array");
-        }
-        std::vector<RealType> params;
-        putContent(params, xmlCoefs);
-        int indx(0);
-        if (params.size() == N * M)
-          for (int igs = 0; igs < N; igs++)
-          {
-            for (int ib = 0; ib < M; ib++)
-            {
-              allowedOrbs(igs, ib) = params[indx++];
-            }
-          }
-        else
-        {
-          app_error() << "Map size is incorrect. parameters given= " << params.size() << " parameters needed= " << M * N
-                      << std::endl;
-        }
-      }
-      xmlCoefs = xmlCoefs->next;
-    }
-  }
-  else if (same_orbital)
-  {
-    if (M < N)
-    {
-      app_error() << "Map size is incorrect. N<M for same orbital mapping" << std::endl;
-    }
-    allowedOrbs = 0;
-    for (int igs = 0; igs < N; igs++)
-      allowedOrbs(igs, igs) = 1;
-  }
-  else if (lapl)
-  {
-    for (int igs = 0; igs < N; igs++)
-      for (int ib = 0; ib < (M / N); ib++)
-        allowedOrbs(igs, igs + N * ib) = 1;
-  }
-  else
-  {
-    for (int igs = 0; igs < N; igs++)
-    {
-      for (int ib = 0; ib < M; ib++)
-      {
-        allowedOrbs(igs, ib) = 1;
-      }
-    }
-  }
-  // Now, look for coefficients element
-  xmlNodePtr xmlCoefs = node->xmlChildrenNode;
-  while (xmlCoefs != NULL)
-  {
-    std::string cname((const char*)xmlCoefs->name);
-    if (cname == "coefficients")
-    {
-      std::string type("0"), id("0");
-      int state = -1;
-      int asize(-1);
-      OhmmsAttributeSet cAttrib;
-      cAttrib.add(id, "id");
-      cAttrib.add(type, "type");
-      cAttrib.add(asize, "size");
-      cAttrib.add(state, "state");
-      cAttrib.put(xmlCoefs);
-      if (state == -1)
-      {
-        app_error() << "You must specify the \"state\" attribute in <linearopt>'s <coefficient>.\n";
-      }
-      if (type != "Array")
-      {
-        // app_error() << "Unknown correlation type " + type + " in OptimizableSPOSet." + "Resetting to \"Array\"\n";
-        xmlNewProp(xmlCoefs, (const xmlChar*)"type", (const xmlChar*)"Array");
-      }
-      std::vector<RealType> params;
-      putContent(params, xmlCoefs);
-      app_log() << "Coefficients for state" << state << ":\n";
-      // std::cerr << "params.size() = " << params.size() << std::endl;
-      //If params is missized resize and initialize to zero.
-      if ((params.size() != asize) && (asize > 0))
-      {
-        std::vector<RealType> t_params(params);
-        params.resize(asize, 0.0);
-        int ipm(0);
-        if (t_params.size() > 0)
-        {
-          for (int ib = 0; ib < M; ib++)
-          {
-#if defined(QMC_COMPLEX)
-            if ((allowedOrbs(state, ib).real() > 0) && (std::abs(t_params[ipm]) > thr))
-            {
-              params[ib] = t_params[ipm++];
-#else
-            if ((allowedOrbs(state, ib) > 0) && (std::abs(t_params[ipm]) > thr))
-            {
-              params[ib] = t_params[ipm++];
-#endif
-            }
-            else
-            {
-              params[ib] = 0;
-            }
-          }
-        }
-      }
-      else
-      {
-        if (params.size())
-          asize = params.size();
-      }
-      //    for (int i=0; i< params.size(); i++) {
-      ValueMatrix_t& cref(*C);
-      int indx = 0;
-      for (int i = 0; i < M; i++)
-      {
-        std::stringstream sstr;
-#ifndef QMC_COMPLEX
-        // 	  ParamPointers.push_back(&(C(state,i)));
-        // 	  ParamIndex.push_back(TinyVector<int,2>(state,i));
-        sstr << id << "_" << indx;
-        if (allowedOrbs(state, i))
-        {
-          ParamPointers.push_back(&(cref(state, i)));
-          ParamIndex.push_back(TinyVector<int, 2>(state, i));
-          if ((indx < asize) && (std::abs(params[indx]) > thr))
-            cref(state, i) = params[indx];
-          else
-            cref(state, i) = 0.0;
-          myVars.insert(sstr.str(), cref(state, i), true, optimize::SPO_P);
-          indx++;
-        }
-        else
-        {
-          cref(state, i) = 0.0;
-          //        myVars.insert(sstr.str(),C(state,i),false,optimize::LINEAR_P);
-        }
-#else
-        //ParamPointers.push_back(&(C(state,i).real()));
-        //ParamPointers.push_back(&(C(state,i).imag()));
-        ParamPointers.push_back(reinterpret_cast<RealType*>(cref[state] + i));
-        ParamPointers.push_back(reinterpret_cast<RealType*>(cref[state] + i) + 1);
-        ParamIndex.push_back(TinyVector<int, 2>(state, i));
-        ParamIndex.push_back(TinyVector<int, 2>(state, i));
-        sstr << id << "_" << 2 * i + 0;
-        myVars.insert(sstr.str(), cref(state, i).real(), true, optimize::SPO_P);
-        sstr << id << "_" << 2 * i + 1;
-        myVars.insert(sstr.str(), cref(state, i).imag(), true, optimize::SPO_P);
-#endif
-      }
-      for (int i = 0; i < params.size(); i++)
-      {
-        char buf[100];
-        snprintf(buf, 100, "  %12.5f\n", params[i]);
-        app_log() << buf;
-      }
-    }
-    xmlCoefs = xmlCoefs->next;
-  }
-  return true;
-}
-
-
-void OptimizableSPOSet::resetTargetParticleSet(ParticleSet& P)
-{
-  GSOrbitals->resetTargetParticleSet(P);
-  if (BasisOrbitals)
-    BasisOrbitals->resetTargetParticleSet(P);
-}
-
-void OptimizableSPOSet::setOrbitalSetSize(int norbs) { OrbitalSetSize = norbs; }
-
-
-void OptimizableSPOSet::checkInVariables(opt_variables_type& active) { active.insertFrom(myVars); }
-
-void OptimizableSPOSet::checkOutVariables(const opt_variables_type& active) { myVars.getIndex(active); }
-
-
-void OptimizableSPOSet::resetParameters(const opt_variables_type& active)
-{
-  for (int i = 0; i < ParamPointers.size(); i++)
-  {
-    int loc=myVars.where(i);
-    if (loc>=0)
-      if(std::abs(active[loc])>thr) {
-        myVars[i]=active[loc];
-        //FIXME the use of "std::real" needs to be removed after the orbital optimization code is made compatible with complex parameters
-        *(ParamPointers[i]) = std::real( myVars[i] );
-      }
-  }
-}
-
-// Obsolete
-void OptimizableSPOSet::evaluateDerivatives(ParticleSet& P,
-                                            int iat,
-                                            const opt_variables_type& active,
-                                            ValueMatrix_t& d_phi,
-                                            ValueMatrix_t& d_lapl_phi)
-{
-  for (int i = 0; i < d_phi.size(); i++)
-    for (int j = 0; j < N; j++)
-      d_phi[i][j] = d_lapl_phi[i][j] = ValueType();
-  // Evaluate basis states
-  if (BasisOrbitals)
-  {
-    BasisOrbitals->evaluate(P, iat, BasisVal);
-    std::vector<TinyVector<int, 2>>::iterator iter;
-    std::vector<TinyVector<int, 2>>& act = ActiveBasis[iat];
-    for (iter = act.begin(); iter != act.end(); iter++)
-    {
-      int elem  = (*iter)[0];
-      int param = (*iter)[1];
-      int loc   = myVars.where(param);
-      if (loc >= 0)
-        ;
-    }
-  }
-  else
-  {}
-}
-
-void OptimizableSPOSet::evaluate(const ParticleSet& P, int iat, ValueVector_t& psi)
-{
-  GSOrbitals->evaluate(P, iat, GSVal);
-  if (BasisOrbitals)
-  {
-    BasisOrbitals->evaluate(P, iat, BasisVal);
-    BLAS::gemv_trans(N, M, C->data(), &(BasisVal[0]), &(psi[0]));
-  }
-  else
-    BLAS::gemv_trans(N, M, C->data(), &(GSVal[N]), &(psi[0]));
-  //       for (int i=0; i<N; i++) {
-  // 	psi[i] = 0.0;
-  // 	for (int j=0; j<M; j++)
-  // 	  psi[i] += C(i,j)*GSVal[N+j];
-  //       }
-  for (int i = 0; i < N; i++)
-    psi[i] += GSVal[i];
-}
-
-void OptimizableSPOSet::evaluate(const ParticleSet& P,
-                                 int iat,
-                                 ValueVector_t& psi,
-                                 GradVector_t& dpsi,
-                                 ValueVector_t& d2psi)
-{
-  GSOrbitals->evaluate(P, iat, GSVal, GSGrad, GSLapl);
-  const ValueMatrix_t& cref(*C);
-  if (BasisOrbitals)
-  {
-    BasisOrbitals->evaluate(P, iat, BasisVal, BasisGrad, BasisLapl);
-    BLAS::gemv_trans(N, M, C->data(), &(BasisVal[0]), &(psi[0]));
-    BLAS::gemv_trans(N, M, C->data(), &(BasisLapl[0]), &(d2psi[0]));
-    for (int iorb = 0; iorb < N; iorb++)
-    {
-      psi[iorb] += GSVal[iorb];
-      dpsi[iorb] = GSGrad[iorb];
-      d2psi[iorb] += GSLapl[iorb];
-      for (int ibasis = 0; ibasis < M; ibasis++)
-        for (int dim = 0; dim < OHMMS_DIM; dim++)
-          dpsi[iorb][dim] += cref(iorb, ibasis) * BasisGrad[ibasis][dim];
-    }
-  }
-  else
-  {
-    for (int iorb = 0; iorb < N; iorb++)
-    {
-      psi[iorb]   = GSVal[iorb];
-      dpsi[iorb]  = GSGrad[iorb];
-      d2psi[iorb] = GSLapl[iorb];
-      for (int ibasis = 0; ibasis < M; ibasis++)
-      {
-        psi[iorb] += cref(iorb, ibasis) * GSVal[N + ibasis];
-        for (int dim = 0; dim < OHMMS_DIM; dim++)
-          dpsi[iorb][dim] += cref(iorb, ibasis) * GSGrad[N + ibasis][dim];
-        d2psi[iorb] += cref(iorb, ibasis) * GSLapl[N + ibasis];
-      }
-    }
-    // BLAS::gemv_trans (N, M, C.data(), &(GSVal[N]),  &(psi[0]));
-    // BLAS::gemv_trans (N, M, C.data(), &(GSLapl[N]), &(d2psi[0]));
-  }
-  // for (int i=0; i<N; i++) {
-  //   psi[i]  += GSVal[i];
-  //   d2psi[i] += GSLapl[i];
-  // }
-}
-
-void OptimizableSPOSet::evaluateBasis(const ParticleSet& P,
-                                      int first,
-                                      int last,
-                                      ValueMatrix_t& basis_val,
-                                      GradMatrix_t& basis_grad,
-                                      ValueMatrix_t& basis_lapl)
-{
-  if (BasisOrbitals)
-    BasisOrbitals->evaluate_notranspose(P, first, last, basis_val, basis_grad, basis_lapl);
-  else
-  {
-    for (int iat = first; iat < last; iat++)
-    {
-      GSOrbitals->evaluate(P, iat, GSVal, GSGrad, GSLapl);
-      for (int i = 0; i < M; i++)
-      {
-        basis_val(iat - first, i)  = GSVal[N + i];
-        basis_grad(iat - first, i) = GSGrad[N + i];
-        basis_lapl(iat - first, i) = GSLapl[N + i];
-      }
-    }
-  }
-}
-
-void OptimizableSPOSet::copyParamsFromMatrix(const opt_variables_type& active,
-                                             const Matrix<RealType>& mat,
-                                             std::vector<RealType>& destVec)
-{
-  for (int ip = 0; ip < myVars.size(); ip++)
-  {
-    int loc = myVars.where(ip);
-    if (loc >= 0)
-    {
-      TinyVector<int, 2> idx = ParamIndex[ip];
-      destVec[loc] += derivScale * mat(idx[0], idx[1]);
-    }
-  }
-}
-
-void OptimizableSPOSet::copyParamsFromMatrix(const opt_variables_type& active,
-                                             const Matrix<ComplexType>& mat,
-                                             std::vector<RealType>& destVec)
-{
-  for (int ip = 0; ip < myVars.size(); ip += 2)
-  {
-    int loc = myVars.where(ip);
-    if (loc >= 0)
-    {
-      TinyVector<int, 2> idx = ParamIndex[ip];
-      assert(ParamIndex[ip + 1][0] == idx[0] && ParamIndex[ip + 1][1] == idx[1]);
-      destVec[loc] += derivScale * mat(idx[0], idx[1]).real();
-      loc = myVars.where(ip + 1);
-      assert(loc >= 0);
-      destVec[loc] += derivScale * mat(idx[0], idx[1]).imag();
-    }
-  }
-}
-
-
-void OptimizableSPOSet::evaluate_notranspose(const ParticleSet& P,
-                                             int first,
-                                             int last,
-                                             ValueMatrix_t& logdet,
-                                             GradMatrix_t& dlogdet,
-                                             ValueMatrix_t& d2logdet)
-{
-  //     std::cerr << "GSValMatrix.size =(" << GSValMatrix.size(0) << ", " << GSValMatrix.size(1) << ")\n";
-  //     std::cerr << "GSGradMatrix.size =(" << GSGradMatrix.size(0) << ", " << GSGradMatrix.size(1) << ")\n";
-  //     std::cerr << "GSLaplMatrix.size =(" << GSLaplMatrix.size(0) << ", " << GSLaplMatrix.size(1) << ")\n";
-  //     std::cerr << "first=" << first << "  last=" << last << std::endl;
-  GSOrbitals->evaluate_notranspose(P, first, last, GSValMatrix, GSGradMatrix, GSLaplMatrix);
-  if (BasisOrbitals)
-  {
-    BasisOrbitals->evaluate_notranspose(P, first, last, BasisValMatrix, BasisGradMatrix, BasisLaplMatrix);
-    //Note to Ken:
-    //Use Numerics/MatrixOperators.h
-    //for C=AB MatrixOperators::product(C,BasisValMatrix,logdet);
-    //for C=AB^t MatrixOperators::ABt(C,BasisValMatrix,logdet);
-    BLAS::gemm('T', 'N', N, N, M, 1.0, C->data(), M, BasisValMatrix.data(), M, 0.0, logdet.data(), N);
-    logdet += GSValMatrix;
-    BLAS::gemm('T', 'N', N, N, M, 1.0, C->data(), M, BasisLaplMatrix.data(), M, 0.0, d2logdet.data(), N);
-    d2logdet += GSLaplMatrix;
-    // Gradient part.
-    for (int dim = 0; dim < OHMMS_DIM; dim++)
-    {
-      for (int i = 0; i < M; i++)
-        for (int j = 0; j < N; j++)
-          GradTmpSrc(i, j) = BasisGradMatrix(i, j)[dim];
-      BLAS::gemm('T', 'N', N, N, M, 1.0, C->data(), M, GradTmpSrc.data(), M, 0.0, GradTmpDest.data(), N);
-      for (int i = 0; i < N; i++)
-        for (int j = 0; j < N; j++)
-          dlogdet(i, j)[dim] = GradTmpDest(i, j) + GSGradMatrix(i, j)[dim];
-    }
-  }
-  else
-  {
-    // HACK HACK HACK
-    //       BLAS::gemm ('T', 'N', N, N, M, 1.0, C.data(),
-    // 		  M, GSValMatrix.data()+N, M+N, 0.0, logdet.data(), N);
-    //       for (int i=0; i<N; i++)
-    // 	for (int j=0; j<N; j++)
-    // 	  logdet(i,j) += GSValMatrix(i,j);
-    for (int iel = 0; iel < N; iel++)
-      for (int iorb = 0; iorb < N; iorb++)
-      {
-        logdet(iel, iorb)              = GSValMatrix(iel, iorb);
-        const ValueType* restrict cptr = (*C)[iorb];
-        for (int ibasis = 0; ibasis < M; ibasis++)
-          logdet(iel, iorb) += cptr[ibasis] * GSValMatrix(iel, N + ibasis);
-      }
-    //      logdet += GSValMatrix;
-    // BLAS::gemm ('T', 'N', N, N, M, 1.0, C.data(),
-    // 		  M, GSLaplMatrix.data()+N, M+N, 0.0, d2logdet.data(), N);
-    for (int iel = 0; iel < N; iel++)
-      for (int iorb = 0; iorb < N; iorb++)
-      {
-        const ValueType* restrict cptr = (*C)[iorb];
-        d2logdet(iel, iorb)            = GSLaplMatrix(iel, iorb);
-        for (int ibasis = 0; ibasis < M; ibasis++)
-          d2logdet(iel, iorb) += cptr[ibasis] * GSLaplMatrix(iel, N + ibasis);
-      }
-    //d2logdet += GSLaplMatrix;
-    // Gradient part.
-    for (int i = 0; i < N; i++)
-      for (int iorb = 0; iorb < N; iorb++)
-      {
-        const ValueType* restrict cptr = (*C)[iorb];
-        dlogdet(i, iorb)               = GSGradMatrix(i, iorb);
-        for (int n = 0; n < M; n++)
-          dlogdet(i, iorb) += cptr[n] * GSGradMatrix(i, N + n);
-      }
-    // for (int dim=0; dim<OHMMS_DIM; dim++) {
-    // 	for (int i=0; i<M; i++)
-    // 	  for (int j=0; j<N; j++)
-    // 	    GradTmpSrc(i,j) = GSGradMatrix(i,j+N)[dim];
-    // 	BLAS::gemm ('T', 'N', N, N, M, 1.0, C.data(), M,
-    // 		     GradTmpSrc.data(), M, 0.0, GradTmpDest.data(), N);
-    // 	for (int i=0; i<N; i++)
-    // 	  for (int j=0; j<N; j++)
-    // 	    dlogdet(i,j)[dim] = GradTmpDest(i,j) + GSGradMatrix(i,j)[dim];
-    //}
-  }
-}
-
-SPOSet* OptimizableSPOSet::makeClone() const
-{
-  SPOSet *gs, *basis(0);
-  OptimizableSPOSet* clone;
-  gs = GSOrbitals->makeClone();
-  if (BasisOrbitals)
-    basis = BasisOrbitals->makeClone();
-  clone             = new OptimizableSPOSet(N, gs, basis);
-  clone->IsCloned   = true;
-  clone->C          = C; //just pointer
-  clone->thr        = thr;
-  clone->myVars     = myVars;
-  clone->derivScale = derivScale;
-#if 0
-  clone->ParamPointers.clear();
-  clone->ParamIndex=ParamIndex;
-  for(int i=0; i<ParamIndex.size() ; i++)
-  {
-#ifndef QMC_COMPLEX
-    clone->ParamPointers.push_back(&(clone->C(ParamIndex[i][0],ParamIndex[i][1])));
-#else
-    int ci=ParamIndex[i][0];
-    int cj=ParamIndex[i][1];
-    //clone->ParamPointers.push_back(&(clone->C(ParamIndex[i][0],ParamIndex[i][1]).real()));
-    //clone->ParamPointers.push_back(&(clone->C(ParamIndex[i][0],ParamIndex[i][1]).imag()));
-    clone->ParamPointers.push_back(reinterpret_cast<RealType*>(clone->C[ci]+cj));
-    clone->ParamPointers.push_back(reinterpret_cast<RealType*>(clone->C[ci]+cj)+1);
-#endif
-  }
-#endif
-  //   for (int i=0; i< N; i++) {
-  //   for (int j=0; j< M; j++) {
-  //     std::stringstream sstr;
-  // #ifndef QMC_COMPLEX
-  //     clone->ParamPointers.push_back(&(clone->C(i,j)));
-  //     clone->ParamIndex.push_back(TinyVector<int,2>(i,j));
-  // #else
-  //     clone->ParamPointers.push_back(&(clone->C(i,j).real()));
-  //     clone->ParamPointers.push_back(&(clone->C(i,j).imag()));
-  //     clone->ParamIndex.push_back(TinyVector<int,2>(i,j));
-  //     clone->ParamIndex.push_back(TinyVector<int,2>(i,j));
-  // #endif
-  //   }
-  //   }
-  return clone;
-}
-
-
-} // namespace qmcplusplus
diff --git a/src/QMCWaveFunctions/OptimizableSPOSet.h b/src/QMCWaveFunctions/OptimizableSPOSet.h
deleted file mode 100644
index 0e076e402..000000000
--- a/src/QMCWaveFunctions/OptimizableSPOSet.h
+++ /dev/null
@@ -1,226 +0,0 @@
-//////////////////////////////////////////////////////////////////////////////////////
-// This file is distributed under the University of Illinois/NCSA Open Source License.
-// See LICENSE file in top directory for details.
-//
-// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
-//
-// File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
-//                    Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
-//                    Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
-//                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
-//                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
-//
-// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
-//////////////////////////////////////////////////////////////////////////////////////
-
-
-#ifndef OPTIMIZABLE_SPO_SET_H
-#define OPTIMIZABLE_SPO_SET_H
-
-#include "QMCWaveFunctions/SPOSet.h"
-#include "Numerics/OptimizableFunctorBase.h"
-#include "Optimize/VariableSet.h"
-
-namespace qmcplusplus
-{
-class OptimizableSPOSet : public SPOSet
-{
-protected:
-  typedef optimize::VariableSet opt_variables_type;
-  ///typedef for name-value lists
-  typedef optimize::VariableSet::variable_map_type variable_map_type;
-
-  // Number of occupied states, number of basis states
-  int N, M;
-  ///number of basis
-  IndexType BasisSetSize;
-  /** pointer to matrix containing the coefficients
-   *
-   * makeClone makes a shallow copy
-   */
-  ValueMatrix_t* C;
-  ///if true, do not clean up
-  bool IsCloned;
-
-  RealType derivScale;
-  RealType thr;
-
-  // If BasisOrbitals==NULL, only GSOrbitals is used and it's evaluate
-  // functions should return N+M orbitals.
-  SPOSet *GSOrbitals, *BasisOrbitals;
-
-  // The first index is the orbital to be optimized, the second is the
-  // basis element
-  //    ValueMatrix_t OptCoefs;
-
-  // This maps the parameter indices to elements in the C
-  // matrix.  Using pointers allows the use of complex C
-  // while maintaining real optimizable parameters.
-  std::vector<RealType*> ParamPointers;
-  // Maps the parameter index in myVars to an index in the C array
-  std::vector<TinyVector<int, 2>> ParamIndex;
-
-  void addParameter(std::string id, int iorb, int basis);
-
-  ValueVector_t GSVal, BasisVal, GSLapl, BasisLapl;
-  GradVector_t GSGrad, BasisGrad;
-
-  ValueMatrix_t GSValMatrix, BasisValMatrix, GSLaplMatrix, BasisLaplMatrix, GradTmpSrc, GradTmpDest;
-  GradMatrix_t GSGradMatrix, BasisGradMatrix;
-
-  // Cache the positions to avoid recomputing GSVal, BasisVal, etc.
-  // if unnecessary.
-  std::vector<PosType> CachedPos;
-
-public:
-  ///set of variables to be optimized;  These are mapped to the
-  ///C matrix.  Update:  Moved to SPOSet
-  // opt_variables_type myVars;
-
-  // For each occupied orbital, this lists which of the M
-  // basis functions are used to construct the optimal orbital.
-  // This must be initialized by the SPOSet builder.  If it is not
-  // initialized, we assume that C is fully dense.
-  // The first element of the TinyVector is the basis element
-  // number.  The second is the corresponding parameter number.
-  std::vector<std::vector<TinyVector<int, 2>>> ActiveBasis;
-
-
-  OptimizableSPOSet() : N(0), M(0), derivScale(10.0), thr(0.0), GSOrbitals(0), BasisOrbitals(0)
-  {
-    Optimizable = true;
-    className   = "OptimizableSPOSet";
-  }
-
-  OptimizableSPOSet(int num_orbs, SPOSet* gsOrbs, SPOSet* basisOrbs = 0)
-      : GSOrbitals(gsOrbs), BasisOrbitals(basisOrbs), derivScale(10.0)
-  {
-    N = num_orbs;
-    setOrbitalSetSize(N);
-    if (BasisOrbitals)
-      M = BasisOrbitals->getOrbitalSetSize();
-    else
-      M = GSOrbitals->getOrbitalSetSize() - N;
-    resize(N, M);
-    Optimizable = true;
-  }
-
-
-  //    bool put(xmlNodePtr cur, SPOPool_t &spo_pool);
-  void resetTargetParticleSet(ParticleSet& P);
-  void setOrbitalSetSize(int norbs);
-
-  void set_active_basis(std::vector<std::vector<int>>& active)
-  {
-    //ActiveBasis = active;
-  }
-
-  // Read coefficients, or create the XML element.
-  bool put(xmlNodePtr node, SPOPool_t& spo_pool);
-
-  // This stores new orbital coefficients int C
-  void resetParameters(const opt_variables_type& optvars);
-
-  // Evaluate the derivative of the optimized orbitals with
-  // respect to the parameters
-  void evaluateDerivatives(ParticleSet& P,
-                           int iat,
-                           const opt_variables_type& active,
-                           ValueMatrix_t& d_phi,
-                           ValueMatrix_t& d_lapl_phi);
-
-  void checkInVariables(opt_variables_type& active);
-  void checkOutVariables(const opt_variables_type& active);
-
-
-  // Evaluate routines.  These call GSOrbitals->evaluate and possibly
-  // BasisOrbitals->evaluate, then does the matrix product with
-  // C.
-  void evaluate(const ParticleSet& P, int iat, ValueVector_t& psi);
-  void evaluate(const ParticleSet& P, int iat, ValueVector_t& psi, GradVector_t& dpsi, ValueVector_t& d2psi);
-  void evaluate_notranspose(const ParticleSet& P,
-                            int first,
-                            int last,
-                            ValueMatrix_t& logdet,
-                            GradMatrix_t& dlogdet,
-                            ValueMatrix_t& d2logdet);
-  void evaluate_notranspose(const ParticleSet& P,
-                            int first,
-                            int last,
-                            ValueMatrix_t& logdet,
-                            GradMatrix_t& dlogdet,
-                            HessMatrix_t& grad_grad_logdet)
-  {
-    APP_ABORT("Need specialization of OptimizableOrbitalSet::evaluate_notranspose() for grad_grad_logdet. \n");
-  }
-
-  void evaluate_notranspose(const ParticleSet& P,
-                            int first,
-                            int last,
-                            ValueMatrix_t& logdet,
-                            GradMatrix_t& dlogdet,
-                            HessMatrix_t& grad_grad_logdet,
-                            GGGMatrix_t& grad_grad_grad_logdet)
-  {
-    APP_ABORT("Need specialization of OptimazableOrbitalSet::evaluate_notranspose() for grad_grad_grad_logdet. \n");
-  }
-
-  void evaluateBasis(const ParticleSet& P,
-                     int first,
-                     int last,
-                     ValueMatrix_t& basis_val,
-                     GradMatrix_t& basis_grad,
-                     ValueMatrix_t& basis_lapl);
-  void copyParamsFromMatrix(const opt_variables_type& active,
-                            const Matrix<RealType>& mat,
-                            std::vector<RealType>& destVec);
-  void copyParamsFromMatrix(const opt_variables_type& active,
-                            const Matrix<ComplexType>& mat,
-                            std::vector<RealType>& destVec);
-
-
-  void resize(int n, int m)
-  {
-    N = n;
-    M = m;
-    if (BasisOrbitals)
-    {
-      GSVal.resize(N);
-      GSGrad.resize(N);
-      GSLapl.resize(N);
-      BasisVal.resize(M);
-      BasisGrad.resize(M);
-      BasisLapl.resize(M);
-      GSValMatrix.resize(N, N);
-      GSGradMatrix.resize(N, N);
-      GSLaplMatrix.resize(N, N);
-      BasisValMatrix.resize(M, N);
-      BasisGradMatrix.resize(M, N);
-      BasisLaplMatrix.resize(M, N);
-    }
-    else
-    {
-      GSVal.resize(N + M);
-      GSGrad.resize(N + M);
-      GSLapl.resize(N + M);
-      GSValMatrix.resize(N, N + M);
-      GSGradMatrix.resize(N, N + M);
-      GSLaplMatrix.resize(N, N + M);
-    }
-    GradTmpSrc.resize(M, N);
-    GradTmpDest.resize(N, N);
-
-    if (C == nullptr)
-      C = new ValueMatrix_t(N, M);
-    else
-      C->resize(N, M);
-    ActiveBasis.resize(N);
-    BasisSetSize = M;
-  }
-
-  // Make a copy of myself
-  SPOSet* makeClone() const;
-};
-} // namespace qmcplusplus
-
-#endif
diff --git a/src/QMCWaveFunctions/SPOSet.cpp b/src/QMCWaveFunctions/SPOSet.cpp
index 7b69f56c6..8728f66b4 100644
--- a/src/QMCWaveFunctions/SPOSet.cpp
+++ b/src/QMCWaveFunctions/SPOSet.cpp
@@ -26,11 +26,10 @@
 
 namespace qmcplusplus
 {
-SPOSet::SPOSet()
+SPOSet::SPOSet(bool ion_deriv, bool optimizable)
     : OrbitalSetSize(0),
-      Optimizable(false),
-      ionDerivs(false),
-      builder_index(-1)
+      ionDerivs(ion_deriv),
+      Optimizable(optimizable)
 #if !defined(ENABLE_SOA)
       ,
       Identity(false),
diff --git a/src/QMCWaveFunctions/SPOSet.h b/src/QMCWaveFunctions/SPOSet.h
index f6490cf19..43d6c857c 100644
--- a/src/QMCWaveFunctions/SPOSet.h
+++ b/src/QMCWaveFunctions/SPOSet.h
@@ -62,16 +62,6 @@ public:
   typedef ParticleSet::Walker_t Walker_t;
   typedef std::map<std::string, SPOSet*> SPOPool_t;
 
-  ///index in the builder list of sposets
-  int builder_index;
-  ///true if SPO is optimizable
-  bool Optimizable;
-  ///number of Single-particle orbitals
-  IndexType OrbitalSetSize;
-  /// Optimizable variables
-  opt_variables_type myVars;
-  ///name of the class
-  std::string className;
   /** name of the object
    *
    * Several user classes can own SPOSet and use objectName as counter
@@ -96,7 +86,7 @@ public:
 #endif
 
   /** constructor */
-  SPOSet();
+  SPOSet(bool ion_deriv = false, bool optimizable = false);
 
   /** destructor
    *
@@ -110,6 +100,9 @@ public:
 #endif
   }
 
+  // accessor function to Optimizable
+  inline bool isOptimizable() const { return Optimizable; }
+
   /** return the size of the orbital set
    * Ye: this needs to be replaced by getOrbitalSetSize();
    */
@@ -472,8 +465,18 @@ protected:
   bool putFromXML(xmlNodePtr coeff_ptr);
   bool putFromH5(const std::string& fname, xmlNodePtr coeff_ptr);
 #endif
+
+protected:
   ///true, if the derived class has non-zero ionic derivatives.
-  bool ionDerivs;
+  const bool ionDerivs;
+  ///true if SPO is optimizable
+  const bool Optimizable;
+  ///number of Single-particle orbitals
+  IndexType OrbitalSetSize;
+  /// Optimizable variables
+  opt_variables_type myVars;
+  ///name of the class
+  std::string className;
 };
 
 #if defined(ENABLE_SMARTPOINTER)
diff --git a/src/QMCWaveFunctions/SPOSetBuilder.cpp b/src/QMCWaveFunctions/SPOSetBuilder.cpp
index 28e12caf6..f98b0232f 100644
--- a/src/QMCWaveFunctions/SPOSetBuilder.cpp
+++ b/src/QMCWaveFunctions/SPOSetBuilder.cpp
@@ -53,7 +53,6 @@ SPOSet* SPOSetBuilder::createSPOSet(xmlNodePtr cur)
   if (sposet)
   {
     //sposet->put(cur); //initialize C and other internal containers
-    sposet->builder_index = sposets.size();
     sposets.push_back(sposet);
   }
   else
diff --git a/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp b/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp
index 41e63df0a..4aecc642a 100644
--- a/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp
+++ b/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp
@@ -36,7 +36,6 @@
 #endif
 #endif
 #include "QMCWaveFunctions/CompositeSPOSet.h"
-#include "QMCWaveFunctions/OptimizableSPOBuilder.h"
 #include "Utilities/ProgressReportEngine.h"
 #include "Utilities/IteratorUtility.h"
 #include "OhmmsData/AttributeSet.h"
@@ -181,11 +180,6 @@ SPOSetBuilder* SPOSetBuilderFactory::createSPOSetBuilder(xmlNodePtr rootNode)
     app_log() << "Harmonic Oscillator SPO set" << std::endl;
     bb = new SHOSetBuilder(targetPtcl, myComm);
   }
-  else if (type == "linearopt")
-  {
-    //app_log()<<"Optimizable SPO set"<< std::endl;
-    bb = new OptimizableSPOBuilder(targetPtcl, ptclPool, myComm, rootNode);
-  }
 #if OHMMS_DIM == 3
   else if (type.find("spline") < type.size())
   {
diff --git a/src/QMCWaveFunctions/lcao/CuspCorrection.cpp b/src/QMCWaveFunctions/lcao/CuspCorrection.cpp
index 552e64227..9ca9c52fb 100644
--- a/src/QMCWaveFunctions/lcao/CuspCorrection.cpp
+++ b/src/QMCWaveFunctions/lcao/CuspCorrection.cpp
@@ -163,8 +163,8 @@ void splitPhiEta(int center, const std::vector<bool>& corrCenter, LCAOrbitalSet&
 
   Phi.myBasisSet->queryOrbitalsForSType(correct_this_center, is_s_orbital);
 
-  int nOrbs = Phi.OrbitalSetSize;
-  int bss   = Phi.BasisSetSize;
+  int nOrbs = Phi.getOrbitalSetSize();
+  int bss   = Phi.getBasisSetSize();
 
   for (int i = 0; i < bss; i++)
   {
@@ -191,8 +191,8 @@ void removeSTypeOrbitals(const std::vector<bool>& corrCenter, LCAOrbitalSet& Phi
 
   Phi.myBasisSet->queryOrbitalsForSType(corrCenter, is_s_orbital);
 
-  int nOrbs = Phi.OrbitalSetSize;
-  int bss   = Phi.BasisSetSize;
+  int nOrbs = Phi.getOrbitalSetSize();
+  int bss   = Phi.getBasisSetSize();
 
   for (int i = 0; i < bss; i++)
   {
diff --git a/src/QMCWaveFunctions/lcao/CuspCorrection.h b/src/QMCWaveFunctions/lcao/CuspCorrection.h
index efa7246ec..2199d144b 100644
--- a/src/QMCWaveFunctions/lcao/CuspCorrection.h
+++ b/src/QMCWaveFunctions/lcao/CuspCorrection.h
@@ -114,7 +114,7 @@ public:
       : targetPtcl(targetP), sourcePtcl(sourceP), curOrb(0), curCenter(0)
   {
     Psi1     = Phi;
-    int norb = Psi1->OrbitalSetSize;
+    int norb = Psi1->getOrbitalSetSize();
     val1.resize(norb);
     grad1.resize(norb);
     lap1.resize(norb);
diff --git a/src/QMCWaveFunctions/lcao/CuspCorrectionConstruction.cpp b/src/QMCWaveFunctions/lcao/CuspCorrectionConstruction.cpp
index 0c9377c1a..20ec184a3 100644
--- a/src/QMCWaveFunctions/lcao/CuspCorrectionConstruction.cpp
+++ b/src/QMCWaveFunctions/lcao/CuspCorrectionConstruction.cpp
@@ -33,14 +33,14 @@ void applyCuspCorrection(const Matrix<CuspCorrectionParameters>& info,
 
   ScopedTimer cuspApplyTimerWrapper(cuspApplyTimer);
 
-  LCAOrbitalSet phi = LCAOrbitalSet(lcwc.myBasisSet);
-  phi.setOrbitalSetSize(lcwc.OrbitalSetSize);
-  phi.BasisSetSize = lcwc.BasisSetSize;
+  LCAOrbitalSet phi = LCAOrbitalSet(lcwc.myBasisSet, lcwc.isOptimizable());
+  phi.setOrbitalSetSize(lcwc.getOrbitalSetSize());
+  phi.BasisSetSize = lcwc.getBasisSetSize();
   phi.setIdentity(false);
 
-  LCAOrbitalSet eta = LCAOrbitalSet(lcwc.myBasisSet);
-  eta.setOrbitalSetSize(lcwc.OrbitalSetSize);
-  eta.BasisSetSize = lcwc.BasisSetSize;
+  LCAOrbitalSet eta = LCAOrbitalSet(lcwc.myBasisSet, lcwc.isOptimizable());
+  eta.setOrbitalSetSize(lcwc.getOrbitalSetSize());
+  eta.BasisSetSize = lcwc.getBasisSetSize();
   eta.setIdentity(false);
 
 
@@ -201,14 +201,14 @@ void generateCuspInfo(int orbital_set_size,
 
   ScopedTimer createCuspTimerWrapper(cuspCreateTimer);
 
-  LCAOrbitalSet phi = LCAOrbitalSet(lcwc.myBasisSet);
-  phi.setOrbitalSetSize(lcwc.OrbitalSetSize);
-  phi.BasisSetSize = lcwc.BasisSetSize;
+  LCAOrbitalSet phi = LCAOrbitalSet(lcwc.myBasisSet, lcwc.isOptimizable());
+  phi.setOrbitalSetSize(lcwc.getOrbitalSetSize());
+  phi.BasisSetSize = lcwc.getBasisSetSize();
   phi.setIdentity(false);
 
-  LCAOrbitalSet eta = LCAOrbitalSet(lcwc.myBasisSet);
-  eta.setOrbitalSetSize(lcwc.OrbitalSetSize);
-  eta.BasisSetSize = lcwc.BasisSetSize;
+  LCAOrbitalSet eta = LCAOrbitalSet(lcwc.myBasisSet, lcwc.isOptimizable());
+  eta.setOrbitalSetSize(lcwc.getOrbitalSetSize());
+  eta.BasisSetSize = lcwc.getBasisSetSize();
   eta.setIdentity(false);
 
 
diff --git a/src/QMCWaveFunctions/lcao/LCAOrbitalBuilder.cpp b/src/QMCWaveFunctions/lcao/LCAOrbitalBuilder.cpp
index 167044edb..1d2c01289 100644
--- a/src/QMCWaveFunctions/lcao/LCAOrbitalBuilder.cpp
+++ b/src/QMCWaveFunctions/lcao/LCAOrbitalBuilder.cpp
@@ -443,15 +443,19 @@ SPOSet* LCAOrbitalBuilder::createSPOSetFromXML(xmlNodePtr cur)
 
   if (myBasisSet == nullptr)
     PRE.error("Missing basisset.", true);
+
+  if (optimize == "yes")
+    app_log() << "  SPOSet " << spo_name << " is optimizable\n";
+
   LCAOrbitalSet* lcos = nullptr;
 #if !defined(QMC_COMPLEX)
   LCAOrbitalSetWithCorrection* lcwc = nullptr;
   if (doCuspCorrection)
-    lcos = lcwc = new LCAOrbitalSetWithCorrection(sourcePtcl, targetPtcl, myBasisSet);
+    lcos = lcwc = new LCAOrbitalSetWithCorrection(sourcePtcl, targetPtcl, myBasisSet, optimize == "yes");
   else
-    lcos = new LCAOrbitalSet(myBasisSet);
+    lcos = new LCAOrbitalSet(myBasisSet, optimize == "yes");
 #else
-  lcos = new LCAOrbitalSet(myBasisSet);
+  lcos = new LCAOrbitalSet(myBasisSet, optimize == "yes");
 #endif
   loadMO(*lcos, cur);
 
@@ -464,7 +468,7 @@ SPOSet* LCAOrbitalBuilder::createSPOSetFromXML(xmlNodePtr cur)
     if (id == "")
       id = spo_name;
 
-    int orbital_set_size = lcos->OrbitalSetSize;
+    int orbital_set_size = lcos->getOrbitalSetSize();
     Matrix<CuspCorrectionParameters> info(num_centers, orbital_set_size);
 
     bool valid = false;
@@ -494,12 +498,6 @@ SPOSet* LCAOrbitalBuilder::createSPOSetFromXML(xmlNodePtr cur)
   }
 #endif
 
-  if (optimize == "yes")
-  {
-    lcos->Optimizable = true;
-    app_log() << "  SPOSet " << spo_name << " is optimizable\n";
-  }
-
   return lcos;
 }
 
diff --git a/src/QMCWaveFunctions/lcao/LCAOrbitalSet.cpp b/src/QMCWaveFunctions/lcao/LCAOrbitalSet.cpp
index 9a844d34b..62cfac296 100644
--- a/src/QMCWaveFunctions/lcao/LCAOrbitalSet.cpp
+++ b/src/QMCWaveFunctions/lcao/LCAOrbitalSet.cpp
@@ -16,11 +16,9 @@
 
 namespace qmcplusplus
 {
-LCAOrbitalSet::LCAOrbitalSet(basis_type* bs)
-    : myBasisSet(nullptr), params_supplied(false), BasisSetSize(0), Identity(true)
+LCAOrbitalSet::LCAOrbitalSet(basis_type* bs, bool optimize)
+    : SPOSet(true, optimize), myBasisSet(nullptr), params_supplied(false), BasisSetSize(0), Identity(true)
 {
-  //This SPOSet has an explicit ion dependence, so set this flag.
-  ionDerivs=true;
   if (bs != nullptr)
     setBasisSet(bs);
 }
diff --git a/src/QMCWaveFunctions/lcao/LCAOrbitalSet.h b/src/QMCWaveFunctions/lcao/LCAOrbitalSet.h
index db198675e..7cdfedf78 100644
--- a/src/QMCWaveFunctions/lcao/LCAOrbitalSet.h
+++ b/src/QMCWaveFunctions/lcao/LCAOrbitalSet.h
@@ -73,7 +73,7 @@ public:
   /** constructor
      * @param bs pointer to the BasisSet
      */
-  LCAOrbitalSet(basis_type* bs = nullptr);
+  LCAOrbitalSet(basis_type* bs, bool optimize);
 
   LCAOrbitalSet(const LCAOrbitalSet& in) = default;
 
@@ -109,16 +109,11 @@ public:
                            const size_t NP2,
                            const std::vector<std::vector<int>>& lookup_tbl);
 
-
   void checkInVariables(opt_variables_type& active)
   {
     if (Optimizable)
-    {
       if (myVars.size())
         active.insertFrom(myVars);
-      else
-        Optimizable = false;
-    }
   }
 
   void checkOutVariables(const opt_variables_type& active)
@@ -127,7 +122,6 @@ public:
       myVars.getIndex(active);
   }
 
-
   ///reset
   void resetParameters(const opt_variables_type& active)
   {
diff --git a/src/QMCWaveFunctions/lcao/LCAOrbitalSetWithCorrection.cpp b/src/QMCWaveFunctions/lcao/LCAOrbitalSetWithCorrection.cpp
index 39fe92cd6..f29b362c6 100644
--- a/src/QMCWaveFunctions/lcao/LCAOrbitalSetWithCorrection.cpp
+++ b/src/QMCWaveFunctions/lcao/LCAOrbitalSetWithCorrection.cpp
@@ -15,8 +15,8 @@
 
 namespace qmcplusplus
 {
-LCAOrbitalSetWithCorrection::LCAOrbitalSetWithCorrection(ParticleSet& ions, ParticleSet& els, basis_type* bs)
-    : LCAOrbitalSet(bs), cusp(ions, els)
+LCAOrbitalSetWithCorrection::LCAOrbitalSetWithCorrection(ParticleSet& ions, ParticleSet& els, basis_type* bs, bool optimize)
+    : LCAOrbitalSet(bs, optimize), cusp(ions, els)
 {}
 
 void LCAOrbitalSetWithCorrection::setOrbitalSetSize(int norbs)
diff --git a/src/QMCWaveFunctions/lcao/LCAOrbitalSetWithCorrection.h b/src/QMCWaveFunctions/lcao/LCAOrbitalSetWithCorrection.h
index 2d3432a31..55b76508c 100644
--- a/src/QMCWaveFunctions/lcao/LCAOrbitalSetWithCorrection.h
+++ b/src/QMCWaveFunctions/lcao/LCAOrbitalSetWithCorrection.h
@@ -34,7 +34,7 @@ struct LCAOrbitalSetWithCorrection : public LCAOrbitalSet
      * @param bs pointer to the BasisSet
      * @param rl report level
      */
-  LCAOrbitalSetWithCorrection(ParticleSet& ions, ParticleSet& els, basis_type* bs = nullptr);
+  LCAOrbitalSetWithCorrection(ParticleSet& ions, ParticleSet& els, basis_type* bs, bool optimize);
 
   LCAOrbitalSetWithCorrection(const LCAOrbitalSetWithCorrection& in) = default;
 
diff --git a/src/QMCWaveFunctions/tests/test_cusp_corr.cpp b/src/QMCWaveFunctions/tests/test_cusp_corr.cpp
index b5b0ace28..8a674c12e 100644
--- a/src/QMCWaveFunctions/tests/test_cusp_corr.cpp
+++ b/src/QMCWaveFunctions/tests/test_cusp_corr.cpp
@@ -108,14 +108,14 @@ TEST_CASE("CuspCorrection He", "[wavefunction]")
 
   typedef NGOBuilder::CenteredOrbitalType COT;
   OrbType bs_phi(lcob->myBasisSet);
-  bs_phi.setOrbitalSetSize(lcob->OrbitalSetSize);
+  bs_phi.setOrbitalSetSize(lcob->getOrbitalSetSize());
   bs_phi.BasisSetSize = lcob->BasisSetSize;
   bs_phi.setIdentity(false);
 
   *(bs_phi.C) = *(lcob->C);
 
   OrbType bs_eta(lcob->myBasisSet);
-  bs_eta.setOrbitalSetSize(lcob->OrbitalSetSize);
+  bs_eta.setOrbitalSetSize(lcob->getOrbitalSetSize());
   bs_eta.BasisSetSize = lcob->BasisSetSize;
   bs_eta.setIdentity(false);
   *(bs_eta.C) = *(lcob->C);
diff --git a/src/QMCWaveFunctions/tests/test_soa_cusp_corr.cpp b/src/QMCWaveFunctions/tests/test_soa_cusp_corr.cpp
index 2bd88d073..211ab2bc6 100644
--- a/src/QMCWaveFunctions/tests/test_soa_cusp_corr.cpp
+++ b/src/QMCWaveFunctions/tests/test_soa_cusp_corr.cpp
@@ -129,14 +129,14 @@ TEST_CASE("applyCuspInfo", "[wavefunction]")
   REQUIRE(lcob != NULL);
 
 
-  LCAOrbitalSet phi = LCAOrbitalSet(lcob->myBasisSet);
-  phi.setOrbitalSetSize(lcob->OrbitalSetSize);
-  phi.BasisSetSize = lcob->BasisSetSize;
+  LCAOrbitalSet phi = LCAOrbitalSet(lcob->myBasisSet, lcob->isOptimizable());
+  phi.setOrbitalSetSize(lcob->getOrbitalSetSize());
+  phi.BasisSetSize = lcob->getBasisSetSize();
   phi.setIdentity(false);
 
-  LCAOrbitalSet eta = LCAOrbitalSet(lcob->myBasisSet);
-  eta.setOrbitalSetSize(lcob->OrbitalSetSize);
-  eta.BasisSetSize = lcob->BasisSetSize;
+  LCAOrbitalSet eta = LCAOrbitalSet(lcob->myBasisSet, lcob->isOptimizable());
+  eta.setOrbitalSetSize(lcob->getOrbitalSetSize());
+  eta.BasisSetSize = lcob->getBasisSetSize();
   eta.setIdentity(false);
 
   *(eta.C) = *(lcob->C);
diff --git a/tests/molecules/He_ae/CMakeLists.txt b/tests/molecules/He_ae/CMakeLists.txt
index a1bfba87a..7b0c5d0dc 100644
--- a/tests/molecules/He_ae/CMakeLists.txt
+++ b/tests/molecules/He_ae/CMakeLists.txt
@@ -1,6 +1,8 @@
 # Deterministic test for QMC optimizer 
 
-IF (NOT MIXED_PRECISION)
+IF (MIXED_PRECISION OR (QMC_COMPLEX AND NOT ENABLE_SOA))
+  MESSAGE_VERBOSE("Skipping optimizer deterministic He tests in mixed precision (QMC_MIXED_PRECISION=1) or complex with aos (QMC_COMPLEX=1 and ENABLE_SOA=0)")
+ELSE()
   SIMPLE_RUN_AND_CHECK(
     deterministic-He_ae-opt
     "${CMAKE_SOURCE_DIR}/tests/molecules/He_ae"
@@ -19,7 +21,5 @@ IF (NOT MIXED_PRECISION)
                     TRUE
                     2 He_VMC_SCALARS # VMC
                     )
-ELSE()
-  MESSAGE_VERBOSE("Skipping optimizer deterministic tests in mixed precision (QMC_MIXED_PRECISION=1)")
 ENDIF()
 
diff --git a/tests/test_automation/jenkins_oxygen_cpu.sh b/tests/test_automation/jenkins_oxygen_cpu.sh
index 8158e283b..c5a6e5fd1 100755
--- a/tests/test_automation/jenkins_oxygen_cpu.sh
+++ b/tests/test_automation/jenkins_oxygen_cpu.sh
@@ -29,10 +29,12 @@ mkdir -p /dev/shm/${BUILD_TAG}-build
 cd /dev/shm/${BUILD_TAG}-build
 
 cmake -DQMC_COMPLEX=0 -DQMC_MIXED_PRECISION=0 -DENABLE_SOA=0 -DBUILD_AFQMC=1 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}
@@ -51,10 +53,12 @@ mkdir -p ${BUILD_TAG}-build
 cd ${BUILD_TAG}-build
 
 cmake -DQMC_COMPLEX=0 -DQMC_MIXED_PRECISION=1 -DBUILD_AFQMC=1 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}
@@ -73,10 +77,12 @@ mkdir -p ${BUILD_TAG}-build
 cd ${BUILD_TAG}-build
 
 cmake -DQMC_COMPLEX=1 -DQMC_MIXED_PRECISION=0 -DENABLE_SOA=0 -DBUILD_AFQMC=1 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}
@@ -95,10 +101,12 @@ mkdir -p ${BUILD_TAG}-build
 cd ${BUILD_TAG}-build
 
 cmake -DQMC_COMPLEX=1 -DQMC_MIXED_PRECISION=1 -DBUILD_AFQMC=1 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}
diff --git a/tests/test_automation/jenkins_oxygen_cuda.sh b/tests/test_automation/jenkins_oxygen_cuda.sh
index 7d2f2705d..a6376c6a3 100755
--- a/tests/test_automation/jenkins_oxygen_cuda.sh
+++ b/tests/test_automation/jenkins_oxygen_cuda.sh
@@ -30,10 +30,12 @@ mkdir -p /dev/shm/${BUILD_TAG}-build
 cd /dev/shm/${BUILD_TAG}-build
 
 time cmake -DQMC_COMPLEX=0 -DQMC_MIXED_PRECISION=0 -DENABLE_SOA=0 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" -DENABLE_CUDA=1 /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}
@@ -52,10 +54,12 @@ mkdir -p ${BUILD_TAG}-build
 cd ${BUILD_TAG}-build
 
 time cmake -DQMC_COMPLEX=0 -DQMC_MIXED_PRECISION=1 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" -DENABLE_CUDA=1 /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}
@@ -74,10 +78,12 @@ mkdir -p ${BUILD_TAG}-build
 cd ${BUILD_TAG}-build
 
 time cmake -DQMC_COMPLEX=1 -DQMC_MIXED_PRECISION=0 -DENABLE_SOA=0 -DBUILD_AFQMC=1 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" -DENABLE_CUDA=1 /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}
@@ -96,10 +102,12 @@ mkdir -p ${BUILD_TAG}-build
 cd ${BUILD_TAG}-build
 
 time cmake -DQMC_COMPLEX=1 -DQMC_MIXED_PRECISION=1 -DBUILD_AFQMC=1 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" -DENABLE_CUDA=1 /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}
diff --git a/tests/test_automation/jenkins_oxygen_gpu.sh b/tests/test_automation/jenkins_oxygen_gpu.sh
index b39b02f5f..3d798a84b 100755
--- a/tests/test_automation/jenkins_oxygen_gpu.sh
+++ b/tests/test_automation/jenkins_oxygen_gpu.sh
@@ -30,10 +30,12 @@ mkdir -p /dev/shm/${BUILD_TAG}-build
 cd /dev/shm/${BUILD_TAG}-build
 
 time cmake -DQMC_COMPLEX=0 -DQMC_MIXED_PRECISION=0 -DENABLE_SOA=0 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" -DQMC_CUDA=1 /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}
@@ -52,10 +54,12 @@ mkdir -p ${BUILD_TAG}-build
 cd ${BUILD_TAG}-build
 
 time cmake -DQMC_COMPLEX=0 -DQMC_MIXED_PRECISION=1 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" -DQMC_CUDA=1 /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}
@@ -74,10 +78,12 @@ mkdir -p ${BUILD_TAG}-build
 cd ${BUILD_TAG}-build
 
 time cmake -DQMC_COMPLEX=1 -DQMC_MIXED_PRECISION=0 -DENABLE_SOA=0 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" -DQMC_CUDA=1 /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}
@@ -96,10 +102,12 @@ mkdir -p ${BUILD_TAG}-build
 cd ${BUILD_TAG}-build
 
 time cmake -DQMC_COMPLEX=1 -DQMC_MIXED_PRECISION=1 -DCMAKE_C_COMPILER="mpicc" -DCMAKE_CXX_COMPILER="mpicxx" -DQMC_CUDA=1 /dev/shm/${BUILD_TAG}-src 2>&1 | tee cmake.out
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
 make -j 8
-ctest -L unit --output-on-failure
+if [[ $? -ne 0 ]] ; then exit 1 ; fi
 
+ctest -L unit --output-on-failure
 ret=$?
 if [[ ${ret} -ne 0 ]] ; then
   exit_code=${ret}