qmcpack_input.py

##################################################################
##  (c) Copyright 2015-  by Jaron T. Krogel                     ##
##################################################################


#====================================================================#
#  qmcpack_input.py                                                  #
#    Supports I/O and manipulation of QMCPACK's xml input file.      #
#                                                                    #
#  Content summary:                                                  #
#    QmcpackInput                                                    #
#      SimulationInput class for QMCPACK.                            #
#      Represents the QMCPACK input file as a nested collection of   #
#        objects mirroring the structure of the XML input.           #
#      XML attributes and <parameter/> elements are joined into      #
#        a common keyword representation.                            #
#      Full tranlations of test QMCPACK input files into executable  #
#        Python code can be found at the end of this file for        #
#        reference.                                                  #
#                                                                    #
#    BundledQmcpackInput                                             #
#      Class represents QMCPACK input when provided in 'bundled'     #
#        form, i.e. an input file containing a list of XML input     #
#        files.                                                      #
#      A BundledQmcpackInput object contains many QmcpackInput       #
#        objects.                                                    #
#                                                                    #
#    QmcpackInputTemplate                                            #
#      Class supports Nexus interaction with text input files        #
#        provided by users as 'templates'.                           #
#      Users can mark keywords in a template input file, generate    #
#        variations on that input file through this class, and       #
#        then use the resulting input files in Nexus workflows.      #
#                                                                    #
#    generate_qmcpack_input                                          #
#      User-facing function to compose specific QMCPACK input files. #
#      Though any QMCPACK input file can be composed directly with   #
#        Python code, only a more limited selection can be           #
#        generated directly.                                         #
#      Calls many other functions to generate specific XML elements: #
#        generate_simulationcell                                     #
#        generate_particlesets                                       #
#        generate_sposets                                            #
#        generate_sposet_builder                                     #
#        generate_bspline_builder                                    #
#        generate_heg_builder                                        #
#        generate_determinantset                                     #
#        generate_determinantset_old                                 #
#        generate_hamiltonian                                        #
#        generate_jastrows                                           #
#        generate_jastrow                                            #
#        generate_jastrow1                                           #
#        generate_bspline_jastrow2                                   #
#        generate_pade_jastrow2                                      #
#        generate_jastrow2                                           #
#        generate_jastrow3                                           #
#        generate_kspace_jastrow                                     #
#        generate_opt                                                #
#        generate_opts                                               #
#                                                                    #
#   QIxml, Names                                                     #
#     Class represents a generic XML element.                        #
#     Supports read/write and setting default values.                #
#     Derived classes represent specific elements and contain        #
#       a keyword specification for each element.                    #
#     Support for new XML elements is enabled by creating            #
#       a corresponding class with the allowed keywords, etc.,       #
#       and then adding it to the global 'classes' list below.       #
#     See classes simulation, project, application, random, include, #
#       mcwalkerset, qmcsystem, simulationcell, particleset, group,  #
#       sposet, bspline_builder, heg_builder, composite_builder,     #
#       wavefunction, determinantset, basisset, grid, atomicbasisset,# 
#       basisgroup, radfunc, slaterdeterminant, determinant,         #
#       occupation, multideterminant, detlist, ci, jastrow1,         #
#       jastrow2, jastrow3, correlation, var, coefficients,          #
#       coefficient, hamiltonian, coulomb, constant, pseudopotential,# 
#       pseudo, mpc, localenergy, energydensity, reference_points,   #
#       spacegrid, origin, axis, chiesa, density, nearestneighbors,  #
#       neighbor_trace, dm1b, spindensity, structurefactor, init,    #
#       scalar_traces, array_traces, particle_traces, traces, loop,  #
#       linear, cslinear, vmc, dmc.                                  #
#                                                                    #
#   QIxmlFactory                                                     #
#     Class supports comprehension of XML elements that share the    #
#       same XML tag (e.g. <pairpot/>), but have many different      #
#       and distinct realizations (e.g. coulomb, mpc, etc.).         #
#     See QIxmlFactory instances pairpot, estimator,                 #
#       sposet_builder, jastrow, and qmc.                            #
#                                                                    #
#   collection                                                       #
#     Container class representing an ordered set of plural XML      #
#       elements named according to an XML attribute (e.g. named     #
#       particlesets).                                               #
#     XML elements are listed by name in the collection to allow     #
#       intuitive interactive navigation by the user.                #
#                                                                    #
#   Param (single global instance is param)                          #
#     Handles all reads/writes of attributes and text contents       #
#       of XML elements.                                             #
#     Automatically distinguishes string, int, float, and array      #
#       data of any of these types.                                  #
#                                                                    #
#   Note: those seeking to add new XML elements should be aware      #
#     of the global data structures described below.                 #
#                                                                    #
#   classes                                                          #
#     Global list of QIxml classes.                                  #
#     Each new QIxml class must be added here.                       #
#                                                                    #
#   types                                                            #
#     Global dict of elements that are actually simple types (to be  #
#     interpreted the same way as <parameter/> elements) and also    #
#     factory instances.                                             #
#                                                                    #
#   plurals                                                          #
#     Global obj of elements that can be plural (i.e. multiple can   #
#     be present with the same tag.  A QmcpackInput object will      #
#     contain collection instances containing the related XML        #
#     elements.  Each collection is named according to the keys in   #
#     plurals.                                                       #
#                                                                    #
#   Names.set_expanded_names                                         #
#     Function to specify mappings from all-lowercase names used     #
#     in QmcpackInput objects to the names expected by QMCPACK in    #
#     the actual XML input file.  QMCPACK does not follow any        #
#     consistent naming convention (camel-case, hyphenation, and     #
#     separation via underscore are all present in names) and some   #
#     names are case-sensitive while others are not.  This function  #
#     allows developers to cope with this diversity upon write,      #
#     while preserving a uniform representation (all lowercase) in   #
#     QmcpackInput objects.                                          #
#                                                                    #
#====================================================================#


import os
import inspect
import keyword
import numpy as np
from numpy import fromstring,empty,array,float64,\
    loadtxt,ndarray,dtype,sqrt,pi,arange,exp,eye,\
    ceil,mod,dot,abs,identity,floor,linalg,where,isclose
from io import StringIO
from superstring import string2val
from generic import obj,hidden
from xmlreader import XMLreader,XMLelement
from developer import DevBase,error
from periodic_table import is_element
from structure import Structure,Jellium,get_kpath
from physical_system import PhysicalSystem
from simulation import SimulationInput,SimulationInputTemplate
from pwscf_input import array_to_string as pwscf_array_string
from debug import ci as interact

yesno_dict     = {True:'yes' ,False:'no'}
truefalse_dict = {True:'true',False:'false'}
onezero_dict   = {True:'1'   ,False:'0'}
boolmap={'yes':True,'no':False,'true':True,'false':False,'1':True,'0':False}

def is_int(var):
    try:
        int(var)
        return True
    except ValueError:
        return False
    #end try
#end def is_int

def is_float(var):
    try:
        float(var)
        return True
    except ValueError:
        return False
    #end try
#end def is_float

def is_array(var,type):
    try:
        if isinstance(var,str):
            array(var.split(),type)
        else:
            array(var,type)
        #end if
        return True
    except ValueError:
        return False
    #end try
#end def is_float_array


def attribute_to_value(attr):
    if is_int(attr):
        val = int(attr)
    elif is_float(attr):
        val = float(attr)
    elif is_array(attr,int):
        val = array(attr.split(),int)
        if val.size==9:
            val.shape = 3,3
        #end if
    elif is_array(attr,float):
        val = array(attr.split(),float)
    else:
        val = attr
    #end if
    return val
#end def attribute_to_value



#local write types
def yesno(var):
    return render_bool(var,'yes','no')
#end def yesno

def yesnostr(var):
    if isinstance(var,str):
        return var
    else:
        return yesno(var)
    #end if
#end def yesnostr

def onezero(var):
    return render_bool(var,'1','0')
#end def onezero

def truefalse(var):
    return render_bool(var,'true','false')
#end def onezero

def render_bool(var,T,F):
    if isinstance(var,bool) or var in (1,0):
        if var:
            return T
        else:
            return F
        #end if
    elif var in (T,F):
        return var
    else:
        error('Invalid QMCPACK input encountered.\nUser provided an invalid value of "{}" when yes/no was expected.\nValid options are: "{}", "{}", True, False, 1, 0'.format(var,T,F))

    #end if
#end def render_bool


bool_write_types = set([yesno,onezero,truefalse])




class QIobj(DevBase):

    afqmc_mode = False

    # user settings
    permissive_read  = False
    permissive_write = False
    permissive_init  = False

    @staticmethod
    def settings(
        permissive_read  = False,
        permissive_write = False,
        permissive_init  = False,
        ):
        QIobj.permissive_read  = permissive_read 
        QIobj.permissive_write = permissive_write
        QIobj.permissive_init  = permissive_init 
    #end def settings
#end class QIobj


class meta(obj):
    None
#end class meta


class section(QIobj):
    def __init__(self,*args,**kwargs):
        self.args   = args
        self.kwargs = kwargs
    #end def __init__
#end class section


class collection(hidden):
    def __init__(self,*elements):
        hidden.__init__(self)
        if len(elements)==1 and isinstance(elements[0],list):
            elements = elements[0]
        elif len(elements)==1 and isinstance(elements[0],collection):
            elements = elements[0].__dict__.values()
        #end if
        self.hidden().order = []
        for element in elements:
            self.add(element)
        #end for
    #end def __init__

    def __setitem__(self,name,value):
        #self.error('elements can only be set via the add function')
        self.add(value,key=name)
    #end def __setitem__

    def __delitem__(self,name):
        #self.error('elements can only be deleted via the remove function')
        self.remove(name)
    #end def __delitem__

    def __setattr__(self,name,value):
        #self.error('elements can only be set via the add function')
        self.add(value,key=name)
    #end def __setattr__

    def __delattr__(self,name):
        #self.error('elements can only be deleted via the remove function')
        self.remove(name)
    #end def __delattr__

    def add(self,element,strict=True,key=None):
        if not isinstance(element,QIxml):
            self.error('collection cannot be formed\nadd attempted for non QIxml element\ntype received: {0}'.format(element.__class__.__name__))
        #end if
        keyin = key
        key   = None
        public = self.public()
        identifier = element.identifier
        missing_identifier = False
        if not element.tag in plurals_inv and element.collection_id is None:
            self.error('collection cannot be formed\n  encountered non-plural element\n  element class: {0}\n  element tag: {1}\n  tags allowed in a collection: {2}'.format(element.__class__.__name__,element.tag,sorted(plurals_inv.keys())))
        elif identifier is None:
            key = len(public)
        elif isinstance(identifier,str):
            if identifier in element:
                key = element[identifier]
            else:
                missing_identifier = True
            #end if
        else:
            key = ''
            for ident in identifier:
                if ident in element:
                    key+=element[ident]
                #end if
            #end for
            missing_identifier = key==''
        #end if
        if missing_identifier:
            key = len(public)
        #end if
        if keyin is not None and not isinstance(key,int) and keyin.lower()!=key.lower():
            self.error('attempted to add key with incorrect name\nrequested key: {0}\ncorrect key: {1}'.format(keyin,key))
        #end if
        #if key in public:
        #    self.error('attempted to add duplicate key to collection: {0}\n keys present: {1}'.format(key,sorted(public.keys())))
        ##end if
        public[key] = element
        self.hidden().order.append(key)
        return True
    #end def add

    def remove(self,key):
        public = self.public()
        if key in public:
            del public[key]
            self.hidden().order.remove(key)
        else:
            raise KeyError
        #end if
    #end def remove

    def get_single(self,preference=None):
        if len(self)>0:
            if preference!=None and preference in self:
                return self[preference]
            else:
                return self.list()[0]
            #end if
        else:
            #return self
            return None
        #end if
    #end def get_single

    def list(self):
        lst = []
        for key in self.hidden().order:
            lst.append(self[key])
        #end for
        return lst
    #end def list

    def pairlist(self):
        pairs = []
        for key in self.hidden().order:
            pairs.append((key,self[key]))
        #end for
        return pairs
    #end def pairlist
#end class collection


def make_collection(elements):
    return collection(*elements)
#end def make_collection


class classcollection(QIobj):
    def __init__(self,*classes):
        if len(classes)==1 and isinstance(classes[0],list):
            classes = classes[0]
        #end if
        self.classes = classes
    #end def __init__
#end class classcollection


class QmcpackInputCollections(QIobj):
    def add(self,element):
        if element.tag in plurals_inv:
            cname = plurals_inv[element.tag]
            if not cname in self:
                coll = collection()
                success = coll.add(element,strict=False)
                if success:
                    self[cname] = coll
                #end if
            else:
                self[cname].add(element,strict=False)
            #end if
        #end if
    #end def add

    def get(self,cname,label=None):
        v = None
        if cname in self:
            if label is None:
                v = self[cname]
            elif label in self[cname]:
                v = self[cname][label]
            #end if
        #end if
        return v
    #end def get
#end class QmcpackInputCollections
QIcollections = QmcpackInputCollections()


class Names(QIobj):
    condensed_names = obj()
    expanded_names = None

    rsqmc_expanded_names = None
    afqmc_expanded_names = None

    escape_names = set(keyword.kwlist+['write'])
    escaped_names = list(escape_names)
    for i in range(len(escaped_names)):
        escaped_names[i]+='_'
    #end for
    escaped_names = set(escaped_names)

    @staticmethod
    def set_expanded_names(**kwargs):
        exnames = obj(**kwargs)
        Names.expanded_names = exnames
        Names.rsqmc_expanded_names = exnames
    #end def set_expanded_names

    @staticmethod
    def set_afqmc_expanded_names(**kwargs):
        Names.afqmc_expanded_names = obj(**kwargs)
    #end def set_afqmc_expanded_names

    @staticmethod
    def use_rsqmc_expanded_names():
        Names.expanded_names = Names.rsqmc_expanded_names
    #end def use_rsqmc_expanded_names

    @staticmethod
    def use_afqmc_expanded_names():
        Names.expanded_names = Names.afqmc_expanded_names
    #end def use_afqmc_expanded_names

    def expand_name(self,condensed):
        expanded = condensed
        cname = self.condense_name(condensed)
        if cname in self.escaped_names:
            cname = cname[:-1]
            expanded = cname
        #end if
        if cname in self.expanded_names:
            expanded = self.expanded_names[cname]
        #end if
        return expanded
    #end def expand_name

    def condense_name(self,expanded):
        condensed = expanded
        condensed = condensed.replace('___','_').replace('__','_')
        condensed = condensed.replace('-','_').replace(' ','_')
        condensed = condensed.lower()
        if condensed in self.escape_names:
            condensed += '_'
        #end if
        self.condensed_names[expanded]=condensed
        return condensed
    #end def condense_name

    def condense_names(self,*namelists):
        out = []
        for namelist in namelists:
            exp = obj()
            for expanded in namelist:
                condensed = self.condense_name(expanded)
                exp[condensed]=expanded
            #end for
            out.append(exp)
        #end for
        return out
    #end def condense_names

    def condensed_name_report(self):
        print()
        print('Condensed Name Report:')
        print('----------------------')
        keylist = array(list(self.condensed_names.keys()))
        order = array(list(self.condensed_names.values())).argsort()
        keylist = keylist[order]
        for expanded in keylist: 
            condensed = self.condensed_names[expanded]
            if expanded!=condensed:
                print("    {0:15} = '{1}'".format(condensed,expanded))
            #end if
        #end for
        print()
        print()
    #end def condensed_name_report
#end class Names




class QIxml(Names):

    def init_from_args(self,args):
        print()
        print('In init from args (not implemented).')
        print('Possible reasons for incorrect entry:  ')
        print('  Is xml element {0} meant to be plural?'.format(self.__class__.__name__))
        print('    If so, add it to the plurals object.')
        print()
        print('Arguments received:')
        print(args)
        print()
        self.not_implemented()
    #end def init_from_args



    @classmethod
    def init_class(cls):
        cls.class_set_optional(
            tag         = cls.__name__,
            identifier  = None,
            attributes  = [],
            elements    = [],
            text        = None,
            parameters  = [],
            attribs     = [],
            costs       = [],
            h5tags      = [],
            types       = obj(),
            write_types = obj(),
            attr_types  = None,
            precision   = None,
            defaults    = obj(),
            collection_id = None,
            exp_names   = None,
            )
        for v in ['attributes','elements','parameters','attribs','costs','h5tags']:
            names = cls.class_get(v)
            for i in range(len(names)):
                if names[i] in cls.escape_names:
                    names[i]+='_'
                #end if
            #end for
        #end for
        cls.params = cls.parameters + cls.attribs + cls.costs + cls.h5tags
        cls.plurals_inv = obj()
        for e in cls.elements:
            if e in plurals_inv:
                cls.plurals_inv[e] = plurals_inv[e]
            #end if
        #end for
        cls.plurals = cls.plurals_inv.inverse()
        if cls.exp_names is not None:
            cls.expanded_names = obj(Names.expanded_names,cls.exp_names)
        #end if
    #end def init_class


    def write(self,indent_level=0,pad='   ',first=False):
        param.set_precision(self.get_precision())
        if not QIobj.permissive_write:
            self.check_junk(exit=True)
        #end if
        indent  = indent_level*pad
        ip = indent+pad
        ipp= ip+pad
        expanded_tag = self.expand_name(self.tag)
        c = indent+'<'+expanded_tag
        for a in self.attributes:
            if a in self:
                val = self[a]
                if isinstance(val,str):
                    val = self.expand_name(val)
                #end if
                c += ' '+self.expand_name(a)+'='
                if a in self.write_types:
                    c += '"'+self.write_types[a](val)+'"'
                else:
                    c += '"'+param.write(val)+'"'
                #end if
            #end if
        #end for
        #if first:
        #    c+=' xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://www.mcc.uiuc.edu/qmc/schema/molecu.xsd"'
        ##end if
        #no_contents = len(set(self.keys())-set(self.elements)-set(self.plurals.keys()))==0
        no_contents = len(set(self.keys())-set(self.attributes))==0
        if no_contents:
            c += '/>\n'
        else:
            c += '>\n'
            for v in self.h5tags:
                if v in self:
                    if v in self.write_types:
                        write_type = self.write_types[v]
                    else:
                        write_type = None
                    #end if
                    c += param.write(self[v],name=self.expand_name(v),tag='h5tag',mode='elem',pad=ip,write_type=write_type)
                #end if
            #end for
            for v in self.costs:
                if v in self:
                    c += param.write(self[v],name=self.expand_name(v),tag='cost',mode='elem',pad=ip)
                #end if
            #end for
            for p in self.parameters:
                if p in self:
                    if p in self.write_types:
                        write_type = self.write_types[p]
                    else:
                        write_type = None
                    #end if
                    c += param.write(self[p],name=self.expand_name(p),mode='elem',pad=ip,write_type=write_type)
                #end if
            #end for
            for a in self.attribs:
                if a in self:
                    if a in self.write_types:
                        write_type = self.write_types[a]
                    else:
                        write_type = None
                    #end if
                    c += param.write(self[a],name=self.expand_name(a),tag='attrib',mode='elem',pad=ip,write_type=write_type)
                #end if
            #end for
            elements = self.elements
            if self.afqmc_mode and 'afqmc_order' in self.__class__.__dict__:
                elements = self.afqmc_order
            #end if
            for e in elements:
                if e in self:
                    elem = self[e]
                    if isinstance(elem,QIxml):
                        c += elem.write(indent_level+1)
                    else:
                        begin = '<'+e+'>'                        
                        contents = param.write(elem)
                        end = '</'+e+'>'
                        if contents.strip()=='':
                            c += ip+begin+end+'\n'
                        else:                            
                            c += ip+begin+'\n'
                            c += ipp+contents+'\n'
                            c += ip+end+'\n'
                        #end if
                    #end if
                elif e in plurals_inv and plurals_inv[e] in self:
                    coll = self[plurals_inv[e]]
                    if not isinstance(coll,collection):
                        self.error('write failed\n  element {0} is not a collection\n  contents of element {0}:\n{1}'.format(plurals_inv[e],str(coll)))
                    #end if
                    for instance in coll.list():
                        c += instance.write(indent_level+1)
                    #end for
                #end if
            #end for
            if self.text!=None:
                c = c.rstrip('\n')
                c+=param.write(self[self.text],mode='elem',pad=ip,tag=None,normal_elem=True)
            #end if
            c+=indent+'</'+expanded_tag+'>\n'
        #end if
        param.reset_precision()

        return c
    #end def write


    def __init__(self,*args,**kwargs):
        if Param.metadata==None:
            Param.metadata = meta()
        #end if
        if len(args)==1:
            a = args[0]
            if isinstance(a,XMLelement):
                self.init_from_xml(a)
            elif isinstance(a,section):
                self.init_from_inputs(a.args,a.kwargs)
            elif isinstance(a,self.__class__):
                self.transfer_from(a)
            else:
                self.init_from_inputs(args,kwargs)
            #end if
        else:
            self.init_from_inputs(args,kwargs)
        #end if
        QIcollections.add(self)
    #end def __init__


    def init_from_xml(self,xml):
        al,el = self.condense_names(xml._attributes.keys(),xml._elements.keys())
        xa,sa = set(al.keys()) , set(self.attributes)
        attr = xa & sa
        junk = xa-attr
        junk_elem = []
        for e,ecap in el.items():
            value = xml._elements[ecap]
            if (isinstance(value,list) or isinstance(value,tuple)) and e in self.plurals_inv.keys():
                if e not in types:
                    self.error('input element "{}" is unknown'.format(e))
                #end if
                p = self.plurals_inv[e]
                plist = []
                for instance in value:
                    plist.append(types[e](instance))
                #end for
                self[p] = make_collection(plist)
            elif e in self.elements:
                if e not in types:
                    self.error('input element "{}" is unknown'.format(e))
                #end if
                self[e] = types[e](value)
            elif e in ['parameter','attrib','cost','h5tag']:
                if isinstance(value,XMLelement):
                    value = [value]
                #end if
                for p in value:
                    name = self.condense_name(p.name)
                    if name in self.params:
                        self[name] = param(p)
                    else:
                        junk_elem.append(name)
                    #end if
                #end for
            else:
                junk_elem.append(e)
            #end if
        #end for
        junk = junk | set(junk_elem)
        if QmcpackInput.profile_collection!=None:
            self.collect_profile(xml,al,el,junk)
        #end for
        if not QIobj.permissive_read:
            self.check_junk(junk)
        #end if
        if self.attr_types!=None:
            typed_attr = attr & set(self.attr_types.keys())
            attr -= typed_attr
            for a in typed_attr:
                self[a] = self.attr_types[a](xml._attributes[al[a]])
            #end for
        #end if
        for a in attr:
            if a in self.write_types and self.write_types[a] in bool_write_types:
                aval = xml._attributes[al[a]]
                if aval in boolmap:
                    self[a] = boolmap[aval]
                else:
                    self.error('{0} is not a valid value for boolean attribute {1}\n  valid values are: {2}'.format(aval,a,boolmap.keys()))
                #end if
            else:
                self[a] = attribute_to_value(xml._attributes[al[a]])
            #end if
        #end for
        if self.text!=None:
            self[self.text] = param(xml)
        #end if
    #end def init_from_xml


    def init_from_inputs(self,args,kwargs):
        if len(args)>0:
            if len(args)==1 and isinstance(args[0],self.__class__):
                self.transfer_from(args[0])
            elif len(args)==1 and isinstance(args[0],dict):
                self.init_from_kwargs(args[0])
            else:
                self.init_from_args(args)
            #end if
        #end if
        self.init_from_kwargs(kwargs)
    #end def init_from_inputs


    def init_from_kwargs(self,kwargs):
        ks=[]
        kmap = dict()
        for key,val in kwargs.items():
            ckey = self.condense_name(key)
            ks.append(ckey)
            kmap[ckey] = val
        #end for
        ks = set(ks)
        kwargs = kmap
        h5tags     = ks & set(self.h5tags)
        costs      = ks & set(self.costs)
        parameters = ks & set(self.parameters)
        attribs    = ks & set(self.attribs)
        attr = ks & set(self.attributes)
        elem = ks & set(self.elements)
        plur = ks & set(self.plurals.keys())
        if self.text!=None:
            text = ks & set([self.text])
        else:
            text = set()
        #end if
        if not QIobj.permissive_init:
            junk = ks -attr -elem -plur -h5tags -costs -parameters -attribs -text
            self.check_junk(junk,exit=True)
        #end if

        for v in h5tags:
            self[v] = param(kwargs[v])
        #end for
        for v in costs:
            self[v] = param(kwargs[v])
        #end for
        for v in parameters:
            self[v] = param(kwargs[v])
        #end for
        for v in attribs:
            self[v] = param(kwargs[v])
        #end for
        for a in attr:
            self[a] = kwargs[a]
        #end for
        for e in elem:
            self[e] = types[e](kwargs[e])
        #end for
        for p in plur:
            e = self.plurals[p]
            kwcoll = kwargs[p]
            if isinstance(kwcoll,collection):
                plist = kwcoll.list()
            elif isinstance(kwcoll,(list,tuple)):
                plist = kwcoll
            else:
                self.error('init failed\n  encountered non-list collection')
            #end if
            ilist = []
            for instance in plist:
                ilist.append(types[e](instance))
            #end for
            self[p] = make_collection(ilist)
        #end for
        for t in text:
            self[t] = kwargs[t]
        #end for
    #end def init_from_kwargs


    def incorporate_defaults(self,elements=False,overwrite=False,propagate=True):
        for name,value in self.defaults.items():
            defval=None
            if isinstance(value,classcollection):
                if elements:
                    coll=[]
                    for cls in value.classes:
                        ins = cls()
                        ins.incorporate_defaults()
                        coll.append(ins)
                    #end for
                    defval = make_collection(coll)
                #end if
            elif inspect.isclass(value):
                if elements:
                    defval = value()
                #end if
            elif inspect.isfunction(value):
                defval = value()
            else:
                defval = value
            #end if
            if defval!=None:
                if overwrite or not name in self:
                    self[name] = defval
                #end if
            #end if
        #end for
        if propagate:
            for name,value in self.items():
                if isinstance(value,QIxml):
                    value.incorporate_defaults(elements,overwrite)
                elif isinstance(value,collection):
                    for v in value:
                        if isinstance(v,QIxml):
                            v.incorporate_defaults(elements,overwrite)
                        #end if
                    #end for
                #end if
            #end for
        #end if
    #end def incorporate_defaults                    


    def check_junk(self,junk=None,exit=False):
        if junk is None:
            ks = set(self.keys())
            h5tags     = ks & set(self.h5tags)
            costs      = ks & set(self.costs)
            parameters = ks & set(self.parameters)
            attribs    = ks & set(self.attribs)
            attr = ks & set(self.attributes)
            elem = ks & set(self.elements)
            plur = ks & set(self.plurals.keys())
            if self.text is not None:
                text = ks & set([self.text])
            else:
                text = set()
            #end if
            junk = ks -attr -elem -plur -h5tags -costs -parameters -attribs -text
        #end if
        if len(junk)>0:
            oname = ''
            if self.tag!=self.__class__.__name__:
                oname = ' ('+self.__class__.__name__+')'
            #end if
            msg = '{0}{1} does not have the following attributes/elements:\n'.format(self.tag,oname)
            for jname in junk:
                msg+='    '+jname+'\n'
            #end for
            #if QmcpackInput.profile_collection is None:
            #    self.error(msg,'QmcpackInput',exit=exit,trace=exit)
            ##end if

            print(obj(dict(self.__class__.__dict__)))

            self.error(msg,'QmcpackInput',exit=exit,trace=exit)
        #end if
    #end def check_junk


    def collect_profile(self,xml,al,el,junk):
        attributes = obj(**al)
        parameters = obj()
        elements   = obj()
        for e,ecap in el.items():
            if e=='parameter':
                parameters[e]=ecap
            else:
                elements[e]=ecap
            #end if
        #end for
        profile = obj(
            attributes = attributes,
            parameters = parameters,
            elements   = elements,
            xml        = xml,
            junk       = junk
            )
        #xml = xml.copy()
        xname = xml._name
        if xname[-1:].isdigit():
            xname = xname[:-1]
        elif xname[-2:].isdigit():
            xname = xname[:-2]
        #end if
        #xml._name = xname

        if len(profile.junk)>0:
            print('  '+xname+' (found '+str(junk)+')')
            for sector in 'attributes elements'.split():
                missing = []
                for n in profile.junk:
                    if n in profile[sector]:
                        missing.append(profile[sector][n])
                    #end if
                #end for
                if len(missing)>0:
                    ms= '    '+sector+':'
                    for m in missing:
                        ms+=' '+m
                    #end for
                    print(ms)
                #end if
            #end for
            if 'parameter' in profile.xml:
                params = obj()
                for p in profile.xml.parameter:
                    params[p.name.lower()] = p.text.strip()
                #end for
                missing = []
                for n in profile.junk:
                    if n in params:
                        missing.append(n)
                    #end if
                #end for
                if len(missing)>0:
                    ms= '    parameters:'
                    for m in missing:
                        ms+=' '+m
                    #end for
                    print(ms)
                #end if
            #end if
            if junk!=set(['analysis']) and junk!=set(['ratio']) and junk!=set(['randmo']) and junk!=set(['printeloc', 'source']) and junk!=set(['warmup_steps']) and junk!=set(['sposet_collection']) and junk!=set(['eigensolve', 'atom']) and junk!=set(['maxweight', 'reweightedvariance', 'unreweightedvariance', 'energy', 'exp0', 'stabilizerscale', 'minmethod', 'alloweddifference', 'stepsize', 'beta', 'minwalkers', 'nstabilizers', 'bigchange', 'usebuffer']) and junk!=set(['loop2']) and junk!=set(['random']) and junk!=set(['max_steps']):
                exit()
            #end if
        #end if

        pc = QmcpackInput.profile_collection
        if not xname in pc:
            pc[xname] = obj()
        #end if
        pc[xname].append(profile)
    #end def collect_profile



    def get_single(self,preference):
        return self
    #end def get_single


    def get(self,names,namedict=None,host=False,root=True):
        if namedict is None:
            namedict = {}
        #end if
        if isinstance(names,str):
            names = [names]
        #end if
        if root and not host:
            if self.identifier!=None and self.identifier in self:
                identity = self[self.identifier]
            else:
                identity = None
            #end if
            for name in names:
                if name==self.tag:
                    namedict[name]=self
                elif name==identity:
                    namedict[name]=self
                #end if
            #end for
        #end if
        for name in names:
            loc = None
            if name in self:
                loc = name
            elif name in plurals_inv and plurals_inv[name] in self:
                loc = plurals_inv[name]
            #end if
            name_absent = not name in namedict 
            not_element = False
            if not name_absent:
                not_xml  = not isinstance(namedict[name],QIxml)
                not_coll = not isinstance(namedict[name],collection)
                not_element = not_xml and not_coll
            #end if
            if loc!=None and (name_absent or not_element):
                if host:
                    namedict[name] = self
                else:
                    namedict[name] = self[loc]
                #end if
            #end if
        #end for
        for name,value in self.items():
            if isinstance(value,QIxml):
                value.get(names,namedict,host,root=False)
            elif isinstance(value,collection):
                for n,v in value.items():
                    name_absent = not n in namedict 
                    not_element = False
                    if not name_absent:
                        not_xml  = not isinstance(namedict[n],QIxml)
                        not_coll = not isinstance(namedict[n],collection)
                        not_element = not_xml and not_coll
                    #end if
                    if n in names and (name_absent or not_element):
                        if host:
                            namedict[n] = value
                        else:
                            namedict[n] = v
                        #end if
                    #end if
                    if isinstance(v,QIxml):
                        v.get(names,namedict,host,root=False)
                    #end if
                #end if
            #end if
        #end for
        if root:
            namelist = []
            for name in names:
                if name in namedict:
                    namelist.append(namedict[name])
                else:
                    namelist.append(None)
                #end if
            #end for
            if len(namelist)==1:
                return namelist[0]
            else:
                return namelist
            #end if
        #end if
    #end def get

    def remove(self,*names):
        if len(names)==1 and not isinstance(names[0],str):
            names = names[0]
        #end if
        remove = []
        for name in names:
            attempt = True
            if name in self:
                rname = name
            elif name in plurals_inv and plurals_inv[name] in self:
                rname = plurals_inv[name]
            else:
                attempt = False
            #end if
            if attempt:
                val = self[rname]
                if isinstance(val,QIxml) or isinstance(val,collection):
                    remove.append(rname)
                #end if
            #end if
        #end for
        for name in remove:
            del self[name]
        #end for
        for name,value in self.items():
            if isinstance(value,QIxml):
                value.remove(*names)
            elif isinstance(value,collection):
                for element in value:
                    if isinstance(element,QIxml):
                        element.remove(*names)
                    #end if
                #end for
            #end if
        #end for
    #end def remove


    def assign(self,**kwargs):
        for var,vnew in kwargs.items():
            if var in self:
                val = self[var]
                not_coll = not isinstance(val,collection)
                not_xml  = not isinstance(val,QIxml)
                not_arr  = not isinstance(val,ndarray)
                if not_coll and not_xml and not_arr:
                    self[var] = vnew
                #end if
            #end if
        #end for
        for vname,val in self.items():
            if isinstance(val,QIxml):
                val.assign(**kwargs)
            elif isinstance(val,collection):
                for v in val:
                    if isinstance(v,QIxml):
                        v.assign(**kwargs)
                    #end if
                #end for
            #end if
        #end for
    #end def assign


    def replace(self,*args,**kwargs):
        if len(args)==2 and isinstance(args[0],str) and isinstance(args[1],str):
            vold,vnew = args
            args = [(vold,vnew)]
        #end for
        for valpair in args:
            vold,vnew = valpair
            for var,val in self.items():
                not_coll = not isinstance(val,collection)
                not_xml  = not isinstance(val,QIxml)
                not_arr  = not isinstance(val,ndarray)
                if not_coll and not_xml and not_arr and val==vold:
                    self[var] = vnew
                #end if
            #end for
        #end for
        for var,valpair in kwargs.items():
            vold,vnew = valpair
            if var in self:
                val = self[var]
                if vold==None:
                    self[var] = vnew
                else:
                    not_coll = not isinstance(val,collection)
                    not_xml  = not isinstance(val,QIxml)
                    not_arr  = not isinstance(val,ndarray)
                    if not_coll and not_xml and not_arr and val==vold:
                        self[var] = vnew
                    #end if
                #end if
            #end if
        #end for
        for vname,val in self.items():
            if isinstance(val,QIxml):
                val.replace(*args,**kwargs)
            elif isinstance(val,collection):
                for v in val:
                    if isinstance(v,QIxml):
                        v.replace(*args,**kwargs)
                    #end if
                #end for
            #end if
        #end for
    #end def replace


    def combine(self,other):
        #elemental combine only
        for name,element in other.items():
            plural = isinstance(element,collection)
            single = isinstance(element,QIxml)
            if single or plural:
                elem = []
                single_name = None
                plural_name = None
                if single:
                    elem.append(element)
                    single_name = name
                    if name in plurals_inv:
                        plural_name = plurals_inv[name]
                    #end if
                else:
                    elem.extend(element.values())
                    plural_name = name
                    single_name = plurals[name]
                #end if
                if single_name in self:
                    elem.append(self[single_name])
                    del self[single_name]
                elif plural_name!=None and plural_name in self:
                    elem.append(self[plural_name])
                    del self[plural_name]
                #end if
                if len(elem)==1:
                    self[single_name]=elem[0]
                elif plural_name==None:
                    self.error('attempting to combine non-plural elements: '+single_name)
                else:
                    self[plural_name] = make_collection(elem)
                #end if
            #end if
        #end for
    #end def combine

                    
    def move(self,**elemdests):        
        names = list(elemdests.keys())
        hosts = self.get_host(names)
        dests = self.get(list(elemdests.values()))
        if len(names)==1:
            hosts = [hosts]
            dests = [dests]
        #end if
        for i in range(len(names)):
            name = names[i]
            host = hosts[i]
            dest = dests[i]
            if host!=None and dest!=None and id(host)!=id(dest):
                if not name in host:
                    name = plurals_inv[name]
                #end if
                dest[name] = host[name]
                del host[name]
            #end if
        #end for
    #end def move



    def pluralize(self):
        make_plural = []
        for name,value in self.items():
            if isinstance(value,QIxml):
                if name in plurals_inv:
                    make_plural.append(name)
                #end if
                value.pluralize()
            elif isinstance(value,collection):
                if name in plurals_inv:
                    make_plural.append(name)
                #end if
                for v in value:
                    if isinstance(v,QIxml):
                        v.pluralize()
                    #end if
                #end for
            #end if
        #end for
        for name in make_plural:
            value = self[name]
            del self[name]
            plural_name = plurals_inv[name]
            self[plural_name] = make_collection([value])
        #end for
    #end def pluralize


    def difference(self,other,root=True):
        if root:
            q1 = self.copy()
            q2 = other.copy()
        else:
            q1 = self
            q2 = other
        #end if
        if q1.__class__!=q2.__class__:
            different = True
            diff = None
            d1 = q1
            d2 = q2
        else:
            cls = q1.__class__
            s1 = set(q1.keys())
            s2 = set(q2.keys())
            shared  = s1 & s2
            unique1 = s1 - s2
            unique2 = s2 - s1
            different = len(unique1)>0 or len(unique2)>0
            diff = cls()
            d1 = cls()
            d2 = cls()
            d1.transfer_from(q1,unique1)
            d2.transfer_from(q2,unique2)
            for k in shared:
                value1 = q1[k]
                value2 = q2[k]
                is_coll1 = isinstance(value1,collection)
                is_coll2 = isinstance(value2,collection)
                is_qxml1 = isinstance(value1,QIxml)
                is_qxml2 = isinstance(value2,QIxml)
                if is_coll1!=is_coll2 or is_qxml1!=is_qxml2:
                    self.error('values for '+k+' are of inconsistent types\n  difference could not be taken')
                #end if
                if is_qxml1 and is_qxml2:
                    kdifferent,kdiff,kd1,kd2 = value1.difference(value2,root=False)
                elif is_coll1 and is_coll2:
                    ks1 = set(value1.keys())
                    ks2 = set(value2.keys())
                    kshared  = ks1 & ks2
                    kunique1 = ks1 - ks2
                    kunique2 = ks2 - ks1
                    kdifferent = len(kunique1)>0 or len(kunique2)>0
                    kd1 = collection()
                    kd2 = collection()
                    kd1.transfer_from(value1,kunique1)
                    kd2.transfer_from(value2,kunique2)
                    kdiff = collection()
                    for kk in kshared:
                        v1 = value1[kk]
                        v2 = value2[kk]
                        if isinstance(v1,QIxml) and isinstance(v2,QIxml):
                            kkdifferent,kkdiff,kkd1,kkd2 = v1.difference(v2,root=False)
                            kdifferent = kdifferent or kkdifferent
                            if kkdiff!=None:
                                kdiff[kk]=kkdiff
                            #end if
                            if kkd1!=None:
                                kd1[kk]=kkd1
                            #end if
                            if kkd2!=None:
                                kd2[kk]=kkd2
                            #end if
                        #end if
                    #end for
                else:
                    if isinstance(value1,ndarray):
                        a1 = value1.ravel()
                    else:
                        a1 = array([value1])
                    #end if
                    if isinstance(value2,ndarray):
                        a2 = value2.ravel()
                    else:
                        a2 = array([value2])
                    #end if
                    if len(a1)!=len(a2):
                        kdifferent = True
                    elif len(a1)==0:
                        kdifferent = False
                    elif (isinstance(a1[0],float) or isinstance(a2[0],float)) and not  (isinstance(a1[0],str)   or isinstance(a2[0],str)):
                        kdifferent = abs(a1-a2).max()/max(1e-99,abs(a1).max(),abs(a2).max()) > 1e-6
                    else:
                        kdifferent = not (a1==a2).all()
                    #end if
                    if kdifferent:
                        kdiff = (value1,value2)
                        kd1   = value1
                        kd2   = value2
                    else:
                        kdiff = None
                        kd1   = None
                        kd2   = None
                    #end if
                #end if
                different = different or kdifferent
                if kdiff!=None:
                    diff[k] = kdiff
                #end if
                if kd1!=None:
                    d1[k] = kd1
                #end if
                if kd2!=None:
                    d2[k] = kd2  
                #end if
            #end for
        #end if
        if root:
            if diff!=None:
                diff.remove_empty()
            #end if
            d1.remove_empty()
            d2.remove_empty()
        #end if 
        return different,diff,d1,d2
    #end def difference

    def remove_empty(self):
        names = list(self.keys())
        for name in names:
            value = self[name]
            if isinstance(value,QIxml):
                value.remove_empty()
                if len(value)==0:
                    del self[name]
                #end if
            elif isinstance(value,collection):
                ns = list(value.keys())
                for n in ns:
                    v = value[n]
                    if isinstance(v,QIxml):
                        v.remove_empty()
                        if len(v)==0:
                            del value[n]
                        #end if
                    #end if
                #end for
                if len(value)==0:
                    del self[name]
                #end if
            #end if
        #end for
    #end def remove_empty

    def get_host(self,names):
        return self.get(names,host=True)
    #end def get_host

    def get_precision(self):
        return self.__class__.class_get('precision')
    #end def get_precision
#end class QIxml



class QIxmlFactory(Names):
    def __init__(self,name,types,typekey='',typeindex=-1,typekey2='',default=None):
        self.name = name
        self.types = types
        self.typekey = typekey
        self.typekey2 = typekey2
        self.typeindex = typeindex
        self.default = default
    #end def __init__

    def __call__(self,*args,**kwargs):
        #emulate QIxml.__init__
        #get the value of the typekey
        a  = args
        kw = kwargs
        found_type = False
        if len(args)>0:
            v = args[0]
            if isinstance(v,XMLelement):
                kw = v._attributes
            elif isinstance(v,section):
                a  = v.args
                kw = v.kwargs
            elif isinstance(v,tuple(self.types.values())):
                found_type = True
                type = v.__class__.__name__
            #end if
        #end if
        if not found_type:
            if self.typekey in kw.keys():
                type = kw[self.typekey]
            elif self.typekey2 in kw.keys():
                type = kw[self.typekey2]
            elif self.default!=None:
                type = self.default
            elif self.typeindex==-1:
                self.error('QMCPACK input file is misformatted\ncannot identify type for <{0}/> element\nwith contents:\n{1}\nplease find the XML element matching this description in the input file to identify the problem\nmost likely, it is missing attributes "{2}" or "{3}"'.format(self.name,str(v).rstrip(),self.typekey,self.typekey2))
            else:
                type = a[self.typeindex]
            #end if
        #end if
        type = self.condense_name(type)
        if type in self.types:
            return self.types[type](*args,**kwargs)
        else:
            msg = self.name+' factory is not aware of the following subtype:\n'
            msg+= '    '+type+'\n'
            self.error(msg,exit=False,trace=False)
        #end if
    #end def __call__

    def init_class(self):
        None # this is for compatibility with QIxml only (do not overwrite)
    #end def init_class
#end class QIxmlFactory



class Param(Names):        
    metadata = None

    def __init__(self):
        self.reset_precision()
    #end def __init__

    def reset_precision(self):
        self.precision   = None
        self.prec_format = None
    #end def reset_precision

    def set_precision(self,precision):
        if precision is None:
            self.reset_precision()
        elif not isinstance(precision,str):
            self.error('attempted to set precision with non-string: {0}'.format(precision))
        else:
            self.precision   = precision
            self.prec_format = '{0:'+precision+'}'
        #end if
    #end def set_precision

    def __call__(self,*args,**kwargs):
        if len(args)==0:
            self.error('no arguments provided, should have received one XMLelement')
        elif not isinstance(args[0],XMLelement):
            return args[0]
            #self.error('first argument is not an XMLelement')
        #end if
        return self.read(args[0])
    #end def __call__

    def read(self,xml):
        val = ''
        attr = set(xml._attributes.keys())
        other_attr = attr-set(['name'])
        if 'name' in attr and len(other_attr)>0:
            oa = obj()
            for a in other_attr:
                oa[a] = xml._attributes[a]
            #end for
            self.metadata[xml.name] = oa
        #end if
        if 'text' in xml:
            token = xml.text.split('\n',1)[0].split(None,1)[0]
            try:
                if is_int(token):
                    val = loadtxt(StringIO(xml.text),int)
                elif is_float(token):
                    val = loadtxt(StringIO(xml.text),float)
                else:
                    val = array(xml.text.split())
                #end if
            except:
                if is_int(token):
                    val = array(xml.text.split(),dtype=int)
                elif is_float(token):
                    val = array(xml.text.split(),dtype=float)
                else:
                    val = array(xml.text.split())
                #end if
            #end try
            if val.size==1:
                val = val.ravel()[0]
            #end if
        #end if
        return val
    #end def read


    def write(self,value,mode='attr',tag='parameter',name=None,pad='   ',write_type=None,normal_elem=False):
        c = ''
        attr_mode = mode=='attr'
        elem_mode = mode=='elem'
        if not attr_mode and not elem_mode:
            self.error(mode+' is not a valid mode.  Options are attr,elem.')
        #end if
        if isinstance(value,list) or isinstance(value,tuple):
            value = array(value)
        #end if
        if attr_mode:
            if isinstance(value,ndarray):
                arr = value.ravel()
                for v in arr:
                    c+=self.write_val(v)+' '
                #end for
                c=c[:-1]
            else:
                c = self.write_val(value)
            #end if
        elif elem_mode:
            c+=pad
            is_array = isinstance(value,ndarray)
            is_single = not (is_array and value.size>1)
            if tag!=None:
                if is_single:
                    max_len = 20
                    rem_len = max(0,max_len-len(name))
                else:
                    rem_len = 0
                #end if
                other=''
                if name in self.metadata:
                    for a,v in self.metadata[name].items():
                        other +=' '+self.expand_name(a)+'="'+self.write_val(v)+'"'
                    #end for
                #end if
                c+='<'+tag+' name="'+name+'"'+other+rem_len*' '+'>'
                pp = pad+'   '
            else:
                pp = pad
            #end if
            if is_array:
                if normal_elem:
                    c+='\n'
                #end if
                if tag!=None:
                    c+='\n'
                #end if
                ndim = len(value.shape)
                if ndim==1:
                    line_len = 70
                    if tag!=None:
                        c+=pp
                    #end if
                    line = ''
                    for v in value:
                        line+=self.write_val(v)+' '
                        if len(line)>line_len:
                            c+=line+'\n'
                            line = ''
                        #end if
                    #end for
                    if len(line)>0:
                        c+=line
                    #end if
                    c=c[:-1]+'\n'
                elif ndim==2:
                    nrows,ncols = value.shape
                    fmt=pp
                    if value.dtype == dtype(float):
                        if self.precision is None:
                            vfmt = ':16.8f' # must have 8 digits of post decimal accuracy to meet qmcpack tolerance standards
                            #vfmt = ':16.8e'
                        else:
                            vfmt = ': '+self.precision
                        #end if
                    else:
                        vfmt = ''
                    #end if
                    for nc in range(ncols):
                        fmt+='{'+str(nc)+vfmt+'}  '
                    #end for
                    fmt = fmt[:-2]+'\n'
                    for nr in range(nrows):
                        c+=fmt.format(*value[nr])
                    #end for
                else:
                    self.error('only 1 and 2 dimensional arrays are supported for xml formatting.\n  Received '+ndim+' dimensional array.')
                #end if
            else:
                if write_type!=None:
                    val = write_type(value)
                else:
                    val = value
                #end if
                #c += '    '+str(val)
                c += '    {0:<10}'.format(self.write_val(val))
            #end if
            if tag!=None:
                c+=pad+'</'+tag+'>\n'
            #end if
        #end if
        return c
    #end def write
            

    def write_val(self,val):
        if self.precision!=None and isinstance(val,float):
            return self.prec_format.format(val)
        else:
            return str(val)
        #end if
    #end def write_val

    def init_class(self):
        None
    #end def init_class
#end class Param
param = Param()




class simulation(QIxml):
    #            afqmc
    attributes = ['method']
    #            rsqmc
    elements   = ['project','random','include','qmcsystem','particleset',
                  'wavefunction','hamiltonian','estimators','init','traces',
                  'qmc','loop','mcwalkerset','cmc']+\
                  ['afqmcinfo','walkerset','propagator','execute'] # afqmc
    afqmc_order = ['project','random','afqmcinfo','hamiltonian',
                   'wavefunction','walkerset','propagator','execute']
    write_types = obj(random=yesno)
#end class simulation


class project(QIxml):
    attributes = ['id','series']
    parameters = ['driver_version','maxcpusecs','max_seconds']
    elements   = ['application','host','date','user']
#end class project

class application(QIxml):
    attributes = ['name','role','class','version']
#end class application

class host(QIxml):
    text = 'value'
#end class host

class date(QIxml):
    text = 'value'
#end class date

class user(QIxml):
    text = 'value'
#end class user

class random(QIxml):
    attributes = ['seed','parallel']
    write_types= obj(parallel=truefalse)
#end class random

class include(QIxml):
    attributes = ['href']
#end def include

class mcwalkerset(QIxml):
    attributes = ['fileroot','version','collected','node','nprocs','href','target','file','walkers']
    write_types = obj(collected=yesno)
#end class mcwalkerset

class qmcsystem(QIxml):
    attributes = ['dim'] #,'wavefunction','hamiltonian']  # breaks QmcpackInput
    elements = ['simulationcell','particleset','wavefunction','hamiltonian','random','init','mcwalkerset','estimators']
#end class qmcsystem



class simulationcell(QIxml):
    attributes = ['name','tilematrix']
    parameters = ['lattice','reciprocal','bconds','lr_dim_cutoff','lr_tol','lr_handler','rs','nparticles','scale','uc_grid']
#end class simulationcell

class particleset(QIxml):
    attributes = ['name','size','random','random_source','randomsrc','charge','source']
    elements   = ['group','simulationcell']
    attribs    = ['ionid','position']
    write_types= obj(random=yesno)
    identifier = 'name'
#end class particleset

class group(QIxml):
    attributes = ['name','size','mass'] # mass attr and param, bad bad bad!!!
    parameters = ['charge','valence','atomicnumber','mass','lmax',
                  'cutoff_radius','spline_radius','spline_npoints'] 
    attribs    = ['position']
    identifier = 'name'
#end class group



class sposet(QIxml):
    attributes = ['basisset','type','name','group','size',
                  'index_min','index_max','energy_min','energy_max',
                  'spindataset','cuspinfo','sort','gpu','href','twistnum',
                  'gs_sposet','basis_sposet','same_k','frequency','mass',
                  'source','version','precision','tilematrix',
                  'meshfactor']
    elements   = ['occupation','coefficient','coefficients']
    text       = 'spos'
    identifier = 'name'
#end class sposet

class bspline_builder(QIxml):
    tag         = 'sposet_builder'
    identifier  = 'type'
    attributes  = ['type','href','sort','tilematrix','twistnum','twist','source',
                   'version','meshfactor','gpu','transform','precision','truncate',
                   'lr_dim_cutoff','shell','randomize','key','buffer','rmax_core','dilation','tag','hybridrep','gpusharing']
    elements    = ['sposet']
    write_types = obj(gpu=yesno,sort=onezero,transform=yesno,truncate=yesno,randomize=truefalse,hybridrep=yesno,gpusharing=yesno)
#end class bspline_builder

class heg_builder(QIxml):
    tag        = 'sposet_builder'
    identifier = 'type'
    attributes = ['type','twist']
    elements   = ['sposet']
#end class heg_builder

class molecular_orbital_builder(QIxml):
    tag = 'sposet_builder'
    identifier = 'type'
    attributes = ['name','type','transform','source','cuspcorrection']
    elements   = ['basisset','sposet'] 
#end class molecular_orbital_builder

class composite_builder(QIxml):
    tag = 'sposet_builder'
    identifier = 'type'
    attributes = ['type']
    elements   = ['sposet']
#end class composite_builder

sposet_builder = QIxmlFactory(
    name    = 'sposet_builder',
    types   = dict(bspline=bspline_builder,
                   einspline=bspline_builder,
                   heg=heg_builder,
                   composite=composite_builder,
                   molecularorbital = molecular_orbital_builder),
    typekey = 'type'
    )



class wavefunction(QIxml):
    #            rsqmc                        afqmc
    attributes = ['name','target','id','ref']+['info','type']
    #            afqmc
    parameters = ['filetype','filename','cutoff']
    elements   = ['sposet_builder','determinantset','jastrow','override_variational_parameters']
    identifier = 'name','id'
#end class wavefunction

class determinantset(QIxml):
    attributes = ['type','href','sort','tilematrix','twistnum','twist','source','version','meshfactor','gpu','transform','precision','truncate','lr_dim_cutoff','shell','randomize','key','rmax_core','dilation','name','cuspcorrection','tiling','usegrid','meshspacing','shell2','src','buffer','bconds','keyword','hybridrep','pbcimages','gpusharing']
    elements   = ['basisset','sposet','slaterdeterminant','multideterminant','spline','backflow','cubicgrid']
    h5tags     = ['twistindex','twistangle','rcut']
    write_types = obj(gpu=yesno,sort=onezero,transform=yesno,truncate=yesno,randomize=truefalse,cuspcorrection=yesno,usegrid=yesno,gpusharing=yesno)
#end class determinantset

class spline(QIxml):
    attributes = ['method']
    elements   = ['grid']
#end class spline

class cubicgrid(QIxml):
    attributes = ['method']
    elements   = ['grid']
#end class cubicgrid

class basisset(QIxml):
    attributes = ['ecut','name','ref','type','source','transform','key']
    elements   = ['grid','atomicbasisset']
    write_types = obj(transform=yesno)
#end class basisset

class grid(QIxml):
    attributes = ['dir','npts','closed','type','ri','rf','rc','step']
    #identifier = 'dir'
#end class grid

class atomicbasisset(QIxml):
    attributes = ['type','elementtype','expandylm','href','normalized','name','angular']
    elements   = ['grid','basisgroup']
    identifier = 'elementtype'
    write_types= obj(#expandylm=yesno,
                     normalized=yesno)
#end class atomicbasisset

class basisgroup(QIxml):
    attributes = ['rid','ds','n','l','m','zeta','type','s','imin','source']
    parameters = ['b']
    elements   = ['radfunc']
    #identifier = 'rid'
#end class basisgroup

class radfunc(QIxml):
    attributes = ['exponent','node','contraction','id','type']
    precision  = '16.12e'
#end class radfunc

class slaterdeterminant(QIxml):
    attributes = ['optimize','delay_rank','gpu','matrix_inverter']
    elements   = ['determinant']
    write_types = obj(optimize=yesno,gpu=yesno)
#end class slaterdeterminant

class determinant(QIxml):
    attributes = ['id','group','sposet','size','ref','spin','href','orbitals','spindataset','name','cuspinfo','debug']
    elements   = ['occupation','coefficient']
    identifier = 'id'
    write_types = obj(debug=yesno)
#end class determinant

class occupation(QIxml):
    attributes = ['mode','spindataset','size','pairs','format']
    text       = 'contents'
#end class occupation

class multideterminant(QIxml):
    attributes = ['optimize','spo_up','spo_dn']
    elements   = ['detlist']
#end class multideterminant

class detlist(QIxml):
    attributes = ['size','type','nca','ncb','nea','neb','nstates','cutoff','href']
    elements   = ['ci','csf']
#end class detlist

class ci(QIxml):
    attributes = ['id','coeff','qc_coeff','alpha','beta']
    #identifier = 'id'
    attr_types = obj(alpha=str,beta=str)
    precision  = '16.12e'
#end class ci

class csf(QIxml):
    attributes = ['id','exctlvl','coeff','coeff_real','coeff_imag','qchem_coeff','occ']
    elements   = ['det']
    attr_types = obj(occ=str)
#end class csf

class det(QIxml):
    attributes = ['id','coeff','alpha','beta']
    attr_types = obj(alpha=str,beta=str)
#end class det

class backflow(QIxml):
    attributes = ['optimize']
    elements   = ['transformation']
    write_types = obj(optimize=yesno)
#end class backflow

class transformation(QIxml):
    attributes = ['name','type','function','source']
    elements   = ['correlation']
    identifier = 'name'
#end class transformation

class jastrow1(QIxml):
    tag = 'jastrow'
    attributes = ['type','name','function','source','print','spin','transform']
    elements   = ['correlation','distancetable','grid']
    identifier = 'name'
    write_types = obj(print=yesno,spin=yesno,transform=yesno)
#end class jastrow1

class jastrow2(QIxml):
    tag = 'jastrow'
    attributes = ['type','name','function','print','spin','init','kc','transform','source','optimize']
    elements   = ['correlation','distancetable','basisset','grid','basisgroup']
    parameters = ['b','longrange']
    identifier = 'name'
    write_types = obj(print=yesno,transform=yesno,optimize=yesno)
#end class jastrow2

class jastrow3(QIxml):
    tag = 'jastrow'
    attributes = ['type','name','function','print','source']
    elements   = ['correlation']
    identifier = 'name'
    write_types = obj(print=yesno)
#end class jastrow3

class kspace_jastrow(QIxml):
    tag = 'jastrow'
    attributes = ['type','name','source']
    elements   = ['correlation']
    identifier = 'name'
    write_types = obj(optimize=yesno)
#end class kspace_jastrow

class rpa_jastrow(QIxml):
    tag = 'jastrow'
    attributes = ['type','name','source','function','kc']
    parameters = ['longrange']
    identifier = 'name'
    write_types = obj(longrange=yesno)
#end class rpa_jastrow

class correlation(QIxml):
    attributes = ['elementtype','speciesa','speciesb','size','ispecies','especies',
                  'especies1','especies2','isize','esize','rcut','cusp','pairtype',
                  'kc','type','symmetry','cutoff','spindependent','dimension','init',
                  'species']
    parameters = ['a','b','c','d']
    elements   = ['coefficients','var','coefficient']
    identifier = 'speciesa','speciesb','elementtype','especies1','especies2','ispecies'
    write_types = obj(init=yesno)
#end class correlation

class var(QIxml):
    attributes = ['id','name','optimize']
    text       = 'value'
    identifier = 'id'
    write_types=obj(optimize=yesno)
#end class var

class coefficients(QIxml):
    attributes = ['id','type','optimize','state','size','cusp','rcut']
    text       = 'coeff'
    write_types= obj(optimize=yesno)
    exp_names  = obj(array='Array')
#end class coefficients

class coefficient(QIxml):  # this is bad!!! coefficients/coefficient
    attributes = ['id','type','size','dataset','spindataset']
    text       = 'coeff'
    precision  = '16.12e'
#end class coefficient

class distancetable(QIxml):
    attributes = ['source','target']
#end class distancetable

jastrow = QIxmlFactory(
    name = 'jastrow',
    types   = dict(one_body=jastrow1,two_body=jastrow2,jastrow1=jastrow1,jastrow2=jastrow2,eei=jastrow3,jastrow3=jastrow3,kspace=kspace_jastrow,kspace_jastrow=kspace_jastrow,rpa=rpa_jastrow,rpa_jastrow=rpa_jastrow),
    typekey = 'type'
    )

class override_variational_parameters(QIxml):
    attributes = ['href']
#end class override_variational_parameters



class estimators(QIxml):
    elements = ['estimator']
#end class estimators

class hamiltonian(QIxml):
    #            rsqmc                              afqmc
    attributes = ['name','type','target','default']+['info']
    #            afqmc
    parameters = ['filetype','filename']
    elements   = ['pairpot','constant','estimator']
    identifier = 'name'
#end class hamiltonian

class coulomb(QIxml):
    tag = 'pairpot'
    attributes  = ['type','name','source','target','physical','forces']
    write_types = obj(physical=yesno)
    identifier  = 'name'
#end class coulomb

class constant(QIxml):
    attributes = ['type','name','source','target','forces']
    write_types= obj(forces=yesno)
#end class constant

class pseudopotential(QIxml):
    tag = 'pairpot'
    attributes = ['type','name','source','wavefunction','format','target','forces','dla']
    elements   = ['pseudo']
    write_types= obj(forces=yesno,dla=yesno)
    identifier = 'name'
#end class pseudopotential

class pseudo(QIxml):
    attributes = ['elementtype','href','format','cutoff','lmax','nrule','l_local']
    elements   = ['header','local','grid']
    identifier = 'elementtype'
#end class pseudo

class mpc(QIxml):
    tag='pairpot'
    attributes=['type','name','source','target','ecut','physical']
    write_types = obj(physical=yesno)    
    identifier='name'
#end class mpc

class cpp(QIxml):
    tag = 'pairpot'
    attributes = ['type','name','source','target']
    elements   = ['element']
    identifier = 'name'
#end class cpp

class element(QIxml):
    attributes = ['name','alpha','rb']
#end class element

pairpot = QIxmlFactory(
    name  = 'pairpot',
    types = dict(coulomb=coulomb,pseudo=pseudopotential,
                 pseudopotential=pseudopotential,mpc=mpc,
                 cpp=cpp),
    typekey = 'type'
    )


class header(QIxml):
    attributes = ['symbol','atomic-number','zval']
#end class header

class local(QIxml):
    elements = ['grid']
#end class local


class localenergy(QIxml):
    tag = 'estimator'
    attributes = ['name','hdf5']
    write_types= obj(hdf5=yesno)
    identifier = 'name'
#end class localenergy

class energydensity(QIxml):
    tag = 'estimator'
    attributes  = ['type','name','dynamic','static','ion_points']
    elements    = ['reference_points','spacegrid']
    identifier  = 'name'
    write_types = obj(ion_points=yesno)
#end class energydensity

class reference_points(QIxml):
    attributes = ['coord']
    text = 'points'
#end class reference_points

class spacegrid(QIxml):
    attributes = ['coord','min_part','max_part']
    elements   = ['origin','axis']
#end class spacegrid

class origin(QIxml):
    attributes = ['p1','p2']
#end class origin

class axis(QIxml):
    attributes = ['p1','p2','scale','label','grid']
    identifier = 'label'
#end class axis

class chiesa(QIxml):
    tag = 'estimator'
    attributes = ['name','type','source','psi','wavefunction','target']
    identifier = 'name'
#end class chiesa

class density(QIxml):
    tag = 'estimator'
    attributes = ['name','type','delta','x_min','x_max','y_min','y_max','z_min','z_max']
    identifier = 'type'
#end class density

class nearestneighbors(QIxml):
    tag = 'estimator'
    attributes = ['type']
    elements   = ['neighbor_trace']
    identifier = 'type'
#end class nearestneighbors

class neighbor_trace(QIxml):
    attributes = ['count','neighbors','centers']
    identifier = 'neighbors','centers'
#end class neighbor_trace

class spindensity(QIxml):
    tag = 'estimator'
    attributes  = ['type','name','report']
    parameters  = ['dr','grid','cell','center','corner','voronoi','test_moves']
    write_types = obj(report=yesno)
    identifier  = 'name'
#end class spindensity

class structurefactor(QIxml):
    tag = 'estimator'
    attributes  = ['type','name','report']
    write_types = obj(report=yesno)
    identifier  = 'name'
#end class structurefactor

class force(QIxml):
    tag = 'estimator'
    attributes = ['type','name','mode','source','species','target','addionion']
    parameters = ['rcut','nbasis','weightexp']
    identifier = 'name'
    write_types= obj(addionion=yesno)
#end class force

class forwardwalking(QIxml):
    tag = 'estimator'
    attributes = ['type','blocksize']
    elements   = ['observable']
    identifier = 'name'
#end class forwardwalking

class pressure(QIxml):
    tag = 'estimator'
    attributes = ['type','potential','etype','function']
    parameters = ['kc']
    identifier = 'type'
#end class pressure

class dmccorrection(QIxml):
    tag = 'estimator'
    attributes = ['type','blocksize','max','frequency']
    elements   = ['observable']
    identifier = 'type'
#end class dmccorrection

class nofk(QIxml):
    tag = 'estimator'
    attributes = ['type','name','wavefunction']
    identifier = 'name'
#end class nofk

class mpc_est(QIxml):
    tag = 'estimator'
    attributes = ['type','name','physical']
    write_types = obj(physical=yesno)
    identifier = 'name'
#end class mpc_est

class sk(QIxml):
    tag = 'estimator'
    attributes = ['name','type','hdf5']
    identifier = 'name'
    write_types = obj(hdf5=yesno)
#end class sk

class skall(QIxml):
    tag = 'estimator'
    attributes = ['name','type','hdf5','source','target','writeionion']
    identifier = 'name'
    write_types = obj(hdf5=yesno,writeionion=yesno)
#end class skall

class gofr(QIxml):
    tag = 'estimator'
    attributes = ['type','name','num_bin','rmax','source']
    identifier = 'name'
#end class gofr

class flux(QIxml):
    tag = 'estimator'
    attributes = ['type','name']
    identifier = 'name'
#end class flux

class momentum(QIxml): # legacy
    tag = 'estimator'
    attributes = ['type','name','grid','samples','hdf5','wavefunction','kmax','kmax0','kmax1','kmax2']
    identifier = 'name'
    write_types = obj(hdf5=yesno)
#end class momentum

class momentumdistribution(QIxml): # batched
    tag = 'estimator'
    attributes = ['type','name','grid','samples','hdf5','wavefunction','kmax','kmax0','kmax1','kmax2']
    identifier = 'name'
    write_types = obj(hdf5=yesno)
#end class momentumdistribution

class dm1b(QIxml): # legacy
    tag         = 'estimator'
    identifier  = 'type'
    attributes  = ['type','name','reuse']#reuse is a temporary dummy keyword
    parameters  = ['energy_matrix','basis_size','integrator','points','scale','basis','evaluator','center','check_overlap','check_derivatives','acceptance_ratio','rstats','normalized','volume_normed']
    write_types = obj(energy_matrix=yesno,check_overlap=yesno,check_derivatives=yesno,acceptance_ratio=yesno,rstats=yesno,normalized=yesno,volume_normed=yesno)
#end class dm1b

class onebodydensitymatrices(QIxml): # batched
    tag         = 'estimator'
    identifier  = 'type'
    attributes  = ['type','name','reuse']#reuse is a temporary dummy keyword
    parameters  = ['energy_matrix','basis_size','integrator','points','scale','basis','evaluator','center','check_overlap','check_derivatives','acceptance_ratio','rstats','normalized','volume_normed']
    write_types = obj(energy_matrix=yesno,check_overlap=yesno,check_derivatives=yesno,acceptance_ratio=yesno,rstats=yesno,normalized=yesno,volume_normed=yesno)
#end class onebodydensitymatrices


# afqmc estimators
class back_propagation(QIxml):
    tag = 'estimator'
    attributes = ['name']
    parameters = ['naverages','block_size','ortho','nsteps']
    elements   = ['onerdm']
    identifier = 'name'
#end class back_propagation

estimator = QIxmlFactory(
    name  = 'estimator',
    types = dict(localenergy            = localenergy,
                 energydensity          = energydensity,
                 chiesa                 = chiesa,
                 density                = density,
                 nearestneighbors       = nearestneighbors,
                 dm1b                   = dm1b,
                 spindensity            = spindensity,
                 structurefactor        = structurefactor,
                 force                  = force,
                 forwardwalking         = forwardwalking,
                 pressure               = pressure,
                 dmccorrection          = dmccorrection,
                 nofk                   = nofk,
                 mpc                    = mpc_est,
                 sk                     = sk,
                 skall                  = skall,
                 gofr                   = gofr,
                 flux                   = flux,
                 momentum               = momentum,
                 momentumdistribution   = momentumdistribution,
                 onebodydensitymatrices = onebodydensitymatrices,
                 # afqmc estimators   
                 back_propagation       = back_propagation,
                 ),
    typekey  = 'type',
    typekey2 = 'name'
    )


class observable(QIxml):
    attributes = ['name','max','frequency']
    identifier = 'name'
#end class observable



class init(QIxml):
    attributes = ['source','target']
#end class


class scalar_traces(QIxml):
    attributes  = ['defaults']
    text        = 'quantities'
    write_types = obj(defaults=yesno)
#end class scalar_traces

class array_traces(QIxml):
    attributes  = ['defaults']
    text        = 'quantities'
    write_types = obj(defaults=yesno)
#end class array_traces

class particle_traces(QIxml): # legacy
    attributes  = ['defaults']
    text        = 'quantities'
    write_types = obj(defaults=yesno)
#end class particle_traces

class traces(QIxml):
    attributes = ['write','throttle','format','verbose','scalar','array',
                  'scalar_defaults','array_defaults',
                  'particle','particle_defaults']
    elements = ['scalar_traces','array_traces','particle_traces']
    write_types = obj(write_=yesno,verbose=yesno,scalar=yesno,array=yesno,
                      scalar_defaults=yesno,array_defaults=yesno,
                      particle=yesno,particle_defaults=yesno)
#end class traces


class record(QIxml):
    attributes = ['name','stride']
#end class record


class loop(QIxml):
    collection_id = 'qmc'
    attributes = ['max']
    elements = ['qmc','init']
    def unroll(self):
        calculations=[]
        calcs = []
        if 'qmc' in self:
            calcs = [self.qmc]
        elif 'calculations' in self:
            calcs = self.calculations
        #end if
        for n in range(self.max):
            for i in range(len(calcs)):
                calculations.append(calcs[i].copy())
            #end for
        #end for
        return make_collection(calculations)
    #end def unroll
#end class loop


class optimize(QIxml):
    text = 'parameters'
#end class optimize

class cg_optimizer(QIxml):
    tag        = 'optimizer'
    attributes = ['method']
    parameters = ['max_steps','tolerance','stepsize','friction','epsilon',
                  'xybisect','verbose','max_linemin','tolerance_g','length_cg',
                  'rich','xypolish','gfactor']
#end class cg_optimizer

class flex_optimizer(QIxml):
    tag        = 'optimizer'
    attributes = ['method']
    parameters = ['max_steps','tolerance','stepsize','epsilon',
                  'xybisect','verbose','max_linemin','tolerance_g','length_cg',
                  'rich','xypolish','gfactor']
#end class flex_optimizer



optimizer = QIxmlFactory(
    name    = 'optimizer',
    types   = dict(cg=cg_optimizer,flexopt=flex_optimizer),
    typekey = 'method',
    )



class optimize_qmc(QIxml):
    collection_id = 'qmc'
    tag = 'qmc'
    attributes = ['method','move','renew','completed','checkpoint','gpu']
    parameters = ['blocks','steps','timestep','walkers','minwalkers','useweight',
                  'power','correlation','maxweight','usedrift','min_walkers',
                  'minke','samples','warmupsteps','minweight','warmupblocks',
                  'maxdispl','tau','tolerance','stepsize','epsilon',
                  'en_ref','usebuffer','substeps','stepsbetweensamples',
                  'samplesperthread','max_steps','nonlocalpp']
    elements = ['optimize','optimizer','estimator']
    costs    = ['energy','variance','difference','weight','unreweightedvariance','reweightedvariance']
    write_types = obj(renew=yesno,completed=yesno)
#end class optimize_qmc

class linear(QIxml):
    collection_id = 'qmc'
    tag = 'qmc'
    attributes = ['method','move','profiling','kdelay', # batched
                  'checkpoint','gpu','trace']           # legacy - batched
    elements   = ['estimator']
    parameters = ['total_walkers','walkers_per_rank','crowds','opt_num_crowds',   # batched
                  'walkers','warmupsteps','blocks','steps','substeps','timestep', # who knows
                  'usedrift','stepsbetweensamples','samples','minmethod',
                  'minwalkers','maxweight','nonlocalpp','use_nonlocalpp_deriv',
                  'usebuffer','alloweddifference','gevmethod','beta','exp0',
                  'bigchange','stepsize','stabilizerscale','nstabilizers',
                  'max_its','cgsteps','eigcg','stabilizermethod',
                  'rnwarmupsteps','walkersperthread','minke','gradtol','alpha',
                  'tries','min_walkers','samplesperthread',
                  'shift_i','shift_s','max_relative_change','max_param_change',
                  'chase_lowest','chase_closest','block_lm','nblocks','nolds',
                  'nkept','max_seconds'
                  ]
    costs      = ['energy','unreweightedvariance','reweightedvariance','variance','difference']
    write_types = obj(gpu=yesno,usedrift=yesno,nonlocalpp=yesno,usebuffer=yesno,use_nonlocalpp_deriv=yesno,chase_lowest=yesno,chase_closest=yesno,block_lm=yesno)
#end class linear

class cslinear(QIxml):
    collection_id = 'qmc'
    tag = 'qmc'
    attributes = ['method','move','checkpoint','gpu','trace']
    elements   = ['estimator']
    parameters = ['walkers','warmupsteps','blocks','steps','substeps','timestep',
                  'usedrift','stepsbetweensamples','samples','minmethod',
                  'minwalkers','maxweight','nonlocalpp','usebuffer',
                  'alloweddifference','gevmethod','beta','exp0','bigchange',
                  'stepsize','stabilizerscale','nstabilizers','max_its',
                  'stabilizermethod','cswarmupsteps','alpha_error','gevsplit',
                  'beta_error','use_nonlocalpp_deriv']
    costs      = ['energy','unreweightedvariance','reweightedvariance']
    write_types = obj(gpu=yesno,usedrift=yesno,nonlocalpp=yesno,use_nonlocalpp_deriv=yesno,usebuffer=yesno)
#end class cslinear

class vmc(QIxml):
    collection_id = 'qmc'
    tag = 'qmc'
    attributes = ['method','move','profiling','kdelay',         # batched
                  'multiple','warp','gpu','checkpoint','trace', # legacy - batched
                  'target','completed','id'] 
    elements   = ['estimator', # batched
                  'record']    # legacy - batched
    parameters = ['total_walkers','walkers_per_rank','crowds','warmupsteps',         # batched
                  'blocks','steps','substeps','timestep','maxcpusecs','rewind',
                  'storeconfigs','checkproperties','recordconfigs','current',
                  'stepsbetweensamples','samplesperthread','samples','usedrift',
                  'walkers','nonlocalpp','tau','walkersperthread','reconfiguration', # legacy - batched
                  'dmcwalkersperthread','current','ratio','firststep',
                  'minimumtargetwalkers','max_seconds']
    write_types = obj(usedrift=yesno,profiling=yesno,                   # batched
                      gpu=yesno,nonlocalpp=yesno,reconfiguration=yesno, # legacy - batched
                      ratio=yesno,completed=yesno)
#end class vmc

class dmc(QIxml):
    collection_id = 'qmc'
    tag = 'qmc'
    attributes = ['method','move','profiling','kdelay',         # batched
                  'gpu','multiple','warp','checkpoint','trace', # legacy - batched
                  'target','completed','id','continue']
    elements   = ['estimator']
    parameters = ['total_walkers','walkers_per_rank','crowds','warmupsteps',
                  'crowd_serialize_walkers',            # batched
                  'blocks','steps','substeps','timestep','maxcpusecs','rewind',
                  'storeconfigs','checkproperties','recordconfigs','current',
                  'stepsbetweensamples','samplesperthread','samples','reconfiguration',
                  'nonlocalmoves','maxage','alpha','gamma','reserve','use_nonblocking',
                  'branching_cutoff_scheme','feedback','sigmabound',
                  'walkers','nonlocalmove','pop_control','targetwalkers',               # legacy - batched
                  'minimumtargetwalkers','energybound','feedback','recordwalkers',
                  'fastgrad','popcontrol','branchinterval','usedrift','storeconfigs',
                  'en_ref','tau','alpha','gamma','max_branch','killnode','swap_walkers',
                  'swap_trigger','branching_cutoff_scheme','l2_diffusion','maxage',
                  'max_seconds']
    write_types = obj(usedrift=yesno,profiling=yesno,reconfiguration=yesno,
                      crowd_serialize_walkers=yesno,    # batched
                      nonlocalmoves=yesnostr,use_nonblocking=yesno,
                      gpu=yesno,fastgrad=yesno,completed=yesno,killnode=yesno, # legacy - batched
                      swap_walkers=yesno,l2_diffusion=yesno)
#end class dmc

class rmc(QIxml):
    collection_id = 'qmc'
    tag = 'qmc'
    attributes = ['method','multiple','target','observables','target','warp']
    parameters = ['blocks','steps','chains','cuts','bounce','clone','walkers','timestep','trunclength','maxtouch','mass','collect']
    elements = ['qmcsystem']
    write_types = obj(collect=yesno)
#end class rmc

class vmc_batch(QIxml):
    # Do not assume all of the parameters below are supported.
    # These were simply copied over from legacy drivers because the 
    # batched driver compatible inputs have yet not been listed anywhere. 
    collection_id = 'qmc'
    tag = 'qmc'
    attributes = ['method','move','profiling','kdelay','checkpoint']
    elements   = ['estimator']
    parameters = ['total_walkers','walkers_per_rank','crowds','warmupsteps','blocks','steps','substeps','timestep','maxcpusecs','rewind','storeconfigs','checkproperties','recordconfigs','current','stepsbetweensamples','samplesperthread','samples','usedrift']
    write_types = obj(usedrift=yesno,profiling=yesno)
#end class vmc_batch

class dmc_batch(QIxml):
    # Do not assume all of the parameters below are supported.
    # These were simply copied over from legacy drivers because the 
    # batched driver compatible inputs have yet not been listed anywhere. 
    collection_id = 'qmc'
    tag = 'qmc'
    attributes = ['method','move','profiling','kdelay','checkpoint']
    elements   = ['estimator']
    parameters = ['total_walkers','walkers_per_rank','crowd_serialize_walkers','crowds','warmupsteps','blocks','steps','substeps','timestep','maxcpusecs','rewind','storeconfigs','checkproperties','recordconfigs','current','stepsbetweensamples','samplesperthread','samples','reconfiguration','nonlocalmoves','maxage','alpha','gamma','reserve','use_nonblocking','branching_cutoff_scheme','feedback','sigmabound']
    write_types = obj(usedrift=yesno,profiling=yesno,reconfiguration=yesno,nonlocalmoves=yesnostr,use_nonblocking=yesno, crowd_serialize_walkers=yesno)
#end class dmc_batch

class linear_batch(QIxml):
    # Do not assume all of the parameters below are supported.
    # These were simply copied over from legacy drivers because the 
    # batched driver compatible inputs have yet not been listed anywhere. 
    collection_id = 'qmc'
    tag = 'qmc'
    attributes = ['method','move','profiling','kdelay']
    elements   = ['estimator']
    parameters = ['walkers','warmupsteps','blocks','steps','substeps','timestep',
                  'usedrift','stepsbetweensamples','samples','minmethod',
                  'minwalkers','maxweight','nonlocalpp','use_nonlocalpp_deriv',
                  'usebuffer','alloweddifference','gevmethod','beta','exp0',
                  'bigchange','stepsize','stabilizerscale','nstabilizers',
                  'max_its','cgsteps','eigcg','stabilizermethod',
                  'rnwarmupsteps','walkersperthread','minke','gradtol','alpha',
                  'tries','min_walkers','samplesperthread',
                  'shift_i','shift_s','max_relative_change','max_param_change',
                  'chase_lowest','chase_closest','block_lm','nblocks','nolds',
                  'nkept',
                  'crowds','opt_num_crowds'
                  ]
    costs      = ['energy','unreweightedvariance','reweightedvariance','variance','difference']
    write_types = obj(usedrift=yesno,nonlocalpp=yesno,usebuffer=yesno,use_nonlocalpp_deriv=yesno,chase_lowest=yesno,chase_closest=yesno,block_lm=yesno)
#end class linear_batch

class wftest(QIxml):
    collection_id = 'qmc'
    tag = 'qmc'
    attributes = ['method','checkpoint', 'gpu', 'move', 'multiple', 'warp']
    parameters = ['ratio','walkers','clone','source','hamiltonianpbyp','orbitalutility','printeloc','basic','virtual_move']
    #elements   = ['printeloc','source']
    write_types = obj(ratio=yesno,clone=yesno,hamiltonianpbyp=yesno,orbitalutility=yesno,printeloc=yesno,basic=yesno,virtual_move=yesno)
#end class wftest    

class setparams(QIxml):
    collection_id = 'qmc'
    tag = 'qmc'
    attributes = ['method','move','checkpoint','gpu']
    parameters = ['alpha','blocks','warmupsteps','stepsbetweensamples','timestep','samples','usedrift']
    elements   = ['estimator']
#end class setparams    

qmc = QIxmlFactory(
    name = 'qmc',
    types   = dict(linear=linear,cslinear=cslinear,vmc=vmc,dmc=dmc,loop=loop,optimize=optimize_qmc,wftest=wftest,rmc=rmc,setparams=setparams,vmc_batch=vmc_batch,dmc_batch=dmc_batch,linear_batch=linear_batch),
    typekey = 'method',
    default = 'loop'
    )



class cmc(QIxml):
    attributes = ['method','target']
#end class cmc



# afqmc elements

class afqmcinfo(QIxml):
    attributes = ['name']
    parameters = ['nmo','naea','naeb']
#end class afqmcinfo

class walkerset(QIxml):
    attributes = ['name','type']
    parameters = ['walker_type']
#end class walkerset

class propagator(QIxml):
    attributes  = ['name','info']
    parameters  = ['hybrid']
    write_types = obj(hybrid=yesno)
#end class propagator

class execute(QIxml):
    attributes = ['info','ham','wfn','wset','prop']
    parameters = ['ncores','nwalkers','blocks','steps','timestep']
    elements   = ['estimator']
#end class execute

class onerdm(QIxml):
    None
#end class onerdm





class gen(QIxml):
    attributes = []
    elements   = []
#end class gen


classes = [   #standard classes
    simulation,project,application,random,qmcsystem,simulationcell,particleset,
    group,hamiltonian,constant,pseudopotential,coulomb,pseudo,mpc,chiesa,density,
    localenergy,energydensity,spacegrid,origin,axis,wavefunction,
    determinantset,slaterdeterminant,basisset,grid,determinant,occupation,
    jastrow1,jastrow2,jastrow3,
    correlation,coefficients,loop,linear,cslinear,vmc,dmc,vmc_batch,dmc_batch,linear_batch,
    atomicbasisset,basisgroup,init,var,traces,scalar_traces,particle_traces,array_traces,
    reference_points,nearestneighbors,neighbor_trace,dm1b,
    coefficient,radfunc,spindensity,structurefactor,
    sposet,bspline_builder,composite_builder,heg_builder,include,
    multideterminant,detlist,ci,mcwalkerset,csf,det,
    optimize,cg_optimizer,flex_optimizer,optimize_qmc,wftest,kspace_jastrow,
    header,local,force,forwardwalking,observable,record,rmc,pressure,dmccorrection,
    nofk,mpc_est,flux,distancetable,cpp,element,spline,setparams,
    backflow,transformation,cubicgrid,molecular_orbital_builder,cmc,sk,skall,gofr,
    host,date,user,rpa_jastrow,momentum,override_variational_parameters,
    momentumdistribution,onebodydensitymatrices,estimators,
    # afqmc classes
    afqmcinfo,walkerset,propagator,execute,back_propagation,onerdm
    ]
types = dict( #simple types and factories
    #host           = param,
    #date           = param,
    #user           = param,
    pairpot        = pairpot,
    estimator      = estimator,
    sposet_builder = sposet_builder,
    jastrow        = jastrow,
    qmc            = qmc,
    optimizer      = optimizer,
    )
plurals = obj(
    particlesets    = 'particleset',
    groups          = 'group',    
    hamiltonians    = 'hamiltonian', 
    pairpots        = 'pairpot',
    pseudos         = 'pseudo',
    estimators      = 'estimator',
    spacegrids      = 'spacegrid',
    axes            = 'axis',
    wavefunctions   = 'wavefunction',
    grids           = 'grid',
    determinants    = 'determinant',
    correlations    = 'correlation',
    jastrows        = 'jastrow',
    basisgroups     = 'basisgroup',
    calculations    = 'qmc',
    vars            = 'var',
    neighbor_traces = 'neighbor_trace',
    sposet_builders = 'sposet_builder',
    sposets         = 'sposet',
    radfuncs        = 'radfunc',
    #qmcsystems      = 'qmcsystem',  # not a good idea
    atomicbasissets = 'atomicbasisset',
    cis             = 'ci',
    csfs            = 'csf',
    dets            = 'det',
    observables     = 'observable',
    optimizes       = 'optimize',
    #coefficientss   = 'coefficients', # bad plurality of qmcpack
    constants       = 'constant',
    mcwalkersets    = 'mcwalkerset',
    transformations = 'transformation',
    )
plurals_inv = plurals.inverse()
plural_names = set(plurals.keys())
single_names = set(plurals.values())
Names.set_expanded_names(
    elementtype      = 'elementType',
    energydensity    = 'EnergyDensity',
    gevmethod        = 'GEVMethod',
    localenergy      = 'LocalEnergy',
    lr_dim_cutoff    = 'LR_dim_cutoff',
    lr_tol           = 'LR_tol',
    lr_handler       = 'LR_handler',
    minmethod        = 'MinMethod',
    one_body         = 'One-Body',
    speciesa         = 'speciesA',
    speciesb         = 'speciesB',
    substeps         = 'subSteps',
    two_body         = 'Two-Body',
    usedrift         = 'useDrift',
    maxweight        = 'maxWeight',
    warmupsteps      = 'warmupSteps',
    twistindex       = 'twistIndex',
    twistangle       = 'twistAngle',
    usebuffer        = 'useBuffer',
    mpc              = 'MPC',
    kecorr           = 'KEcorr',
    ionion           = 'IonIon',
    elecelec         = 'ElecElec',
    pseudopot        = 'PseudoPot',
    posarray         = 'posArray',
    #array            = 'Array',  # handle separately, namespace collision
    atomicbasisset   = 'atomicBasisSet',
    basisgroup       = 'basisGroup',
    expandylm        = 'expandYlm',
    mo               = 'MO',
    numerical        = 'Numerical',
    nearestneighbors = 'NearestNeighbors',
    cuspcorrection   = 'cuspCorrection',
    cuspinfo         = 'cuspInfo',
    exctlvl          = 'exctLvl',
    pairtype         = 'pairType',
    printeloc        = 'printEloc',
    spindependent    = 'spinDependent',
    l_local          = 'l-local',
    pbcimages        = 'PBCimages',
    dla              = 'DLA',
    l2_diffusion     = 'L2_diffusion',
    maxage           = 'MaxAge',
    sigmabound       = 'sigmaBound',
    )
# afqmc names
Names.set_afqmc_expanded_names(
    afqmcinfo        = 'AFQMCInfo',
    nmo              = 'NMO',
    naea             = 'NAEA',
    naeb             = 'NAEB',
    hamiltonian      = 'Hamiltonian',
    wavefunction     = 'Wavefunction',
    walkerset        = 'WalkerSet',
    propagator       = 'Propagator',
    onerdm           = 'OneRDM',
    nwalkers         = 'nWalkers',
    estimator        = 'Estimator',
    )
for c in classes:
    c.init_class()
    types[c.__name__] = c
#end for


#set default values
simulation.defaults.set(
    project      = project,
    qmcsystem    = qmcsystem,
    calculations = lambda:list()
    )
project.defaults.set(
    series=0,
    application = application
    )
application.defaults.set(
    name='qmcpack',role='molecu',class_='serial',version='1.0'
    )
#simulationcell.defaults.set(
#    bconds = 'p p p',lr_dim_cutoff=15
#    )
wavefunction.defaults.set(
    name='psi0'
    )
#determinantset.defaults.set(
#    type='einspline',tilematrix=lambda:eye(3,dtype=int),meshfactor=1.,gpu=False,precision='double'
#    )
#occupation.defaults.set(
#    mode='ground',spindataset=0
#    )
jastrow1.defaults.set(
    name='J1',type='one-body',function='bspline',print=True,source='ion0',
    correlation=correlation
    )
jastrow2.defaults.set(
    name='J2',type='two-body',function='bspline',print=True,
    correlation=correlation
    )
jastrow3.defaults.set(
    name='J3',type='eeI',function='polynomial',print=True,source='ion0',
    correlation=correlation
    )
correlation.defaults.set(
    coefficients=coefficients
    )
coefficients.defaults.set(
    type='Array'
    )
#hamiltonian.defaults.set(
#    name='h0',type='generic',target='e',
#    constant = constant,
#    pairpots = classcollection(coulomb,pseudopotential,mpc),
#    estimators = classcollection(chiesa),
#    )
#coulomb.defaults.set(
#    name='ElecElec',type='coulomb',source='e',target='e'
#    )
#constant.defaults.set(
#    name='IonIon',type='coulomb',source='ion0',target='ion0'
#    )
#pseudopotential.defaults.set(
#    name='PseudoPot',type='pseudo',source='ion0',wavefunction='psi0',format='xml'
#    )
#mpc.defaults.set(
#    name='MPC',type='MPC',ecut=60.0,source='e',target='e',physical=False
#    )
localenergy.defaults.set(
    name='LocalEnergy',hdf5=True
    )
#chiesa.defaults.set(
#    name='KEcorr',type='chiesa',source='e',psi='psi0'
#    )
#energydensity.defaults.set(
#    type='EnergyDensity',name='EDvoronoi',dynamic='e',static='ion0',
#    spacegrid = spacegrid
#    )
#spacegrid.defaults.set(
#    coord='voronoi'
#    )
density.defaults.set(
    type='density',name='Density'
    )
spindensity.defaults.set(
    type='spindensity',name='SpinDensity'
    )
skall.defaults.set(
    type='skall',name='skall',source='ion0',target='e',hdf5=True
    )
force.defaults.set(
    type='Force',name='force'
    )
pressure.defaults.set(
    type='Pressure'
    )
momentum.defaults.set(
    type='momentum'
    )
momentumdistribution.defaults.set(
    type='MomentumDistribution',name='nofk'
    )
dm1b.defaults.set(
    type = 'dm1b',name='DensityMatrices',energy_matrix=False,
    evaluator='matrix',
    )
onebodydensitymatrices.defaults.set(
    type = 'OneBodyDensityMatrices',name='DensityMatrices',energy_matrix=False,
    evaluator='matrix',
    )


linear.defaults.set(
     method = 'linear',move='pbyp',checkpoint=-1,
     #estimators = classcollection(localenergy)
#  #jtk
#    method='linear',move='pbyp',checkpoint=-1,gpu=True,
#    energy=0, reweightedvariance=0, unreweightedvariance=0,
#    warmupsteps       = 20,
#    usedrift          = True,
#    timestep          = .5,
#    minmethod         ='rescale',
#    stepsize          = .5,
#    beta              = 0.05,
#    alloweddifference = 1e-8,
#    bigchange         = 1.1,
#    cgsteps           = 3,
#    eigcg             = 1,
#    exp0              = -6,
#    maxweight         = 1e9,
#    minwalkers        = .5,
#    nstabilizers      = 10,
#    stabilizerscale   = .5,
#    usebuffer         = True,
    )
cslinear.defaults.set(
    method='cslinear', move='pbyp', checkpoint=-1,
    #estimators = classcollection(localenergy)
  #jtk
    #method='cslinear',move='pbyp',checkpoint=-1,gpu=True,
    #energy=0,reweightedvariance=0,unreweightedvariance=1.,
    #warmupsteps=5,steps=2,usedrift=True,timestep=.5,
    #minmethod='quartic',gevmethod='mixed',exp0=-15,
    #nstabilizers=5,stabilizerscale=3,stepsize=.35,
    #alloweddifference=1e-5,beta=.05,bigchange=5.,
    #estimators=classcollection(localenergy)
  #lschulen
    #method='cslinear', move='pbyp', checkpoint=-1, gpu=True,
    #energy=0, reweightedvariance=0, unreweightedvariance=0,
    #warmupsteps       = 20,
    ##steps             = 5,
    #usedrift          = True,
    #timestep          = .8,
    #nonlocalpp        = False,
    #minmethod         = 'rescale',
    #stepsize          = .4,
    #beta              = .05,
    #gevmethod         = 'mixed',
    #alloweddifference = 1e-4,
    #bigchange         = 9.,
    #exp0              = -16,
    #max_its           = 1,
    #maxweight         = 1e9,
    #minwalkers        = .5,
    #nstabilizers      = 3,
    #stabilizerscale   = 1,
    #usebuffer         = False,
    #estimators = classcollection(localenergy)
  #jmm
    #method='cslinear', move='pbyp', checkpoint=-1, gpu=True,
    #energy=0, reweightedvariance=0, unreweightedvariance=0,
    #warmupsteps       = 20,
    #usedrift          = True,
    #timestep          = .5,
    #nonlocalpp        = True,
    #minmethod         = 'quartic',
    #stepsize          = .4,
    #beta              = 0.0,
    #gevmethod         = 'mixed',
    #alloweddifference = 1.0e-4,
    #bigchange         = 9.0,
    #exp0              = -16,
    #max_its           = 1,
    #maxweight         = 1e9,
    #minwalkers        = 0.5,
    #nstabilizers      = 3,
    #stabilizerscale   = 1.0,
    #usebuffer         = True,
    #estimators = classcollection(localenergy)
    )
vmc.defaults.set(
    method='vmc',move='pbyp',
    #walkers     = 1,
    #warmupsteps = 50,
    #substeps    = 3,
    #usedrift    = True,
    #timestep    = .5,
    #estimators = classcollection(localenergy)
    )
dmc.defaults.set(
    method='dmc',move='pbyp',
    #warmupsteps   = 20,
    #timestep      = .01,
    #nonlocalmoves = True,
    #estimators = classcollection(localenergy)
    )
vmc_batch.defaults.set(
    method='vmc_batch',move='pbyp',
    )
dmc_batch.defaults.set(
    method='dmc_batch',move='pbyp',
    )
linear_batch.defaults.set(
    method='linear_batch',move='pbyp',
    )



# afqmc defaults
afqmcinfo.defaults.set(
    name = 'info0',
    )
walkerset.defaults.set(
    name = 'wset0',
    )
propagator.defaults.set(
    name = 'prop0',
    info = 'info0',
    )
execute.defaults.set(
    info = 'info0',
    ham  = 'ham0',
    wfn  = 'wfn0',
    wset = 'wset0',
    prop = 'prop0',
    )
back_propagation.defaults.set(
    name='back_propagation'
    )





def set_rsqmc_mode():
    QIobj.afqmc_mode = False
    Names.use_rsqmc_expanded_names()
#end def set_rsqmc_mode

def set_afqmc_mode():
    QIobj.afqmc_mode = True
    Names.use_afqmc_expanded_names()
#end def set_afqmc_mode




class QmcpackInput(SimulationInput,Names):
    
    profile_collection = None

    opt_methods = set(['opt','linear','cslinear','linear_batch'])

    simulation_type = simulation

    default_metadata = meta(
        lattice    = dict(units='bohr'),
        reciprocal = dict(units='2pi/bohr'),
        ionid      = dict(datatype='stringArray'),
        position   = dict(datatype='posArray', condition=0)
        )

    @staticmethod
    def settings(**kwargs):
        QIobj.settings(**kwargs)
    #end def settings

    def __init__(self,arg0=None,arg1=None):
        Param.metadata = None
        filepath = None
        metadata = None
        element  = None
        if arg0==None and arg1==None:
            None
        elif isinstance(arg0,str) and arg1==None:
            filepath = arg0
        elif isinstance(arg0,QIxml) and arg1==None:
            element = arg0
        elif isinstance(arg0,meta) and isinstance(arg1,QIxml):
            metadata = arg0
            element  = arg1
        else:
            self.error('input arguments of types '+arg0.__class__.__name__+' and '+arg0.__class__.__name__+' cannot be used to initialize QmcpackInput')
        #end if
        if metadata!=None:
            self._metadata = metadata
        else:
            self._metadata = meta()
        #end if
        if filepath!=None:
            self.read(filepath)
        elif element!=None:
            #simulation = arg0
            #self.simulation = self.simulation_type(simulation)
            elem_class = element.__class__
            if elem_class.identifier!=None:
                name = elem_class.identifier
            else:
                name = elem_class.__name__
            #end if
            self[name] = elem_class(element)
        #end if
        Param.metadata = None
        QIcollections.clear()
    #end def __init__

    def is_afqmc_input(self):
        is_afqmc = False
        if 'simulation' in self:
            sim = self.simulation
            is_afqmc = 'method' in sim and sim.method.lower()=='afqmc'
        #end if
        return is_afqmc
    #end def is_afqmc_input

    def get_base(self):
        elem_names = list(self.keys())
        elem_names.remove('_metadata')
        if len(elem_names)>1:
            self.error('qmcpack input cannot have more than one base element\n  You have provided '+str(len(elem_names))+': '+str(elem_names))
        #end if
        return self[elem_names[0]]
    #end def get_base

    def get_basename(self):
        elem_names = list(self.keys())
        elem_names.remove('_metadata')
        if len(elem_names)>1:
            self.error('qmcpack input cannot have more than one base element\n  You have provided '+str(len(elem_names))+': '+str(elem_names))
        #end if
        return elem_names[0]
    #end def get_basename

    def read(self,filepath=None,xml=None):
        if xml!=None or os.path.exists(filepath):
            element_joins=['qmcsystem']
            element_aliases=dict(loop='qmc')
            xml = XMLreader(filepath,element_joins,element_aliases,warn=False,xml=xml).obj
            xml.condense()
            self._metadata = meta() #store parameter/attrib attribute metadata
            Param.metadata = self._metadata
            if 'simulation' in xml:
                self.simulation = simulation(xml.simulation)
            else:
                #try to determine the type
                elements = []
                keys = []
                error = False
                for key,value in xml.items():
                    if isinstance(key,str) and key[0]!='_':
                        if key in types:
                            elements.append(types[key](value))
                            keys.append(key)
                        else:
                            self.error('element '+key+' is not a recognized type',exit=False)
                            error = True
                        #end if
                    #end if
                #end for
                if error:
                    self.error('cannot read input xml file')
                #end if
                if len(elements)==0:
                    self.error('no valid elements were found for input xml file')
                #end if
                for i in range(len(elements)):
                    elem = elements[i]
                    key  = keys[i]
                    if isinstance(elem,QIxml):
                        if elem.identifier!=None:
                            name = elem.identifier
                        else:
                            name = elem.tag
                        #end if
                    else:
                        name = key
                    #end if
                    self[name] = elem
                #end for
            #end if
            Param.metadata = None
        else:
            self.error('the filepath you provided does not exist.\n  Input filepath: '+filepath)
        #end if
        return self
    #end def read


    def write_text(self,filepath=None):
        set_rsqmc_mode()
        if self.is_afqmc_input():
            set_afqmc_mode()
        #end if
        c = ''
        header = '''<?xml version="1.0"?>
'''
        c+= header
        if len(self._metadata)==0:
            Param.metadata = self.default_metadata
        else:
            Param.metadata = self._metadata
        #end if
        base = self.get_base()
        c+=base.write(first=True)
        Param.metadata = None
        set_rsqmc_mode()
        return c
    #end def write_text


    def unroll_calculations(self,modify=True):
        qmc = []
        sim = self.simulation
        if 'calculations' in sim:
            calcs = sim.calculations
        elif 'qmc' in sim:
            calcs = [sim.qmc]
        else:
            calcs = []
        #end if
        for i in range(len(calcs)):
            c = calcs[i]
            if isinstance(c,loop):
                qmc.extend(c.unroll())
            else:
                qmc.append(c)
            #end if
        #end for
        qmc = make_collection(qmc)
        if modify:
            self.simulation.calculations = qmc
        #end if
        return qmc
    #end def unroll_calculations

    def get(self,*names):
        base = self.get_base()
        return base.get(names)
    #end def get

    def remove(self,*names):
        base = self.get_base()
        base.remove(*names)
    #end def remove
    
    def assign(self,**kwargs):
        base = self.get_base()
        base.assign(**kwargs)
    #end def assign
    
    def replace(self,*args,**kwargs):# input is list of keyword=(oldval,newval)
        base = self.get_base()
        base.replace(*args,**kwargs)
    #end def replace

    def move(self,**elemdests):
        base = self.get_base()
        base.move(**elemdests)
    #end def move
            

    def get_host(self,names):
        base = self.get_base()
        return base.get_host(names)
    #end if

    def incorporate_defaults(self,elements=False,overwrite=False,propagate=False):
        base = self.get_base()
        base.incorporate_defaults(elements,overwrite,propagate)
    #end def incorporate_defaults

    def pluralize(self):
        base = self.get_base()
        base.pluralize()
    #end def pluralize

    def standard_placements(self):
        self.move(particleset='qmcsystem',wavefunction='qmcsystem',hamiltonian='qmcsystem')
    #end def standard_placements

    def difference(self,other):
        s1 = self.copy()
        s2 = other.copy()
        b1 = s1.get_basename()
        b2 = s2.get_basename()
        q1 = s1[b1]
        q2 = s2[b2]
        if b1!=b2:
            different = True
            d1 = q1
            d2 = q2
            diff = None
        else:
            s1.standard_placements()
            s2.standard_placements()
            s1.pluralize()
            s2.pluralize()
            different,diff,d1,d2 = q1.difference(q2,root=False)
        #end if
        if diff!=None:
            diff.remove_empty()
        #end if
        d1.remove_empty()
        d2.remove_empty()
        return different,diff,d1,d2
    #end def difference

    def remove_empty(self):
        base = self.get_base()
        base.remove_empty()
    #end def remove_empty

    def read_xml(self,filepath=None,xml=None):
        if os.path.exists(filepath):
            element_joins=['qmcsystem']
            element_aliases=dict(loop='qmc')
            if xml is None:
                xml = XMLreader(filepath,element_joins,element_aliases,warn=False).obj
            else:
                xml = XMLreader(None,element_joins,element_aliases,warn=False,xml=xml).obj
            #end if
            xml.condense()
        else:
            self.error('the filepath you provided does not exist.\n  Input filepath: '+filepath)
        #end if
        return xml
    #end def read_xml

    def include_xml(self,xmlfile,replace=True,exists=True):
        xml = self.read_xml(xmlfile)
        Param.metadata = self._metadata
        for name,exml in xml.items():
            if not name.startswith('_'):
                qxml = types[name](exml)
                qname = qxml.tag
                host = self.get_host(qname)
                if host==None and exists:
                    self.error('host xml section for '+qname+' not found','QmcpackInput')
                #end if
                if qname in host:
                    section_name = qname
                elif qname in plurals_inv and plurals_inv[qname] in host:
                    section_name = plurals_inv[qname]
                else:
                    section_name = None
                #end if
                if replace:
                    if section_name!=None:
                        del host[section_name]
                    #end if
                    host[qname] = qxml
                else:
                    if section_name==None:
                        host[qname] = qxml
                    else:
                        section = host[section_name]
                        if isinstance(section,collection):
                            section[qxml.identifier] = qxml
                        elif section_name in plurals_inv:
                            coll = collection()
                            coll[section.identifier] = section
                            coll[qxml.identifier]    = qxml
                            del host[section_name]
                            host[plurals_inv[section_name]] = coll
                        else:
                            section.combine(qxml)
                        #end if
                    #end if
                #end if
            #end if
        #end for
        Param.metadata = None
    #end def include_xml

    # This include functionality is currently not being used
    # The rationale is essentially this:
    #   -Having includes explicitly represented in the input file object
    #    makes it very difficult to search for various components
    #    i.e. where is the particleset? the wavefunction? a particular determinant?
    #   -Difficulty in locating components makes it difficult to modify them
    #   -Includes necessarily introduce greater variability in input file structure
    #    and it is difficult to ensure every possible form is preserved each and 
    #    every time a modification is made
    #   -The only time it is undesirable to incorporate the contents of an
    #    include directly into the input file object is if the data is large
    #    e.g. for an xml wavefunction or pseudopotential.
    #    In these cases, an external file should be provided that contains
    #    only the large object in question (pseudo or wavefunction).
    #    This is already done for pseudopotentials and should be done for 
    #    wavefunctions, e.g. multideterminants.
    #    Until that time, wavefunctions will be explicitly read into the full
    #    input file.
    def add_include(self,element_type,href,placement='on'):
        # check the element type
        elems = ['cell','ptcl','wfs','ham']
        emap  = obj(
            simulationcell = 'cell',
            particleset    = 'ptcl',
            wavefunction   = 'wfs',
            hamiltonian    = 'ham'
            )
        if not element_type in elems:
            self.error('cannot add include for element of type {0}\n  valid element types are {1}'.format(element_type,elems))
        #end if
        # check the requested placement
        placements = ('before','on','after')
        if not placement in placements:
            self.error('cannot add include for element with placement {0}\n  valid placements are {1}'.format(placement,list(placements)))
        #end if
        # check that the base element is a simulation
        base = self.get_base()
        if not isinstance(base,simulation):
            self.error('an include can only be added to simulation\n  attempted to add to {0}'.format(base.__class__.__name__))
        #end if
        # gather a list of current qmcsystems
        if 'qmcsystem' in base:
            qslist = [(0,base.qmcsystem)]
            del base.qmcsystem
        elif 'qmcsystems' in base:
            qslist = base.qmcsystems.pairlist()
            del base.qmcsystems
        else:
            qslist = []
        #end if
        # organize the elements of the qmcsystems
        cur_elems = obj()
        for elem in elems:
            for place in placements:
                cur_elems[elem,place] = None
            #end for
        #end for
        for qskey,qs in qslist:
            if isinstance(qs,include):
                inc = qs
                ekey = qskey.split('_')[1]
                if not ekey in elems:
                    self.error('encountered invalid element key: {0}\n  valid keys are: {1}'.format(ekey,elems))
                #end if
                if cur_elems[ekey,'on'] is None:
                    cur_elems[ekey,'before'] = ekey,inc
                else:
                    cur_elems[ekey,'after' ] = ekey,inc
                #end if
            elif not isinstance(qs,qmcsystem):
                self.error('expected qmcsystem element, got {0}'.format(qs.__class__.__name__))
            else:
                for elem in qmcsystem.elements:
                    elem_plural = elem+'s'
                    name  = None
                    if elem in qs:
                        name = elem
                    elif elem_plural in qs:
                        name = elem_plural
                    #end if
                    if name!=None:
                        cur_elems[emap[elem],'on'] = name,qs[name]
                        del qs[name]
                    #end if
                #end for
                residue = list(qs.keys())
                if len(residue)>0:
                    self.error('extra keys found in qmcsystem: {0}'.format(sorted(residue)))
                #end if
            #end if
        #end for
        for elem in elems:
            pbef = cur_elems[elem,'before']
            pon  = cur_elems[elem,'on'    ]
            paft = cur_elems[elem,'after' ] 
            if pon is None:
                if not pbef is None and paft is None:
                    cur_elems[elem,'on'    ] = pbef
                    cur_elems[elem,'before'] = None
                elif not paft is None and pbef is None:
                    cur_elems[elem,'on'    ] = paft
                    cur_elems[elem,'after' ] = None
                #end if
            #end if
        #end for
        # insert the new include
        inc_name  = 'include_'+element_type
        inc_value = include(href=href)
        cur_elems[element_type,placement] = inc_name,inc_value
        # create a collection of qmcsystems
        qmcsystems = collection()
        qskey = ''
        qs    = qmcsystem()
        for elem in elems:
            for place in placements:
                cur_elem = cur_elems[elem,place]
                if cur_elem!=None:
                    name,value = cur_elem
                    if isinstance(value,include):
                        if len(qskey)>0:
                            qmcsystems.add(qs,key=qskey)
                            qskey = ''
                            qs    = qmcsystem()
                        #end if
                        qmcsystems.add(value,key=name)
                    else:
                        qskey += elem[0]
                        qs[name] = value
                    #end if
                #end if
            #end for
        #end for
        if len(qskey)>0:
            qmcsystems.add(qs,key=qskey)
        #end if
        # attach the collection to the input file
        base.qmcsystems = qmcsystems
    #end def add_include


    def get_output_info(self,*requests):
        project = self.simulation.project
        prefix = project.id
        series = project.series

        qmc = []
        calctypes = set()
        outfiles = []

        if not self.is_afqmc_input():
            qmc_ur = self.unroll_calculations(modify=False)
            n=0
            for qo in qmc_ur:
                q = obj()
                q.prefix = prefix
                q.series = series+n
                n+=1
                method = qo.method
                if method in self.opt_methods:
                    q.type = 'opt'
                else:
                    q.type = method
                #end if
                calctypes.add(q.type)
                q.method = method
                fprefix = prefix+'.s'+str(q.series).zfill(3)+'.'
                files = obj()
                files.scalar = fprefix+'scalar.dat'
                files.stat   = fprefix+'stat.h5'
                # apparently this one is no longer generated by default as of r5756
                #files.config = fprefix+'storeConfig.h5' 
                if q.type=='opt':
                    files.opt = fprefix+'opt.xml'
                elif q.type=='dmc':
                    files.dmc = fprefix+'dmc.dat'
                #end if
                outfiles.extend(files.values())
                q.files = files
                qmc.append(q)
            #end for
        else:
            q = obj()
            q.prefix = prefix
            q.series = series
            q.type   = 'afqmc'
            q.method = 'afqmc'
            calctypes.add(q.type)
            fprefix = prefix+'.s'+str(q.series).zfill(3)+'.'
            files = obj()
            files.scalar = fprefix+'scalar.dat'
            outfiles.extend(files.values())
            q.files = files
            qmc.append(q)
        #end if

        res = dict(qmc=qmc,calctypes=calctypes,outfiles=outfiles)

        values = []
        for req in requests:
            if req in res:
                values.append(res[req])
            else:
                self.error(req+' is not a valid output info request')
            #end if
        #end for
        if len(values)==1:
            return values[0]
        else:
            return values
        #end if
    #end def get_output_info


    def generate_jastrows(self,size=None,j1func='bspline',j1size=8,j2func='bspline',j2size=8):
        if size!=None:
            j1size = size
            j2size = size
        #end if

        #self.remove('jastrow')
        lattice,particlesets,wavefunction = self.get('lattice','particleset','wavefunction')
        no_lattice = lattice==None
        no_particleset = particlesets==None
        no_wavefunction = wavefunction==None
        if no_lattice:
            self.error('a simulationcell lattice must be present to generate jastrows',exit=False)
        #end if
        if no_particleset:
            self.error('a particleset must be present to generate jastrows',exit=False)
        #end if
        if no_wavefunction:
            self.error('a wavefunction must be present to generate jastrows',exit=False)
        #end if
        if no_lattice or no_particleset or no_wavefunction:
            self.error('jastrows cannot be generated')
        #end if
        if isinstance(particlesets,QIxml):
            particlesets = make_collection([particlesets])
        #end if
        if not 'e' in particlesets:
            self.error('electron particleset (e) not found\n particlesets: '+str(particlesets.keys()))
        #end if


        jastrows = collection()

        cell = Structure(lattice)
        volume = cell.volume()
        rcut   = cell.rmin()

        #use the rpa jastrow for electrons (modeled after Luke's tool)
        size = j2size
        e = particlesets.e
        nelectrons = 0
        for g in e.groups:
            nelectrons += g.size
        #end for
        density = nelectrons/volume
        wp = sqrt(4*pi*density)
        dr = rcut/size
        r = .02 + dr*arange(size)
        uuc = .5/(wp*r)*(1.-exp(-r*sqrt(wp/2)))*exp(-(2*r/rcut)**2)
        udc = .5/(wp*r)*(1.-exp(-r*sqrt(wp)))*exp(-(2*r/rcut)**2)
        jastrows.J2 = jastrow2(
            name = 'J2',type='Two-Body',function=j2func,print='yes',
            correlations = collection(
                uu = correlation(speciesA='u',speciesB='u',size=size,
                                 coefficients=section(id='uu',type='Array',coeff=uuc)),
                ud = correlation(speciesA='u',speciesB='d',size=size,
                                 coefficients=section(id='ud',type='Array',coeff=udc))
                )
            )

        #generate electron-ion jastrows, if ions present
        ions = []
        for name in particlesets.keys():
            if name=='i' or name.startswith('ion'):
                ions.append(name)
            #end if
        #end for
        if len(ions)>0:
            size = j1size
            j1 = []
            for ion in ions:
                i = particlesets[ion]
                if 'group' in i:
                    groups = [i.group]
                else:
                    groups = i.groups
                #end if
                corr = []
                for g in groups:
                    elem = g.name
                    c=correlation(
                        elementtype=elem,
                        cusp=0.,
                        size=size,
                        coefficients=section(
                            id='e'+elem,
                            type='Array',
                            coeff=size*[0]
                            )
                        )
                    corr.append(c)
                #end for
                j=jastrow1(
                    name='J1_'+ion,
                    type='One-Body',
                    function=j1func,
                    source=ion,
                    print='yes',
                    correlations = corr
                    )
                j1.append(j)
            #end for
            if len(j1)==1:
                j1[0].name='J1'
            #end if
            for j in j1:
                jastrows[j.name]=j
            #end for
        #end if

        if 'J2' in wavefunction.jastrows:
            J2 = wavefunction.jastrows.J2
            if 'function' in J2 and J2.function.lower()=='bspline':
                c = wavefunction.jastrows.J2.correlations
                ctot = abs(array(c.uu.coefficients.coeff)).sum() + abs(array(c.ud.coefficients.coeff)).sum()
                if ctot < 1e-3:
                    wavefunction.jastrows.J2 = jastrows.J2
                #end if
            #end if
        #end if

        #only add the jastrows if ones of the same type 
        # (one-body,two-body,etc) are not already present
        for jastrow in jastrows:
            jtype = jastrow.type.lower().replace('-','_')
            has_jtype = False
            for wjastrow in wavefunction.jastrows:
                wjtype = wjastrow.type.lower().replace('-','_')
                has_jtype = has_jtype or wjtype==jtype
            #end for
            if not has_jtype:
                wavefunction.jastrows[jastrow.name] = jastrow
            #end if
        #end for
    #end def generate_jastrows


    def incorporate_system(self,system):
        self.warn('incorporate_system may or may not work\n  please check the qmcpack input produced\n  if it is wrong, please contact the developer')
        system = system.copy()
        system.check_folded_system()
        system.change_units('B')
        #system.structure.group_atoms()
        system.structure.order_by_species()
        particles  = system.particles
        structure  = system.structure
        net_charge = system.net_charge
        net_spin   = system.net_spin

        qs,sc,ham,ps = self.get('qmcsystem','simulationcell','hamiltonian','particleset')

        old_eps_name = None
        old_ips_name = None
        if ps!=None:
            if isinstance(ps,particleset):
                ps = make_collection([ps])
            #end if
            for pname,pset in ps.items():
                g0name = list(pset.groups.keys())[0]
                g0 = pset.groups[g0name]
                if abs(-1-g0.charge)<1e-2:
                    old_eps_name = pname
                elif 'ionid' in pset:
                    old_ips_name = pname
                #end if
            #end for
        #end if
        del ps
        self.remove('particleset')
        if qs==None:
            qs = qmcsystem()
            qs.incorporate_defaults(elements=False,propagate=False)
            self.simulation.qmcsystem = qs
        #end if
        if sc==None:
            sc = simulationcell()
            sc.incorporate_defaults(elements=False,propagate=False)
            qs.simulationcell = sc
        #end if
        if ham==None:
            ham = hamiltonian()
            ham.incorporate_defaults(elements=False,propagate=False)
            qs.hamiltonian = ham
        elif isinstance(ham,collection):
            if 'h0' in ham:
                ham = ham.h0
            elif len(ham)==1:
                ham = ham.list()[0]
            else:
                self.error('cannot find hamiltonian for system incorporation')
            #end if
        #end if

        elem = structure.elem
        pos  = structure.pos

        if len(structure.axes)>0: #exclude systems with open boundaries
            #setting the 'lattice' (cell axes) requires some delicate care
            #  qmcpack will fail if this is even 1e-10 off of what is in 
            #  the wavefunction hdf5 file from pwscf
            if structure.folded_structure!=None:
                fs = structure.folded_structure
                axes = array(pwscf_array_string(fs.axes).split(),dtype=float)
                axes.shape = fs.axes.shape
                axes = dot(structure.tmatrix,axes)
                if abs(axes-structure.axes).sum()>1e-5:
                    self.error('supercell axes do not match tiled version of folded cell axes\n  you may have changed one set of axes (super/folded) and not the other\n  folded cell axes:\n'+str(fs.axes)+'\n  supercell axes:\n'+str(structure.axes)+'\n  folded axes tiled:\n'+str(axes))
                #end if
            else:
                axes = array(pwscf_array_string(structure.axes).split(),dtype=float)
                axes.shape = structure.axes.shape
            #end if
            structure.adjust_axes(axes)

            sc.lattice = axes
        #end if    

        elns = particles.get_electrons()
        ions = particles.get_ions()
        eup  = elns.up_electron
        edn  = elns.down_electron

        particlesets = []
        eps = particleset(
            name='e',random=True,
            groups = [
                group(name='u',charge=-1,mass=eup.mass,size=eup.count),
                group(name='d',charge=-1,mass=edn.mass,size=edn.count)
                ]
            )
        particlesets.append(eps)
        if len(ions)>0:
            if sc!=None and 'bconds' in sc and tuple(sc.bconds)!=('p','p','p'):
                eps.randomsrc = 'ion0'  
            #end if
            ips = particleset(
                name='ion0',
                )
            groups = []
            ham.pluralize()
            pseudos = ham.get('pseudo')
            if pseudos==None:
                pp = ham.get('PseudoPot')
                if pp!=None:
                    pseudos = collection()
                    pp.pseudos = pseudos
                #end if
            #end if
            for ion in ions:
                gpos = pos[elem==ion.name]
                g = group(
                    name         = ion.name,
                    charge       = ion.charge,
                    valence      = ion.charge,
                    atomicnumber = ion.protons,
                    mass         = ion.mass,
                    position     = gpos,
                    size         = len(gpos)
                    )
                groups.append(g)
                if pseudos!=None and not ion.name in pseudos:
                    pseudos[ion.name] = pseudo(elementtype=ion.name,href='MISSING.xml')
                #end if
            #end for
            ips.groups = make_collection(groups)
            particlesets.append(ips)
        #end if
        qs.particlesets = make_collection(particlesets)

        if old_eps_name!=None:
            self.replace(old_eps_name,'e')
        #end if
        if old_ips_name!=None and len(ions)>0:
            self.replace(old_ips_name,'ion0')
        #end if
            
        udet,ddet = self.get('updet','downdet')

        if udet!=None:
            udet.size = elns.up_electron.count
        #end if
        if ddet!=None:
            ddet.size = elns.down_electron.count
        #end if

        if abs(net_spin) > 1e-1:
            if ddet!=None:
                if 'occupation' in ddet:
                    ddet.occupation.spindataset = 1
                else:
                    ss = self.get('sposets')
                    ss[ddet.sposet].spindataset = 1
                #end if
            #end if
        #end if
    #end def incorporate_system
        
    def get_electron_particle_set(self):

        input = self.copy()
        input.pluralize()
        return input.get('particlesets').e

    #end def get_electron_particle_set

    def return_system(self,structure_only=False):
        input = self.copy()
        input.pluralize()
        axes,ps,H = input.get('lattice','particlesets','hamiltonian')

        if ps is None:
            return None
        #end if

        # find electrons and ions
        have_ions = True
        have_jellium = False
        ions = None
        elns = None
        ion_list = []
        for name,p in ps.items():
            if 'ionid' in p:
                ion_list.append(p)
            elif name.startswith('e'):
                elns = p
            #end if
        #end for
        if len(ion_list)==0: #try to identify ions by positive charged groups
            for name,p in ps.items():
                if 'groups' in p:
                    for g in p.groups:
                        if 'charge' in g and g.charge>0:
                            ion_list.append(p)
                            break
                        #end if
                    #end for
                #end if
            #end for
        #end if
        if len(ion_list)==1:
            ions = ion_list[0]
        elif len(ion_list)>1:
            self.error('ability to handle multiple ion particlesets has not been implemented')
        #end if
        if ions is None and elns!=None and 'groups' in elns:
            simcell = input.get('simulationcell')
            if simcell!=None and 'rs' in simcell:
                have_ions = False
                have_jellium = True
            elif not 'pairpots' in H:
                have_ions = False
            #end if
        #end if
        if elns==None:
            self.error('could not find electron particleset')
        #end if
        if ions==None and have_ions:
            self.error('could not find ion particleset')
        #end if

        #compute spin and electron charge
        net_spin   = 0
        eln_charge = 0
        for spin,eln in elns.groups.items():
            if spin[0]=='u':
                net_spin+=eln.size
            elif spin[0]=='d':
                net_spin-=eln.size
            #end if
            eln_charge += eln.charge*eln.size
        #end if

        #get structure and ion charge
        structure  = None
        ion_charge = 0
        valency    = dict()
        if have_ions:
            elem = None
            if 'ionid' in ions:
                if isinstance(ions.ionid,str):
                    elem = [ions.ionid]
                else:
                    elem = list(ions.ionid)
                #end if
                pos  = ions.position
            elif 'size' in ions and ions.size==1:
                elem = [ions.groups.list()[0].name]
                pos  = [[0,0,0]]
            elif 'groups' in ions:
                elem = []
                pos  = []
                for group in ions.groups:
                    if 'position' in group:
                        nions = group.size
                        elem.extend(nions*[group.name])
                        if group.size==1:
                            pos.extend([list(group.position)])
                        else:
                            pos.extend(list(group.position))
                        #end if
                    #end if
                #end for
                if len(elem)==0:
                    elem = None
                    pos  = None
                else:
                    elem  = array(elem)
                    pos   = array(pos)
                    order = elem.argsort()
                    elem  = elem[order]
                    pos   = pos[order]
                #end if
            #end if
            if elem is None:
                self.error('could not read ions from ion particleset')
            #end if
            if axes is None:
                center = (0,0,0)
            else:
                md = input._metadata
                if 'position' in md and 'condition' in md['position'] and md['position']['condition']==1:
                    pos = dot(pos,axes)
                #end if
                center = axes.sum(0)/2
            #end if

            # pos must be a 2D array, shape (N,3)
            # reshape single atom case, shape (3,) as shape (1,3)
            pos = np.asarray(pos)
            if len(pos.flatten())==3:
                pos.shape = (1,3)
            #end if

            structure = Structure(axes=axes,elem=elem,pos=pos,center=center,units='B')

            for name,element in ions.groups.items():
                if 'charge' in element:
                    valence = element.charge
                elif 'valence' in element:
                    valence = element.valence
                elif 'atomic_number' in element:
                    valence = element.atomic_number
                else:
                    self.error('could not identify valency of '+name)
                #end if
                valency[name] = valence
                count = list(elem).count(name)
                ion_charge += valence*count
            #end for
        elif have_jellium:
            structure  = Jellium(rs=simcell.rs,background_charge=-eln_charge)
            ion_charge = structure.background_charge
        #end if

        net_charge = ion_charge + eln_charge

        system = PhysicalSystem(structure,net_charge,net_spin,**valency) 
        
        if structure_only:
            return structure
        else:
            return system
        #end if
    #end def return_system


    def get_ion_particlesets(self):
        ions = obj()
        ps = self.get('particlesets')
        #try to identify ions by positive charged groups
        for name,p in ps.items():
            if name.startswith('ion') or name.startswith('atom'):
                ions[name] = p
            elif 'groups' in p:
                for g in p.groups:
                    if 'charge' in g and g.charge>0:
                        ions[name] = p
                        break
                    #end if
                #end for
            #end if
        #end for
        return ions
    #end def get_ion_particlesets


    def get_pp_files(self):
        pp_files = []
        h = self.get('hamiltonian')
        if h != None:
            pp = None
            if 'pairpots' in h:
                for pairpot in h.pairpots:
                    if 'type' in pairpot and pairpot.type=='pseudo':
                        pp = pairpot
                    #end if
                #end for
            elif 'pairpot' in h and 'type' in h.pairpot and h.pairpot.type=='pseudo':
                pp = h.pairpot
            #end if
            if pp!=None:
                if 'pseudo' in pp and 'href' in pp.pseudo:
                    pp_files.append(pp.pseudo.href)
                elif 'pseudos' in pp:
                    for pseudo in pp.pseudos:
                        if 'href' in pseudo:
                            pp_files.append(pseudo.href)
                        #end if
                    #end for
                #end if
            #end if
        #end if
        return pp_files
    #end def get_pp_files


    def remove_physical_system(self):
        qs = self.simulation.qmcsystem
        if 'simulationcell' in qs:
            del qs.simulationcell
        #end if
        if 'particlesets' in qs:
            del qs.particlesets
        #end if
        for name in qs.keys():
            if isinstance(qs[name],particleset):
                del qs[name]
            #end if
        #end for
        self.replace('ion0','i')
    #end def remove_physical_system

        
    def cusp_correction(self):
        cc = False
        if not self.is_afqmc_input():
            ds = self.get('determinantset')
            cc_var = ds!=None and 'cuspcorrection' in ds and ds.cuspcorrection==True
            cc_run = len(self.simulation.calculations)==0
            cc = cc_var and cc_run
        #end if
        return cc
    #end def cusp_correction


    def get_qmc(self,series):
        qmc = None
        calcs        = self.get('calculations')
        series_start = self.get('series')
        if calcs!=None:
            if series_start is None:
                qmc = calcs[series]
            else:
                qmc = calcs[series-series_start]
            #end if
        #end if
        return qmc
    #end def get_qmc


    def bundle(self,inputs,filenames):
        return BundledQmcpackInput(inputs,filenames)
    #end def bundle

    
    def trace(self,quantity,values):
        return TracedQmcpackInput(quantity,values,self)
    #end def trace

    
    def twist_average(self,twistnums):
        return self.trace('twistnum',twistnums)
    #end def twist_average
#end class QmcpackInput



# base class for bundled qmcpack input
#  not used on its own
class BundledQmcpackInput(SimulationInput):
    
    def __init__(self,inputs,filenames):
        self.inputs = obj()
        for input in inputs:
            self.inputs.append(input)
        #end for
        self.filenames = filenames
    #end def __init__


    def get_output_info(self,*requests):
        outfiles = []

        for index,input in self.inputs.items():
            outfs = input.get_output_info('outfiles')
            infile = self.filenames[index]
            outfile= infile.rsplit('.',1)[0]+'.g'+str(index).zfill(3)+'.qmc'
            outfiles.append(infile)
            outfiles.append(outfile)
            twfile = infile.rsplit('.',3)[0]+'.twist_info.dat'
            outfiles.append(twfile)
            for outf in outfs:
                prefix,rest = outf.split('.',1)
                outfiles.append(prefix+'.g'+str(index).zfill(3)+'.'+rest)
            #end for
        #end for

        values = []
        for req in requests:
            if req=='outfiles':
                values.append(outfiles)
            else:
                values.append(None)
            #end if
        #end for
        if len(values)==1:
            return values[0]
        else:
            return values
        #end if
    #end def get_output_info

        
    def generate_filenames(self,infile):
        self.not_implemented()
    #end def generate_filenames
        

    def write(self,filepath=None):
        if filepath!=None and not 'filenames' in self:
            infile = os.path.split(filepath)[1]
            if not infile.endswith('.xml'):
                infile+='.xml'
            #end if
            self.generate_filenames(infile)
        #end if
        if filepath==None:
            c = ''
            for i in range(len(self.inputs)):
                c += self.filenames[i]+'\n'
            #end for
            return c
        else:
            path,file  = os.path.split(filepath)
            #if file!=self.filenames[-1]:
            #    self.error('main filenames do not match\n  internal: '+self.filenames[-1]+'\n  inputted: '+file)
            ##end if
            c = ''
            for i in range(len(self.inputs)):
                input = self.inputs[i]
                bfile = self.filenames[i]
                c += bfile+'\n'
                bfilepath = os.path.join(path,bfile)
                input.write(bfilepath)
            #end for
            fobj = open(filepath,'w')
            fobj.write(c)
            fobj.close()
        #end if
    #end def write
#end class BundledQmcpackInput



class TracedQmcpackInput(BundledQmcpackInput):
    def __init__(self,quantity=None,values=None,input=None):
        self.quantities = obj()
        self.variables = obj()
        self.inputs = obj()
        self.filenames = None
        if quantity!=None and values!=None and input!=None:
            self.bundle_inputs(quantity,values,input)
        #end if
    #end def __init__

    def bundle_inputs(self,quantity,values,input):
        range = len(self.inputs),len(self.inputs)+len(values)
        self.quantities.append(obj(quantity=quantity,range=range))
        for value in values:
            inp = input.copy()
            qhost = inp.get_host(quantity)                               
            print(qhost)
            if qhost!=None:
                qhost[quantity] = value
            else:
                self.error('quantity '+quantity+' was not found in '+input.__class__.__name__)
            #end if
            self.variables.append(obj(quantity=quantity,value=value))
            self.inputs.append(inp)
        #end for
    #end def bundle_inputs


    def generate_filenames(self,infile):
        prefix,ext = infile.split('.',1)
        if not ext.endswith('xml'):
            ext+='.xml'
        #end if
        self.filenames = []
        for i in range(len(self.variables)):
            var = self.variables[i]
            q = var.quantity
            v = var.value
            bfile = prefix+'.g'+str(i).zfill(3)+'.'+q+'_'+str(v)+'.'+ext
            self.filenames.append(bfile)
        #end if
        self.filenames.append(prefix+'.in')
    #end def generate_filenames
#end class TracedQmcpackInput




class QmcpackInputTemplate(SimulationInputTemplate):
    def preprocess(self,contents,filepath=None):
        if filepath!=None:
            basepath,filename = os.path.split(filepath)
            c = contents
            contents=''
            for line in c.splitlines():
                if '<include' in line and '/>' in line:
                    tokens = line.replace('<include','').replace('/>','').split()
                    for token in tokens:
                        if token.startswith('href'):
                            include_file = token.replace('href','').replace('=','').replace('"','').strip()
                            include_path = os.path.join(basepath,include_file)
                            if os.path.exists(include_path):
                                icont = open(include_path,'r').read()+'\n'
                                line = ''
                                for iline in icont.splitlines():
                                    if not '<?' in iline:
                                        line+=iline+'\n'
                                    #end if
                                #end for
                            #end if
                        #end if
                    #end for
                #end if
                contents+=line+'\n'
            #end for
        #end if
        return contents
    #end def preprocess


    def get_output_info(self,*args,**kwargs):
        # just pretend
        return []
    #end def get_output_info
#end class QmcpackInputTemplate




def generate_simulationcell(bconds='ppp',lr_dim_cutoff=15,lr_tol=None,lr_handler=None,system=None):
    bconds = tuple(bconds)
    sc = simulationcell(bconds=bconds)
    periodic = 'p' in bconds
    axes_valid = system!=None and len(system.structure.axes)>0
    if periodic:
        sc.lr_dim_cutoff = lr_dim_cutoff
        if lr_tol is not None:
            sc.lr_tol = lr_tol
        #end if
        if lr_handler is not None:
            sc.lr_handler = lr_handler
        #end if
        if not axes_valid:
            QmcpackInput.class_error('invalid axes in generate_simulationcell\nargument system must be provided\naxes of the structure must have non-zero dimension')
        #end if
    #end if
    if axes_valid:
        system.check_folded_system()
        system.change_units('B')
        structure = system.structure
        if isinstance(structure,Jellium):
            sc.rs         = structure.rs()
            sc.nparticles = system.particles.count_electrons()
        else:
            #setting the 'lattice' (cell axes) requires some delicate care
            #  qmcpack will fail if this is even 1e-10 off of what is in 
            #  the wavefunction hdf5 file from pwscf
            if structure.folded_structure!=None:
                fs = structure.folded_structure
                axes = array(pwscf_array_string(fs.axes).split(),dtype=float)
                axes.shape = fs.axes.shape
                axes = dot(structure.tmatrix,axes)
                if abs(axes-structure.axes).sum()>1e-5:
                    QmcpackInput.class_error('in generate_simulationcell\nsupercell axes do not match tiled version of folded cell axes\nyou may have changed one set of axes (super/folded) and not the other\nfolded cell axes:\n'+str(fs.axes)+'\nsupercell axes:\n'+str(structure.axes)+'\nfolded axes tiled:\n'+str(axes))
                #end if
            else:
                axes = array(pwscf_array_string(structure.axes).split(),dtype=float)
                axes.shape = structure.axes.shape
            #end if
            structure.adjust_axes(axes)

            sc.lattice = axes
        #end if
    #end if    
    return sc
#end def generate_simulationcell


def generate_particlesets(electrons   = 'e',
                          ions        = 'ion0',
                          up          = 'u',
                          down        = 'd',
                          system      = None,
                          randomsrc   = True,
                          hybrid_rcut = None,
                          hybrid_lmax = None,
                          ):
    if system is None:
        QmcpackInput.class_error('generate_particlesets argument system must not be None')
    #end if

    ename = electrons
    iname = ions
    uname = up
    dname = down

    del electrons
    del ions
    del up
    del down

    system.check_folded_system()
    system.change_units('B')

    particles  = system.particles
    structure  = system.structure
    net_charge = system.net_charge
    net_spin   = system.net_spin

    elns = particles.get_electrons()
    ions = particles.get_ions()
    eup  = elns.up_electron
    edn  = elns.down_electron

    particleset_groups = []
    if eup.count > 0:
        particleset_groups.append(group(name=uname,charge=-1,mass=eup.mass,size=eup.count))
    #end if
    if edn.count > 0:
        particleset_groups.append(group(name=dname,charge=-1,mass=edn.mass,size=edn.count))
    #end if

    particlesets = []
    eps = particleset(
        name   = ename,
        random = True,
        groups = particleset_groups,
        )
    particlesets.append(eps)
    if len(ions)>0:
        # maintain consistent order
        ion_species,ion_counts = structure.order_by_species()
        elem = structure.elem
        pos  = structure.pos
        if randomsrc:
            eps.randomsrc = iname
        #end if
        ips = particleset(name=iname)
        # handle hybrid rep
        hybridrep = hybrid_rcut is not None or hybrid_lmax is not None
        if hybridrep:
            hybrid_vars = (
                ('hybrid_rcut',hybrid_rcut),
                ('hybrid_lmax',hybrid_lmax),                
                )
            for hvar,hval in hybrid_vars:
                if not isinstance(hval,obj):
                    QmcpackInput.class_error('generate_particlesets argument "{0}" must be of type obj\nyou provided type: {1}\nwith value: {2}'.format(hvar,hval.__class__.__name__,hval))
                #end if
                if set(hval.keys())!=set(ion_species):
                    QmcpackInput.class_error('generate_particsets argument "{0}" is incorrect\none entry must be present for each atomic species\natomic species present in the simulation: {1}\nvalues provided for the following species: {2}'.format(hvar,sorted(ion_species),sorted(hval.keys())))
                #end if
            #end for
        #end if
        # make groups
        groups = []
        for ion_spec in ion_species:
            ion = ions[ion_spec]
            gpos = pos[elem==ion.name]
            g = group(
                name         = ion.name,
                charge       = ion.charge,
                valence      = ion.charge,
                atomicnumber = ion.protons,
                mass         = ion.mass,
                position     = gpos,
                size         = len(gpos)
                )
            if hybridrep:
                g.lmax           = hybrid_lmax[ion_spec]
                g.cutoff_radius  = hybrid_rcut[ion_spec]
            #end if
            groups.append(g)
        #end for
        ips.groups = make_collection(groups)
        particlesets.append(ips)
    #end if
    particlesets = make_collection(particlesets)
    return particlesets
#end def generate_particlesets


def generate_sposets(type           = None,
                     occupation     = None,
                     spin_polarized = False,
                     nup            = None,
                     ndown          = None,
                     spo_up         = 'spo_u',
                     spo_down       = 'spo_d',
                     system         = None,
                     sposets        = None,
                     spindatasets   = False):
    ndn = ndown
    if type is None:
        QmcpackInput.class_error('cannot generate sposets\n  type of sposet not specified')
    #end if
    if sposets!=None:
        for spo in sposets:
            spo.type = type
        #end for
    elif occupation=='slater_ground':
        have_counts = not (nup is None or ndown is None)
        if system is None and not have_counts:
            QmcpackInput.class_error('cannot generate sposets in occupation mode {0}\n  arguments nup & ndown or system must be given to generate_sposets'.format(occupation))
        elif not have_counts:
            elns = system.particles.get_electrons()
            nup  = elns.up_electron.count
            ndn  = elns.down_electron.count
        #end if
        if not spin_polarized:
            if nup==ndn:
                sposets = [sposet(type=type,name='spo_ud',spindataset=0,size=nup)]
            else:
                sposets = [sposet(type=type,name=spo_up,  spindataset=0,size=nup),
                           sposet(type=type,name=spo_down,spindataset=0,size=ndn)]
            #end if
        else:
            sposets_list = []
            if nup > 0:
                sposets_list.append(sposet(type=type,name=spo_up,  spindataset=0,size=nup))
            #end if
            if ndn > 0:
                sposets_list.append(sposet(type=type,name=spo_down,spindataset=1,size=ndn))
            #end if
            sposets = sposets_list
        #end if
        if not spindatasets:
            for spo in sposets:
                del spo.spindataset
            #end for
        #end if
    else:
        QmcpackInput.class_error('cannot generate sposets in occupation mode {0}\n  generate_sposets currently supports the following occupation modes:\n  slater_ground'.format(occupation))
    #end if
    return make_collection(sposets)
#end def generate_sposets


def generate_sposet_builder(type,*args,**kwargs):
    if type=='bspline' or type=='einspline':
        return generate_bspline_builder(type,*args,**kwargs)
    elif type=='heg':
        return generate_heg_builder(*args,**kwargs)
    else:
        QmcpackInput.class_error('cannot generate sposet_builder\n  sposet_builder of type {0} is unrecognized'.format(type))
    #end if
#end def generate_sposet_builder


def generate_bspline_builder(type           = 'bspline',
                             meshfactor     = 1.0,
                             precision      = 'float',
                             twistnum       = None, 
                             twist          = None,
                             sort           = None,
                             version        = '0.10',
                             truncate       = False,
                             buffer         = None,
                             spin_polarized = False,
                             hybridrep      = None,
                             href           = 'MISSING.h5',
                             ions           = 'ion0',
                             spo_up         = 'spo_u',
                             spo_down       = 'spo_d',
                             sposets        = None,
                             system         = None,
                             orbitals_cpu   = None,
                             gpusharing     = None,
                             ):
    tilematrix = identity(3,dtype=int)
    if system!=None:
        tilematrix = system.structure.tilematrix()
    #end if
    bsb = bspline_builder(
        type       = type,
        meshfactor = meshfactor,
        precision  = precision,
        tilematrix = tilematrix,
        href       = href,
        version    = version,
        truncate   = truncate,
        source     = ions,
        sposets    = generate_sposets(
            type           = type,
            occupation     = 'slater_ground',
            spin_polarized = spin_polarized,
            system         = system,
            sposets        = sposets,
            spindatasets   = True
            )
        )
    if sort is not None:
        bsb.sort = sort
    #end if
    if truncate and buffer!=None:
        bsb.buffer = buffer
    #end if
    if hybridrep is not None:
        bsb.hybridrep = hybridrep
    #end if
    if twist is not None:
        bsb.twistnum = system.structure.select_twist(twist)
    elif twistnum is not None:
        bsb.twistnum = twistnum
    elif len(system.structure.kpoints)==1:
        bsb.twistnum = 0
    else:
        bsb.twistnum = None
    #end if
    if orbitals_cpu is not None and orbitals_cpu:
        bsb.gpu = False
    #end if
    if gpusharing is not None:
        bsb.gpusharing = gpusharing
    #end if
    return bsb
#end def generate_bspline_builder


def generate_heg_builder(twist          = None,
                         spin_polarized = False,
                         spo_up         = 'spo_u',
                         spo_down       = 'spo_d',
                         sposets        = None,
                         system         = None
                         ):
    type = 'heg'
    hb = heg_builder(
        type    = type,
        sposets = generate_sposets(
            type           = type,
            occupation     = 'slater_ground',
            spin_polarized = spin_polarized,
            system         = system,
            sposets        = sposets
            )
        )
    if twist!=None:
        hb.twist = tuple(twist)
    #end if
    return hb
#end def generate_heg_builder


def partition_sposets(sposet_builder,partition,partition_meshfactors=None):
    ssb = sposet_builder
    spos_in =ssb.sposets
    del ssb.sposets
    if isinstance(partition,(dict,obj)):
        partition_indices  = sorted(partition.keys())
        partition_contents = partition
    else:
        partition_indices  = list(partition)
        partition_contents = None
    #end if
    if partition_meshfactors is not None:
        if partition_contents is None:
            partition_contents = obj()
            for p in partition_indices:
                partition_contents[p] = obj()
            #end for
        #end if
        for p,mf in zip(partition_indices,partition_meshfactors):
            partition_contents[p].meshfactor = mf
        #end for
    #end if
    # partition each spo in the builder and create a corresponding composite spo
    comp_spos = []
    part_spos = []
    for spo in spos_in.list():
        part_spo_names = []
        part_ranges = partition_indices+[spo.size]
        for i in range(len(partition_indices)):
            index_min = part_ranges[i]
            index_max = part_ranges[i+1]
            if index_min>spo.size:
                break
            elif index_max>spo.size:
                index_max = spo.size
            #end if
            part_spo_name = spo.name+'_'+str(index_min)
            part_spo = sposet(**spo)
            part_spo.name = part_spo_name
            if index_min==0:
                part_spo.size = index_max
            else:
                part_spo.index_min = index_min
                part_spo.index_max = index_max
                del part_spo.size
            #end if
            if partition_contents is not None:
                part_spo.set(**partition_contents[index_min])
            #end if
            part_spos.append(part_spo)
            part_spo_names.append(part_spo_name)
        #end for
        comp_spo = sposet(
            name = spo.name,
            size = spo.size,
            spos = part_spo_names,
            )
        comp_spos.append(comp_spo)
    #end for
        
    ssb.sposets = make_collection(part_spos)

    cssb = composite_builder(
        type = 'composite',
        sposets = make_collection(comp_spos),
        )

    return [ssb,cssb]
#end def partition_sposets


def generate_determinantset(up             = 'u',
                            down           = 'd',
                            spo_up         = 'spo_u',
                            spo_down       = 'spo_d',
                            spin_polarized = False,
                            delay_rank     = None,
                            matrix_inv_cpu = None,
                            system         = None
                            ):
    if system is None:
        QmcpackInput.class_error('generate_determinantset argument system must not be None')
    #end if
    elns = system.particles.get_electrons()
    nup  = elns.up_electron.count
    ndn  = elns.down_electron.count
    if not spin_polarized and nup==ndn:
        spo_u = 'spo_ud'
        spo_d = 'spo_ud'
    else:
        spo_u = spo_up
        spo_d = spo_down
    #end if
    determinants_list = []
    if nup > 0:
        determinants_list.append(
            determinant(
                id     = 'updet',
                group  = up,
                sposet = spo_u,
                size   = nup
                )
        )
    #end if
    if ndn > 0:
        determinants_list.append(
            determinant(
                id     = 'downdet',
                group  = down,
                sposet = spo_d,
                size   = ndn
                )
        )
    #end if
    dset = determinantset(
        slaterdeterminant = slaterdeterminant(
            determinants = collection(*determinants_list)
            )
        )
    if delay_rank is not None:
        dset.slaterdeterminant.delay_rank = delay_rank
    #end if
    if matrix_inv_cpu is not None and matrix_inv_cpu:
        dset.slaterdeterminant.matrix_inverter = 'host'
    #end if
    return dset
#end def generate_determinantset


def check_excitation_type(excitation):

    # Possible spin channels or spin states
    exc_spins = obj(
        up      = 1, # 'up'
        down    = 2, # 'down'
        singlet = 3, # 'singlet'
        triplet = 4, # 'triplet'
        )
    # Possible orbital excitation types
    exc_types = obj(
        band    = 1, # '0 45 3 46'   # Type 1
        energy  = 2, # '-215 +216'   # Type 2
        kpoint  = 3, # 'L vb F cb'   # Type 3
        lowest  = 4, # 'lowest'      # Type 4
        )

    exc_spin = None
    exc_type = None

    # Check that 'excitation' is correctly formated
    format_failed = False
    # Extract elements form excitation
    if not isinstance(excitation,(tuple,list)) or len(excitation) != 2:
        format_failed = True
    else:
        exc1,exc2 = excitation
        if not isinstance(exc1,str) or not isinstance(exc2,str):
            format_failed = True
        #end if
    #end if

    # Check first element
    if not format_failed:
        if exc1.lower() not in ('up','down','singlet','triplet'):
            format_failed = True
        else:
            exc_spin = exc_spins[exc1.lower()]
        #end if
    #end if

    # Check second element
    if not format_failed:
        if any(substr in exc2.lower() for substr in ('vb','cb','lowest')):
            if exc2.lower()=='lowest':
                exc_type = exc_types.lowest
            elif len(exc2.split())!=4:
                format_failed = True
            else:
                exc_type = exc_types.kpoint
            #end if
        else:
            tmp = None
            try:
                tmp = array(exc2.split(),dtype=int)
            except:
                format_failed = True
            #end try
            if not tmp is None:
                if len(tmp)==4:
                    # '0 45 3 46'
                    if not tmp[0]>=0 or not tmp[1]>=0 or not tmp[2]>=0 or not tmp[3]>=0:
                        format_failed = True
                    #end if
                    exc_type = exc_types.band
                elif len(tmp)==2:
                    # '-215 +216'
                    if not tmp[0]<0 or not tmp[1]>0:
                        format_failed = True
                    #end if
                    exc_type = exc_types.energy
                else:
                    format_failed = True
                #end if
            #end if
        #end if
    #end if
    
    if format_failed:

        msg  = 'excitation must be a tuple or list with with two elements.\n'
        msg += 'The first element must be either "up", "down", "singlet", or "triplet"\n'
        msg += 'and the second element must be a band format (e.g. "0 45 3 46"),\n'
        msg += 'energy format (e.g. "-215 +216"), kpoint format (e.g. "L vb F cb"),\n'
        msg += 'or lowest format (e.g. "lowest").\n'
        msg += 'You Provided: {0}'
        msg = msg.format(excitation)

        QmcpackInput.class_error(msg)

    #end if

    return exc_spin,exc_type,exc_spins,exc_types,exc1,exc2
#end def check_excitation_type


def generate_determinantset_old(type           = 'bspline',
                                meshfactor     = 1.0,
                                precision      = 'float',
                                twistnum       = None, 
                                twist          = None,
                                spin_polarized = False,
                                hybridrep      = None,
                                source         = 'ion0',
                                href           = 'MISSING.h5',
                                excitation     = None,
                                delay_rank     = None,
                                gpusharing     = None,
                                system         = None
                                ):
    if system is None:
        QmcpackInput.class_error('generate_determinantset argument system must not be None')
    #end if
    elns = system.particles.get_electrons()
    down_spin = 0
    if spin_polarized:
        down_spin=1
    #end if
    tilematrix = identity(3,dtype=int)
    if system!=None:
        tilematrix = system.structure.tilematrix()
    #end if
    nup = elns.up_electron.count
    ndn = elns.down_electron.count
    determinants_list = []
    if nup > 0:
        determinants_list.append(
            determinant(
                id   = 'updet',
                size = nup,
                occupation=section(mode='ground',spindataset=0)
                ),
        )
    #end if
    if ndn > 0:
        determinants_list.append(
            determinant(
                id   = 'downdet',
                size = ndn,
                occupation=section(mode='ground',spindataset=down_spin)
                )
        )
    #end if
    dset = determinantset(
        type       = type,
        meshfactor = meshfactor,
        precision  = precision,
        tilematrix = tilematrix,
        href       = href,
        source     = source,
        slaterdeterminant = slaterdeterminant(
            determinants = collection(*determinants_list)
            )
        )
    if twist is not None:
        dset.twistnum = system.structure.select_twist(twist)
    elif twistnum is not None:
        dset.twistnum = twistnum
    elif len(system.structure.kpoints)==1:
        dset.twistnum = 0
    else:
        dset.twistnum = None
    #end if
    if hybridrep is not None:
        if hybridrep=='yes' or hybridrep=='no':
            dset.hybridrep = hybridrep
        else:
            dset.hybridrep = yesno_dict[hybridrep]
        #end if
    #end if
    if delay_rank is not None:
        dset.slaterdeterminant.delay_rank = delay_rank
    #end if
    if gpusharing is not None:
        dset.gpusharing = gpusharing
    #end if
    if excitation is not None:

        exc_spin,exc_type,exc_spins,exc_types,exc1,exc2 = check_excitation_type(excitation)

        if exc_spin==exc_spins.up:
            sdet = dset.get('updet')
        elif exc_spin==exc_spins.down:
            sdet = dset.get('downdet')
        elif exc_spin in (exc_spins.singlet,exc_spins.triplet):

            # Are there an equal number of up and down electrons?
            # If no, then exit. Currently, singlet and triplet 
            # excitations are assumed to have ms = 0.
            if elns.down_electron.count != elns.up_electron.count:
                QmcpackInput.class_error('The \'singlet\' and \'triplet\' excitation types currently assume number of up and down electrons is the same for the reference ground state. Otherwise, one should use \'up\' or \'down\' types.\nFor your system: Nup={} and Ndown={}.\nWe plan to expand to additional cases in the future.'.format(elns.up_electron.count,elns.down_electron.count))
            #end if

            coeff_sign = ''
            if exc_spin==exc_spins.triplet:
                coeff_sign = '-'
            #end if

            if down_spin:
                sposet_list = [sposet(name            = 'spo_u',
                                      spindataset     = 0,
                                      size            = elns.up_electron.count+1,
                                      occupation      = section(mode='ground'),
                                      coefficient     = section(size=90,spindataset=0),
                                      spos            = ''
                                     ),
                               sposet(name            = 'spo_d',
                                      spindataset     = 1,
                                      size            = elns.up_electron.count+1,
                                      occupation      = section(mode='ground'),
                                      coefficient     = section(spindataset=1),
                                      spos            = ''
                                     )]
            else:
                sposet_list = [sposet(name            = 'spo_ud',
                                      spindataset     = 0,
                                      size            = elns.up_electron.count+1,
                                      occupation      = section(mode='ground'),
                                      coefficient     = section(spindataset=0),
                                      spos            = ''
                                     )]
            #end if

            dset = determinantset(
                type       = type,
                meshfactor = meshfactor,
                precision  = precision,
                tilematrix = tilematrix,
                twistnum   = twistnum,
                href       = href,
                source     = source,
                sposets    = sposet_list,
                multideterminant = multideterminant(
                    optimize = 'no',
                    spo_up='spo_u' if down_spin else 'spo_ud',
                    spo_dn='spo_d' if down_spin else 'spo_ud',
                    detlist = detlist(
                        size = '1',
                        type = 'CSF',
                        nca  = '0',
                        ncb  = '0',
                        nea = elns.up_electron.count,
                        neb = elns.down_electron.count,
                        cutoff = '0.001',
                        csf = csf(
                            id          = 'CSF_0',
                            exctLvl     = '1',
                            coeff       = '1.0',
                            coeff_real  = '1.0',
                            qchem_coeff = '1.0',
                            dets = collection(
                                det(
                                    id='csf_00',
                                    coeff='0.70710678118654752440',
                                    ),
                                det(
                                    id='csf_01',
                                    coeff=coeff_sign+'0.70710678118654752440',
                                    ),
                                )
                            )
                        ),
                    )
                )
            
            if exc_type in (exc_types.energy,exc_types.lowest):

                nup = elns.up_electron.count 
                if exc_type==exc_types.lowest:
                    exc_orbs = [nup,nup+1]
                else:
                    # assume excitation of form '-216 +217' or '-216 217'
                    exc_orbs = array(exc2.split(),dtype=int)
                    exc_orbs[0] *= -1
                #end if

                for sp in dset.sposets:
                    sp.size=exc_orbs[1]
                #end for

                dset.multideterminant.detlist.nstates = exc_orbs[1]

                dset.multideterminant.detlist.csf.occ = '2'*nup+'0'*(exc_orbs[1]-nup-1)+'1'
                dset.multideterminant.detlist.csf.occ = dset.multideterminant.detlist.csf.occ[:exc_orbs[0]-1]+'1'+dset.multideterminant.detlist.csf.occ[exc_orbs[0]:]

                dset.multideterminant.detlist.csf.dets[0].alpha = '1'*(exc_orbs[0]-1)+'0'+'1'*(nup-exc_orbs[0])+'0'*(exc_orbs[1]-nup-1)+'1'
                dset.multideterminant.detlist.csf.dets[0].beta = '1'*nup+'0'*(exc_orbs[1]-nup)

                dset.multideterminant.detlist.csf.dets[1].alpha = '1'*nup+'0'*(exc_orbs[1]-nup)
                dset.multideterminant.detlist.csf.dets[1].beta = '1'*(exc_orbs[0]-1)+'0'+'1'*(nup-exc_orbs[0])+'0'*(exc_orbs[1]-nup-1)+'1'

            elif exc_type == exc_types.kpoint: 
                QmcpackInput.class_error('{} excitation is not yet available for kpoint type'.format(exc1))
            else: 
                QmcpackInput.class_error('{} excitation is not yet available for band type'.format(exc1))
            #end if

            return dset

        #end if

        occ = sdet.occupation
        occ.pairs    = 1
        occ.mode     = 'excited'
        occ.contents = '\n'+exc2+'\n'
        # add new input format
        if exc_type == exc_types.kpoint:
            # assume excitation of form 'gamma vb k cb' or 'gamma vb-1 k cb+1'
            excitation = exc2.upper().split(' ')
            if len(excitation) == 4:
                k_1, band_1, k_2, band_2 = excitation
            else:
                QmcpackInput.class_error('excitation with vb-cb band format works only with special k-points')
            #end if
            
            vb = int(sdet.size / abs(linalg.det(tilematrix))) -1  # Separate for each spin channel
            cb = vb+1
            # Convert band_1, band_2 to band indexes
            bands = [band_1, band_2]
            for bnum, b in enumerate(bands):
                b = b.lower()
                if 'cb' in b:
                    if '-' in b:
                        b = b.split('-')
                        bands[bnum] = cb - int(b[1])
                    elif '+' in b:
                        b = b.split('+')
                        bands[bnum] = cb + int(b[1])
                    else:
                        bands[bnum] = cb
                    #end if
                elif 'vb' in b:
                    if '-' in b:
                        b = b.split('-')
                        bands[bnum] = vb - int(b[1])
                    elif '+' in b:
                        b = b.split('+')
                        bands[bnum] = vb + int(b[1])
                    else:
                        bands[bnum] = vb
                    #end if
                else:
                    QmcpackInput.class_error('{0} in excitation has the wrong formatting'.format(b))
                #end if
            #end for
            band_1, band_2 = bands
            
            # Convert k_1 k_2 to wavevector indexes
            structure = system.structure.get_smallest().copy()
            structure.change_units('A')
            kpath       = get_kpath(structure=structure)
            kpath_label = array(kpath['explicit_kpoints_labels'])
            kpath_rel   = kpath['explicit_kpoints_rel']
            
            k1_in = k_1
            k2_in = k_2
            if k_1 in kpath_label and k_2 in kpath_label:   
                k_1 = kpath_rel[where(kpath_label == k_1)][0]
                k_2 = kpath_rel[where(kpath_label == k_2)][0]

                #kpts = nscf.input.k_points.kpoints
                kpts = structure.kpoints_unit()
                found_k1 = False
                found_k2 = False
                for knum, k in enumerate(kpts):
                    if isclose(k_1, k).all():
                        k_1 = knum
                        found_k1 = True
                    #end if
                    if isclose(k_2, k).all():
                        k_2 = knum
                        found_k2 = True
                    #end if
                #end for
                if not found_k1 or not found_k2:
                    QmcpackInput.class_error('Requested special kpoint is not in the tiled cell\nRequested "{}", present={}\nRequested "{}", present={}\nAvailable kpoints: {}'.format(k1_in,found_k1,k2_in,found_k2,sorted(set(kpath_label))))
                #end if
            else:
                QmcpackInput.class_error('Excitation wavevectors are not found in the kpath\nlabels requested: {} {}\nlabels present: {}'.format(k_1,k_2,sorted(set(kpath_label))))
            #end if

            #Write everything in band (ti,bi) format
            occ.contents = '\n'+str(k_1)+' '+str(band_1)+' '+str(k_2)+' '+str(band_2)+'\n'
            occ.format = 'band'
            
        elif exc_type == exc_types.energy:
            # assume excitation of form '-216 +217'
            occ.format = 'energy'
        elif exc_type == exc_types.lowest: # Type 4
            occ.format = 'energy'
            if exc_spin == exc_spins.up:
                nel = elns.up_electron.count 
            else:
                nel = elns.down_electron.count 
            #end if
            excitation = '-{} +{}'.format(nel,nel+1) 
            occ.contents = '\n'+excitation+'\n'
        else: #Type 1
            # assume excitation of form '6 36 6 37'
            occ.format   = 'band'
        #end if
    #end if
    return dset
#end def generate_determinantset_old


def generate_hamiltonian(name         = 'h0',
                         type         = 'generic',
                         electrons    = 'e',
                         ions         = 'ion0',
                         wavefunction = 'psi0',
                         pseudos      = None,
                         dla          = None,
                         format       = 'xml',
                         estimators   = None,
                         system       = None,
                         interactions = 'default',
                         ):
    if system is None:
        QmcpackInput.class_error('generate_hamiltonian argument system must not be None')
    #end if

    ename   = electrons
    iname   = ions
    wfname  = wavefunction
    ppfiles = pseudos
    del electrons
    del ions
    del pseudos
    del wavefunction

    particles = system.particles
    if particles.count_electrons()==0:
        QmcpackInput.class_error('cannot generate hamiltonian, no electrons present')
    #end if

    pairpots = []
    if interactions!=None:
        pairpots.append(coulomb(name='ElecElec',type='coulomb',source=ename,target=ename))
        if particles.count_ions()>0:
            pairpots.append(coulomb(name='IonIon',type='coulomb',source=iname,target=iname))
            ions = particles.get_ions()
            if not system.pseudized:
                pairpots.append(coulomb(name='ElecIon',type='coulomb',source=iname,target=ename))
            else:
                if ppfiles is None or len(ppfiles)==0:
                    QmcpackInput.class_error('cannot generate hamiltonian\n  system is pseudized, but no pseudopotentials have been provided\n  please provide pseudopotential files via the pseudos keyword')
                #end if
                if isinstance(ppfiles,list):
                    pplist = ppfiles
                    ppfiles = obj()
                    for pppath in pplist:
                        if '/' in pppath:
                            ppfile = pppath.split('/')[-1]
                        else:
                            ppfile = pppath
                        #end if
                        element = ppfile.split('.')[0]
                        if len(element)>2:
                            element = element[0:2]
                        #end if
                        ppfiles[element] = pppath
                    #end for
                #end if
                pseudos = collection()
                for ion in ions:
                    label = ion.name
                    iselem,symbol = is_element(ion.name,symbol=True)
                    if label in ppfiles:
                        ppfile = ppfiles[label]
                    elif symbol in ppfiles:
                        ppfile = ppfiles[symbol]
                    else:
                        QmcpackInput.class_error('pseudos provided to generate_hamiltonian are incomplete\n  a pseudopotential for ion of type {0} is missing\n  pseudos provided:\n{1}'.format(ion.name,str(ppfiles)))
                    #end if
                    pseudos.add(pseudo(elementtype=label,href=ppfile))
                #end for
                pp = pseudopotential(name='PseudoPot',type='pseudo',source=iname,wavefunction=wfname,format=format,pseudos=pseudos)
                if dla is not None:
                    pp.dla = dla
                #end if
                pairpots.append(pp)
            #end if
        #end if
    #end if

    ests = []
    if estimators is not None:
        for estimator in estimators:
            if isinstance(estimator,QIxml):
                estimator = estimator.copy()
            #end if
            est=estimator
            if isinstance(estimator,str):
                estname = estimator.lower().replace(' ','_').replace('-','_').replace('__','_')
                if estname=='mpc':
                    pairpots.append(mpc(name='MPC',type='MPC',ecut=60.0,source=ename,target=ename,physical=False))
                    est = None
                elif estname=='chiesa':
                    est = chiesa(name='KEcorr',type='chiesa',source=ename,psi=wfname)
                elif estname=='localenergy':
                    est = localenergy(name='LocalEnergy')
                elif estname=='skall':
                    est = skall(name='SkAll',type='skall',source=iname,target=ename,hdf5=True)
                elif estname=='energydensity':
                    est = energydensity(
                        type='EnergyDensity',name='EDvoronoi',dynamic=ename,static=iname,
                        spacegrid = spacegrid(coord='voronoi')
                        )
                elif estname=='pressure':
                    est = pressure(type='Pressure')
                else:
                    QmcpackInput.class_error('estimator '+estimator+' has not yet been enabled in generate_basic_input')
                #end if
            elif not isinstance(estimator,QIxml):
                QmcpackInput.class_error('generate_hamiltonian received an invalid estimator\n  an estimator must either be a name or a QIxml object\n  inputted estimator type: {0}\n  inputted estimator contents: {1}'.format(estimator.__class__.__name__,estimator))
            elif isinstance(estimator,energydensity):
                est.set_optional(
                    type = 'EnergyDensity',
                    dynamic = ename,
                    static  = iname,
                    )
            elif isinstance(estimator,dm1b):
                est = process_dm1b_estimator(estimator,wfname)
            #end if
            if est!=None:
                ests.append(est)
            #end if
        #end for
    #end if
    estimators = ests

    hmltn = hamiltonian(
        name   = name,
        type   = type,
        target = ename
        )

    if len(pairpots)>0:
        hmltn.pairpots = make_collection(pairpots)
    #end if

    if len(estimators)>0:
        hmltn.estimators = make_collection(estimators)
    #end if

    return hmltn
#end def generate_hamiltonian


def generate_estimators_batched(estimators,
                                electrons    = 'e',
                                ions         = 'ion0',
                                wavefunction = 'psi0',
                                ):
    assert len(estimators)>0
    ename  = electrons
    iname  = ions
    wfname = wavefunction
    del electrons
    del ions

    ests = []
    for estimator in estimators:
        if isinstance(estimator,QIxml):
            estimator = estimator.copy()
        #end if
        est = estimator
        if isinstance(estimator,str):
            estname = estimator.lower().replace(' ','_').replace('-','_').replace('__','_')
            #if estname=='chiesa':
            #    est = chiesa(name='KEcorr',type='chiesa',source=ename,psi=wfname)
            #else:
            QmcpackInput.class_error('estimator '+estimator+' has not yet been enabled in generate_estimators')
            ##end if
        elif not isinstance(estimator,QIxml):
                QmcpackInput.class_error('generate_estimators received an invalid estimator\n  an estimator must either be a name or a QIxml object\n  inputted estimator type: {0}\n  inputted estimator contents: {1}'.format(estimator.__class__.__name__,estimator))
        elif isinstance(estimator,momentum):
            estimator.type = 'MomentumDistribution'
        elif isinstance(estimator,onebodydensitymatrices):
            est = process_dm1b_estimator(estimator,wfname)
        #end if
        if est is not None:
            ests.append(est)
        #end if
    #end for
    estimators = make_collection(ests)
    return estimators
#end def generate_estimators_batched


def process_dm1b_estimator(dm,wfname):
    reuse = False
    if 'reuse' in dm:
        reuse = bool(dm.reuse)
        del dm.reuse
    #end if
    basis = []
    builder = None
    maxed = False
    if reuse and 'basis' in dm and isinstance(dm.basis,sposet):
        spo = dm.basis
        # get sposet size
        if 'size' in dm.basis:
            size = spo.size
            del spo.size
        elif 'index_max' in dm.basis:
            size = spo.index_max
            del spo.index_max
        else:
            QmcpackInput.class_error('cannot generate estimator dm1b\n  basis sposet provided does not have a "size" attribute')
        #end if
        try:
            # get sposet from wavefunction
            wf = QIcollections.get('wavefunctions',wfname)
            dets = wf.get('determinant')
            det  = dets.get_single()
            if 'sposet' in det:
                rsponame = det.sposet
            else:
                rsponame = det.id
            #end if
            builders = QIcollections.get('sposet_builders')
            if builders is None:
                builders = [wf.sposet_builders.bspline]
            #end if
            rspo = None
            for bld in builders:
                if rsponame in bld.sposets:
                    builder = bld
                    rspo    = bld.sposets[rsponame]
                    break
                #end if
            #end for
            basis.append(rsponame)
            # adjust current sposet
            spo.index_min = rspo.size
            spo.index_max = size
            maxed = rspo.size>=size
        except Exception as e:
            msg = 'cannot generate estimator dm1b\n  '
            if wf is None:
                QmcpackInput.class_error(msg+'wavefunction {0} not found'.format(wfname))
            elif dets is None or det is None:
                QmcpackInput.class_error(msg+'determinant not found')
            elif builders is None:
                QmcpackInput.class_error(msg+'sposet_builders not found')
            elif rspo is None:
                QmcpackInput.class_error(msg+'sposet {0} not found'.format(rsponame))
            else:
                QmcpackInput.class_error(msg+'cause of failure could not be determined\n  see the following error message:\n{0}'.format(e))
            #end if
        #end if
    #end if
    # put the basis sposet in the appropriate builder
    if isinstance(dm.basis,sposet) and not maxed:
        spo = dm.basis
        del dm.basis
        if not 'type' in spo:
            QmcpackInput.class_error('cannot generate estimator dm1b\n  basis sposet provided does not have a "type" attribute')
        #end if
        if not 'name' in spo:
            spo.name = 'spo_dm'
        #end if
        builders = QIcollections.get('sposet_builders')
        if not spo.type in builders:
            bld = generate_sposet_builder(spo.type,sposets=[spo])
            builders.add(bld)
        else:
            bld = builders[spo.type]
            bld.sposets.add(spo)
        #end if
        basis.append(spo.name)
    #end if
    dm.basis = basis
    dm.incorporate_defaults(elements=False,overwrite=False,propagate=False)
    return dm
#end def process_dm1b_estimator



def generate_jastrows(jastrows,system=None,return_list=False,check_ions=False):
    jin = []
    have_ions = True
    if check_ions and system!=None:
        have_ions = system.particles.count_ions()>0
    #end if
    if isinstance(jastrows,str):
        jorders = set(jastrows.replace('generate',''))
        if '1' in jorders and have_ions:
            jterm = generate_jastrow('J1','bspline',8,system=system)
        #end if
        if '2' in jorders:
            jterm = generate_jastrow('J2','bspline',8,system=system)
        #end if
        if '3' in jorders and have_ions:
            jterm = generate_jastrow('J3','polynomial',3,3,4.0,system=system)
        #end if
        if 'k' in jorders:
            kcut = max(system.rpa_kf())
            nksh = system.structure.count_kshells(kcut)
            jterm = generate_kspace_jastrow(kc1=0, kc2=kcut, nk1=0, nk2=nksh)
        #end if
        jin.append(jterm)
        if len(jin)==0:
            QmcpackInput.class_error('jastrow generation requested but no orders specified (1,2,and/or 3)')
        #end if
    else:
        jset = set(['J1','J2','J3'])
        for jastrow in jastrows:
            if isinstance(jastrow,QIxml):
                jin.append(jastrow)
            elif isinstance(jastrow,dict) or isinstance(jastrow,obj):
                jdict = dict(**jastrow)
                if not 'type' in jastrow:
                    QmcpackInput.class_error("could not determine jastrow type from input\n  field 'type' must be 'J1', 'J2', or 'J3'\n  object you provided: "+str(jastrow))
                #end if
                jtype = jdict['type']
                if not jtype in jset:
                    QmcpackInput.class_error("invalid jastrow type provided\n  field 'type' must be 'J1', 'J2', or 'J3'\n  object you provided: "+str(jdict))
                #end if
                del jdict['type']
                if 'system' in jdict:
                    jsys = jdict['system']
                    del jdict['system']
                else:
                    jsys = system
                #end if
                jterm = generate_jastrow(jtype,system=jsys,**jdict)
                if jterm is not None:
                    jin.append(jterm)
                #end if
                del jtype
                del jsys
            elif jastrow[0] in jset:
                jin.append(generate_jastrow(jastrow,system=system))
            else:
                QmcpackInput.class_error('starting jastrow unrecognized:\n  '+str(jastrow))
            #end if
        #end for
    #end if
    if return_list:
        return jin
    else:
        wf = wavefunction(jastrows=jin)
        wf.pluralize()
        return wf.jastrows
    #end if
#end def generate_jastrows



def generate_jastrows_alt(
    J1           = False,
    J2           = False,
    J3           = False,
    J1_size      = None,
    J1_rcut      = None,
    J1_dr        = 0.5,
    J2_size      = None,
    J2_rcut      = None,
    J2_dr        = 0.5,
    J2_init      = 'zero',
    J3_isize     = 3,
    J3_esize     = 3,
    J3_rcut      = 5.0,
    J1_rcut_open = 5.0,
    J2_rcut_open = 10.0,
    system       = None,
    ):
    if system is None:
        QmcpackInput.class_error('input variable "system" is required to generate jastrows','generate_jastrows_alt')
    elif system.structure.units!='B':
        system = system.copy()
        system.structure.change_units('B')
    #end if

    openbc = system.structure.is_open()

    natoms = system.particles.count_ions()
    nelec  = system.particles.count_electrons()

    jastrows = []
    J2 |= J3
    J1 |= J2
    if not J1 and not J2 and not J3:
        J1 = True
        J2 = True
    #end if
    rwigner = None
    if J1:
        if natoms<1:
            QmcpackInput.class_error('One-body Jastrow (J1) requested, but no atoms are present','generate_jastrows_alt')
        #end if
        if J1_rcut is None:
            if openbc:
                J1_rcut = J1_rcut_open
            else:
                if rwigner is None:
                    rwigner = system.structure.rwigner(1)
                #end if
                J1_rcut = rwigner 
            #end if
        #end if
        if J1_size is None:
            J1_size = int(ceil(J1_rcut/J1_dr))
        #end if
        J = generate_jastrow('J1','bspline',J1_size,J1_rcut,system=system)
        jastrows.append(J)
    #end if
    if J2:
        if nelec<2:
            QmcpackInput.class_error('Two-body Jastrow (J2) requested, but not enough electrons are present.\nElectrons required: 2 or more\nElectrons present: {}'.format(nelec),'generate_jastrows_alt')
        #end if
        if J2_rcut is None:
            if openbc:
                J2_rcut = J2_rcut_open
            else:
                if rwigner is None:
                    rwigner = system.structure.rwigner(1)
                #end if
                J2_rcut = rwigner
            #end if
        #end if
        if J2_size is None:
            J2_size = int(ceil(J2_rcut/J2_dr))
        #end if
        J = generate_jastrow('J2','bspline',J2_size,J2_rcut,init=J2_init,system=system)
        jastrows.append(J)
    #end if
    if J3:
        if natoms<1 or nelec<2:
            QmcpackInput.class_error('Three-body Jastrow (J3) requested, but not enough particles are present.\nAtoms required: 1 or more\nElectrons required: 2 or more\nAtoms present: {}\nElectrons present: {}'.format(natoms,nelec),'generate_jastrows_alt')
        #end if
        if not openbc:
            if rwigner is None:
                rwigner = system.structure.rwigner(1)
            #end if
            J3_rcut = min(J3_rcut,rwigner)
        #end if
        J = generate_jastrow('J3','polynomial',J3_esize,J3_isize,J3_rcut,system=system)
        jastrows.append(J)
    #end if

    return jastrows
#end def generate_jastrows_alt


def generate_jastrow(descriptor,*args,**kwargs):
    keywords = set(['function','size','rcut','elements','coeff','cusp','ename',
                    'iname','spins','density','Buu','Bud','system','isize','esize','init'])
    if not 'system' in kwargs:
        kwargs['system'] = None
    #end if
    system = kwargs['system']
    del kwargs['system']
    if system!=None:
        system.change_units('B')
    #end if
    if isinstance(descriptor,str):
        descriptor = [descriptor]
    #end if
    ikw=0
    for i in range(len(descriptor)):
        if descriptor[i] in keywords:
            break
        #end if
        ikw += 1
    #end for
    dargs = descriptor[1:ikw]
    if len(dargs)>0:
        args = dargs
    #end if
    for i in range(ikw,len(descriptor),2):
        d = descriptor[i]
        if isinstance(d,str):
            if d in keywords:
                kwargs[d] = descriptor[i+1]
            else:
                QmcpackInput.class_error('keyword {0} is unrecognized\n  valid options are: {1}'.format(d,str(keywords)),'generate_jastrow')
            #end if
        #end if
    #end for
    kwargs['system'] = system
    jtype = descriptor[0]
    if jtype=='J1':
        jastrow = generate_jastrow1(*args,**kwargs)
    elif jtype=='J2':
        jastrow = generate_jastrow2(*args,**kwargs)
    elif jtype=='J3':
        jastrow = generate_jastrow3(*args,**kwargs)
    else:
        QmcpackInput.class_error('jastrow type unrecognized: '+jtype)
    #end if
    return jastrow
#end def generate_jastrow



def generate_jastrow1(function='bspline',size=8,rcut=None,coeff=None,cusp=0.,ename='e',iname='ion0',elements=None,system=None,**elemargs):
    noelements = elements is None
    nosystem   = system is None
    noelemargs = len(elemargs)==0
    isopen     = False
    isperiodic = False
    rwigner = 1e99
    if noelements and nosystem and noelemargs:
        QmcpackInput.class_error('must specify elements or system','generate_jastrow1')
    #end if
    if noelements:
        elements = []
    #end if
    if not nosystem:
        elements.extend(list(set(system.structure.elem)))
        isopen     = system.structure.is_open()
        isperiodic = system.structure.is_periodic()
        if not isopen and isperiodic:
            rwigner = system.structure.rwigner()
        #end if
    #end if
    if not noelemargs:
        elements.extend(elemargs.keys())
    #end if
    # remove duplicate elements
    eset = set()
    elements = [ e for e in elements if e not in eset and not eset.add(e) ]     
    corrs = []
    for i in range(len(elements)):
        element = elements[i]
        if cusp == 'Z':
            QmcpackInput.class_error('need to implement Z cusp','generate_jastrow1')
        else:
            lcusp  = cusp
        #end if
        lrcut  = rcut
        lcoeff = size*[0]
        if coeff!=None:
            if element in coeff:
                lcoeff = coeff[element]
            else:
                lcoeff = coeff[i]
            #end if
        #end if
        if element in elemargs:
            v = elemargs[element]
            if 'cusp' in v:
                lcusp = v['cusp']
            #end if
            if 'rcut' in v:
                lrcut = v['rcut']
            #end if
            if 'size' in v and not 'coeff' in v:
                lcoeff = v['size']*[0]
            #end if
            if 'coeff' in v:
                lcoeff = v['coeff']
            #end if
        #end if
        corr = correlation(
            elementtype = element,
            size        = len(lcoeff),
            cusp        = cusp,
            coefficients=section(
                id    = ename+element,
                type  = 'Array',
                coeff = lcoeff
                )
            )            
        if lrcut!=None:
            if isperiodic and lrcut>rwigner:
                QmcpackInput.class_error('rcut must not be greater than the simulation cell wigner radius\nyou provided: {0}\nwigner radius: {1}'.format(lrcut,rwigner),'generate_jastrow1')
                
            corr.rcut = lrcut
        elif isopen:
            QmcpackInput.class_error('rcut must be provided for an open system','generate_jastrow1')
        elif isperiodic:
            corr.rcut = rwigner
        #end if
        corrs.append(corr)
    #end for
    j1 = jastrow1(
        name         = 'J1',
        type         = 'One-Body',
        function     = function,
        source       = iname,
        print       = True,
        correlations = corrs
        )
    return j1
#end def generate_jastrow1



def generate_bspline_jastrow2(size=8,rcut=None,coeff=None,spins=('u','d'),density=None,system=None,init='rpa'):
    if coeff is None and system is None and (init=='rpa' and density is None or rcut is None):
        QmcpackInput.class_error('rcut and density or system must be specified','generate_bspline_jastrow2')
    #end if
    isopen      = False
    isperiodic  = False
    allperiodic = False
    rwigner     = 1e99
    if system!=None:
        isopen      = system.structure.is_open()
        isperiodic  = system.structure.is_periodic()
        allperiodic = system.structure.all_periodic()
        if not isopen and isperiodic:
            rwigner = system.structure.rwigner()
        #end if
        volume = system.structure.volume()
        if isopen: 
            if rcut is None:
                QmcpackInput.class_error('rcut must be provided for an open system','generate_bspline_jastrow2')
            #end if
            if init=='rpa':
                init = 'zero'
            #end if
        else:
            if rcut is None and isperiodic:
                rcut = rwigner
            #end if
            nelectrons = system.particles.count_electrons()
            density = nelectrons/volume
        #end if
    elif init=='rpa':
        init = 'zero'
    #end if
    if coeff is None:
        if init=='rpa':
            if not allperiodic:
                QmcpackInput.class_error('rpa initialization can only be used for fully periodic systems','generate_bspline_jastrow2')
            #end if
            wp = sqrt(4*pi*density)
            dr = rcut/size
            r = .02 + dr*arange(size)
            uuc = .5/(wp*r)*(1.-exp(-r*sqrt(wp/2)))*exp(-(2*r/rcut)**2)
            udc = .5/(wp*r)*(1.-exp(-r*sqrt(wp))  )*exp(-(2*r/rcut)**2)
            coeff = [uuc,udc]
        elif init=='zero' or init==0:
            coeff = [size*[0],size*[0]]
        else:
            QmcpackInput.class_error(str(init)+' is not a valid value for parameter init\n  valid options are: rpa, zero','generate_bspline_jastrow2')
        #end if
    elif len(coeff)!=2:
        QmcpackInput.class_error('must provide 2 sets of coefficients (uu,ud)','generate_bspline_jastrow2')
    #end if
    size = len(coeff[0])
    uname,dname = spins
    uuname = uname+uname
    udname = uname+dname
    corrs = [
        correlation(speciesA=uname,speciesB=uname,size=size,
                    coefficients=section(id=uuname,type='Array',coeff=coeff[0])),
        correlation(speciesA=uname,speciesB=dname,size=size,
                    coefficients=section(id=udname,type='Array',coeff=coeff[1]))
        ]
    if rcut!=None:
        if isperiodic and rcut>rwigner:
            QmcpackInput.class_error('rcut must not be greater than the simulation cell wigner radius\nyou provided: {0}\nwigner radius: {1}'.format(rcut,rwigner),'generate_jastrow2')
        #end if
        for corr in corrs:
            corr.rcut=rcut
        #end for
    #end if
    j2 = jastrow2(
        name = 'J2',type='Two-Body',function='bspline',print=True,
        correlations = corrs
        )
    return j2
#end def generate_bspline_jastrow2


def generate_pade_jastrow2(Buu=None,Bud=None,spins=('u','d'),system=None):
    if Buu is None:
        Buu = 2.0
    #end if
    if Bud is None:
        Bud = float(Buu)
    #end if
    uname,dname = spins
    uuname = uname+uname
    udname = uname+dname
    cuu = var(id=uuname+'_b',name='B',value=Buu)
    cud = var(id=udname+'_b',name='B',value=Bud)
    corrs = [
        correlation(speciesA=uname,speciesB=uname,
                    vars=[cuu]),
        correlation(speciesA=uname,speciesB=dname,
                    vars=[cud])
        ]
    j2 = jastrow2(
        name = 'J2',type='Two-Body',function='pade',
        correlations = corrs
        )
    return j2
#end def generate_pade_jastrow2



def generate_jastrow2(function='bspline',*args,**kwargs):
    if not 'spins' in kwargs:
        kwargs['spins'] = ('u','d')
    #end if
    spins = kwargs['spins']
    if not isinstance(spins,tuple) and not isinstance(spins,list):
        QmcpackInput.class_error('spins must be a list or tuple of u/d spin names\n  you provided: '+str(spins))
    #end if
    if len(spins)!=2:
        QmcpackInput.class_error('name for up and down spins must be specified\n  you provided: '+str(spins))
    #end if
    if not isinstance(function,str):
        QmcpackInput.class_error('function must be a string\n  you provided: '+str(function),'generate_jastrow2')
    #end if
    if function=='bspline':
        j2 = generate_bspline_jastrow2(*args,**kwargs)
    elif function=='pade':
        j2 = generate_pade_jastrow2(*args,**kwargs)
    else:
        QmcpackInput.class_error('function is invalid\n  you provided: {0}\n  valid options are: bspline or pade'.format(function),'generate_jastrow2')
    #end if
    if 'system' in kwargs and kwargs['system'] is not None:
        nup,ndn = kwargs['system'].particles.electron_counts()
        if nup<2:
            del j2.correlations.uu
        #end if
        if nup<1 or ndn<1:
            del j2.correlations.ud
        #end if
        if len(j2.correlations)==0:
            j2=None
        #end if
    #end if
    return j2
#end def generate_jastrow2



def generate_jastrow3(function='polynomial',esize=3,isize=3,rcut=4.,coeff=None,iname='ion0',spins=('u','d'),elements=None,system=None):
    if elements is None and system is None:
        QmcpackInput.class_error('must specify elements or system','generate_jastrow3')
    elif elements is None:
        elements = list(set(system.structure.elem))
    #end if
    if coeff!=None:
        QmcpackInput.class_error('handling coeff is not yet implemented for generate jastrow3')
    #end if
    if len(spins)!=2:
        QmcpackInput.class_error('must specify name for up and down spins\n  provided: '+str(spins),'generate_jastrow3')
    #end if
    if rcut is None:
        QmcpackInput.class_error('must specify rcut','generate_jastrow3')
    #end if
    if system!=None and system.structure.is_periodic():
        rwigner = system.structure.rwigner()
        if rcut>rwigner:
            QmcpackInput.class_error('rcut must not be greater than the simulation cell wigner radius\nyou provided: {0}\nwigner radius: {1}'.format(rcut,rwigner),'generate_jastrow3')
        #end if
    #end if
    uname,dname = spins
    uuname = uname+uname
    udname = uname+dname
    corrs=[]
    for element in elements:
        corrs.append(
            correlation(
                especies1=uname,especies2=uname,ispecies=element,esize=esize,
                isize=isize,rcut=rcut,
                coefficients=section(id=uuname+element,type='Array',optimize=True))
            )
        corrs.append(
            correlation(
                especies1=uname,especies2=dname,ispecies=element,esize=esize,
                isize=isize,rcut=rcut,
                coefficients=section(id=udname+element,type='Array',optimize=True))
            )
    #end for
    jastrow = jastrow3(
        name = 'J3',type='eeI',function=function,print=True,source=iname,
        correlations = corrs
        )
    return jastrow
#end def generate_jastrow3


def generate_kspace_jastrow(kc1=0, kc2=0, nk1=0, nk2=0,
  symm1='isotropic', symm2='isotropic', coeff1=None, coeff2=None):
  """Generate <jastrow type="kSpace">

  Parameters
  ----------
    kc1 : float, optional
      kcut for one-body Jastrow, default 0
    kc2 : float, optional
      kcut for two-body Jastrow, default 0
    nk1 : int, optional
      number of coefficients for one-body Jastrow, default 0
    nk2 : int, optional
      number of coefficients for two-body Jastrow, default 0
    symm1 : str, optional
      one of ['crystal', 'isotropic', 'none'], default 'isotropic'
    symm2 : str, optional
      one of ['crystal', 'isotropic', 'none'], default is 'isotropic'
    coeff1 : list, optional
      one-body Jastrow coefficients, default None
    coeff2 : list, optional
      list, optional two-body Jastrow coefficients, default None
  Returns
  -------
    jk: QIxml
      kspace_jastrow qmcpack_input element
  """

  if coeff1 is None: coeff1 = [0]*nk1
  if coeff2 is None: coeff2 = [0]*nk2
  if len(coeff1) != nk1:
    QmcpackInput.class_error('coeff1 mismatch', 'generate_kspace_jastrow')
  #end if
  if len(coeff2) != nk2:
    QmcpackInput.class_error('coeff2 mismatch', 'generate_kspace_jastrow')
  #end if

  corr1 = correlation(
    type = 'One-Body',
    symmetry = symm1,
    kc = kc1,
    coefficients = section(
      id = 'cG1', type = 'Array',
      coeff = coeff1
    )
  )
  corr2 = correlation(
    type = 'Two-Body',
    symmetry = symm2,
    kc = kc2,
    coefficients = section(
      id = 'cG2', type = 'Array',
      coeff = coeff2
     )
  )
  jk = kspace_jastrow(
    type = 'kSpace',
    name = 'Jk',
    source = 'ion0',
    correlations = collection([corr1, corr2])
  )
  return jk
# end def generate_kspace_jastrow


def count_jastrow_params(jastrows):
    if isinstance(jastrows,QIxml):
        jastrows = [jastrows]
    #end if
    params = 0
    for jastrow in jastrows:
        name = jastrow.name
        if 'type' in jastrow:
            type = jastrow.type.lower()
        else:
            type = ''
        #end if
        jastrow.pluralize()
        if name=='J1' or type=='one-body':
            for correlation in jastrow.correlations:
                params += correlation.size
            #end for
        elif name=='J2' or type=='two-body':
            for correlation in jastrow.correlations:
                params += correlation.size
            #end for
        elif name=='J3' or type=='eeI':
            for correlation in jastrow.correlations:
                params += correlation.esize
                params += correlation.isize
            #end for
        #end if
    #end for
    return params
#end def count_jastrow_params


def generate_energydensity(
        name      = None,
        dynamic   = None,
        static    = None,
        coord     = None,
        grid      = None,
        scale     = None,
        ion_grids = None,
        system    = None,
        ):
    if dynamic is None:
        dynamic = 'e'
    #end if
    if static is None:
        static = 'ion0'
    #end if
    refp = None
    sg = []
    if coord is None:
        QmcpackInput.class_error('coord must be provided','generate_energydensity')
    elif coord=='voronoi':
        if name is None:
            name = 'EDvoronoi'
        #end if
        sg.append(spacegrid(coord=coord))
    elif coord=='cartesian':
        if name is None:
            name = 'EDcell'
        #end if
        if grid is None:
            QmcpackInput.class_error('grid must be provided for cartesian coordinates','generate_energydensity')
        #end if
        axes = [
            axis(p1='a1',scale='.5',label='x'),
            axis(p1='a2',scale='.5',label='y'),
            axis(p1='a3',scale='.5',label='z'),
            ]
        n=0
        for ax in axes:
           ax.grid = '-1 ({0}) 1'.format(grid[n])
           n+=1
        #end for
        sg.append(spacegrid(coord=coord,origin=origin(p1='zero'),axes=axes))
    elif coord=='spherical':
        if name is None:
            name = 'EDatom'
        #end if
        if ion_grids is None:
            QmcpackInput.class_error('ion_grids must be provided for spherical coordinates','generate_energydensity')
        #end if
        refp = reference_points(coord='cartesian',points='\nr1 1 0 0\nr2 0 1 0\nr3 0 0 1\n')
        if system is None:
            i=1
            for scale,g1,g2,g3 in ion_grids:
                grid = g1,g2,g3
                axes = [
                    axis(p1='r1',scale=scale,label='r'),
                    axis(p1='r2',scale=scale,label='phi'),
                    axis(p1='r3',scale=scale,label='theta'),
                    ]
                n=0
                for ax in axes:
                    ax.grid = '0 ({0}) 1'.format(grid[n])
                    n+=1
                #end for
                sg.append(spacegrid(coord=coord,origin=origin(p1=static+str(i)),axes=axes))
                i+=1
            #end for
        else:
            ig = ion_grids
            ion_grids = obj()
            for e,s,g1,g2,g3 in ig:
                ion_grids[e] = s,(g1,g2,g3)
            #end for
            missing = set(ion_grids.keys())
            i=1
            for e in system.structure.elem:
                if e in ion_grids:
                    scale,grid = ion_grids[e]
                    axes = [
                        axis(p1='r1',scale=scale,label='r'),
                        axis(p1='r2',scale=scale,label='phi'),
                        axis(p1='r3',scale=scale,label='theta'),
                        ]
                    n=0
                    for ax in axes:
                        ax.grid = '0 ({0}) 1'.format(grid[n])
                        n+=1
                    #end for
                    sg.append(spacegrid(coord=coord,origin=origin(p1=static+str(i)),axes=axes))
                    if e in missing:
                        missing.remove(e)
                    #end if
                #end if
                i+=1
            #end for
            if len(missing)>0:
                QmcpackInput.class_error('ion species not found for spherical grid\nspecies not found: {0}\nspecies present: {1}'.format(sorted(missing),sorted(set(list(system.structure.elem)))),'generate_energydensity')
            #end if
        #end if
    else:
        QmcpackInput.class_error('unsupported coord type\ncoord type provided: {0}\nsupported coord types: voronoi, cartesian, spherical'.format(coord),'generate_energydensity')
    #end if
    ed = energydensity(
        type       = 'EnergyDensity',
        name       = name,
        dynamic    = dynamic,
        static     = static,
        spacegrids = sg,
        )
    if refp is not None:
        ed.reference_points = refp
    #end if
    return ed
#end def generate_energydensity


opt_map = dict(linear=linear,cslinear=cslinear,linear_batch=linear_batch)
def generate_opt(method,
                 repeat           = 1,
                 energy           = None,
                 rw_variance      = None,
                 urw_variance     = None,
                 params           = None,
                 jastrows         = None,
                 processes        = None,
                 walkers_per_proc = None,
                 threads          = None,
                 blocks           = 2000,
                 #steps            = 5,
                 decorr           = 10,
                 min_walkers      = None, #use e.g. 128 for gpu's
                 timestep         = .5,
                 nonlocalpp       = False,
                 sample_factor    = 1.0):
    if not method in opt_map:
        QmcpackInput.class_error('section cannot be generated for optimization method '+method)
    #end if
    if energy is None and rw_variance is None and urw_variance is None:
        QmcpackInput.class_error('at least one cost parameter must be specified\n options are: energy, rw_variance, urw_variance')
    #end if
    if params is None and jastrows is None:
        QmcpackInput.class_error('must provide either number of opt parameters (params) or a list of jastrow objects (jastrows)')
    #end if
    if processes is None:
        QmcpackInput.class_error('must specify total number of processes')
    elif walkers_per_proc is None and threads is None:
        QmcpackInput.class_error('must specify walkers_per_proc or threads')
    #end if

    if params is None:
        params = count_jastrow_params(jastrows)
    #end if
    samples = max(100000,100*params**2)
    samples = int(round(sample_factor*samples))
    samples_per_proc = int(round(float(samples)/processes))

    if walkers_per_proc is None:
        walkers = 1
        walkers_per_proc = threads
    else:
        walkers = walkers_per_proc
    #end if
    tot_walkers = processes*walkers_per_proc
    if min_walkers!=None:
        tot_walkers = max(min_walkers,tot_walkers)
        walkers = int(ceil(float(tot_walkers)/processes-.001))
        if threads!=None and mod(walkers,threads)!=0:
            walkers = threads*int(ceil(float(walkers)/threads-.001))
        #end if
    #end if
    #blocks = int(ceil(float(decorr*samples)/(steps*tot_walkers)))
    blocks = min(blocks,samples_per_proc*decorr)

    opt = opt_map[method]()
 
    opt.set(
        walkers    = walkers,
        blocks     = blocks,
        #steps      = steps,
        samples    = samples,
        substeps   = decorr,
        timestep   = timestep,
        nonlocalpp = nonlocalpp,
        stepsbetweensamples = 1
        )
    if energy!=None:
        opt.energy = energy
    #end if
    if rw_variance!=None:
        opt.reweightedvariance = rw_variance
    #end if
    if urw_variance!=None:
        opt.unreweightedvariance = urw_variance
    #end if
    
    opt.incorporate_defaults(elements=True)

    if repeat>1:
        opt = loop(max=repeat,qmc=opt)
    #end if

    return opt
#end def generate_opt


def generate_opts(opt_reqs,**kwargs):
    opts = []
    for opt_req in opt_reqs:
        opts.append(generate_opt(*opt_req,**kwargs))
    #end for
    return opts
#end def generate_opts




# legacy driver defaults

opt_legacy_defaults = obj(
    method          = 'linear',
    minmethod       = 'quartic',
    cost            = 'variance',
    cycles          = 12,
    var_cycles      = 0,
    var_samples     = None,
    init_cycles     = 0,
    init_samples    = None,
    init_minwalkers = 1e-4,
    )

shared_opt_legacy_defaults = obj(
    samples              = 204800,
    nonlocalpp           = True,
    use_nonlocalpp_deriv = True,
    warmupsteps          = 300,                
    blocks               = 100,                
    steps                = 1,                  
    substeps             = 10,                 
    timestep             = 0.3,
    usedrift             = False,  
    max_seconds          = None,
    )

linear_quartic_legacy_defaults = obj(
    minwalkers        = 0.3,
    usebuffer         = True,
    exp0              = -6,
    bigchange         = 10.0,
    alloweddifference = 1e-04,
    stepsize          = 0.15,
    nstabilizers      = 1,
    **shared_opt_legacy_defaults
    )
linear_oneshift_legacy_defaults = obj(
    minwalkers = 0.5,
    #shift_i    = 0.01,
    #shift_s    = 1.00,
    **shared_opt_legacy_defaults
    )
linear_adaptive_legacy_defaults = obj(
    minwalkers          = 0.3,
    max_relative_change = 10.0,
    max_param_change    = 0.3,
    shift_i             = 0.01,
    shift_s             = 1.00,
    **shared_opt_legacy_defaults
    )

opt_method_legacy_defaults = obj({
    ('linear'  ,'quartic' ) : linear_quartic_legacy_defaults,
    ('linear'  ,'rescale' ) : linear_quartic_legacy_defaults,
    ('linear'  ,'linemin' ) : linear_quartic_legacy_defaults,
    ('cslinear','quartic' ) : linear_quartic_legacy_defaults,
    ('cslinear','rescale' ) : linear_quartic_legacy_defaults,
    ('cslinear','linemin' ) : linear_quartic_legacy_defaults,
    ('linear'  ,'adaptive') : linear_adaptive_legacy_defaults,
    ('linear'  ,'oneshift') : linear_oneshift_legacy_defaults,
    ('linear'  ,'oneshiftonly') : linear_oneshift_legacy_defaults,
    })
del shared_opt_legacy_defaults
del linear_quartic_legacy_defaults
del linear_oneshift_legacy_defaults
del linear_adaptive_legacy_defaults

allowed_opt_method_legacy_inputs = set(linear.attributes+linear.parameters
                                       +cslinear.attributes+cslinear.parameters)

vmc_legacy_defaults = obj(
    walkers     = 1,
    warmupsteps = 50,
    blocks      = 800,
    steps       = 10,
    substeps    = 3,
    timestep    = 0.3,
    checkpoint  = -1,
    usedrift    = None,
    max_seconds = None,
    )
vmc_test_legacy_defaults = obj(
    warmupsteps = 10,
    blocks      = 20,
    steps       =  4,
    ).set_optional(**vmc_legacy_defaults)
vmc_noJ_legacy_defaults = obj(
    warmupsteps = 200,
    blocks      = 800,
    steps       = 100,
    ).set_optional(**vmc_legacy_defaults)

dmc_legacy_defaults = obj(
    warmupsteps             = 20,
    blocks                  = 200,
    steps                   = 10,
    timestep                = 0.01,
    checkpoint              = -1,
    vmc_samples             = 2048,
    vmc_samplesperthread    = None, 
    vmc_walkers             = None,
    vmc_warmupsteps         = 30,
    vmc_blocks              = 40,
    vmc_steps               = 10,
    vmc_substeps            = 3,
    vmc_timestep            = 0.3,
    vmc_usedrift            = None,
    vmc_checkpoint          = -1,
    eq_dmc                  = False,
    eq_warmupsteps          = 20,
    eq_blocks               = 20,
    eq_steps                = 5,
    eq_timestep             = 0.02,
    eq_checkpoint           = -1,
    ntimesteps              = 1,
    timestep_factor         = 0.5,    
    nonlocalmoves           = None,
    branching_cutoff_scheme = None,
    maxage                  = None,
    feedback                = None,
    sigmabound              = None,
    max_seconds             = None,
    )
dmc_test_legacy_defaults = obj(
    vmc_warmupsteps = 10,
    vmc_blocks      = 20,
    vmc_steps       =  4,
    eq_warmupsteps  =  2,
    eq_blocks       =  5,
    eq_steps        =  2,
    warmupsteps     =  2,
    blocks          = 10,
    steps           =  2,
    ).set_optional(**dmc_legacy_defaults)
dmc_noJ_legacy_defaults = obj(
    warmupsteps     =  40,
    blocks          = 400,
    steps           =  20,
    ).set_optional(**dmc_legacy_defaults)


# batched driver defaults

opt_batched_defaults = obj(
    method          = 'linear',
    minmethod       = 'quartic',
    cost            = 'variance',
    cycles          = 12,
    var_cycles      = 0,
    var_samples     = None,
    init_cycles     = 0,
    init_samples    = None,
    init_minwalkers = 1e-4,
    )

shared_opt_batched_defaults = obj(
    samples              = 204800,
    nonlocalpp           = True,
    use_nonlocalpp_deriv = True,
    warmupsteps          = 300,                
    blocks               = 100,                
    steps                = 1,                  
    substeps             = 10,                 
    timestep             = 0.3,
    usedrift             = False,  
    )

linear_quartic_batched_defaults = obj(
    minwalkers        = 0.3,
    usebuffer         = True,
    exp0              = -6,
    bigchange         = 10.0,
    alloweddifference = 1e-04,
    stepsize          = 0.15,
    nstabilizers      = 1,
    **shared_opt_batched_defaults
    )
linear_oneshift_batched_defaults = obj(
    minwalkers = 0.5,
    #shift_i    = 0.01,
    #shift_s    = 1.00,
    **shared_opt_batched_defaults
    )
linear_adaptive_batched_defaults = obj(
    minwalkers          = 0.3,
    max_relative_change = 10.0,
    max_param_change    = 0.3,
    shift_i             = 0.01,
    shift_s             = 1.00,
    **shared_opt_batched_defaults
    )

opt_method_batched_defaults = obj({
    ('linear'  ,'quartic' ) : linear_quartic_batched_defaults,
    ('linear'  ,'rescale' ) : linear_quartic_batched_defaults,
    ('linear'  ,'linemin' ) : linear_quartic_batched_defaults,
    ('cslinear','quartic' ) : linear_quartic_batched_defaults,
    ('cslinear','rescale' ) : linear_quartic_batched_defaults,
    ('cslinear','linemin' ) : linear_quartic_batched_defaults,
    ('linear'  ,'adaptive') : linear_adaptive_batched_defaults,
    ('linear'  ,'oneshift') : linear_oneshift_batched_defaults,
    ('linear'  ,'oneshiftonly') : linear_oneshift_batched_defaults,
    })
del shared_opt_batched_defaults
del linear_quartic_batched_defaults
del linear_oneshift_batched_defaults
del linear_adaptive_batched_defaults

allowed_opt_method_batched_inputs = set(linear.attributes+linear.parameters
                                       +cslinear.attributes+cslinear.parameters)

vmc_batched_defaults = obj(
    total_walkers    = None,
    walkers_per_rank = None,
    warmupsteps      = 50,
    blocks           = 800,
    steps            = 10,
    substeps         = 3,
    timestep         = 0.3,
    usedrift         = False,
    checkpoint       = None,
    maxcpusecs       = None,
    crowds           = None,
    )
vmc_test_batched_defaults = obj(
    warmupsteps = 10,
    blocks      = 20,
    steps       =  4,
    ).set_optional(**vmc_batched_defaults)
vmc_noJ_batched_defaults = obj(
    warmupsteps = 200,
    blocks      = 800,
    steps       = 100,
    ).set_optional(**vmc_batched_defaults)

dmc_batched_defaults = obj(
    total_walkers           = None,
    walkers_per_rank        = None,
    warmupsteps             = 20,
    blocks                  = 200,
    steps                   = 10,
    substeps                = None,
    timestep                = 0.01,
    checkpoint              = None,
    vmc_warmupsteps         = 30,
    vmc_blocks              = 40,
    vmc_steps               = 10,
    vmc_substeps            = 3,
    vmc_timestep            = 0.3,
    vmc_usedrift            = False,
    vmc_checkpoint          = None,
    eq_dmc                  = False,
    eq_warmupsteps          = 20,
    eq_blocks               = 20,
    eq_steps                = 5,
    eq_timestep             = 0.02,
    eq_checkpoint           = None,
    ntimesteps              = 1,
    timestep_factor         = 0.5,    
    nonlocalmoves           = None,
    branching_cutoff_scheme = None,
    crowd_serialize_walkers = None,
    crowds                  = None,
    reconfiguration         = None,
    maxage                  = None,
    feedback                = None,
    sigmabound              = None,
    )
dmc_test_batched_defaults = obj(
    vmc_warmupsteps = 10,
    vmc_blocks      = 20,
    vmc_steps       =  4,
    eq_warmupsteps  =  2,
    eq_blocks       =  5,
    eq_steps        =  2,
    warmupsteps     =  2,
    blocks          = 10,
    steps           =  2,
    ).set_optional(**dmc_batched_defaults)
dmc_noJ_batched_defaults = obj(
    warmupsteps     =  40,
    blocks          = 400,
    steps           =  20,
    ).set_optional(**dmc_batched_defaults)



# collected defaults for opt methods

opt_method_defaults = obj(
    legacy  = opt_method_legacy_defaults,
    batched = opt_method_batched_defaults,
    )


# collected defaults for all drivers

qmc_defaults = obj(
    legacy = obj(
        opt      = opt_legacy_defaults,
        vmc      = vmc_legacy_defaults,
        vmc_test = vmc_test_legacy_defaults,
        vmc_noJ  = vmc_noJ_legacy_defaults,
        dmc      = dmc_legacy_defaults,
        dmc_test = dmc_test_legacy_defaults,
        dmc_noJ  = dmc_noJ_legacy_defaults,
        ),
    batched = obj(
        opt      = opt_batched_defaults,
        vmc      = vmc_batched_defaults,
        vmc_test = vmc_test_batched_defaults,
        vmc_noJ  = vmc_noJ_batched_defaults,
        dmc      = dmc_batched_defaults,
        dmc_test = dmc_test_batched_defaults,
        dmc_noJ  = dmc_noJ_batched_defaults,
        ),
    )

del opt_legacy_defaults
del vmc_legacy_defaults
del vmc_test_legacy_defaults
del vmc_noJ_legacy_defaults
del dmc_legacy_defaults
del dmc_test_legacy_defaults
del dmc_noJ_legacy_defaults

del opt_batched_defaults
del vmc_batched_defaults
del vmc_test_batched_defaults
del vmc_noJ_batched_defaults
del dmc_batched_defaults
del dmc_test_batched_defaults
del dmc_noJ_batched_defaults



def generate_opt_calculations(driver,**kwargs):
    if driver=='legacy':
        calcs = generate_legacy_opt_calculations(**kwargs)
    elif driver=='batched':
        calcs = generate_batched_opt_calculations(**kwargs)
    else:
        error('Cannot generate calculations for unrecognized driver.\nUnrecognized driver: {}'.format(driver))
    #end if
    return calcs
#end def generate_opt_calculations


def generate_vmc_calculations(driver,**kwargs):
    if driver=='legacy':
        calcs = generate_legacy_vmc_calculations(**kwargs)
    elif driver=='batched':
        calcs = generate_batched_vmc_calculations(**kwargs)
    else:
        error('Cannot generate calculations for unrecognized driver.\nUnrecognized driver: {}'.format(driver))
    #end if
    return calcs
#end def generate_vmc_calculations


def generate_dmc_calculations(driver,**kwargs):
    if driver=='legacy':
        calcs = generate_legacy_dmc_calculations(**kwargs)
    elif driver=='batched':
        calcs = generate_batched_dmc_calculations(**kwargs)
    else:
        error('Cannot generate calculations for unrecognized driver.\nUnrecognized driver: {}'.format(driver))
    #end if
    return calcs
#end def generate_dmc_calculations



def generate_legacy_opt_calculations(
    method     ,
    cost       ,
    cycles     ,
    var_cycles ,
    var_samples,
    init_cycles,
    init_samples,
    init_minwalkers,
    loc        = 'generate_opt_calculations',
    **opt_inputs
    ):

    methods = obj(linear=linear,cslinear=cslinear)
    if method not in methods:
        error('invalid optimization method requested\ninvalid method: {0}\nvalid options are: {1}'.format(method,sorted(methods.keys())),loc)
    #end if
    opt = methods[method]

    opt_inputs = obj(opt_inputs)
    invalid = set(opt_inputs.keys())-allowed_opt_method_legacy_inputs
    oneshift = False
    if len(invalid)>0:
        error('invalid optimization inputs provided\ninvalid inputs: {}\nvalid options are: {}'.format(sorted(invalid),sorted(allowed_opt_method_legacy_inputs)))
    #end if
    for k in list(opt_inputs.keys()):
        if opt_inputs[k] is None:
            del opt_inputs[k]
        #end if
    #end for
    if 'minmethod' in opt_inputs and opt_inputs.minmethod.lower().startswith('oneshift'):
        opt_inputs.minmethod = 'OneShiftOnly'
        oneshift = True
    #end if

    if cost=='variance':
        cost = (0.0,1.0,0.0)
    elif cost=='energy':
        cost = (1.0,0.0,0.0)
    elif isinstance(cost,(tuple,list)) and (len(cost)==2 or len(cost)==3):
        if len(cost)==2:
            cost = (cost[0],0.0,cost[1])
        #end if
    else:
        error('invalid optimization cost function encountered\ninvalid cost fuction: {0}\nvalid options are: variance, energy, (0.95,0.05), etc'.format(cost),loc)
    #end if
    opt_calcs = []
    if var_cycles>0:
        vmin_opt = opt(
            energy               = 0.0,
            unreweightedvariance = 1.0,
            reweightedvariance   = 0.0,
            **opt_inputs
            )
        if var_samples is not None:
            vmin_opt.samples = var_samples
        #end if
        opt_calcs.append(loop(max=var_cycles,qmc=vmin_opt))
    #end if
    if init_cycles>0:
        init_opt = opt(**opt_inputs)
        if init_samples is not None:
            init_opt.samples = init_samples
        #end if
        init_opt.minwalkers = init_minwalkers
        if not oneshift:
            init_opt.energy               = cost[0]
            init_opt.unreweightedvariance = cost[1]
            init_opt.reweightedvariance   = cost[2]
        #end if
        opt_calcs.append(loop(max=init_cycles,qmc=init_opt))
    #end if

    cost_opt = opt(**opt_inputs)
    if not oneshift:
        cost_opt.energy               = cost[0]
        cost_opt.unreweightedvariance = cost[1]
        cost_opt.reweightedvariance   = cost[2]
    #end if

    opt_calcs.append(loop(max=cycles,qmc=cost_opt))
    return opt_calcs
#end def generate_legacy_opt_calculations



def generate_legacy_vmc_calculations(
    walkers    ,
    warmupsteps,
    blocks     ,
    steps      ,
    substeps   ,
    timestep   ,
    checkpoint ,
    usedrift   ,
    max_seconds,
    loc        = 'generate_vmc_calculations',
    ):

    vmc_calc = vmc(
        walkers     = walkers,
        warmupsteps = warmupsteps,
        blocks      = blocks,
        steps       = steps,
        substeps    = substeps,
        timestep    = timestep,
        checkpoint  = checkpoint,
        )

    if usedrift is not None:
        vmc_calc.usedrift = usedrift
    #end if
    if max_seconds is not None:
        vmc_calc.max_seconds = max_seconds
    #end if

    vmc_calcs = [vmc_calc]

    return vmc_calcs
#end def generate_legacy_vmc_calculations



def generate_legacy_dmc_calculations(
    warmupsteps            ,
    blocks                 ,
    steps                  ,
    timestep               ,
    checkpoint             ,
    vmc_samples            ,
    vmc_samplesperthread   , 
    vmc_walkers            ,
    vmc_warmupsteps        ,
    vmc_blocks             ,
    vmc_steps              ,
    vmc_substeps           ,
    vmc_timestep           ,
    vmc_usedrift           ,
    vmc_checkpoint         ,
    eq_dmc                 ,
    eq_warmupsteps         ,
    eq_blocks              ,
    eq_steps               ,
    eq_timestep            ,
    eq_checkpoint          ,
    ntimesteps             ,
    timestep_factor        ,    
    nonlocalmoves          ,
    branching_cutoff_scheme,
    maxage                 ,
    feedback               ,
    sigmabound             ,
    max_seconds            ,
    loc                 = 'generate_dmc_calculations',
    ):

    if vmc_samples is None and vmc_samplesperthread is None and vmc_walkers is None:
        error('vmc samples (dmc walkers) not specified\nplease provide one of the following keywords: vmc_samples, vmc_samplesperthread, vmc_walkers',loc)
    #end if
    if vmc_walkers is None:
        vmc_walkers = 1
    #end if

    vmc_calc = vmc(
        walkers     = vmc_walkers,
        warmupsteps = vmc_warmupsteps,
        blocks      = vmc_blocks,
        steps       = vmc_steps,
        substeps    = vmc_substeps,
        timestep    = vmc_timestep,
        checkpoint  = vmc_checkpoint,
        )
    if vmc_samplesperthread is not None:
        vmc_calc.samplesperthread = vmc_samplesperthread
    elif vmc_samples is not None:
        vmc_calc.samples = vmc_samples
    #end if
    if vmc_usedrift is not None:
        vmc_calc.usedrift = vmc_usedrift
    #end if
    if max_seconds is not None:
        vmc_calc.max_seconds = max_seconds
    #end if

    dmc_calcs = [vmc_calc]
    if eq_dmc:
        dmc_calcs.append(
            dmc(
                warmupsteps   = eq_warmupsteps,
                blocks        = eq_blocks,
                steps         = eq_steps,
                timestep      = eq_timestep,
                checkpoint    = eq_checkpoint,
                )
            )
    #end if
    tfac = 1.0
    for n in range(ntimesteps):
        sfac = 1.0/tfac
        dmc_calcs.append(
            dmc(
                warmupsteps   = int(sfac*warmupsteps),
                blocks        = blocks,
                steps         = int(sfac*steps),
                timestep      = tfac*timestep,
                checkpoint    = checkpoint,
                )
            )
        tfac *= timestep_factor
    #end for

    optional_dmc_inputs = obj(
        nonlocalmoves           = nonlocalmoves,
        branching_cutoff_scheme = branching_cutoff_scheme,
        maxage                  = maxage    ,
        feedback                = feedback  ,
        sigmabound              = sigmabound,
        max_seconds             = max_seconds,
        )
    for calc in dmc_calcs:
        if isinstance(calc,dmc):
            for k,v in optional_dmc_inputs.items():
                if v is not None:
                    calc[k] = v
                #end if
            #end for
        #end if
    #end for
    
    return dmc_calcs
#end def generate_legacy_dmc_calculations



def generate_batched_opt_calculations(
    method     ,
    cost       ,
    cycles     ,
    var_cycles ,
    var_samples,
    init_cycles,
    init_samples,
    init_minwalkers,
    loc        = 'generate_opt_calculations',
    **opt_inputs
    ):

    methods = obj(linear=linear)
    if method not in methods:
        error('invalid optimization method requested\ninvalid method: {0}\nvalid options are: {1}'.format(method,sorted(methods.keys())),loc)
    #end if
    opt = methods[method]

    opt_inputs = obj(opt_inputs)
    invalid = set(opt_inputs.keys())-allowed_opt_method_batched_inputs
    oneshift = False
    if len(invalid)>0:
        error('invalid optimization inputs provided\ninvalid inputs: {}\nvalid options are: {}'.format(sorted(invalid),sorted(allowed_opt_method_batched_inputs)))
    #end if
    if 'minmethod' in opt_inputs and opt_inputs.minmethod.lower().startswith('oneshift'):
        opt_inputs.minmethod = 'OneShiftOnly'
        oneshift = True
    #end if

    if cost=='variance':
        cost = (0.0,1.0,0.0)
    elif cost=='energy':
        cost = (1.0,0.0,0.0)
    elif isinstance(cost,(tuple,list)) and (len(cost)==2 or len(cost)==3):
        if len(cost)==2:
            cost = (cost[0],0.0,cost[1])
        #end if
    else:
        error('invalid optimization cost function encountered\ninvalid cost fuction: {0}\nvalid options are: variance, energy, (0.95,0.05), etc'.format(cost),loc)
    #end if
    opt_calcs = []
    if var_cycles>0:
        vmin_opt = opt(
            energy               = 0.0,
            unreweightedvariance = 1.0,
            reweightedvariance   = 0.0,
            **opt_inputs
            )
        if var_samples is not None:
            vmin_opt.samples = var_samples
        #end if
        opt_calcs.append(loop(max=var_cycles,qmc=vmin_opt))
    #end if
    if init_cycles>0:
        init_opt = opt(**opt_inputs)
        if init_samples is not None:
            init_opt.samples = init_samples
        #end if
        init_opt.minwalkers = init_minwalkers
        if not oneshift:
            init_opt.energy               = cost[0]
            init_opt.unreweightedvariance = cost[1]
            init_opt.reweightedvariance   = cost[2]
        #end if
        opt_calcs.append(loop(max=init_cycles,qmc=init_opt))
    #end if

    cost_opt = opt(**opt_inputs)
    if not oneshift:
        cost_opt.energy               = cost[0]
        cost_opt.unreweightedvariance = cost[1]
        cost_opt.reweightedvariance   = cost[2]
    #end if

    opt_calcs.append(loop(max=cycles,qmc=cost_opt))
    return opt_calcs
#end def generate_batched_opt_calculations



def generate_batched_vmc_calculations(
    total_walkers    ,
    walkers_per_rank ,     
    warmupsteps      ,
    blocks           ,
    steps            ,
    substeps         ,
    timestep         ,
    usedrift         ,
    checkpoint       ,
    maxcpusecs       ,
    crowds           ,
    loc              = 'generate_vmc_calculations',
    ):
    
    if total_walkers is not None and walkers_per_rank is not None:
        error('Only one of "total_walkers" and "walkers_per_rank" may be provided.',loc)
    #end if

    vmc_inputs = obj(
        warmupsteps = warmupsteps,
        blocks      = blocks,
        steps       = steps,
        substeps    = substeps,
        timestep    = timestep,
        usedrift    = usedrift,
        )
    optional_vmc_inputs = obj(
        total_walkers = total_walkers,
        walkers_per_rank = walkers_per_rank,
        #checkpoint       = checkpoint, # no checkpointing support yet
        maxcpusecs       = maxcpusecs,
        crowds           = crowds,
        )
    for name,value in optional_vmc_inputs.items():
        if value is not None:
            vmc_inputs[name] = value
        #end if
    #end for

    vmc_calcs = [vmc(**vmc_inputs)]

    return vmc_calcs
#end def generate_batched_vmc_calculations



def generate_batched_dmc_calculations(
    total_walkers          ,
    walkers_per_rank       ,
    warmupsteps            ,
    blocks                 ,
    steps                  ,
    substeps               ,
    timestep               ,
    checkpoint             ,
    vmc_warmupsteps        ,
    vmc_blocks             ,
    vmc_steps              ,
    vmc_substeps           ,
    vmc_timestep           ,
    vmc_usedrift           ,
    vmc_checkpoint         ,
    eq_dmc                 ,
    eq_warmupsteps         ,
    eq_blocks              ,
    eq_steps               ,
    eq_timestep            ,
    eq_checkpoint          ,
    ntimesteps             ,
    timestep_factor        ,    
    nonlocalmoves          ,
    branching_cutoff_scheme,
    crowd_serialize_walkers,
    crowds                 ,
    reconfiguration        ,
    maxage                 ,
    feedback               ,
    sigmabound             ,
    loc                 = 'generate_dmc_calculations',
    ):

    if total_walkers is None and walkers_per_rank is None:
        error('DMC walker count not specified via "total_walkers" or "walkers_per_rank".\nPlease provide at least one of these.\n\nWarning: use care in the selection of these parameters.\nPerformance critically depends on the walker count and the batched QMCPACK \ndrivers make no effort to prevent substantial under-utilization.',loc)
    elif total_walkers is not None and walkers_per_rank is not None:
        error('Only one of "total_walkers" and "walkers_per_rank" may be provided.',loc)
    #end if

    vmc_inputs = obj(
        warmupsteps = vmc_warmupsteps,
        blocks      = vmc_blocks,
        steps       = vmc_steps,
        substeps    = vmc_substeps,
        timestep    = vmc_timestep,
        usedrift    = vmc_usedrift,
        )
    optional_vmc_inputs = obj(
        total_walkers    = total_walkers,
        walkers_per_rank = walkers_per_rank,
        crowds           = crowds,
        #checkpoint       = vmc_checkpoint, # not supported yet
        )
    for name,value in optional_vmc_inputs.items():
        if value is not None:
            vmc_inputs[name] = value
        #end if
    #end for

    vmc_calc = vmc(**vmc_inputs)

    dmc_calcs = [vmc_calc]
    if eq_dmc:
        dmc_calcs.append(
            dmc(
                warmupsteps   = eq_warmupsteps,
                blocks        = eq_blocks,
                steps         = eq_steps,
                timestep      = eq_timestep,
                #checkpoint    = eq_checkpoint, # not supported yet
                )
            )
    #end if
    tfac = 1.0
    for n in range(ntimesteps):
        sfac = 1.0/tfac
        dmc_calcs.append(
            dmc(
                warmupsteps   = int(sfac*warmupsteps),
                blocks        = blocks,
                steps         = int(sfac*steps),
                timestep      = tfac*timestep,
                #checkpoint    = checkpoint,  # not supported yet
                )
            )
        tfac *= timestep_factor
    #end for

    optional_dmc_inputs = obj(
        total_walkers           = total_walkers,
        walkers_per_rank        = walkers_per_rank,
        substeps                = substeps,
        nonlocalmoves           = nonlocalmoves,
        branching_cutoff_scheme = branching_cutoff_scheme,
        crowd_serialize_walkers = crowd_serialize_walkers,
        crowds                  = crowds,
        reconfiguration         = reconfiguration,
        maxage                  = maxage,
        feedback                = feedback,
        sigmabound              = sigmabound,
        )
    for calc in dmc_calcs:
        if isinstance(calc,dmc):
            for name,value in optional_dmc_inputs.items():
                if value is not None:
                    calc[name] = value
                #end if
            #end for
        #end if
    #end for
    
    return dmc_calcs
#end def generate_batched_dmc_calculations



def generate_qmcpack_input(**kwargs):
    QIcollections.clear()
    system = kwargs.get('system',None)
    if isinstance(system,PhysicalSystem):
        system.update_particles()
    #end if
    selector = kwargs.pop('input_type','basic')
    if selector=='basic':
        inp = generate_basic_input(**kwargs)
    elif selector=='basic_afqmc':
        inp = generate_basic_afqmc_input(**kwargs)
    elif selector=='opt_jastrow':
        inp = generate_opt_jastrow_input(**kwargs)
    else:
        QmcpackInput.class_error('selection '+str(selector)+' has not been implemented for qmcpack input generation')
    #end if
    return inp
#end def generate_qmcpack_input



def read_jastrows(filepath):
    qi = QmcpackInput(filepath)
    qi.pluralize()
    jastrows = qi.get('jastrows')
    return jastrows
#end def read_jastrows



gen_basic_input_defaults = obj(
    id             = 'qmc',            
    series         = 0,                
    purpose        = '',     
    maxcpusecs     = None,
    max_seconds    = None,
    seed           = None,             
    bconds         = None,             
    truncate       = False,            
    buffer         = None,             
    lr_dim_cutoff  = 15,               
    lr_tol         = None,               
    lr_handler     = None,               
    remove_cell    = False,            
    randomsrc      = True,            
    meshfactor     = 1.0,              
    orbspline      = None,             
    precision      = 'float',          
    twistnum       = None,             
    twist          = None,             
    spin_polarized = None,             
    partition      = None,             
    partition_mf   = None,             
    hybridrep      = None,             
    hybrid_rcut    = None,             
    hybrid_lmax    = None,             
    orbitals_h5    = 'MISSING.h5',     
    run_path       = None,
    check_paths    = True,
    excitation     = None,             
    system         = 'missing',        
    pseudos        = None,
    dla            = None,
    delay_rank     = None,
    jastrows       = 'generateJ12',    
    interactions   = 'all',            
    corrections    = 'default',        
    observables    = None,             
    estimators     = None,             
    traces         = None,             
    calculations   = None,             
    det_format     = 'new',            
    J1             = False,            
    J2             = False,            
    J3             = False,            
    J1_size        = None,             
    J1_rcut        = None,             
    J1_dr          = 0.5,              
    J2_size        = None,             
    J2_rcut        = None,             
    J2_dr          = 0.5,              
    J2_init        = 'zero',           
    J3_isize       = 3,                
    J3_esize       = 3,                
    J3_rcut        = 5.0,              
    J1_rcut_open   = 5.0,              
    J2_rcut_open   = 10.0,
    driver         = 'legacy', # legacy,batched
    # batched driver inputs
    orbitals_cpu   = None,     # place/evaluate orbitals on cpu if on gpu
    matrix_inv_cpu = None,     # evaluate matrix inverse on cpu if on gpu
    # legacy cuda inputs
    gpusharing     = None,
    qmc            = None,     # opt,vmc,vmc_test,dmc,dmc_test
    )

def generate_basic_input(**kwargs):
    # capture inputs
    kw = obj(kwargs)
    # apply general defaults
    kw.set_optional(**gen_basic_input_defaults)
    valid = set(gen_basic_input_defaults.keys())
    # apply method specific defaults
    if kw.qmc is not None:
        if kw.driver not in qmc_defaults:
            QmcpackInput.class_error('Invalid input for argument "driver".\nInvalid input: {}\nValid options are: {}'.format(kw.driver,sorted(qmc_defaults.keys())),'generate_qmcpack_input')
        #end if
        qmc_driver_defaults = qmc_defaults[kw.driver]
        if kw.qmc not in qmc_driver_defaults:
            QmcpackInput.class_error('Invalid input for argument "qmc".\nInvalid input: {}\nValid options are: {}'.format(kw.qmc,sorted(qmc_driver_defaults.keys())),'generate_qmcpack_input')
        #end if
        qmc_keys = ['driver']
        kw.set_optional(**qmc_driver_defaults[kw.qmc])
        qmc_keys += list(qmc_driver_defaults[kw.qmc].keys())
        if kw.qmc=='opt':
            opt_method_driver_defaults = opt_method_defaults[kw.driver]
            key = (kw.method,kw.minmethod.lower())
            if key not in opt_method_driver_defaults:
                QmcpackInput.class_error('invalid input for arguments "method,minmethod".\nInvalid input: {}\nValid options are: {}'.format(key,sorted(opt_method_driver_defaults.keys())),'generate_qmcpack_input')
            #end if
            kw.set_optional(**opt_method_driver_defaults[key])
            qmc_keys += list(opt_method_driver_defaults[key].keys())
            del key
        #end if
        valid |= set(qmc_keys)
    #end if
    # screen for invalid keywords
    invalid_kwargs = set(kw.keys())-valid
    if len(invalid_kwargs)>0:
        QmcpackInput.class_error('invalid input parameters encountered.\nInvalid input parameters: {0}\nValid options are: {1}'.format(sorted(invalid_kwargs),sorted(valid)),'generate_qmcpack_input')
    #end if

    batched = kw.driver=='batched'
    legacy  = kw.driver=='legacy'

    if kw.system=='missing':
        QmcpackInput.class_error('argument "system" is missing.\nIf you really do not want particlesets to be generated, set system to None.','generate_qmcpack_input')
    #end if
    if kw.bconds is None:
        if kw.system is not None:
            s = kw.system.structure
            kw.bconds = s.bconds
            if len(kw.bconds)==0 or not s.has_axes():
                kw.bconds = 'nnn'
            #end if
        else:
            kw.bconds = 'ppp'
        #end if
    #end if
    if kw.corrections=='default' and tuple(kw.bconds)==tuple('ppp'):
        if not batched:
            kw.corrections = ['mpc','chiesa']
        else:
            kw.corrections = ['mpc']
        #end if
    elif isinstance(kw.corrections,(list,tuple)):
        None
    else:
        kw.corrections = []
    #end if
    if kw.observables is None:
        #observables = ['localenergy']
        kw.observables = []
    #end if
    if kw.estimators is None:
        kw.estimators = []
    #end if
    kw.estimators = kw.estimators + kw.observables + kw.corrections
    if kw.calculations is None:
        kw.calculations = []
    #end if
    if kw.spin_polarized is None:
        kw.spin_polarized = kw.system.net_spin>0
    #end if
    if kw.partition is not None:
        kw.det_format = 'new'
    #end if
    if kw.hybrid_rcut is not None or kw.hybrid_lmax is not None:
        kw.hybridrep = True
    #end if

    metadata = QmcpackInput.default_metadata.copy()

    proj = project(
        id             = kw.id,
        series         = kw.series,
        application    = application(),
        driver_version = kw.driver,
        )
    if batched:
        if kw.maxcpusecs is not None:
            proj.maxcpusecs = kw.maxcpusecs
        #end if
        if kw.max_seconds is not None:
            proj.max_seconds = kw.max_seconds
        #end if
    #end if

    simcell = generate_simulationcell(
        bconds        = kw.bconds,
        lr_dim_cutoff = kw.lr_dim_cutoff,
        lr_tol        = kw.lr_tol,
        lr_handler    = kw.lr_handler,
        system        = kw.system,
        )

    if kw.system is not None:
        kw.system.structure.set_bconds(kw.bconds)
        particlesets = generate_particlesets(
            system      = kw.system,
            randomsrc   = kw.randomsrc,
            hybrid_rcut = kw.hybrid_rcut,
            hybrid_lmax = kw.hybrid_lmax,
            )
    #end if


    if kw.det_format=='new':
        if kw.excitation is not None:
            QmcpackInput.class_error('user provided "excitation" input argument with new style determinant format.\nPlease add det_format="old" and try again','generate_qmcpack_input')
        #end if
        if kw.system is not None and isinstance(kw.system.structure,Jellium):
            ssb = generate_sposet_builder(
                type           = 'heg',
                twist          = kw.twist,
                spin_polarized = kw.spin_polarized,
                system         = kw.system,
                )
        else:
            if kw.orbspline is None:
                kw.orbspline = 'bspline'
            #end if
            if kw.orbitals_h5!='MISSING.h5':
                orbfile_exists = os.path.exists(kw.orbitals_h5)
                if kw.check_paths and not orbfile_exists:
                    QmcpackInput.class_error('user provided "orbitals_h5" path does not exist\nPath provided: {}\nTo disable this check, set check_paths=False'.format(kw.orbitals_h5),'generate_qmcpack_input')
                #end if
                if kw.run_path is not None:
                    kw.orbitals_h5 = os.path.relpath(kw.orbitals_h5,kw.run_path)
                #end if
            #end if
            ssb = generate_sposet_builder(
                type           = kw.orbspline,
                twist          = kw.twist,
                twistnum       = kw.twistnum,
                meshfactor     = kw.meshfactor,
                precision      = kw.precision,
                truncate       = kw.truncate,
                buffer         = kw.buffer,
                hybridrep      = kw.hybridrep,
                href           = kw.orbitals_h5,
                spin_polarized = kw.spin_polarized,
                system         = kw.system,
                orbitals_cpu   = kw.orbitals_cpu,
                gpusharing     = kw.gpusharing,
                )
        #end if
        if kw.partition is None:
            spobuilders = [ssb]
        else:
            spobuilders = partition_sposets(
                sposet_builder = ssb,
                partition      = kw.partition,
                partition_meshfactors = kw.partition_mf,
                )
        #end if

        dset = generate_determinantset(
            spin_polarized = kw.spin_polarized,
            delay_rank     = kw.delay_rank,
            matrix_inv_cpu = kw.matrix_inv_cpu,
            system         = kw.system,
            )
    elif kw.det_format=='old':
        spobuilders = None
        if kw.orbspline is None:
            kw.orbspline = 'einspline'
        #end if
        dset = generate_determinantset_old(
            type           = kw.orbspline,
            twistnum       = kw.twistnum,
            meshfactor     = kw.meshfactor,
            precision      = kw.precision,
            hybridrep      = kw.hybridrep,
            href           = kw.orbitals_h5,
            spin_polarized = kw.spin_polarized,
            excitation     = kw.excitation,
            delay_rank     = kw.delay_rank,
            gpusharing     = kw.gpusharing,
            system         = kw.system,
            )
    else:
        QmcpackInput.class_error('argument "det_format" is invalid.\nReceived: {0}\nValid options are: new, old'.format(det_format),'generate_qmcpack_input')
    #end if


    wfn = wavefunction(        
        name           = 'psi0',
        target         = 'e',
        determinantset = dset,
        )

    if isinstance(kw.jastrows,str) and kw.jastrows.endswith('.xml'):
        if not os.path.exists(kw.jastrows):
            QmcpackInput.class_error('user provided "jastrows" file path does not exist\nFile path provided: {}'.format(kw.jastrows),'generate_qmcpack_input')
        #end if
        jastrows = read_jastrows(kw.jastrows)
        if jastrows is None:
            QmcpackInput.class_error('no jastrows found at user provided "jastrows" file.\nFile path provided: {}'.format(kw.jastrows),'generate_qmcpack_input')
        #end if
        kw.jastrows = jastrows
    elif kw.J1 or kw.J2 or kw.J3:
        kw.jastrows = generate_jastrows_alt(
            J1           = kw.J1          ,
            J2           = kw.J2          ,
            J3           = kw.J3          ,
            J1_size      = kw.J1_size     ,
            J1_rcut      = kw.J1_rcut     ,
            J1_dr        = kw.J1_dr       ,
            J2_size      = kw.J2_size     ,
            J2_rcut      = kw.J2_rcut     ,
            J2_dr        = kw.J2_dr       ,
            J2_init      = kw.J2_init     ,
            J3_isize     = kw.J3_isize    ,
            J3_esize     = kw.J3_esize    ,
            J3_rcut      = kw.J3_rcut     ,
            J1_rcut_open = kw.J1_rcut_open,
            J2_rcut_open = kw.J2_rcut_open,
            system       = kw.system      ,
            )
    #end if
    if kw.jastrows is not None:
        wfn.jastrows = generate_jastrows(kw.jastrows,kw.system,check_ions=True)
    #end if

    h_estimators = kw.estimators
    d_estimators = None
    if batched:
        h_estimators = []
        d_estimators = []
        if isinstance(kw.estimators,list):
            for est in kw.estimators:
                if isinstance(est,str) and est.lower()=='mpc':
                    h_estimators.append(est)
                else:
                    d_estimators.append(est)
                #end if
            #end for
        #end if
    #end if

    hmltn = generate_hamiltonian(
        system       = kw.system,
        pseudos      = kw.pseudos,
        dla          = kw.dla,
        interactions = kw.interactions,
        estimators   = h_estimators,
        )

    if spobuilders is not None:
        wfn.sposet_builders = make_collection(spobuilders)
    #end if

    qmcsys = qmcsystem(
        simulationcell  = simcell,
        wavefunction    = wfn,
        hamiltonian     = hmltn,
        )

    if kw.system is not None:
        qmcsys.particlesets = particlesets
    #end if

    sim = simulation(
        project   = proj,
        qmcsystem = qmcsys,
        )

    if kw.seed is not None:
        sim.random = random(seed=kw.seed)
    #end if

    if kw.traces is not None:
        sim.traces = kw.traces
    #end if

    if len(kw.calculations)==0 and kw.qmc is not None:
        qmc_inputs = kw.obj(*qmc_keys)
        if kw.qmc=='opt':
            kw.calculations = generate_opt_calculations(**qmc_inputs)
        elif 'vmc' in kw.qmc:
            kw.calculations = generate_vmc_calculations(**qmc_inputs)
        elif 'dmc' in kw.qmc:
            kw.calculations = generate_dmc_calculations(**qmc_inputs)
        #end if
    #end if
    if batched and d_estimators is not None and len(d_estimators)>0:
        estimators = generate_estimators_batched(d_estimators)
        for calc in kw.calculations:
            if isinstance(calc,loop):
                calc = calc.qmc
            #end if
            calc.estimators = estimators.copy()
        #end for
    #end if
    sim.calculations = make_collection(kw.calculations).copy()

    qi = QmcpackInput(metadata,sim)

    qi.incorporate_defaults(elements=False,overwrite=False,propagate=True)

    if kw.remove_cell:
        qi.remove_physical_system()
    #end if

    for calc in sim.calculations:
        if isinstance(calc,loop):
            calc = calc.qmc
        #end if
        if isinstance(calc,(linear,cslinear,linear_batch)) and 'nonlocalpp' not in calc:
            calc.nonlocalpp           = True
            calc.use_nonlocalpp_deriv = True
        #end if
    #end for

    return qi
#end def generate_basic_input



gen_basic_afqmc_input_defaults = obj(
    id          = 'qmc',            
    series      = 0,   
    seed        = None,
    nmo         = None,
    naea        = None,
    naeb        = None,
    ham_file    = None,
    wfn_file    = None,
    wfn_type    = 'NOMSD',
    cutoff      = 1e-8,
    wset_type   = 'shared',
    walker_type = 'CLOSED',
    hybrid      = True,
    ncores      = 1,
    nwalkers    = 10,
    blocks      = 10000,
    steps       = 10,
    timestep    = 0.005,
    estimators  = None,
    info_name   = 'info0',
    ham_name    = 'ham0',
    wfn_name    = 'wfn0',
    wset_name   = 'wset0',
    prop_name   = 'prop0',
    system      = None,
    run_path    = None,
    )

def generate_basic_afqmc_input(**kwargs):
    # capture inputs
    kw = obj(kwargs)
    gen_info = obj()
    for k,v in kw.items():
        if not isinstance(v,obj):
            gen_info[k] = v
        #end if
    #end for
    # apply general defaults
    kw.set_optional(**gen_basic_afqmc_input_defaults)
    valid = set(gen_basic_afqmc_input_defaults.keys())
    # screen for invalid keywords
    invalid_kwargs = set(kw.keys())-valid
    if len(invalid_kwargs)>0:
        QmcpackInput.class_error('invalid input parameters encountered\ninvalid input parameters: {0}\nvalid options are: {1}'.format(sorted(invalid_kwargs),sorted(valid)),'generate_qmcpack_input')
    #end if

    metadata = meta(
        generation_info = gen_info.copy(),
        )

    sim = simulation(
        method = 'afqmc',
        )

    sim.project = project(
        id     = kw.id,
        series = kw.series,
        )

    if kw.seed is not None:
        sim.random = random(seed=kw.seed)
    #end if

    info = afqmcinfo(
        name = kw.info_name,
        )
    if kw.nmo is not None:
        info.nmo = kw.nmo
    #end if
    if kw.naea is not None:
        info.naea = kw.naea
    #end if
    if kw.naeb is not None:
        info.naeb = kw.naeb
    #end if
    sim.afqmcinfo = info

    if kw.ham_file is None and kw.wfn_file is not None:
        kw.ham_file = kw.wfn_file
    elif kw.ham_file is not None and kw.wfn_file is None:
        kw.wfn_file = kw.ham_file
    elif kw.ham_file is None and kw.wfn_file is None:
        kw.ham_file = 'MISSING.h5'
        kw.wfn_file = 'MISSING.h5'
    #end if
    def get_filetype(filename,loc):
        if filename.endswith('.h5'):
            filetype = 'hdf5'
        else:
            QmcpackInput.class_error('Type of {} file "{}" is unrecognized.\n The following file extensions are allowed: .h5'.format(loc,filename))
        #end if
        return filetype
    #end def get_filetype
    
    ham = hamiltonian(
        name     = kw.ham_name,
        info     = info.name,
        filetype = get_filetype(kw.ham_file,'hamiltonian'),
        filename = kw.ham_file,
        )
    sim.hamiltonian = ham

    wfn = wavefunction(
        name     = kw.wfn_name,
        info     = info.name,
        filetype = get_filetype(kw.wfn_file,'wavefunction'),
        filename = kw.wfn_file,
        )
    if kw.wfn_type is not None:
        wfn.type = kw.wfn_type
    #end if
    if kw.cutoff is not None:
        wfn.cutoff = kw.cutoff
    #end if
    sim.wavefunction = wfn

    wset = walkerset(
        name        = kw.wset_name,
        )
    if kw.wset_type is not None:
        wset.type = kw.wset_type
    #end if
    if kw.walker_type is not None:
        wset.walker_type = kw.walker_type
    #end if
    sim.walkerset = wset

    prop = propagator(
        name = kw.prop_name,
        info = info.name,
        )
    if kw.hybrid is not None:
        prop.hybrid = kw.hybrid
    #end if
    sim.propagator = prop

    exe = execute(
        info = info.name,
        ham  = ham.name,
        wfn  = wfn.name,
        wset = wset.name,
        prop = prop.name,
        )
    for k in execute.parameters:
        if k in kw and kw[k] is not None:
            exe[k] = kw[k]
        #end if
    #end for
    estimators = []
    valid_estimators = (back_propagation,)
    if kw.estimators is not None:
        for est in kw.estimators:
            invalid = False
            if isinstance(est,QIxml):
                est = est.copy()
            else:
                invalid = True
            #end if
            invalid |= not isinstance(est,valid_estimators)
            if invalid:
                valid_names = [e.__class__.__name__ for e in valid_estimators]
                QmcpackInput.class_error('invalid estimator input encountered\nexpected one of the following: {}\ninputted type: {}\ninputted value: {}'.format(valid_names,est.__class__.__name__,est))
            #end if
            est.incorporate_defaults()
            estimators.append(est)
        #end for
    #end if
    if len(estimators)>0:
        exe.estimators = make_collection(estimators)
    #end if
    sim.execute = exe
    
    qi = QmcpackInput(metadata,sim)

    return qi
#end def generate_basic_afqmc_input



def generate_opt_jastrow_input(id  = 'qmc',
                               series           = 0,
                               purpose          = '',
                               seed             = None,
                               bconds           = None,
                               remove_cell      = False,
                               meshfactor       = 1.0,
                               precision        = 'float',
                               twistnum         = None, 
                               twist            = None,
                               spin_polarized   = False,
                               orbitals_h5      = 'MISSING.h5',
                               system           = None,
                               pseudos          = None,
                               jastrows         = 'generateJ12',
                               corrections      = None,
                               observables      = None,
                               processes        = None,
                               walkers_per_proc = None,
                               threads          = None,
                               decorr           = 10,
                               min_walkers      = None, #use e.g. 128 for gpu's
                               timestep         = 0.5,
                               nonlocalpp       = False,
                               sample_factor    = 1.0,
                               opt_calcs        = None,
                               det_format       = 'new'):
    jastrows = generate_jastrows(jastrows,system)

    if opt_calcs is None:
        opt_calcs = [
            ('linear', 4,  0,  0, 1.0),
            ('linear', 4, .8, .2,   0)
            ]
    #end if
    opts = []
    for opt_calc in opt_calcs:
        if isinstance(opt_calc,QIxml):
            opts.append(opt_calc)
        elif len(opt_calc)==5:
            if opt_calc[0] in opt_map:
                opts.append(
                    generate_opt(
                        *opt_calc,
                         jastrows         = jastrows,
                         processes        = processes,
                         walkers_per_proc = walkers_per_proc,
                         threads          = threads,
                         decorr           = decorr,
                         min_walkers      = min_walkers,
                         timestep         = timestep,
                         nonlocalpp       = nonlocalpp,
                         sample_factor    = sample_factor
                         )
                    )
            else:
                QmcpackInput.class_error('optimization method '+opt_calc[0]+' has not yet been implemented')
            #end if
        else:
            QmcpackInput.class_error('optimization calculation is ill formatted\n  opt calc provided: \n'+str(opt_calc))
        #end if
    #end if

    input = generate_basic_input(
        id             = id             ,
        series         = series         ,
        purpose        = purpose        ,
        seed           = seed           ,
        bconds         = bconds         ,
        remove_cell    = remove_cell    ,
        meshfactor     = meshfactor     ,
        precision      = precision      ,
        twistnum       = twistnum       ,
        twist          = twist          ,
        spin_polarized = spin_polarized ,
        orbitals_h5    = orbitals_h5    ,
        system         = system         ,
        pseudos        = pseudos        ,
        jastrows       = jastrows       ,
        corrections    = corrections    ,
        observables    = observables    ,
        calculations   = opts           ,
        det_format     = det_format     ,
        )

    return input
#end def generate_opt_jastrow_input







if __name__=='__main__':

    filepath = './example_input_files/c_boron/qmcpack.in.xml'

    element_joins=['qmcsystem']
    element_aliases=dict(loop='qmc')
    xml = XMLreader(filepath,element_joins,element_aliases,warn=False).obj
    xml.condense()

    qi = QmcpackInput()
    qi.read(filepath)


    s = qi.simulation
    q = s.qmcsystem
    c = s.calculations
    h = q.hamiltonian
    p = q.particlesets
    w = q.wavefunction
    j = w.jastrows
    co= j.J1.correlations.B.coefficients


    qi.write('./output/qmcpack.in.xml')

    #qi.condensed_name_report()
    #exit()


    test_ret_system    = 1
    test_gen_input     = 0
    test_difference    = 0
    test_moves         = 0
    test_defaults      = 0
    test_substitution  = 0
    test_generation    = 0


    if test_ret_system:
        from structure import generate_structure
        from physical_system import PhysicalSystem
        
        system = PhysicalSystem(
            structure = generate_structure('diamond','fcc','Ge',(2,2,2),scale=5.639,units='A'),
            net_charge = 1,
            net_spin   = 1,
            Ge = 4
        )

        gi = generate_qmcpack_input('basic',system=system)
        
        rsys = gi.return_system()

        print(rsys)

    #end if


    if test_gen_input:
        from structure import generate_structure
        from physical_system import PhysicalSystem
        
        system = PhysicalSystem(
            structure = generate_structure('diamond','fcc','Ge',(2,2,2),scale=5.639,units='A'),
            net_charge = 1,
            net_spin   = 1,
            Ge = 4
        )

        gi = generate_qmcpack_input('basic',system=system)
        
        print(gi)

        print(gi.write())
    #end if



    if test_difference:
        tstep = QmcpackInput('./example_input_files/luke_tutorial/diamond-dmcTsteps.xml')
        opt   = QmcpackInput('./example_input_files/luke_tutorial/opt-diamond.xml')

        different,diff,d1,d2 = tstep.difference(tstep)
        different,diff,d1,d2 = tstep.difference(opt)
        
    #end if



    if test_moves:
        print(50*'=')
        sim = qi.simulation
        print(repr(sim))
        print(repr(sim.qmcsystem))
        print(50*'=')
        qi.move(particleset='simulation')
        print(repr(sim))
        print(repr(sim.qmcsystem))
        print(50*'=')
        qi.standard_placements()
        print(repr(sim))
        print(repr(sim.qmcsystem))

        qi.pluralize()
    #end if



    if test_defaults:
        q=QmcpackInput(
            simulation(
                qmcsystem=section(
                    simulationcell = section(),
                    wavefunction = section(),
                    hamiltonian = section()             
                    ),
                calculations = [
                    cslinear(),
                    vmc(),
                    dmc()
                    ]
                )            
            )

        #q.simulation = simulation()

        q.incorporate_defaults(elements=True)

        print(q)
    #end if


    if test_substitution:
        q = qi.copy()

        q.remove('simulationcell','particleset','wavefunction')
        q.write('./output/qmcpack.remove.xml')
        q.include_xml('./example_input_files/energy_density/Si.ptcl.xml',replace=False)
        q.include_xml('./example_input_files/energy_density/Si.wfs.xml',replace=False)
        q.write('./output/qmcpack.replace.xml')

        qnj = QmcpackInput()
        qnj.read('./example_input_files/jastrowless/opt_jastrow.in.xml')

        qnj.generate_jastrows(size=6)
        qnj.write('./output/jastrow_gen.in.xml')

    #end if
    


    if test_generation:

        q=QmcpackInput(
            meta(
                lattice    = {'units':'bohr'},
                reciprocal = {'units':'2pi/bohr'},
                ionid      = {'datatype':'stringArray'},
                position   = {'datatype':'posArray', 'condition':0}
                ),
            simulation(
                project = section(
                    id='C16B',
                    series = 0,
                    application = section(
                        name = 'qmcpack',
                        role = 'molecu',
                        class_ = 'serial',
                        version = .2
                        ),
                    host = 'kraken',
                    date = '3 May 2012',
                    user = 'jtkrogel'
                    ),
                random = section(seed=13),
                qmcsystem = section(
                    simulationcell = section(
                        name = 'global',
                        lattice = array([[1,1,0],[1,0,1],[0,1,1]]),
                        reciprocal = array([[1,1,-1],[1,-1,1],[-1,1,1]]),
                        bconds = 'p p p',
                        LR_dim_cutoff = 15            
                        ),
                    particlesets = [
                        particleset(
                            name = 'ion0',
                            size = 32,
                            groups=[
                                group(
                                    name='C',
                                    charge=4.
                                    ),
                                group(
                                    name='B',
                                    charge = 3.
                                    )
                                ],
                            ionid = ['B','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C',
                                     'B','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C'],
                            position = array([
                                    [ 0.00, 0.00, 0.00],[ 1.68, 1.68, 1.68],[ 3.37, 3.37, 0.00],
                                    [ 5.05, 5.05, 1.68],[ 3.37, 0.00, 3.37],[ 5.05, 1.68, 5.05],
                                    [ 6.74, 3.37, 3.37],[ 8.42, 5.05, 5.05],[ 0.00, 3.37, 3.37],
                                    [ 1.68, 5.05, 5.05],[ 3.37, 6.74, 3.37],[ 5.05, 8.42, 5.05],
                                    [ 3.37, 3.37, 6.74],[ 5.05, 5.05, 8.42],[ 6.74, 6.74, 6.74],
                                    [ 8.42, 8.42, 8.42],[ 6.74, 6.74, 0.00],[ 8.42, 8.42, 1.68],
                                    [10.11,10.11, 0.00],[11.79,11.79, 1.68],[10.11, 6.74, 3.37],
                                    [11.79, 8.42, 5.05],[13.48,10.11, 3.37],[15.16,11.79, 5.05],
                                    [ 6.74,10.11, 3.37],[ 8.42,11.79, 5.05],[10.11,13.48, 3.37],
                                    [11.79,15.16, 5.05],[10.11,10.11, 6.74],[11.79,11.79, 8.42],
                                    [13.48,13.48, 6.74],[15.16,15.16, 8.42]])
                            ),
                        particleset(
                            name='e',
                            random = 'yes',
                            random_source = 'ion0',
                            groups=[
                                group(
                                    name='u',
                                    size=64,
                                    charge=-1
                                    ),
                                group(
                                    name='d',
                                    size=63,
                                    charge=-1
                                    )                    
                                ]
                            ),
                        ],
                    hamiltonians = [
                        hamiltonian(
                            name='h0',
                            type='generic',
                            target='e',
                            pairpots=[
                                pairpot(
                                    type = 'coulomb',
                                    name = 'ElecElec',
                                    source = 'e',
                                    target = 'e'
                                    ),
                                pairpot(
                                    type = 'pseudo',
                                    name = 'PseudoPot',
                                    source = 'ion0',
                                    wavefunction='psi0',
                                    format='xml',
                                    pseudos = [
                                        pseudo(
                                            elementtype='B',
                                            href='B.pp.xml'
                                            ),
                                        pseudo(
                                            elementtype='C',
                                            href='C.pp.xml'
                                            )
                                        ]
                                    )
                                ],
                            constant = section(
                                type='coulomb',
                                name='IonIon',
                                source='ion0',
                                target='ion0'
                                ),
                            estimators = [
                                estimator(
                                    type='energydensity',
                                    name='edvoronoi',
                                    dynamic='e',
                                    static='ion0',
                                    spacegrid = section(
                                        coord = 'voronoi'
                                        )
                                    ),
                                energydensity(
                                    name='edchempot',
                                    dynamic='e',
                                    static='ion0',
                                    spacegrid=spacegrid(
                                        coord='voronoi',
                                        min_part=-4,
                                        max_part=5
                                        )
                                    ),
                                estimator(
                                    type='energydensity',
                                    name='edcell',
                                    dynamic='e',
                                    static='ion0',
                                    spacegrid = section(
                                        coord = 'cartesian',
                                        origin = section(p1='zero'),
                                        axes   = (
                                            axis(label='x',p1='a1',scale=.5,grid='-1 (192) 1'),
                                            axis(label='y',p1='a2',scale=.5,grid='-1 (1) 1'),
                                            axis(label='z',p1='a3',scale=.5,grid='-1 (1) 1')
                                            )
        #                                axes   = collection(
        #                                    x = section(p1='a1',scale=.5,grid='-1 (192) 1'),
        #                                    y = section(p1='a2',scale=.5,grid='-1 (1) 1'),
        #                                    z = section(p1='a3',scale=.5,grid='-1 (1) 1')
        #                                    )
                                        )
                                    )
                                ]
                            )
                        ],
                    wavefunction = section(
                        name = 'psi0',
                        target = 'e',
                        determinantset = section(
                            type='bspline',
                            href='Si.pwscf.h5',
                            sort = 1,
                            tilematrix = array([[1,0,0],[0,1,0],[0,0,1]]),
                            twistnum = 0,
                            source = 'ion0',
                            slaterdeterminant = section(
                                determinants=[
                                    determinant(
                                        id='updet',
                                        size=64,
                                        occupation = section(
                                            mode='ground',
                                            spindataset=0
                                            )
                                        ),
                                    determinant(
                                        id='downdet',
                                        size=63,
                                        occupation = section(
                                            mode='ground',
                                            spindataset=1
                                            )
                                        )
                                    ]
                                ),
                            ),
                        jastrows = [
                            jastrow(
                                type='two-body',
                                name='J2',
                                function='bspline',
                                print='yes',
                                correlations = [
                                    correlation(
                                        speciesA='u',
                                        speciesB='u',
                                        size=6,
                                        rcut=3.9,
                                        coefficients = section(
                                            id='uu',
                                            type='Array',
                                            coeff=[0,0,0,0,0,0]
                                            )
                                        ),
                                    correlation(
                                        speciesA='u',
                                        speciesB='d',
                                        size=6,
                                        rcut=3.9,
                                        coefficients = section(
                                            id='ud',
                                            type='Array',
                                            coeff=[0,0,0,0,0,0]
                                            )
                                        )
                                    ]
                                ),
                            jastrow(
                                type='one-body',
                                name='J1',
                                function='bspline',
                                source='ion0',
                                print='yes',
                                correlations = [
                                    correlation(
                                        elementtype='C',
                                        size=6,
                                        rcut=3.9,
                                        coefficients = section(
                                            id='eC',
                                            type='Array',
                                            coeff=[0,0,0,0,0,0]
                                            )
                                        ),
                                    correlation(
                                        elementtype='B',
                                        size=6,
                                        rcut=3.9,
                                        coefficients = section(
                                            id='eB',
                                            type='Array',
                                            coeff=[0,0,0,0,0,0]
                                            )
                                        )
                                    ]
                                )
                            ]
                        ),

                    ),
                calculations=[
                    loop(max=4,
                         qmc=qmc(
                            method='cslinear',
                            move='pbyp',
                            checkpoint=-1,
                            gpu='no',
                            blocks = 3125,
                            warmupsteps = 5,
                            steps = 2,
                            samples = 80000,
                            timestep = .5,
                            usedrift = 'yes',
                            minmethod = 'rescale',
                            gevmethod = 'mixed',
                            exp0=-15,
                            nstabilizers = 5,
                            stabilizerscale = 3,
                            stepsize=.35,
                            alloweddifference=1e-5,
                            beta = .05,
                            bigchange = 5.,
                            energy = 0.,
                            unreweightedvariance = 0.,
                            reweightedvariance = 0.,
                            estimators=[
                                estimator(
                                    name='LocalEnergy',
                                    hdf5='no'
                                    )
                                ]
                            )
                        ),
                    qmc(
                        method = 'vmc',
                        multiple = 'no',
                        warp = 'no',
                        move = 'pbyp',
                        walkers = 1,
                        blocks = 2,
                        steps = 500,
                        substeps = 3,
                        timestep = .5,
                        usedrift = 'yes',
                        estimators=[
                            estimator(
                                name='LocalEnergy',
                                hdf5='yes'
                                )
                            ]
                        ),
                    qmc(
                        method='dmc',
                        move='pbyp',
                        walkers = 72,
                        blocks = 2,
                        steps = 50,
                        timestep = .01,
                        nonlocalmove = 'yes',
                        estimators=[
                            estimator(
                                name='LocalEnergy',
                                hdf5='no'
                                )
                            ]            
                        )
                    ]
                )
            )

        q.write('./output/gen.in.xml')





        #something broke this, check later
        exit()
        qs=QmcpackInput(
            simulation = section(
                project = section(
                    id='C16B',series = 0,
                    application = section(name='qmcpack',role='molecu',class_='serial',version=.2),
                    host='kraken',date='3 May 2012',user='jtkrogel'
                    ),
                random = section(seed=13),
                qmcsystem = section(
                    simulationcell = section(
                        name='global',bconds='p p p',lr_dim_cutoff=15,
                        lattice    = [[1,1,0] ,[1,0,1] ,[0,1,1]],
                        reciprocal = [[1,1,-1],[1,-1,1],[-1,1,1]],
                        ),
                    particlesets = collection(
                        ion0=particleset(
                            size=32,
                            groups=collection(
                                C = group(charge=4.),
                                B = group(charge=3.)),
                            ionid = ('B','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C',
                                     'B','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C'),
                            position = [[ 0.00, 0.00, 0.00],[ 1.68, 1.68, 1.68],[ 3.37, 3.37, 0.00],
                                        [ 5.05, 5.05, 1.68],[ 3.37, 0.00, 3.37],[ 5.05, 1.68, 5.05],
                                        [ 6.74, 3.37, 3.37],[ 8.42, 5.05, 5.05],[ 0.00, 3.37, 3.37],
                                        [ 1.68, 5.05, 5.05],[ 3.37, 6.74, 3.37],[ 5.05, 8.42, 5.05],
                                        [ 3.37, 3.37, 6.74],[ 5.05, 5.05, 8.42],[ 6.74, 6.74, 6.74],
                                        [ 8.42, 8.42, 8.42],[ 6.74, 6.74, 0.00],[ 8.42, 8.42, 1.68],
                                        [10.11,10.11, 0.00],[11.79,11.79, 1.68],[10.11, 6.74, 3.37],
                                        [11.79, 8.42, 5.05],[13.48,10.11, 3.37],[15.16,11.79, 5.05],
                                        [ 6.74,10.11, 3.37],[ 8.42,11.79, 5.05],[10.11,13.48, 3.37],
                                        [11.79,15.16, 5.05],[10.11,10.11, 6.74],[11.79,11.79, 8.42],
                                        [13.48,13.48, 6.74],[15.16,15.16, 8.42]]
                            ),
                        e=particleset(
                            random='yes',random_source='ion0',
                            groups = collection(
                                u=group(size=64,charge=-1),
                                d=group(size=63,charge=-1))
                            ),
                        ),
                    hamiltonian = section(
                        name='h0',type='generic',target='e',
                        pairpots=collection(
                            ElecElec = coulomb(name='ElecElec',source='e',target='e'),
                            PseudoPot = pseudopotential(
                                source='ion0',wavefunction='psi0',format='xml',
                                pseudos = collection(
                                    B = pseudo(href='B.pp.xml'),
                                    C = pseudo(href='C.pp.xml'))
                                )
                            ),
                        constant = section(type='coulomb',name='IonIon',source='ion0',target='ion0'),
                        estimators = collection(
                            edvoronoi = energydensity(
                                dynamic='e',static='ion0',spacegrid=section(coord ='voronoi')
                                ),
                            edchempot = energydensity(
                                dynamic='e',static='ion0',
                                spacegrid=section(coord='voronoi',min_part=-4,max_part=5)
                                ),
                            edcell = energydensity(
                                dynamic='e',static='ion0',
                                spacegrid = section(
                                    coord = 'cartesian',
                                    origin = section(p1='zero'),
                                    axes = collection(
                                        x = axis(p1='a1',scale=.5,grid='-1 (192) 1'),
                                        y = axis(p1='a2',scale=.5,grid='-1 (1) 1'),
                                        z = axis(p1='a3',scale=.5,grid='-1 (1) 1'))
                                    )
                                )
                            )
                        ),
                    wavefunction = section(
                        name = 'psi0',target = 'e',
                        determinantset = section(
                            type='bspline',href='Si.pwscf.h5',sort=1,twistnum=0,source='ion0',
                            tilematrix=(1,0,0,0,1,0,0,0,1),
                            slaterdeterminant = section(
                                determinants=collection(
                                    updet = determinant(
                                        size=64,
                                        occupation=section(mode='ground',spindataset=0)
                                        ),
                                    downdet = determinant(
                                        size=63,
                                        occupation = section(mode='ground',spindataset=1))
                                    )
                                ),
                            ),
                        jastrows = collection(
                            J2=jastrow2(
                                function='bspline',print='yes',
                                correlations = collection(
                                    uu=correlation(
                                        speciesA='u',speciesB='u',size=6,rcut=3.9,
                                        coefficients = section(id='uu',type='Array',coeff=[0,0,0,0,0,0])
                                        ),
                                    ud=correlation(
                                        speciesA='u',speciesB='d',size=6,rcut=3.9,
                                        coefficients = section(id='ud',type='Array',coeff=[0,0,0,0,0,0])
                                        )
                                    )
                                ),
                            J1=jastrow1(
                                function='bspline',source='ion0',print='yes',
                                correlations = collection(
                                    C=correlation(
                                        size=6,rcut=3.9,
                                        coefficients = section(
                                            id='eC',type='Array',coeff=[0,0,0,0,0,0])
                                        ),
                                    B=correlation(
                                        size=6,rcut=3.9,
                                        coefficients = section(id='eB',type='Array',coeff=[0,0,0,0,0,0])
                                        )
                                    )
                                )
                            )
                        ),
                    ),
                calculations=(
                    loop(max=4,
                         qmc=cslinear(
                            move='pbyp',checkpoint=-1,gpu='no',
                            blocks      = 3125,
                            warmupsteps = 5,
                            steps       = 2,
                            samples     = 80000,
                            timestep    = .5,
                            usedrift    = 'yes',
                            minmethod   = 'rescale',
                            gevmethod   = 'mixed',
                            exp0              = -15,
                            nstabilizers      =  5,
                            stabilizerscale   =  3,
                            stepsize          =  .35,
                            alloweddifference = 1e-5,
                            beta              = .05,
                            bigchange         = 5.,
                            energy               = 0.,
                            unreweightedvariance = 0.,
                            reweightedvariance   = 0.,
                            estimator = localenergy(hdf5='no')
                            )
                        ),
                    vmc(multiple='no',warp='no',move='pbyp',
                        walkers  =  1,
                        blocks   =  2,
                        steps    = 500,
                        substeps =  3,
                        timestep = .5,
                        usedrift = 'yes',
                        estimator = localenergy(hdf5='no')
                        ),
                    dmc(move='pbyp',
                        walkers  =  72,
                        blocks   =   2,
                        steps    =  50,
                        timestep = .01,
                        nonlocalmove = 'yes',
                        estimator = localenergy(hdf5='yes')
                        )
                    )
                )
            )

        qs.write('./output/simple.in.xml')


        est = qs.simulation.qmcsystem.hamiltonian.estimators
        sg = est.edcell.spacegrid
        print(repr(est))

        exit()



















        q=QmcpackInput()
        q.simulation = section(
            project = section('C16B',0,
                application = section('qmcpack','molecu','serial',.2),
                host = 'kraken',
                date = '3 May 2012',
                user = 'jtkrogel'
                ),
            random = (13),
            qmcsystem = section(
                simulationcell = section(
                    units = 'bohr',
                    lattice = array([[1,1,0],[1,0,1],[0,1,1]]),
                    bconds = 'p p p',
                    LR_dim_cutoff = 15            
                    ),
                particlesets = [
                    particleset('ion0', ('C',4), ('B',3),
                        ionid = ['B','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C',
                                 'B','C','C','C','C','C','C','C','C','C','C','C','C','C','C','C'],
                        position = array([
                                [ 0.00, 0.00, 0.00],[ 1.68, 1.68, 1.68],[ 3.37, 3.37, 0.00],
                                [ 5.05, 5.05, 1.68],[ 3.37, 0.00, 3.37],[ 5.05, 1.68, 5.05],
                                [ 6.74, 3.37, 3.37],[ 8.42, 5.05, 5.05],[ 0.00, 3.37, 3.37],
                                [ 1.68, 5.05, 5.05],[ 3.37, 6.74, 3.37],[ 5.05, 8.42, 5.05],
                                [ 3.37, 3.37, 6.74],[ 5.05, 5.05, 8.42],[ 6.74, 6.74, 6.74],
                                [ 8.42, 8.42, 8.42],[ 6.74, 6.74, 0.00],[ 8.42, 8.42, 1.68],
                                [10.11,10.11, 0.00],[11.79,11.79, 1.68],[10.11, 6.74, 3.37],
                                [11.79, 8.42, 5.05],[13.48,10.11, 3.37],[15.16,11.79, 5.05],
                                [ 6.74,10.11, 3.37],[ 8.42,11.79, 5.05],[10.11,13.48, 3.37],
                                [11.79,15.16, 5.05],[10.11,10.11, 6.74],[11.79,11.79, 8.42],
                                [13.48,13.48, 6.74],[15.16,15.16, 8.42]])
                        ),
                    particleset('e', ('u',-1,64), ('d',-1,63), random_source = 'ion0'),
                    ],
                hamiltonian = section('h0','e',
                    pairpots=[
                        coulomb('ElecElec','e','e'),
                        pseudopotential('PseudoPot','ion0','psi0',('B','B.pp.xml'),('C','C.pp.xml')),
                        coulomb('IonIon','ion0','ion0'),
                        ],
                    estimators = [
                        energydensity('edvoronoi','e','ion0','voronoi',-4,5),
                        energydensity('edcell','e','ion0',
                            spacegrid('cartesian',
                                origin = 'zero',
                                x = ('a1',.5,'-1 (192) 1'),
                                y = ('a2',.5,'-1 (1) 1'),
                                z = ('a3',.5,'-1 (1) 1')
                                )
                            )
                        ]
                    ),
                wavefunction = section('psi0','e',
                    determinantset = section('bspline','Si.pwscf.h5','ion0',
                        sort = 1,
                        tilematrix = array([[1,0,0],[0,1,0],[0,0,1]]),
                        twistnum = 0,
                        slaterdeterminant = [
                            determinant('updet',64,'ground',0),
                            determinant('downdet',63,'ground',1)
                            ],
                        jastrows = [
                            twobody('J2','bspline',
                                    ('u','u',3.9,[0,0,0,0,0,0]),
                                    ('u','d',3.9,[0,0,0,0,0,0])),
                            onebody('J1','bspline','ion0',
                                    ('C',3.9,[0,0,0,0,0,0]),                            
                                    ('B',3.9,[0,0,0,0,0,0]))
                            ]
                        )
                    )
                ),
            calculations=[
                loop(4,
                    cslinear(
                        blocks = 3125,
                        warmupsteps = 5,
                        steps = 2,
                        samples = 80000,
                        timestep = .5,
                        minmethod = 'rescale',
                        gevmethod = 'mixed',
                        exp0=-15,
                        nstabilizers = 5,
                        stabilizerscale = 3,
                        stepsize=.35,
                        alloweddifference=1e-5,
                        beta = .05,
                        bigchange = 5.,
                        energy = 0.,
                        unreweightedvariance = 0.,
                        reweightedvariance = 0.,
                        estimator = localenergy(hdf5='no') 
                        )
                    ),
                vmc(
                    blocks = 2,
                    steps = 500,
                    substeps = 3,
                    timestep = .5,
                    estimator = localenergy(hdf5='yes') 
                    ),
                dmc(
                    walkers = 72,
                    blocks = 2,
                    steps = 50,
                    timestep = .01,
                    nonlocalmove = 'yes',
                    estimator = localenergy(hdf5='no') 
                    )
                ]
            )
    #end if
#end if