scnc
/
qmcpack
mirror of https://github.com/QMCPACK/qmcpack.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Major updates to Nexus: see full text of message 

Changes made to render core of Nexus more robust, a couple of bug fixes, 
several additions to functionality, and several files relegated to "legacy".


=================
new functionality
=================
 -added support for templated input for qmcpack 
    generation is no longer required
    should also work for all supported codes now
    also allows driving of generic simulation codes lacking detailed nexus interface

 -added basic interfaces to drive CRYSTAL and properties
  
 -added support to convert to CASINO pseudopotential format for semilocal PPs
    can convert from GAMESS/CRYSTAL to QMCPACK/CASINO
    also convert between QMCPACK and CASINO
  
 -improved handling of gamma point calculations for pwscf
    better autodetection of gamma point
    explicit user override to forbid "KPOINTS gamma"

 -added support for head node/local job submission 
    useful for lightweight serial tools or scripts in workflows
  
 -added "execute" function to Simulation
    will later enable use generation/running of simulations directly within user script
    will allow for dynamic workflows, rather than dynamic execution of static workflows

 -added enum class to developer.py

 -added general formatted read to TextFile class in fileio.py


=========
bug fixes
=========
 -dependents of failed simulations now guaranteed not to run following interruption and restart of user script

 -run directory creation now guaranteed to occur prior to simulation dependency processing


================
added robustness
================
 -restructured global base class (obj) in generic.py
    remove dependence on abilities.py classes
    streamline and expand base class functions
    move to "new-style" Python classes (inherit from object)

 -added nexuscore "namespace" and NexusCore base class 
    better control of core vs. noncore information propagation
    less chance of downstream breakage (e.g. for derived Simulation classes)
    classes inheriting from NexusCore:
      Settings
      ProjectManager
      Machine
      Job
      Simulation
      SimulationImage
      SimulationInput
      SimulationAnalyzer
  
 -added separate SimulationImage class for save/load of Simulation state
    ensures only central information is saved/loaded
    prevents interaction with derived Simulation classes
  
 -improved organization for Settings class
  
 -removed debug log cruft from ProjectManager
  
 -improved definition of SimulationTemplateInput
    now uses string.Template rather than custom code

 -replaced dynamic type detection functionality from types module with inspect

 -simplified UnitConverter class
  
 -made periodic_table process symbols consistently
  

============
housekeeping
============
 -removed unused files
    converters.py  auxilliary.py
  
 -relegated several superseded/less used files to new legacy directory
    abilities.py         differences.py           qmcpack_pp_validation.py
    cascade.py           eos_fit.py               qmcpack_variations.py
    density_analyzer.py  generator.py             qmc.py
    dft.py               qmcpack_calculations.py
  
 -consolidated orbital converters for qmcpack into qmcpack_converters.py
    replaces convert4qmc.py, pw2qmcpack.py, wfconvert.py
  
 -consolidated numerical functionality into numerics.py
    extended_numpy.py is removed along with ill-advised "from numpy import *"
  
 -consolidated files associated with QmcpackAnalyzer
    plotter.py qaobject.py spacegrid.py
  
 -renamed project.py as nexus.py, removing vacuous nexus.py
    also renamed project_base.py as nexus_base.py 

 -separated processing of basis sets to basissets.py

 -removed outdated BundledQmcpack class


git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@6616 e5b18d87-469d-4833-9cc0-8cdfa06e9491
This commit is contained in:
Jaron Krogel 2015-11-11 19:39:04 +00:00
parent 8db35a353c
commit 461c3a63e9
63 changed files with 5698 additions and 6154 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
nexus/library/atomscf.py
View File

@ -1329,7 +1329,7 @@ class AtomscfInput(SimulationInput):
#end def __init__
def read_contents(self,text):
def read_text(self,text,filepath=None):
raw_lines = text.splitlines()
lines = []
for rline in raw_lines:
@ -1447,10 +1447,10 @@ class AtomscfInput(SimulationInput):
del self.lines
del self.n
del self.cur_section
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self,filepath=None):
self.text = ''
self.cur_section = None
self.text += '{0}\n'.format(self.lim)
@ -1546,7 +1546,7 @@ class AtomscfInput(SimulationInput):
del self.text
del self.cur_section
return text
#end def write_contents
#end def write_text
def nonzero(self,name):
@ -1688,21 +1688,21 @@ class AtomscfInput(SimulationInput):
if len(kwargs)>0:
kw = obj(**kwargs)
element = kw.delete_option('element')
state = kw.delete_option('state')
jcoeff = kw.delete_option('jcoeff')
kcoeff = kw.delete_option('kcoeff')
coeff = kw.delete_option('coeff')
pseudo = kw.delete_option('pseudo')
ppformat = kw.delete_option('ppformat')
basis = kw.delete_option('basis')
guess = kw.delete_option('guess')
element = kw.delete_optional('element')
state = kw.delete_optional('state')
jcoeff = kw.delete_optional('jcoeff')
kcoeff = kw.delete_optional('kcoeff')
coeff = kw.delete_optional('coeff')
pseudo = kw.delete_optional('pseudo')
ppformat = kw.delete_optional('ppformat')
basis = kw.delete_optional('basis')
guess = kw.delete_optional('guess')
use_none = kw.delete_option('use_none' ,False)
use_defaults = kw.delete_option('use_defaults',True)
use_state = kw.delete_option('use_state' ,True)
use_pseudo = kw.delete_option('use_pseudo' ,True)
use_basis = kw.delete_option('use_basis' ,True)
use_none = kw.delete_optional('use_none' ,False)
use_defaults = kw.delete_optional('use_defaults',True)
use_state = kw.delete_optional('use_state' ,True)
use_pseudo = kw.delete_optional('use_pseudo' ,True)
use_basis = kw.delete_optional('use_basis' ,True)
if use_none:
use_defaults = False
use_state = False

8
nexus/library/auxiliary.py
View File

@ -1,8 +0,0 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################

3
nexus/library/bundle.py
View File

@ -46,9 +46,6 @@ class SimulationBundle(Simulation):
generic_identifier = 'bundle'
image_directory = 'bundle'
preserve = Simulation.preserve & set(['sims'])
def __init__(self,*sims,**kwargs):
if len(sims)==1 and isinstance(sims[0],list):
sims = sims[0]

238
nexus/library/convert4qmc.py
View File

@ -1,238 +0,0 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# convert4qmc.py #
# Interface to convert4qmc tool (associated with QMCPACK). #
# #
# Content summary: #
# Convert4qmcInput #
# Class representing command line interface of convert4qmc. #
# #
# Convert4qmcAnalyzer #
# Placeholder class for output analysis. #
# #
# Convert4qmc #
# Class representing conver4qmc instance. #
# #
# generate_convert4qmc_input #
# Function to generate arbitrary convert4qmc input. #
# #
# generate_convert4qmc #
# Function to generate Convert4qmc simulation object. #
# #
#====================================================================#
import os
from generic import obj
from simulation import Simulation,SimulationInput,SimulationAnalyzer
from debug import *
class Convert4qmcInput(SimulationInput):
def __init__(self,
app_name = 'convert4qmc',
prefix = None,
gamess_ascii = None,
ci = None,
read_initial_guess = None,
natural_orbitals = None,
threshold = None,
zero_ci = False,
add_3body_J = False
):
self.prefix = prefix
self.app_name = app_name
self.gamess_ascii = gamess_ascii
self.ci = ci
self.read_initial_guess = read_initial_guess
self.natural_orbitals = natural_orbitals
self.threshold = threshold
self.zero_ci = zero_ci
self.add_3body_J = add_3body_J
#end def __init__
def set_app_name(self,app_name):
self.app_name = app_name
#end def set_app_name
def app_command(self):
c = self.app_name
if self.prefix!=None:
c += ' -prefix '+self.prefix
#end if
if self.gamess_ascii!=None:
c += ' -gamessAscii '+self.gamess_ascii
#end if
if self.ci!=None:
c += ' -ci '+self.ci
#end if
if self.threshold!=None:
c += ' -threshold '+str(self.threshold)
#end if
if self.zero_ci:
c += ' -zeroCi'
#end if
if self.read_initial_guess!=None:
c += ' -readInitialGuess '+str(self.read_initial_guess)
#end if
if self.natural_orbitals!=None:
c += ' -NaturalOrbitals '+str(self.natural_orbitals)
#end if
if self.add_3body_J:
c += ' -add3BodyJ'
#end if
return c
#end def app_command
def read(self,filepath):
None
#end def read
def write_contents(self):
return self.app_command()
#end def write_contents
def output_files(self):
prefix = 'sample'
if self.prefix!=None:
prefix = self.prefix
#end if
wfn_file = prefix+'.Gaussian-G2.xml'
ptcl_file = prefix+'.Gaussian-G2.ptcl.xml'
return wfn_file,ptcl_file
#end def output_files
#end class Convert4qmcInput
def generate_convert4qmc_input(**kwargs):
return Convert4qmcInput(**kwargs)
#end def generate_convert4qmc_input
class Convert4qmcAnalyzer(SimulationAnalyzer):
def __init__(self,arg0):
if isinstance(arg0,Simulation):
self.infile = arg0.infile
else:
self.infile = arg0
#end if
#end def __init__
def analyze(self):
None
#end def analyze
#end class Convert4qmcAnalyzer
class Convert4qmc(Simulation):
input_type = Convert4qmcInput
analyzer_type = Convert4qmcAnalyzer
generic_identifier = 'convert4qmc'
application = 'convert4qmc'
application_properties = set(['serial'])
application_results = set(['orbitals'])
def set_app_name(self,app_name):
self.app_name = app_name
self.input.set_app_name(app_name)
#end def set_app_name
def propagate_identifier(self):
None
#self.input.prefix = self.identifier
#end def propagate_identifier
def check_result(self,result_name,sim):
calculating_result = False
if result_name=='orbitals':
calculating_result = True
else:
calculating_result = False
self.error('ability to check for result '+result_name+' has not been implemented')
#end if
return calculating_result
#end def check_result
def get_result(self,result_name,sim):
result = obj()
input = self.input
if result_name=='orbitals':
wfn_file,ptcl_file = input.output_files()
result.location = os.path.join(self.locdir,wfn_file)
else:
self.error('ability to get result '+result_name+' has not been implemented')
#end if
return result
#end def get_result
def incorporate_result(self,result_name,result,sim):
if result_name=='orbitals':
orbpath = os.path.relpath(result.location,self.locdir)
if result.scftyp=='mcscf':
self.input.gamess_ascii = orbpath
self.input.ci = orbpath
elif result.scftyp=='none': # cisd, etc
self.input.gamess_ascii = orbpath
self.input.ci = orbpath
if result.mos>0:
self.input.read_initial_guess = result.mos
elif result.norbitals>0:
self.input.read_initial_guess = result.norbitals
#end if
else:
self.input.gamess_ascii = orbpath
#end if
self.job.app_command = self.input.app_command()
else:
self.error('ability to incorporate result '+result_name+' has not been implemented')
#end if
#end def incorporate_result
def check_sim_status(self):
output = open(os.path.join(self.locdir,self.outfile),'r').read()
#errors = open(os.path.join(self.locdir,self.errfile),'r').read()
success = 'QMCGaussianParserBase::dump' in output
for filename in self.input.output_files():
success &= os.path.exists(os.path.join(self.locdir,filename))
#end for
self.failed = not success
self.finished = self.job.finished
#end def check_sim_status
def get_output_files(self):
output_files = []
return output_files
#end def get_output_files
def app_command(self):
return self.input.app_command()
#end def app_command
#end class Convert4qmc
def generate_convert4qmc(**kwargs):
sim_args,inp_args = Simulation.separate_inputs(kwargs)
if 'identifier' in sim_args and not 'prefix' in inp_args:
inp_args.prefix = sim_args.identifier
#end if
if not 'input' in sim_args:
sim_args.input = generate_convert4qmc_input(**inp_args)
#end if
convert4qmc = Convert4qmc(**sim_args)
return convert4qmc
#end def generate_convert4qmc

14
nexus/library/converters.py
View File

@ -1,14 +0,0 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# converters.py #
# Collection of converters for QMCPACK. #
# #
#====================================================================#
from pw2qmcpack import Pw2qmcpack,Pw2qmcpackInput,Pw2qmcpackAnalyzer
from wfconvert import Wfconvert,WfconvertInput,WfconvertAnalyzer

178
nexus/library/crystal.py
View File

@ -3,9 +3,12 @@
##################################################################
import os
from generic import obj
from periodic_table import pt as ptable
from pseudopotential import GaussianPP
from nexus_base import nexus_noncore
from simulation import Simulation
from developer import error
@ -170,4 +173,179 @@ def write_hamiltonian(theory = None,
#end def write_hamiltonian
class CrystalSim: None
class PropertiesSim: None
from simulation import GenericSimulation
class Crystal_1(GenericSimulation,CrystalSim):
generic_identifier = 'crystal'
application = 'crystal'
infile_extension = '.d12'
application_results = set(['orbitals'])
def check_result(self,result_name,sim):
return result_name=='orbitals'
#end def check_result
def get_result(self,result_name,sim):
result = obj() # its up to the other application how to handle crystal
if result_name!='orbitals':
self.error('ability to get result '+result_name+' has not been implemented')
#end if
return result
#end def get_result
def incorporate_result(self,result_name,result,sim):
self.error('ability to incorporate result '+result_name+' has not been implemented')
#end def incorporate_result
def app_command(self):
return self.app_name+'<'+self.infile
#end def app_command
#end class Crystal_1
from simulation import generate_simulation,input_template_dev
def gen_crystal_molecule_sim(
identifier = 'crystal',
path = None,
job = None,
title = 'atom',
system = None,
pseudos = None,
occupations = None,
formats = None,
theory = None,
levshift = None,
maxcycle = None,
):
required = obj(path=path,system=system,pseudos=pseudos,occupations=occupations,formats=formats,theory=theory,job=job)
for k,v in required.iteritems():
if v is None:
error(k+' is a required input','gen_crystal_molecule_sim')
#end if
#end for
pseudos_in = pseudos
pseudos = []
for pp in pseudos_in:
pseudos.append(os.path.join(nexus_noncore.pseudo_dir,pp))
#end for
crys_input = input_template_dev()
crys_input.read_text(
write_geometry(
title = title,
bcond = 'molecule',
system = system
)+
write_basis(
pseudos = pseudos,
occupations = occupations,
formats = formats,
)+
write_hamiltonian(
theory = theory,
system = system,
levshift = levshift,
maxcycle = maxcycle
)
)
#crys = generate_simulation(
# identifier = identifier,
# path = path,
# job = job(cores=1,serial=True,app_command='crystal<{0}.d12>&{0}.out&'.format(identifier)),
# input = crys_input,
# infile = '{0}.d12'.format(identifier),
# outfile = '{0}.out'.format(identifier),
# )
crys = Crystal_1(
identifier = identifier,
path = path,
job = job,
input = crys_input,
)
return crys
#end def gen_crystal_molecule_sim
class Properties_1(GenericSimulation,PropertiesSim): # same implementation, not 'is a'
generic_identifier = 'properties'
application = 'properties'
infile_extension = '.d3'
application_results = set(['orbitals'])
def check_result(self,result_name,sim):
return result_name=='orbitals'
#end def check_result
def get_result(self,result_name,sim):
result = obj() # its up to the other application how to handle crystal
if result_name!='orbitals':
self.error('ability to get result '+result_name+' has not been implemented')
#end if
return result
#end def get_result
def incorporate_result(self,result_name,result,sim):
if result_name!='orbitals' or not isinstance(sim,Crystal_1):
self.error('ability to incorporate result '+result_name+' has not been implemented')
#end if
#end def incorporate_result
def app_command(self):
return self.app_name+'<'+self.infile
#end def app_command
#end class Properties_1
molecule_text = '''
NEWK
1 0
CRYAPI_OUT
END
'''
periodic_text = '''
NEWK
0 0
1 0
CRYAPI_OUT
END
'''
def gen_properties(**kwargs):
if 'systype' not in kwargs:
error('systype is a required input','gen_properties')
#end if
systype = kwargs['systype']
del kwargs['systype']
if systype=='molecule_qmc':
text = molecule_text
elif systype=='periodic_qmc':
text = periodic_text
else:
error('invalid systype encountered\nsystype provided: {0}\nvalid options are: molecule_qmc, periodic_qmc'.format(systype))
#end if
sim_args,inp_args = Simulation.separate_inputs(kwargs)
if len(inp_args)>0:
error('invalid arguments encountered\ninvalid keywords: {0}'.format(sorted(inp_args.keys())),'gen_properties')
#end if
if not 'input' in sim_args:
sim_args.input = input_template_dev(text=text.strip())
#end if
prop = Properties_1(**sim_args)
return prop
#end def gen_properties

9
nexus/library/debug.py
View File

@ -23,8 +23,15 @@
import code
import inspect
def ci(locs,globs):
def ci(locs=None,globs=None):
if locs is None or globs is None:
cur_frame = inspect.currentframe()
caller_frame = cur_frame.f_back
locs = caller_frame.f_locals
globs = caller_frame.f_globals
#end if
code.interact(local=dict(globs,**locs))
#end def ci

166
nexus/library/developer.py
View File

@ -27,151 +27,51 @@
# non-existent entities are used during runtime execution. #
# This enables the maximum amount of Nexus functionality to be #
# accessed given the available modules. #
# #
# #
#====================================================================#
import sys
import traceback
from generic import obj
exit_call = exit
devlog = sys.stdout
def log(*items):
s=''
for item in items:
s+=str(item)+' '
#end for
s+='\n'
devlog.write(s)
#end def log
def warn(message,location,header=None,post_header=' Warning:'):
pad = 4*' '
if location is None:
header = 'warning:'
else:
header = location+' warning:'
#end if
log('\n '+header)
log(pad+message.replace('\n','\n'+pad))
#end def warn
def error(message,location=None,exit=True,trace=True):
pad = 4*' '
if location is None:
header = 'error:'
else:
header = location+' error:'
#end if
log('\n '+header)
log(pad+message.replace('\n','\n'+pad))
if exit:
log(' exiting.\n')
if trace:
traceback.print_stack()
#end if
exit_call()
#end if
#end def error
from generic import obj,object_interface,log,error,warn
from debug import ci,interact
class DevBase(obj):
user_interface_data = obj()
dev_instructions_data = obj()
def not_implemented(self):
#self.error('a base class function has not been implemented','Developer')
self.error('a base class function has not been implemented',trace=True)
#end def not_implemented
@classmethod
def set_user_interface(cls,class_variables=None,class_functions=None,member_variables=None,member_functions=None):
if class_variables is None:
class_variables = []
if class_functions is None:
class_functions = []
if member_variables is None:
member_variables = []
if member_functions is None:
member_functions = []
ui = cls.user_interface_data
ui.class_variables = class_variables
ui.class_functions = class_functions
ui.member_variables= member_variables
ui.member_functions= member_functions
#end def set_user_interface
@classmethod
def set_dev_instruction(cls,situation='writing a derived class',class_variables=None,class_functions=None,member_variables=None,member_functions=None):
if class_variables is None:
class_variables = []
if class_functions is None:
class_functions = []
if member_variables is None:
member_variables = []
if member_functions is None:
member_functions = []
ins = obj()
cls.dev_instructions_data[situation] = ins
ins.class_variables = class_variables
ins.class_functions = class_functions
ins.member_variables= member_variables
ins.member_functions= member_functions
#end def set_dev_instruction
@classmethod
def write_dev_data(cls,s,dev_data,indent=''):
p1 = ' '
p2 = 2*p1
p3 = 3*p1
dpad = indent+p3
for name,value in dev_data.iteritems():
s+=indent+p1+name+'\n'
if isinstance(value,list):
for v in value:
s+=indent+p2+v+'\n'
#end for
else:
for v,description in value.iteritems():
s+=indent+p2+v+'\n'
s+=dpad+description.replace('\n','\n'+dpad)
#end for
#end for
#end for
return s
#end def write_dev_data
@classmethod
def user_interface(cls):
s='User Interface for '+cls.__name__+' Class\n'
s+= (len(s)-1)*'='+'\n'
s+=cls.write_dev_data(cls.user_interface_data)
return s
#end def user_interface
@classmethod
def developer_instructions(cls):
s='Developer Instructions for '+cls.__name__+' Class\n'
s+= (len(s)-1)*'='+'\n'
i1=' '
i2=2*i1
i3=3*i1
for situation,instruction in cls.developer_instructions.iteritems():
s+=i1+'situation: '+situation+'\n'
s+=i2+'define the following variables and functions:'
s+=cls.write_dev_data(instructions,i3)
#end for
return s
#end def developer_instructions
#end class DevBase
class enum(object_interface):
def __init__(self,*keys):
if len(keys)==1 and isinstance(keys[0],(list,tuple)):
keys = keys[0]
#end if
n=0
for key in keys:
self[key] = n
n+=1
#end for
#end def __init__
# override some object interface methods
# specifically forbid modification
def __setitem__(self,name,value):
self._error('attempted to modify immutable enum object')
#end def __setitem__
def __delitem__(self,name):
self._error('attempted to modify immutable enum object')
#end def __delitem__
def clear(self):
self._error('attempted to modify immutable enum object')
#end def clear
def _clear(self):
enum.clear(self)
#end def _clear
#end class enum
class Void:
@ -220,8 +120,10 @@ class Void:
@classmethod
def __instancecheck__(cls,*args,**kwargs):
Void._class_unavailable()
@classmethod
def __subclasscheck__(cls,*args,**kwargs):
Void._class_unavailable()
@classmethod
def __subclasshook__(cls,*args,**kwargs):
Void._class_unavailable()
@ -441,8 +343,6 @@ class Void:
Void._unavailable(self)
def __setstate__(self,*args,**kwargs):
Void._unavailable(self)
def __hash__(self,*args,**kwargs):
Void._unavailable(self)
def __bytes__(self,*args,**kwargs):
Void._unavailable(self)
def __format__(self,*args,**kwargs):

2
nexus/library/execute.py
View File

@ -11,7 +11,7 @@
# execute #
# Execute an arbitrary shell command within the user #
# environment of the local machine and wait for its completion. #
# #
# #
#====================================================================#

699
nexus/library/extended_numpy.py
View File

@ -1,699 +0,0 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# extended_numpy.py #
# Extension of functionality available in numpy plus various #
# numerical functions. To be merged into numerics.py #
# #
# Content summary: #
# ndgrid #
# Function to construct an arbitrary N-dimensional grid. #
# Similar to ndgrid from MATLAB. #
# #
# simstats #
# Compute statistics of N-dimensional Monte Carlo simulation #
# data, including mean, variance, error, and autocorrelation. #
# #
# simplestats #
# Compute error assuming uncorrelated data. #
# #
# equilibration_length #
# Estimate equilibration point of Monte Carlo time series data #
# using a heuristic algorithm. #
# #
# ttest #
# Implementation of Student's T-test #
# #
# surface_normals #
# Compute vectors normal to a parametric surface. #
# #
# simple_surface #
# Create a parametric surface in Cartesian, cylindrical, or #
# spherical coordinates. #
# #
# func_fit #
# Perform a fit to an arbitrary function using an arbitrary #
# cost metric (e.g. least squares). #
# #
# distance_table #
# Calculate all N^2 pair distances for a set of N points. #
# #
# nearest_neighbors #
# Find k nearest neighbors of N points using a fast algorithm. #
# #
# voronoi_neighbors #
# Find nearest neighbors in the Voronoi sense, that is for #
# each point, find the points whose Voronoi polyhedra contact #
# the Voronoi polyhedron of that point. #
# #
# convex_hull #
# Find the convex hull of a set of points in N dimensions. #
# #
#====================================================================#
import sys
from numpy import *
from numpy.linalg import *
from developer import unavailable
try:
from scipy.special import betainc
from scipy.optimize import fmin
from scipy.spatial import KDTree,Delaunay,Voronoi
scipy_unavailable = False
except ImportError:
betainc = unavailable('scipy.special' ,'betainc')
fmin = unavailable('scipy.optimize','fmin')
KDTree,Delaunay,Voronoi = unavailable('scipy.spatial' ,'KDTree','Delaunay','Voronoi')
scipy_unavailable = True
#end try
########################################################################
############ ndgrid
########################################################################
# retrieved from
# http://www.mailinglistarchive.com/html/matplotlib-users@lists.sourceforge.net/2010-05/msg00055.html
#"""
#n-dimensional gridding like Matlab's NDGRID
#
#Typical usage:
#>>> x, y, z = [0, 1], [2, 3, 4], [5, 6, 7, 8]
#>>> X, Y, Z = ndgrid(x, y, z)
#
#See ?ndgrid for details.
#"""
def ndgrid(*args, **kwargs):
"""
n-dimensional gridding like Matlab's NDGRID
The input *args are an arbitrary number of numerical sequences,
e.g. lists, arrays, or tuples.
The i-th dimension of the i-th output argument
has copies of the i-th input argument.
Optional keyword argument:
same_dtype : If False (default), the result is an ndarray.
If True, the result is a lists of ndarrays, possibly with
different dtype. This can save space if some *args
have a smaller dtype than others.
Typical usage:
>>> x, y, z = [0, 1], [2, 3, 4], [5, 6, 7, 8]
>>> X, Y, Z = ndgrid(x, y, z) # unpacking the returned ndarray into X, Y, Z
Each of X, Y, Z has shape [len(v) for v in x, y, z].
>>> X.shape == Y.shape == Z.shape == (2, 3, 4)
True
>>> X
array([[[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]],
[[1, 1, 1, 1],
[1, 1, 1, 1],
[1, 1, 1, 1]]])
>>> Y
array([[[2, 2, 2, 2],
[3, 3, 3, 3],
[4, 4, 4, 4]],
[[2, 2, 2, 2],
[3, 3, 3, 3],
[4, 4, 4, 4]]])
>>> Z
array([[[5, 6, 7, 8],
[5, 6, 7, 8],
[5, 6, 7, 8]],
[[5, 6, 7, 8],
[5, 6, 7, 8],
[5, 6, 7, 8]]])
With an unpacked argument list:
>>> V = [[0, 1], [2, 3, 4]]
>>> ndgrid(*V) # an array of two arrays with shape (2, 3)
array([[[0, 0, 0],
[1, 1, 1]],
[[2, 3, 4],
[2, 3, 4]]])
For input vectors of different data types, same_dtype=False makes ndgrid()
return a list of arrays with the respective dtype.
>>> ndgrid([0, 1], [1.0, 1.1, 1.2], same_dtype=False)
[array([[0, 0, 0], [1, 1, 1]]),
array([[ 1. , 1.1, 1.2], [ 1. , 1.1, 1.2]])]
Default is to return a single array.
>>> ndgrid([0, 1], [1.0, 1.1, 1.2])
array([[[ 0. , 0. , 0. ], [ 1. , 1. , 1. ]],
[[ 1. , 1.1, 1.2], [ 1. , 1.1, 1.2]]])
"""
same_dtype = kwargs.get("same_dtype", True)
V = [array(v) for v in args] # ensure all input vectors are arrays
shape = [len(v) for v in args] # common shape of the outputs
result = []
for i, v in enumerate(V):
# reshape v so it can broadcast to the common shape
# http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html
zero = zeros(shape, dtype=v.dtype)
thisshape = ones_like(shape)
thisshape[i] = shape[i]
result.append(zero + v.reshape(thisshape))
if same_dtype:
return array(result) # converts to a common dtype
else:
return result # keeps separate dtype for each output
#if __name__ == "__main__":
# import doctest
# doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
########################################################################
############ End ndgrid
########################################################################
def simstats(x,dim=None):
shape = x.shape
ndim = len(shape)
if dim==None:
dim=ndim-1
#end if
permute = dim!=ndim-1
reshape = ndim>2
nblocks = shape[dim]
if permute:
r = range(ndim)
r.pop(dim)
r.append(dim)
permutation = tuple(r)
r = range(ndim)
r.pop(ndim-1)
r.insert(dim,ndim-1)
invperm = tuple(r)
x=x.transpose(permutation)
shape = tuple(array(shape)[array(permutation)])
dim = ndim-1
#end if
if reshape:
nvars = prod(shape[0:dim])
x=x.reshape(nvars,nblocks)
rdim=dim
dim=1
else:
nvars = shape[0]
#end if
mean = x.mean(dim)
var = x.var(dim)
N=nblocks
if ndim==1:
i=0
tempC=0.5
kappa=0.0
mtmp=mean
if abs(var)<1e-15:
kappa = 1.0
else:
ovar=1.0/var
while (tempC>0 and i<(N-1)):
kappa=kappa+2.0*tempC
i=i+1
#tempC=corr(i,x,mean,var)
tempC = ovar/(N-i)*sum((x[0:N-i]-mtmp)*(x[i:N]-mtmp))
#end while
if kappa == 0.0:
kappa = 1.0
#end if
#end if
Neff=(N+0.0)/(kappa+0.0)
if (Neff == 0.0):
Neff = 1.0
#end if
error=sqrt(var/Neff)
else:
error = zeros(mean.shape)
kappa = zeros(mean.shape)
for v in xrange(nvars):
i=0
tempC=0.5
kap=0.0
vtmp = var[v]
mtmp = mean[v]
if abs(vtmp)<1e-15:
kap = 1.0
else:
ovar = 1.0/vtmp
while (tempC>0 and i<(N-1)):
i += 1
kap += 2.0*tempC
tempC = ovar/(N-i)*sum((x[v,0:N-i]-mtmp)*(x[v,i:N]-mtmp))
#end while
if kap == 0.0:
kap = 1.0
#end if
#end if
Neff=(N+0.0)/(kap+0.0)
if (Neff == 0.0):
Neff = 1.0
#end if
kappa[v]=kap
error[v]=sqrt(vtmp/Neff)
#end for
#end if
if reshape:
x = x.reshape(shape)
mean = mean.reshape(shape[0:rdim])
var = var.reshape(shape[0:rdim])
error = error.reshape(shape[0:rdim])
kappa = kappa.reshape(shape[0:rdim])
#end if
if permute:
x=x.transpose(invperm)
#end if
return (mean,var,error,kappa)
#end def simstats
def simplestats(x,dim=None):
if dim==None:
dim=len(x.shape)-1
#end if
osqrtN = 1.0/sqrt(1.0*x.shape[dim])
mean = x.mean(dim)
error = x.var(dim)*osqrtN
return (mean,error)
#end def simplestats
def equilibration_length(x,tail=.5,plot=False,xlim=None,bounces=2):
bounces = max(1,bounces)
eqlen = 0
nx = len(x)
xt = x[int((1.-tail)*nx+.5):]
nxt = len(xt)
if nxt<10:
return eqlen
#end if
#mean = xh.mean()
#sigma = sqrt(xh.var())
xs = array(xt)
xs.sort()
mean = xs[int(.5*(nxt-1)+.5)]
sigma = (abs(xs[int((.5-.341)*nxt+.5)]-mean)+abs(xs[int((.5+.341)*nxt+.5)]-mean))/2
crossings = bounces*[0,0]
if abs(x[0]-mean)>sigma:
s = -sign(x[0]-mean)
ncrossings = 0
for i in range(nx):
dist = s*(x[i]-mean)
if dist>sigma and dist<5*sigma:
crossings[ncrossings]=i
s*=-1
ncrossings+=1
if ncrossings==2*bounces:
break
#end if
#end if
#end for
bounce = crossings[-2:]
bounce[1] = max(bounce[1],bounce[0])
#print len(x),crossings,crossings[1]-crossings[0]+1
eqlen = bounce[0]+random.randint(bounce[1]-bounce[0]+1)
#end if
if plot:
xlims = xlim
del plot,xlim
from matplotlib.pyplot import plot,figure,show,xlim
figure()
ix = arange(nx)
plot(ix,x,'b.-')
plot([0,nx],[mean,mean],'k-')
plot([0,nx],[mean+sigma,mean+sigma],'r-')
plot([0,nx],[mean-sigma,mean-sigma],'r-')
plot(ix[crossings],x[crossings],'r.')
plot(ix[bounce],x[bounce],'ro')
plot([ix[eqlen],ix[eqlen]],[x.min(),x.max()],'g-')
plot(ix[eqlen],x[eqlen],'go')
if xlims!=None:
xlim(xlims)
#end if
show()
#end if
return eqlen
#end def equilibration_length
def ttest(m1,e1,n1,m2,e2,n2):
m1 = float(m1)
e1 = float(e1)
m2 = float(m2)
e2 = float(e2)
v1 = e1**2
v2 = e2**2
t = (m1-m2)/sqrt(v1+v2)
nu = (v1+v2)**2/(v1**2/(n1-1)+v2**2/(n2-1))
x = nu/(nu+t**2)
p = 1.-betainc(nu/2,.5,x)
return p
#end def ttest
def surface_normals(x,y,z):
nu,nv = x.shape
normals = empty((nu,nv,3))
mi=nu-1
mj=nv-1
v1 = empty((3,))
v2 = empty((3,))
v3 = empty((3,))
dr = empty((3,))
dr[0] = x[0,0]-x[1,0]
dr[1] = y[0,0]-y[1,0]
dr[2] = z[0,0]-z[1,0]
drtol = 1e-4
for i in xrange(nu):
for j in xrange(nv):
iedge = i==0 or i==mi
jedge = j==0 or j==mj
if iedge:
dr[0] = x[0,j]-x[mi,j]
dr[1] = y[0,j]-y[mi,j]
dr[2] = z[0,j]-z[mi,j]
if norm(dr)<drtol:
im = mi-1
ip = 1
elif i==0:
im=i
ip=i+1
elif i==mi:
im=i-1
ip=i
#end if
else:
im=i-1
ip=i+1
#end if
if jedge:
dr[0] = x[i,0]-x[i,mj]
dr[1] = y[i,0]-y[i,mj]
dr[2] = z[i,0]-z[i,mj]
if norm(dr)<drtol:
jm = mj-1
jp = 1
elif j==0:
jm=j
jp=j+1
elif j==mj:
jm=j-1
jp=j
#end if
else:
jm=j-1
jp=j+1
#end if
v1[0] = x[ip,j]-x[im,j]
v1[1] = y[ip,j]-y[im,j]
v1[2] = z[ip,j]-z[im,j]
v2[0] = x[i,jp]-x[i,jm]
v2[1] = y[i,jp]-y[i,jm]
v2[2] = z[i,jp]-z[i,jm]
v3 = cross(v1,v2)
onorm = 1./norm(v3)
normals[i,j,:]=v3[:]*onorm
#end for
#end for
return normals
#end def surface_normals
simple_surface_coords = [set(['x','y','z']),set(['r','phi','z']),set(['r','phi','theta'])]
simple_surface_min = {'x':-1.00000000001,'y':-1.00000000001,'z':-1.00000000001,'r':-0.00000000001,'phi':-0.00000000001,'theta':-0.00000000001}
def simple_surface(origin,axes,grid):
matched=False
gk = set(grid.keys())
for c in range(3):
if gk==simple_surface_coords[c]:
matched=True
coord=c
#end if
#end for
if not matched:
print 'Error in simple_surface: invalid coordinate system provided'
print ' provided coordinates:',gk
print ' permitted coordinates:'
for c in range(3):
print ' ',simple_surface_coords[c]
#end for
exit
#end if
for k,v in grid.iteritems():
if min(v)<simple_surface_min[k]:
print 'Error in simple surface: '+k+' cannot be less than '+str(simple_surface_min[k])
print ' actual minimum: '+str(min(v))
sys.exit()
#end if
if max(v)>1.00000000001:
print 'Error in simple surface: '+k+' cannot be more than 1'
print ' actual maximum: '+str(max(v))
sys.exit()
#end if
#end if
u=empty((3,))
r=empty((3,))
if coord==0:
xl = grid['x']
yl = grid['y']
zl = grid['z']
dim = (len(xl),len(yl),len(zl))
npoints = prod(dim)
points = empty((npoints,3))
n=0
for i in xrange(dim[0]):
for j in xrange(dim[1]):
for k in xrange(dim[2]):
r[0] = xl[i]
r[1] = yl[j]
r[2] = zl[k]
points[n,:] = dot(axes,r) + origin
n+=1
#end for
#end for
#end for
elif coord==1:
rl = grid['r']
phil = 2.*pi*grid['phi']
zl = grid['z']
dim = (len(rl),len(phil),len(zl))
npoints = prod(dim)
points = empty((npoints,3))
n=0
for i in xrange(dim[0]):
for j in xrange(dim[1]):
for k in xrange(dim[2]):
u[0] = rl[i]
u[1] = phil[j]
u[2] = zl[k]
r[0] = u[0]*cos(u[1])
r[1] = u[0]*sin(u[1])
r[2] = u[2]
points[n,:] = dot(axes,r) + origin
n+=1
#end for
#end for
#end for
elif coord==2:
rl = grid['r']
phil = 2.*pi*grid['phi']
thetal = pi*grid['theta']
dim = (len(rl),len(phil),len(thetal))
if dim[0]==1:
sgn = -1. #this is to 'fix' surface normals
#sgn = 1. #this is to 'fix' surface normals
else:
sgn = 1.
#end if
npoints = prod(dim)
points = empty((npoints,3))
n=0
for i in xrange(dim[0]):
for j in xrange(dim[1]):
for k in xrange(dim[2]):
u[0] = rl[i]
u[1] = phil[j]
u[2] = thetal[k]
r[0] = sgn*u[0]*sin(u[2])*cos(u[1])
r[1] = sgn*u[0]*sin(u[2])*sin(u[1])
r[2] = sgn*u[0]*cos(u[2])
points[n,:] = dot(axes,r) + origin
n+=1
#end for
#end for
#end for
#end if
if min(dim)!=1:
print 'Error in simple_surface: minimum dimension must be 1'
print ' actual minimum dimension:',str(min(dim))
sys.exit()
#end if
dm = []
for d in dim:
if d>1:
dm.append(d)
#end if
#end for
dm=tuple(dm)
x = points[:,0].reshape(dm)
y = points[:,1].reshape(dm)
z = points[:,2].reshape(dm)
return x,y,z
#end def simple_surface
least_squares = lambda p,x,y,f: ((f(p,x)-y)**2).sum()
def func_fit(x,y,fitting_function,p0,minimizer=least_squares):
f = fitting_function
p = fmin(minimizer,p0,args=(x,y,f),maxiter=10000,maxfun=10000)
return p
#end def func_fit
def distance_table(p1,p2,ordering=0):
n1 = len(p1)
n2 = len(p2)
if not isinstance(p1,ndarray):
p1=array(p1)
#end if
if not isinstance(p2,ndarray):
p2=array(p2)
#end if
dt = zeros((n1,n2))
for i1 in xrange(n1):
for i2 in xrange(n2):
dt[i1,i2] = norm(p1[i1]-p2[i2])
#end for
#end for
if ordering==0:
return dt
else:
if ordering==1:
n=n1
elif ordering==2:
n=n2
dt=dt.T
else:
print 'distance_table Error: ordering must be 1 or 2,\n you provided '+str(ordering)+'\nexiting.'
exit()
#end if
order = empty(dt.shape,dtype=int)
for i in xrange(n):
o = dt[i].argsort()
order[i] = o
dt[i,:] = dt[i,o]
#end for
return dt,order
#end if
#end def distance_table
def nearest_neighbors(n,points,qpoints=None,return_distances=False,slow=False):
extra = 0
if qpoints==None:
qpoints=points
if len(points)>1:
extra=1
elif return_distances:
return array([]),array([])
else:
return array([])
#end if
#end if
if n>len(qpoints)-extra:
print 'nearest_neighbors Error: requested more than the total number of neighbors\n maximum is: {0}\n you requested: {1}\nexiting.'.format(len(qpoints)-extra,n)
exit()
#end if
slow = slow or scipy_unavailable
if not slow:
kt = KDTree(points)
dist,ind = kt.query(qpoints,n+extra)
else:
dtable,order = distance_table(points,qpoints,ordering=2)
dist = dtable[:,0:n+extra]
ind = order[:,0:n+extra]
#end if
if extra==0 and n==1 and not slow:
nn = atleast_2d(ind).T
else:
nn = ind[:,extra:]
#end if
if not return_distances:
return nn
else:
return nn,dist
#end if
#end def nearest_neighbors
def voronoi_neighbors(points):
vor = Voronoi(points)
neighbor_pairs = vor.ridge_points
return neighbor_pairs
#end def voronoi_neighbors
def convex_hull(points,dimension=None,tol=None):
if dimension is None:
np,dimension = points.shape
#end if
d1 = dimension+1
tri = Delaunay(points)
all_inds = empty((d1,),dtype=bool)
all_inds[:] = True
verts = []
have_tol = tol!=None
for ni in range(len(tri.neighbors)):
n = tri.neighbors[ni]
ns = list(n)
if -1 in ns:
i = ns.index(-1)
inds = all_inds.copy()
inds[i] = False
v = tri.vertices[ni]
if have_tol:
iv = range(d1)
iv.pop(i)
c = points[v[iv[1]]]
a = points[v[i]]-c
b = points[v[iv[0]]]-c
bn = norm(b)
d = norm(a-dot(a,b)/(bn*bn)*b)
if d<tol:
inds[i]=True
#end if
#end if
verts.extend(v[inds])
#end if
#end for
verts = list(set(verts))
return verts
#end def convex_hull

34
nexus/library/fileio.py
View File

@ -73,7 +73,37 @@ class TextFile(DevBase):
def readtokens(self,s=None):
return self.readline(s).split()
#end def readtokens
def readtokensf(self,s=None,*formats):
if s!=None:
self.seek(s)
#end if
self.mm.readline()
line = self.mm.readline()
stokens = line.split()
all_same = False
if len(formats)==1 and len(stokens)>1:
format = formats[0]
all_same = True
elif len(formats)>len(stokens):
self.error('formatted line read failed\nnumber of tokens and provided number of formats do not match\nline: {0}\nnumber of tokens: {1}\nnumber of formats provided: {2}'.format(line,len(stokens),len(formats)))
#end if
tokens = []
if all_same:
for stoken in stokens:
tokens.append(format(stoken))
#end for
else:
for i in xrange(len(formats)):
tokens.append(formats[i](stokens[i]))
#end for
#end if
if len(tokens)==1:
return tokens[0]
else:
return tokens
#end if
#end def readtokensf
# extended mmap interface below
def close(self):
@ -464,7 +494,7 @@ class XsfFile(DevBase):
i+=1
#end while
if check and not self.is_valid():
self.error('read failed, not a valid xsf file'.format(filepath))
self.error('read failed, not a valid xsf file')
#end if
#end def read_text

20
nexus/library/gamess.py
View File

@ -96,6 +96,7 @@ class Gamess(Simulation):
analyzer = self.load_analyzer_image()
if result_name=='orbitals':
result.location = os.path.join(self.locdir,self.outfile)
result.outfile = result.location
result.vec = None # vec from punch
result.norbitals = 0 # orbital count in punch
result.mos = 0 # orbital count (MO's) from log file
@ -154,6 +155,25 @@ class Gamess(Simulation):
output_files = []
return output_files
#end def get_output_files
def output_filename(self,name):
name = name.upper()
if name not in GamessInput.file_units:
self.error('gamess does not produce a file matching the requested description: {0}'.format(name))
#end if
unit = GamessInput.file_units[name]
filename = '{0}.F{1}'.format(self.identifier,str(unit).zfill(2))
return filename
#end def output_filename
def output_filepath(self,name):
filename = self.output_filename(name)
filepath = os.path.join(self.locdir,filename)
filepath = os.path.abspath(filepath)
return filepath
#end def
#end class Gamess

1
nexus/library/gamess_analyzer.py
View File

@ -185,7 +185,6 @@ class GamessAnalyzer(SimulationAnalyzer):
# try to get the mulliken/lowdin populations in each ao
if 'counts' in self:
ao_populations = obj()
self.read_ao_populations(log,ao_populations)
try:
self.read_ao_populations(log,ao_populations)
except:

29
nexus/library/gamess_input.py
View File

@ -37,6 +37,7 @@ from numpy import array,ndarray,abs
from generic import obj
from periodic_table import pt
from developer import DevBase
from nexus_base import nexus_noncore
from simulation import SimulationInput
from debug import *
@ -386,7 +387,7 @@ class ContrlGroup(KeywordSpecGroup):
'coord' ,'units' ,'nzvar' ,'pp' ,'local' ,'ispher','qmttol',
'maxit' ,'molplt','pltorb','aimpac','friend','nfflvl','nprint',
'nosym' ,'etollz','inttyp','grdtyp','normf' ,'normp' ,'itol' ,
'icut' ,'iskprp','irest' ,'geom' ,'ecp'
'icut' ,'iskprp','irest' ,'geom' ,'ecp' ,'casino'
])
integers = set([
'mplevl','icharg','mult' ,'nzvar' ,'ispher','maxit' ,'nfflvl',
@ -394,7 +395,7 @@ class ContrlGroup(KeywordSpecGroup):
'irest'
])
reals = set(['qmttol' ,'etollz'])
bools = set(['numgrd' ,'molplt','pltorb','aimpac'])
bools = set(['numgrd' ,'molplt','pltorb','aimpac','casino'])
strings = set([
'scftyp','dfttyp','tddft' ,'vbtyp' ,'cityp' ,'cctyp' ,'cimtyp',
'relwfn','runtyp','exetyp','coord' ,'units' ,'pp' ,'local' ,
@ -812,7 +813,7 @@ class GamessInput(SimulationInput,GIbase):
#end def __init__
def read_contents(self,contents):
def read_text(self,contents,filepath=None):
#print 8*'\n'
#print contents
@ -929,7 +930,7 @@ class GamessInput(SimulationInput,GIbase):
#print self
#exit()
#end def read_contents
#end def read_text
def process_line(self,ls):
@ -951,7 +952,7 @@ class GamessInput(SimulationInput,GIbase):
#end def process_line
def write_contents(self):
def write_text(self,filepath=None):
contents = ''
extra_groups = set(self.keys())-set(self.group_order)
if len(extra_groups)>0:
@ -968,7 +969,7 @@ class GamessInput(SimulationInput,GIbase):
#end if
#end for
return contents
#end def write_contents
#end def write_text
def incorporate_system(self,system):
@ -997,7 +998,7 @@ def generate_gamess_input(**kwargs):
ps_inputs = set('descriptor symmetry system pseudos pseudo_bases'.split())
ps_inputs = set('descriptor symmetry system pseudos pseudo_bases bases'.split())
ps_defaults = obj()
for var in ps_inputs:
ps_defaults[var]=None
@ -1032,7 +1033,7 @@ def generate_any_gamess_input(**kwargs):
invalid_names = kwset-GamessInput.all_name_aliases
if len(invalid_names)>0:
GamessInput.class_error('invalid group names or keywords encountered\n invalid names/keywords provided: {0}\n please check if these group names or keywords are actually valid GAMESS inputs\n if so, unsupported groups can be generated by providing the keywords as a single argument:\n generate_gamess_input(\n ...,\n group_name = obj(assign keywords),\n ...,\n )'.format(sorted(invalid_names)),'generate_gamess_input')
GamessInput.class_error('invalid group names or keywords encountered\ninvalid names/keywords provided: {0}\nplease check if these group names or keywords are actually valid GAMESS inputs\nif so, unsupported groups can be generated by providing the keywords as a single argument:\n generate_gamess_input(\n ...,\n group_name = obj(assign keywords),\n ...,\n )'.format(sorted(invalid_names)),'generate_gamess_input')
#end if
gi = GamessInput()
@ -1113,11 +1114,19 @@ def generate_any_gamess_input(**kwargs):
data+='\n'
#end if
if pskw.pseudos is None:
gi.contrl.coord = 'cart'
if pskw.bases!=None:
bss = nexus_noncore.basissets.bases_by_atom(*pskw.bases)
else:
bss = obj()
gi.contrl.coord = 'cart'
#end if
for i in range(len(elem)):
a = elem[i]
Z = pt[a].atomic_number
data+='{0} {1:3.2f} {2:16.8f} {3:16.8f} {4:16.8f}\n'.format(a,Z,*pos[i])
if a in bss:
data+=bss[a].text+'\n\n'
#end if
#end for
else:
gi.contrl.set(
@ -1125,7 +1134,7 @@ def generate_any_gamess_input(**kwargs):
ecp = 'read'
)
ecp = ''
pps = GamessInput.pseudopotentials.pseudos_by_atom(*pskw.pseudos)
pps = nexus_noncore.pseudopotentials.pseudos_by_atom(*pskw.pseudos)
atoms = set()
for i in range(len(elem)):
a = elem[i]

1335
nexus/library/generic.py
View File

File diff suppressed because it is too large Load Diff

4
nexus/library/graph.py
View File

@ -10,7 +10,7 @@
# Content summary: #
# Graph #
# Wrapper class for pydot functionality #
# #
# #
#====================================================================#
@ -168,7 +168,7 @@ class Graph(DevBase):
graphs = graphs[0]
#end if
for graph in graphs:
self.add_graph(graph)
self.add_subgraph(graph)
#end for
#end def add_subgraphs

0
nexus/library/abilities.py → nexus/library/legacy/abilities.py
View File

0
nexus/library/cascade.py → nexus/library/legacy/cascade.py
View File

0
nexus/library/density_analyzer.py → nexus/library/legacy/density_analyzer.py
View File

4
nexus/library/dft.py → nexus/library/legacy/dft.py
View File

@ -13,11 +13,11 @@
from project_base import Pobj
from developer import DevBase
import pwscf
class DFTDescriptor(Pobj):
class DFTDescriptor(DevBase):
#general description of a dft simulation
# information transferred to pwscf or abinit input
None

0
nexus/library/differences.py → nexus/library/legacy/differences.py
View File

0
nexus/library/eos_fit.py → nexus/library/legacy/eos_fit.py
View File

0
nexus/library/generator.py → nexus/library/legacy/generator.py
View File

8
nexus/library/qmc.py → nexus/library/legacy/qmc.py
View File

@ -12,19 +12,19 @@
#====================================================================#
from project_base import Pobj
from developer import DevBase
import qmcpack
class HamiltonianDescriptor(Pobj):
class HamiltonianDescriptor(DevBase):
None
#end class HamiltonianDescriptor
class WavefunctionDescriptor(Pobj):
class WavefunctionDescriptor(DevBase):
None
#end class WavefunctionDescriptor
class QMCDescriptor(Pobj):
class QMCDescriptor(DevBase):
#general description of a qmc simulation
# information transferred to qmcpack or casino input
None

0
nexus/library/qmcpack_calculations.py → nexus/library/legacy/qmcpack_calculations.py
View File

6
nexus/library/qmcpack_pp_validation.py → nexus/library/legacy/qmcpack_pp_validation.py
View File

@ -26,7 +26,7 @@ from project import generate_pw2casino,generate_pw2qmcpack
from project import generate_sqd
from project import generate_qmcpack,generate_jastrow
from project import loop,vmc,dmc,linear,cslinear
from project import settings,ProjectManager,Job,Pobj
from project import settings,ProjectManager,Job,Settings
from qmcpack_calculations import basic_qmc,standard_qmc
@ -2778,8 +2778,8 @@ class ValidatePPs(ValidatePPBase):
test_order = ['atoms','dimers']
def __init__(self,settings=None,atoms=None,dimers=None):
if not isinstance(settings,Pobj):
self.error('input variable settings must be the settings object from the Project Suite')
if not isinstance(settings,Settings):
self.error('input variable settings must be the settings object from Nexus')
#end if
ValidatePPBase.settings = settings
if atoms!=None:

6
nexus/library/qmcpack_variations.py → nexus/library/legacy/qmcpack_variations.py
View File

@ -16,7 +16,7 @@ from datetime import datetime
# jtk library imports
from generic import obj # generic.py
from project import Qmcpack,QmcpackInput,QmcpackAnalyzer # project.py
from project import Job,Pobj
from project import Job,nexus_core
from developer import DevBase # developer.py
from debug import * # debug.py
@ -258,7 +258,7 @@ class QmcpackVariations(DevBase):
input.assign(**assignments)
# add the relative path location of the wavefunction file
runpath = os.path.join(Pobj.local_directory,Pobj.runs,path)
runpath = os.path.join(nexus_core.local_directory,nexus_core.runs,path)
wftmp = os.path.join(source_path,wf_h5file)
wfpath = os.path.relpath(wftmp,runpath)
@ -270,7 +270,7 @@ class QmcpackVariations(DevBase):
# check that the build exists
app_loc = os.path.join(self.build_path,build_dir,self.app_loc)
if not os.path.exists(app_loc) and not self.generate_only:
if not os.path.exists(app_loc) and not nexus_core.generate_only:
print ' Error: no qmcapp at '+app_loc
error = True
#end if

94
nexus/library/machines.py
View File

@ -39,7 +39,7 @@
# including arguments to the simulation executable, #
# the run launcher (aprun/mpirun, etc), and the job submitter #
# (e.g. qsub). #
# #
# #
#====================================================================#
@ -50,7 +50,8 @@ from socket import gethostname
from subprocess import Popen,PIPE
from numpy import array,mod,floor,ceil,round,log
from generic import obj
from project_base import Pobj
from developer import DevBase
from nexus_base import NexusCore,nexus_core
from debug import *
@ -79,7 +80,7 @@ def cpu_count():
class Options(Pobj):
class Options(DevBase):
def __init__(self,**kwargs):
self.add(**kwargs)
#end def __init__
@ -130,7 +131,7 @@ class Options(Pobj):
class Job(Pobj):
class Job(NexusCore):
machine = None #default machine if none is specified in settings
@ -170,15 +171,15 @@ class Job(Pobj):
name = 'jobname',
directory = None,
sub_dir = None,
app_name = None,
app_name = None, # name of/path to application
app_flags = None,
app_command = None,
app_props = None,
app = None,
app = None, # name of/path to application
env = None,
user_env = True,
presub = '',
postsub = '',
user_env = True, # import user environment
presub = '', # shell text executed just prior to submission
postsub = '', # shell text executed just after submission
outfile = None,
errfile = None,
mode = None,
@ -187,9 +188,9 @@ class Job(Pobj):
queue = None,
bundled_jobs = None,
relative = False,
cores = None,
nodes = None,
threads = 1,
cores = None, # number of cores for the job
nodes = None, # number of nodes for the job
threads = 1, # number of openmp threads for the job
ppn = None,
compiler = None,
override = None,
@ -197,7 +198,8 @@ class Job(Pobj):
app_options = None,
run_options = None,
sub_options = None,
serial = False,
serial = False, # run job serially, no mpi
local = False, # run job locally, no queue submission
days = 0,
hours = 0,
minutes = 0,
@ -234,6 +236,7 @@ class Job(Pobj):
self.run_options = Options()
self.sub_options = Options()
self.serial = serial
self.local = local
self.days = days
self.hours = hours
self.minutes = minutes
@ -295,6 +298,11 @@ class Job(Pobj):
#end if
#end if
if self.serial:
self.cores = 1
self.nodes = None
#end if
if machine!=None:
self.machine = machine
#end if
@ -308,6 +316,9 @@ class Job(Pobj):
if not machine.batch_capable:
self.error('running batched/bundled jobs on {0} is either not possible or not yet implemented, sorry.'.format(machine.name))
#end if
for job in bundled_jobs: # Paul Young's fix for jobs
machine.process_job(job) # with differing node counts
#end for
#end if
self.normalize_time()
@ -335,7 +346,7 @@ class Job(Pobj):
#end if
if self.subdir is None:
if machine.local_directory!=None:
self.subdir = os.path.join(machine.local_directory,sim.runs,sim.path)
self.subdir = os.path.join(machine.local_directory,nexus_core.runs,sim.path)
self.abs_subdir = self.subdir
else:
self.subdir = self.directory
@ -439,7 +450,7 @@ class Job(Pobj):
def determine_end_status(self,status):
if not self.generate_only:
if not nexus_core.generate_only:
self.successful = False # not really implemented yet
#errfile = os.path.join(self.directory,self.errfile)
#if os.path.exists(errfile):
@ -569,7 +580,7 @@ class Job(Pobj):
class Machine(Pobj):
class Machine(NexusCore):
machines = obj()
@ -830,13 +841,11 @@ class Workstation(Machine):
#end def write_job_states
def query_queue(self):
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~query queue',n=1)
#self.write_job_states('query queue')
self.validate()
done = []
self.dlog('qq',self.processes.keys(),n=2)
for pid,process in self.processes.iteritems():
if self.generate_only or not self.monitor:
if nexus_core.generate_only or not nexus_core.monitor:
qpid,status = pid,0
else:
qpid,status = os.waitpid(pid,os.WNOHANG)
@ -850,7 +859,7 @@ class Workstation(Machine):
self.running.remove(iid)
self.finished.add(iid)
done.append(pid)
if not self.generate_only:
if not nexus_core.generate_only:
job.out.close()
job.err.close()
#end if
@ -859,14 +868,11 @@ class Workstation(Machine):
for pid in done:
del self.processes[pid]
#end for
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end query queue',n=1)
#end def query_queue
def submit_jobs(self):
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~submit jobs',n=1)
self.dlog('sj',self.waiting,n=2)
cores_used = 0
for process in self.processes:
cores_used += process.job.cores
@ -880,14 +886,13 @@ class Workstation(Machine):
job_req.append(job)
core_req.append(job.cores)
#end for
core_req = array(core_req)
job_req = array(job_req)
order = core_req.argsort()
job_req = job_req[order]
core_req = array(core_req,dtype=int)
job_req = array(job_req ,dtype=object)
order = core_req.argsort()
job_req = job_req[order]
self.dlog('cores available:',cores_available,n=2)
for job in job_req:
if job.cores>self.cores and not self.generate_only:
if job.cores>self.cores and not nexus_core.generate_only:
self.error('job '+str(job.internal_id)+' is too large to run on this machine\n cores requested: '+str(job.cores)+'\n machine cores: '+str(self.cores))
#end if
if job.cores<=cores_available:
@ -902,7 +907,6 @@ class Workstation(Machine):
break
#end if
#end for
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end submit jobs',n=1)
#end def submit_jobs
@ -925,11 +929,11 @@ class Workstation(Machine):
job.status = job.states.running
process = obj()
process.job = job
if self.generate_only:
if nexus_core.generate_only:
self.log(pad+'Would have executed: '+command)
job.system_id = job.internal_id
else:
if self.monitor:
if nexus_core.monitor:
self.log(pad+'Executing: '+command)
job.out = open(job.outfile,'w')
job.err = open(job.errfile,'w')
@ -1246,18 +1250,11 @@ class Supercomputer(Machine):
def query_queue(self):
# print
# print 'in query queue'
self.system_queue.clear()
if self.queue_querier=='qstat':
out,err = Popen('qstat -a',shell=True,stdout=PIPE,stderr=PIPE,close_fds=True).communicate()
lines = out.splitlines()
for line in lines:
# print line
tokens=line.split()
if len(tokens)>0:
if '.' in tokens[0]:
@ -1327,8 +1324,6 @@ class Supercomputer(Machine):
else:
spid = tokens[0]
#endif
#print tokens
#print '----------'
if spid.isdigit() and len(tokens)==6: #if account is empty, only 6 tokens.
pid = int(spid)
@ -1354,14 +1349,9 @@ class Supercomputer(Machine):
else:
self.error('ability to query queue with '+self.queue_querier+' has not yet been implemented')
#end if
# print
# print 'system_queue'
# print self.system_queue
done = []
for pid,process in self.processes.iteritems():
if not pid in self.system_queue or self.system_queue[pid]=='complete' or self.generate_only:
if not pid in self.system_queue or self.system_queue[pid]=='complete' or nexus_core.generate_only:
job = process.job
job.status = job.states.finished
job.finished = True
@ -1374,10 +1364,6 @@ class Supercomputer(Machine):
for pid in done:
del self.processes[pid]
#end for
# print
# print 'processes'
# print self.processes
#end def query_queue
@ -1410,7 +1396,7 @@ class Supercomputer(Machine):
self.error('job submission file was not written prior to submission\n submission file: '+os.path.join(job.directory,job.subfile))
#end if
command = self.sub_launcher+job.sub_options.write()+' '+job.subfile
if self.generate_only:
if nexus_core.generate_only:
self.log(pad+'Would have executed: '+command)
job.status = job.states.running
process = obj()
@ -2074,14 +2060,14 @@ class Amos(Supercomputer):
job.seconds =0
#end if
if job.nodes<base_partition:
self.warn('!!! ATTENTION !!!\n number of nodes in {0} should not be less than {1}\n you requested: {2}'.format(job.queue,self.base_partition,job.nodes))
self.warn('!!! ATTENTION !!!\n number of nodes in {0} should not be less than {1}\n you requested: {2}'.format(job.queue,base_partition,job.nodes))
elif job.nodes>max_partition:
self.warn('!!! ATTENTION !!!\n number of nodes in {0} should not be more than {1}\n you requested: {2}'.format(job.queue,max_partition,job.nodes))
else:
if job.queue != 'verylarge':
partition = log(float(job.nodes)/base_partition)/log(2.)
if abs(partition-int(partition))>1e-6:
self.warn('!!! ATTENTION !!!\n number of nodes on {0} must be {1} times a power of two\n you requested: {2}\n nearby valid node count: {3}'.format(self.name,self.base_partition,job.nodes,self.base_partition*2**int(round(partition))))
self.warn('!!! ATTENTION !!!\n number of nodes on {0} must be {1} times a power of two\n you requested: {2}\n nearby valid node count: {3}'.format(self.name,base_partition,job.nodes,base_partition*2**int(round(partition))))
elif job.nodes != 3072 and job.nodes != 4096:
self.warn('!!! ATTENTION !!!\n number of nodes on {0} must be 3072 or 4096 you requested {1}'.format(self.name,job.nodes))
#end if

2
nexus/library/memory.py
View File

@ -16,7 +16,7 @@
# #
# stacksize #
# Return stack size. #
# #
# #
#====================================================================#

343
nexus/library/nexus.py Normal file → Executable file
View File

@ -5,10 +5,349 @@
#====================================================================#
# nexus.py #
# The sum of all user-facing Nexus components. #
# Gathering place for user-facing Nexus functions. Management #
# of user-provided settings. #
# #
# Content summary: #
# Settings #
# Class to set 'global' Nexus data. #
# #
# settings #
# A single Settings instance users interact with as a function. #
# #
# run_project #
# User interface to the ProjectManager. #
# Runs all simulations generated by the user. #
# #
#====================================================================#
from project import *
import os
from generic import obj
from nexus_base import NexusCore,nexus_core,nexus_noncore,nexus_core_noncore
from machines import Job,job,Machine,Supercomputer,get_machine
from simulation import generate_simulation,input_template,multi_input_template,generate_template_input,generate_multi_template_input
from project_manager import ProjectManager
from structure import Structure,generate_structure,generate_cell,read_structure
from physical_system import PhysicalSystem,generate_physical_system
from pseudopotential import Pseudopotential,Pseudopotentials
from basisset import BasisSets
from bundle import bundle
from opium import Opium , OpiumInput , OpiumAnalyzer
from sqd import Sqd , SqdInput , SqdAnalyzer , generate_sqd_input , generate_sqd, hunds_rule_filling
from pwscf import Pwscf , PwscfInput , PwscfAnalyzer , generate_pwscf_input , generate_pwscf
from gamess import Gamess , GamessInput , GamessAnalyzer , generate_gamess_input , generate_gamess, FormattedGroup
from vasp import Vasp , VaspInput , VaspAnalyzer , generate_vasp_input , generate_vasp
from qmcpack import Qmcpack, QmcpackInput, QmcpackAnalyzer, generate_qmcpack_input, generate_qmcpack
from qmcpack_converters import Pw2qmcpack , Pw2qmcpackInput , Pw2qmcpackAnalyzer , generate_pw2qmcpack_input , generate_pw2qmcpack
from qmcpack_converters import Wfconvert , WfconvertInput , WfconvertAnalyzer , generate_wfconvert_input , generate_wfconvert
from qmcpack_converters import Convert4qmc, Convert4qmcInput, Convert4qmcAnalyzer, generate_convert4qmc_input, generate_convert4qmc
from pw2casino import Pw2casino , Pw2casinoInput , Pw2casinoAnalyzer , generate_pw2casino_input , generate_pw2casino
from qmcpack import loop,linear,cslinear,vmc,dmc
from qmcpack import generate_jastrows,generate_jastrow,generate_jastrow1,generate_jastrow2,generate_jastrow3,generate_opt,generate_opts
from qmcpack import generate_cusp_correction
from debug import *
#set the machine if known, otherwise user will provide
hostmachine = Machine.get_hostname()
if Machine.exists(hostmachine):
Job.machine = hostmachine
ProjectManager.machine = Machine.get(hostmachine)
#end if
class Settings(NexusCore):
singleton = None
machine_vars = set('''
machine account machine_info interactive_cores
machine_mode
'''.split())
core_assign_vars = set('''
status_only generate_only runs results
pseudo_dir sleep local_directory remote_directory
monitor skip_submit load_images stages
verbose debug trace
'''.split())
core_process_vars = set('''
file_locations mode
'''.split())
noncore_assign_vars = set('''
basis_dir
'''.split())
noncore_process_vars = set()
gamess_vars = set('''
ericfmt mcppath
'''.split())
nexus_core_vars = core_assign_vars | core_process_vars
nexus_noncore_vars = noncore_assign_vars | noncore_process_vars
nexus_vars = nexus_core_vars | nexus_noncore_vars
allowed_vars = nexus_vars | machine_vars | gamess_vars
@staticmethod
def kw_set(vars,source=None):
kw = obj()
if source!=None:
for n in vars:
if n in source:
kw[n]=source[n]
del source[n]
#end if
#end for
#end if
return kw
#end def null_kw_set
def __init__(self):
if Settings.singleton is None:
Settings.singleton = self
else:
self.error('attempted to create a second Settings object\n please just use the original')
#end if
#end def __init__
def error(self,message,header='settings',exit=True,trace=True):
NexusCore.error(self,message,header,exit,trace)
#end def error
# sets up Nexus core class behavior and passes information to broader class structure
def __call__(self,**kwargs):
NexusCore.write_splash()
self.log('Applying user settings')
# guard against invalid settings
not_allowed = set(kwargs.keys()) - Settings.allowed_vars
if len(not_allowed)>0:
self.error('unrecognized variables provided.\nYou provided: {0}\nAllowed variables are: {1}'.format(sorted(not_allowed),sorted(Settings.allowed_vars)))
#end if
# assign simple variables
for name in Settings.core_assign_vars:
if name in kwargs:
nexus_core[name] = kwargs[name]
#end if
#end for
# assign simple variables
for name in Settings.noncore_assign_vars:
if name in kwargs:
nexus_noncore[name] = kwargs[name]
#end if
#end for
# extract settings based on keyword groups
kw = Settings.kw_set(Settings.nexus_vars ,kwargs)
mach_kw = Settings.kw_set(Settings.machine_vars,kwargs)
gamess_kw = Settings.kw_set(Settings.gamess_vars ,kwargs)
if len(kwargs)>0:
self.error('some settings keywords have not been accounted for\nleftover keywords: {0}\nthis is a developer error'.format(sorted(kwargs.keys())))
#end if
# copy input settings
self.transfer_from(mach_kw.copy())
self.transfer_from(gamess_kw.copy())
# process machine settings
self.process_machine_settings(mach_kw)
# process nexus core settings
self.process_core_settings(kw)
# process nexus noncore settings
self.process_noncore_settings(kw)
# transfer select core data to the global namespace
nexus_core_noncore.transfer_from(nexus_core,nexus_core_noncore.keys())
nexus_noncore.set(**nexus_core_noncore.copy()) # prevent write to core namespace
# copy final core and noncore settings
self.transfer_from(nexus_core.copy())
self.transfer_from(nexus_noncore.copy())
# process gamess settings
Gamess.settings(**gamess_kw)
return
#end def __call__
def process_machine_settings(self,mset):
if 'machine_info' in mset:
machine_info = mset.machine_info
if isinstance(machine_info,dict) or isinstance(machine_info,obj):
for machine_name,minfo in machine_info.iteritems():
mname = machine_name.lower()
if Machine.exists(mname):
machine = Machine.get(mname)
machine.incorporate_user_info(minfo)
else:
self.error('machine {0} is unknown\n cannot set machine_info'.format(machine_name))
#end if
#end for
else:
self.error('machine_info must be a dict or obj\n you provided type '+machine_info.__class__.__name__)
#end if
#end if
if 'machine' in mset:
machine_name = mset.machine
if not Machine.exists(machine_name):
self.error('machine {0} is unknown'.format(machine_name))
#end if
Job.machine = machine_name
ProjectManager.machine = Machine.get(machine_name)
if 'account' in mset:
account = mset.account
if not isinstance(account,str):
self.error('account for {0} must be a string\n you provided: {1}'.format(machine_name,account))
#end if
ProjectManager.machine.account = account
#end if
if 'machine_mode' in mset:
machine_mode = mset.machine_mode
if machine_mode in Machine.modes:
machine_mode = Machine.modes[machine_mode]
#end if
if machine_mode==Machine.modes.interactive:
if ProjectManager.machine==None:
ProjectManager.class_error('no machine specified for interactive mode')
#end if
if not isinstance(ProjectManager.machine,Supercomputer):
self.error('interactive mode is not supported for machine type '+ProjectManager.machine.__class__.__name__)
#end if
if not 'interactive_cores' in mset:
self.error('interactive mode requested, but interactive_cores not set')
#end if
ProjectManager.machine = ProjectManager.machine.interactive_representation(mset.interactive_cores)
#end if
#end if
#end if
#end def process_machine_settings
def process_core_settings(self,kw):
# process project manager settings
if nexus_core.debug:
nexus_core.verbose = True
#end if
if 'mode' in kw:
if kw.mode in nexus_core.modes:
nexus_core.mode = kw.mode
else:
self.error('invalid mode specified: {0}\nvalid modes are {1}'.format(kw.mode,sorted(nexus_core.modes.keys())))
#end if
#end if
mode = nexus_core.mode
modes = nexus_core.modes
if mode==modes.stages:
stages = nexus_core.stages
elif mode==modes.all:
stages = list(nexus_core.primary_modes)
else:
stages = [kw.mode]
#end if
allowed_stages = set(nexus_core.primary_modes)
if isinstance(stages,str):
stages = [stages]
#end if
if len(stages)==0:
stages = list(nexus_core.primary_modes)
elif 'all' in stages:
stages = list(nexus_core.primary_modes)
else:
forbidden = set(nexus_core.stages)-allowed_stages
if len(forbidden)>0:
self.error('some stages provided are not primary stages.\n You provided '+str(list(forbidden))+'\n Options are '+str(list(allowed_stages)))
#end if
#end if
nexus_core.mode = modes.stages
nexus_core.stages = stages
nexus_core.stages_set = set(nexus_core.stages)
# process simulation settings
if 'local_directory' in kw:
nexus_core.file_locations.append(kw.local_directory)
#end if
if 'file_locations' in kw:
fl = kw.file_locations
if isinstance(fl,str):
nexus_core.file_locations.extend([fl])
else:
nexus_core.file_locations.extend(list(fl))
#end if
#end if
if not 'pseudo_dir' in kw:
nexus_core.pseudopotentials = Pseudopotentials()
else:
pseudo_dir = kw.pseudo_dir
nexus_core.file_locations.append(pseudo_dir)
if not os.path.exists(pseudo_dir):
self.error('pseudo_dir "{0}" does not exist'.format(pseudo_dir),trace=False)
#end if
files = os.listdir(pseudo_dir)
ppfiles = []
for f in files:
pf = os.path.join(pseudo_dir,f)
if os.path.isfile(pf):
ppfiles.append(pf)
#end if
#end for
nexus_core.pseudopotentials = Pseudopotentials(ppfiles)
#end if
#end def process_core_settings
def process_noncore_settings(self,kw):
if not 'basis_dir' in kw:
nexus_noncore.basissets = BasisSets()
else:
basis_dir = kw.basis_dir
nexus_core.file_locations.append(basis_dir)
if not os.path.exists(basis_dir):
self.error('basis_dir "{0}" does not exist'.format(basis_dir),trace=False)
#end if
files = os.listdir(basis_dir)
bsfiles = []
for f in files:
pf = os.path.join(basis_dir,f)
if os.path.isfile(pf):
bsfiles.append(pf)
#end if
#end for
nexus_noncore.basissets = BasisSets(bsfiles)
#end if
#end def process_noncore_settings
#end class Settings
# create settings functor for UI
settings = Settings()
def run_project(*args,**kwargs):
pm = ProjectManager()
pm.add_simulations(*args,**kwargs)
pm.run_project()
#end def run_project

195
nexus/library/nexus_base.py Normal file
View File

@ -0,0 +1,195 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# nexus_base.py #
# Provides base class functionality and access to 'global' data #
# for core Nexus classes. #
# #
# Content summary: #
# NexusCore #
# Base class for core Nexus classes. #
# Data intended to be 'global' among core classes is assigned #
# by the Settings class. #
# #
# nexus_core #
# Namespace to be accessed by core Nexus classes. #
# These are classes that inherit from NexusCore directly. #
# #
# nexus_noncore #
# Namespace to be accessed in read-only fashion by non-core #
# classes. #
# #
#====================================================================#
import os
import traceback
import gc as garbage_collector
from memory import resident
from generic import obj
from developer import DevBase
# Nexus namespaces
# nexus_core: to be used by NexusCore classes only
# nexus_noncore: allows read only access to some nexus_core data to non-core classes
nexus_core = obj()
nexus_noncore = obj()
modes = obj(
none = 0,
setup = 1,
send_files = 2,
submit = 3,
get_output = 4,
analyze = 5,
stages = 6,
all = 7
)
garbage_collector.enable()
nexus_noncore.set(
basis_dir = None,
basissets = None,
)
# core namespace elements that can be accessed by noncore classes
nexus_core_noncore = obj(
pseudo_dir = None, # used by: Settings, VaspInput
pseudopotentials = None, # used by: Simulation, GamessInput
)
nexus_core.set(
status_only = False, # used by: ProjectManager
generate_only = False, # used by: Simulation,Machine
sleep = 3, # used by: ProjectManager
runs = 'runs', # used by: Simulation,Machine
results = 'results', # used by: Simulation
local_directory = './', # used by: Simulation,Machine
remote_directory = './', # used by: Simulation
file_locations = ['./'], # used by: Simulation
monitor = True, # used by: ProjectManager,Simulation,Machine
skip_submit = False, # used by: Simulation
load_images = True, # used by: ProjectManager
modes = modes, # used by: ProjectManager,Simulation
mode = modes.stages, # used by: Simulation
stages_set = set(), # used by: ProjectManager,Simulation
stages = [], # used by: Simulation
primary_modes = ['setup','send_files','submit','get_output','analyze'], # used by: Settings
dependent_modes = set(['submit']), # used by: ProjectManager,Simulation
verbose = True, # used by: NexusCore
debug = False, # used by: NexusCore
trace = False, # used by: NexusCore
indent = ' ', # used by: NexusCore
emulate = False, # unused
**nexus_core_noncore
)
nexus_core_no_process = set('''
status_only generate_only sleep
'''.split())
class NexusCore(DevBase):
# garbage collector
gc = garbage_collector
# mutable/dynamic nexus core data
wrote_something = False # for pretty printing
working_directory = None
wrote_splash = False
@staticmethod
def write_splash():
return # don't do this yet
if not NexusCore.wrote_splash:
splash_text = '''
_____________________________________________________
_ _______ _______
| \ /| | ____ \ |\ /| |\ /| | ____ \
| \ | | | | \/ | \ / | | | | | | | \/
| \ | | | |__ \ \_/ / | | | | | |_____
| |\ \| | | __| | _ | | | | | |_____ |
| | \ | | | / / \ \ | | | | | |
| | \ | | |____/\ | / \ | | |___| | /\____| |
|/ \_| |_______/ |/ \| |_______| \_______|
_____________________________________________________
Main author: Jaron T. Krogel
____________________________
'''
print splash_text
NexusCore.wrote_splash = True
#end if
#end def write_splash
@staticmethod
def write_end_splash():
return # don't do this yet
splash_text = '''
_____________________________________________________
_____________________________________________________
'''
print splash_text
#end def write_end_splash
def mem_usage(self):
return int(resident()/1e6)
#end def mem_usage
def log(self,*texts,**kwargs):
if nexus_core.verbose:
if len(kwargs)>0:
n = kwargs['n']
else:
n=0
#end if
text=''
for t in texts:
text+=str(t)+' '
#end for
pad = n*nexus_core.indent
self._logfile.write(pad+text.replace('\n','\n'+pad)+'\n')
#end if
NexusCore.wrote_something = True
#end def log
def dlog(self,*texts,**kwargs):
if nexus_core.debug:
#self.log('mem_usage',self.mem_usage(),n=5)
self.log(*texts,**kwargs)
#end if
#end def dlog
def tlog(self,*texts,**kwargs):
if nexus_core.trace:
self.log(*texts,**kwargs)
w,s,j,f,g,a=int(self.setup),int(self.submitted),int(self.job.finished),int(self.finished),int(self.got_output),int(self.analyzed)
self.log('w,s,j,f,g,a',w,s,j,f,g,a,n=kwargs['n']+1)
#self.log('dependencies',self.dependencies.keys(),n=kwargs['n']+1)
#self.log('dependents ',self.dependents.keys(),n=kwargs['n']+1)
#end if
#end def tlog
def enter(self,directory,changedir=True,msg=''):
NexusCore.working_directory = os.getcwd()
self.log(' Entering '+directory,msg)
if changedir:
os.chdir(directory)
#end if
pad = ' '
return pad
#end def enter
def leave(self):
os.chdir(NexusCore.working_directory)
#end def leave
#end class NexusCore

698
nexus/library/numerics.py
View File

@ -6,8 +6,8 @@
#====================================================================#
# numerics.py #
# A collection of useful numerical functions, currently #
# including specialized curve fitting and jack-knife statistical #
# analysis. #
# including specialized curve fitting, statistical analysis, #
# and spatial analysis. #
# #
# Content summary: #
# morse_spect_fit #
@ -34,21 +34,73 @@
# Can be used to obtain error bars of fit parameters, #
# eigenvalues, and other statistical results that depend on #
# the input data in a non-linear fashion. #
# #
# #
# ndgrid #
# Function to construct an arbitrary N-dimensional grid. #
# Similar to ndgrid from MATLAB. #
# #
# simstats #
# Compute statistics of N-dimensional Monte Carlo simulation #
# data, including mean, variance, error, and autocorrelation. #
# #
# simplestats #
# Compute error assuming uncorrelated data. #
# #
# equilibration_length #
# Estimate equilibration point of Monte Carlo time series data #
# using a heuristic algorithm. #
# #
# ttest #
# Implementation of Student's T-test #
# #
# surface_normals #
# Compute vectors normal to a parametric surface. #
# #
# simple_surface #
# Create a parametric surface in Cartesian, cylindrical, or #
# spherical coordinates. #
# #
# func_fit #
# Perform a fit to an arbitrary function using an arbitrary #
# cost metric (e.g. least squares). #
# #
# distance_table #
# Calculate all N^2 pair distances for a set of N points. #
# #
# nearest_neighbors #
# Find k nearest neighbors of N points using a fast algorithm. #
# #
# voronoi_neighbors #
# Find nearest neighbors in the Voronoi sense, that is for #
# each point, find the points whose Voronoi polyhedra contact #
# the Voronoi polyhedron of that point. #
# #
# convex_hull #
# Find the convex hull of a set of points in N dimensions. #
# #
#====================================================================#
import types as pytypes
import sys
import inspect
from numpy import array,ndarray,zeros,linspace,pi,exp,sqrt,polyfit,polyval
from numpy import sum,abs,arange,empty,sin,cos,dot,atleast_2d,ogrid
from numpy import ones_like,sign,random,cross,prod
from numpy.linalg import norm
from generic import obj
from developer import unavailable,warn,error
from unit_converter import convert
from periodic_table import pt as ptable
from debug import *
try:
from scipy.special import betainc
from scipy.optimize import fmin
from scipy.spatial import KDTree,Delaunay,Voronoi
scipy_unavailable = False
except ImportError:
fmin = unavailable('scipy.optimize','fmin')
betainc = unavailable('scipy.special' ,'betainc')
fmin = unavailable('scipy.optimize','fmin')
KDTree,Delaunay,Voronoi = unavailable('scipy.spatial' ,'KDTree','Delaunay','Voronoi')
scipy_unavailable = True
#end try
@ -443,7 +495,7 @@ numerics_jackknife = jackknife
# auxfunc is an additional function to get a jackknife sample of a derived quantity
def jackknife_aux(jsamples,auxfunc,args=None,kwargs=None,position=None,capture=None):
# unpack the argument list if compressed
if type(auxfunc)!=pytypes.FunctionType:
if not inspect.isfunction(auxfunc):
if len(auxfunc)==1:
auxfunc = auxfunc[0]
elif len(auxfunc)==2:
@ -536,3 +588,635 @@ def check_jackknife_inputs(args,kwargs,position):
#end if
return argpos,kwargpos,args,kwargs,position
#end def check_jackknife_inputs
########################################################################
############ ndgrid
########################################################################
# retrieved from
# http://www.mailinglistarchive.com/html/matplotlib-users@lists.sourceforge.net/2010-05/msg00055.html
#"""
#n-dimensional gridding like Matlab's NDGRID
#
#Typical usage:
#>>> x, y, z = [0, 1], [2, 3, 4], [5, 6, 7, 8]
#>>> X, Y, Z = ndgrid(x, y, z)
#
#See ?ndgrid for details.
#"""
def ndgrid(*args, **kwargs):
"""
n-dimensional gridding like Matlab's NDGRID
The input *args are an arbitrary number of numerical sequences,
e.g. lists, arrays, or tuples.
The i-th dimension of the i-th output argument
has copies of the i-th input argument.
Optional keyword argument:
same_dtype : If False (default), the result is an ndarray.
If True, the result is a lists of ndarrays, possibly with
different dtype. This can save space if some *args
have a smaller dtype than others.
Typical usage:
>>> x, y, z = [0, 1], [2, 3, 4], [5, 6, 7, 8]
>>> X, Y, Z = ndgrid(x, y, z) # unpacking the returned ndarray into X, Y, Z
Each of X, Y, Z has shape [len(v) for v in x, y, z].
>>> X.shape == Y.shape == Z.shape == (2, 3, 4)
True
>>> X
array([[[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]],
[[1, 1, 1, 1],
[1, 1, 1, 1],
[1, 1, 1, 1]]])
>>> Y
array([[[2, 2, 2, 2],
[3, 3, 3, 3],
[4, 4, 4, 4]],
[[2, 2, 2, 2],
[3, 3, 3, 3],
[4, 4, 4, 4]]])
>>> Z
array([[[5, 6, 7, 8],
[5, 6, 7, 8],
[5, 6, 7, 8]],
[[5, 6, 7, 8],
[5, 6, 7, 8],
[5, 6, 7, 8]]])
With an unpacked argument list:
>>> V = [[0, 1], [2, 3, 4]]
>>> ndgrid(*V) # an array of two arrays with shape (2, 3)
array([[[0, 0, 0],
[1, 1, 1]],
[[2, 3, 4],
[2, 3, 4]]])
For input vectors of different data types, same_dtype=False makes ndgrid()
return a list of arrays with the respective dtype.
>>> ndgrid([0, 1], [1.0, 1.1, 1.2], same_dtype=False)
[array([[0, 0, 0], [1, 1, 1]]),
array([[ 1. , 1.1, 1.2], [ 1. , 1.1, 1.2]])]
Default is to return a single array.
>>> ndgrid([0, 1], [1.0, 1.1, 1.2])
array([[[ 0. , 0. , 0. ], [ 1. , 1. , 1. ]],
[[ 1. , 1.1, 1.2], [ 1. , 1.1, 1.2]]])
"""
same_dtype = kwargs.get("same_dtype", True)
V = [array(v) for v in args] # ensure all input vectors are arrays
shape = [len(v) for v in args] # common shape of the outputs
result = []
for i, v in enumerate(V):
# reshape v so it can broadcast to the common shape
# http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html
zero = zeros(shape, dtype=v.dtype)
thisshape = ones_like(shape)
thisshape[i] = shape[i]
result.append(zero + v.reshape(thisshape))
if same_dtype:
return array(result) # converts to a common dtype
else:
return result # keeps separate dtype for each output
#if __name__ == "__main__":
# import doctest
# doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
########################################################################
############ End ndgrid
########################################################################
def simstats(x,dim=None):
shape = x.shape
ndim = len(shape)
if dim==None:
dim=ndim-1
#end if
permute = dim!=ndim-1
reshape = ndim>2
nblocks = shape[dim]
if permute:
r = range(ndim)
r.pop(dim)
r.append(dim)
permutation = tuple(r)
r = range(ndim)
r.pop(ndim-1)
r.insert(dim,ndim-1)
invperm = tuple(r)
x=x.transpose(permutation)
shape = tuple(array(shape)[array(permutation)])
dim = ndim-1
#end if
if reshape:
nvars = prod(shape[0:dim])
x=x.reshape(nvars,nblocks)
rdim=dim
dim=1
else:
nvars = shape[0]
#end if
mean = x.mean(dim)
var = x.var(dim)
N=nblocks
if ndim==1:
i=0
tempC=0.5
kappa=0.0
mtmp=mean
if abs(var)<1e-15:
kappa = 1.0
else:
ovar=1.0/var
while (tempC>0 and i<(N-1)):
kappa=kappa+2.0*tempC
i=i+1
#tempC=corr(i,x,mean,var)
tempC = ovar/(N-i)*sum((x[0:N-i]-mtmp)*(x[i:N]-mtmp))
#end while
if kappa == 0.0:
kappa = 1.0
#end if
#end if
Neff=(N+0.0)/(kappa+0.0)
if (Neff == 0.0):
Neff = 1.0
#end if
error=sqrt(var/Neff)
else:
error = zeros(mean.shape)
kappa = zeros(mean.shape)
for v in xrange(nvars):
i=0
tempC=0.5
kap=0.0
vtmp = var[v]
mtmp = mean[v]
if abs(vtmp)<1e-15:
kap = 1.0
else:
ovar = 1.0/vtmp
while (tempC>0 and i<(N-1)):
i += 1
kap += 2.0*tempC
tempC = ovar/(N-i)*sum((x[v,0:N-i]-mtmp)*(x[v,i:N]-mtmp))
#end while
if kap == 0.0:
kap = 1.0
#end if
#end if
Neff=(N+0.0)/(kap+0.0)
if (Neff == 0.0):
Neff = 1.0
#end if
kappa[v]=kap
error[v]=sqrt(vtmp/Neff)
#end for
#end if
if reshape:
x = x.reshape(shape)
mean = mean.reshape(shape[0:rdim])
var = var.reshape(shape[0:rdim])
error = error.reshape(shape[0:rdim])
kappa = kappa.reshape(shape[0:rdim])
#end if
if permute:
x=x.transpose(invperm)
#end if
return (mean,var,error,kappa)
#end def simstats
def simplestats(x,dim=None):
if dim==None:
dim=len(x.shape)-1
#end if
osqrtN = 1.0/sqrt(1.0*x.shape[dim])
mean = x.mean(dim)
error = x.var(dim)*osqrtN
return (mean,error)
#end def simplestats
def equilibration_length(x,tail=.5,plot=False,xlim=None,bounces=2):
bounces = max(1,bounces)
eqlen = 0
nx = len(x)
xt = x[int((1.-tail)*nx+.5):]
nxt = len(xt)
if nxt<10:
return eqlen
#end if
#mean = xh.mean()
#sigma = sqrt(xh.var())
xs = array(xt)
xs.sort()
mean = xs[int(.5*(nxt-1)+.5)]
sigma = (abs(xs[int((.5-.341)*nxt+.5)]-mean)+abs(xs[int((.5+.341)*nxt+.5)]-mean))/2
crossings = bounces*[0,0]
if abs(x[0]-mean)>sigma:
s = -sign(x[0]-mean)
ncrossings = 0
for i in range(nx):
dist = s*(x[i]-mean)
if dist>sigma and dist<5*sigma:
crossings[ncrossings]=i
s*=-1
ncrossings+=1
if ncrossings==2*bounces:
break
#end if
#end if
#end for
bounce = crossings[-2:]
bounce[1] = max(bounce[1],bounce[0])
#print len(x),crossings,crossings[1]-crossings[0]+1
eqlen = bounce[0]+random.randint(bounce[1]-bounce[0]+1)
#end if
if plot:
xlims = xlim
del plot,xlim
from matplotlib.pyplot import plot,figure,show,xlim
figure()
ix = arange(nx)
plot(ix,x,'b.-')
plot([0,nx],[mean,mean],'k-')
plot([0,nx],[mean+sigma,mean+sigma],'r-')
plot([0,nx],[mean-sigma,mean-sigma],'r-')
plot(ix[crossings],x[crossings],'r.')
plot(ix[bounce],x[bounce],'ro')
plot([ix[eqlen],ix[eqlen]],[x.min(),x.max()],'g-')
plot(ix[eqlen],x[eqlen],'go')
if xlims!=None:
xlim(xlims)
#end if
show()
#end if
return eqlen
#end def equilibration_length
def ttest(m1,e1,n1,m2,e2,n2):
m1 = float(m1)
e1 = float(e1)
m2 = float(m2)
e2 = float(e2)
v1 = e1**2
v2 = e2**2
t = (m1-m2)/sqrt(v1+v2)
nu = (v1+v2)**2/(v1**2/(n1-1)+v2**2/(n2-1))
x = nu/(nu+t**2)
p = 1.-betainc(nu/2,.5,x)
return p
#end def ttest
def surface_normals(x,y,z):
nu,nv = x.shape
normals = empty((nu,nv,3))
mi=nu-1
mj=nv-1
v1 = empty((3,))
v2 = empty((3,))
v3 = empty((3,))
dr = empty((3,))
dr[0] = x[0,0]-x[1,0]
dr[1] = y[0,0]-y[1,0]
dr[2] = z[0,0]-z[1,0]
drtol = 1e-4
for i in xrange(nu):
for j in xrange(nv):
iedge = i==0 or i==mi
jedge = j==0 or j==mj
if iedge:
dr[0] = x[0,j]-x[mi,j]
dr[1] = y[0,j]-y[mi,j]
dr[2] = z[0,j]-z[mi,j]
if norm(dr)<drtol:
im = mi-1
ip = 1
elif i==0:
im=i
ip=i+1
elif i==mi:
im=i-1
ip=i
#end if
else:
im=i-1
ip=i+1
#end if
if jedge:
dr[0] = x[i,0]-x[i,mj]
dr[1] = y[i,0]-y[i,mj]
dr[2] = z[i,0]-z[i,mj]
if norm(dr)<drtol:
jm = mj-1
jp = 1
elif j==0:
jm=j
jp=j+1
elif j==mj:
jm=j-1
jp=j
#end if
else:
jm=j-1
jp=j+1
#end if
v1[0] = x[ip,j]-x[im,j]
v1[1] = y[ip,j]-y[im,j]
v1[2] = z[ip,j]-z[im,j]
v2[0] = x[i,jp]-x[i,jm]
v2[1] = y[i,jp]-y[i,jm]
v2[2] = z[i,jp]-z[i,jm]
v3 = cross(v1,v2)
onorm = 1./norm(v3)
normals[i,j,:]=v3[:]*onorm
#end for
#end for
return normals
#end def surface_normals
simple_surface_coords = [set(['x','y','z']),set(['r','phi','z']),set(['r','phi','theta'])]
simple_surface_min = {'x':-1.00000000001,'y':-1.00000000001,'z':-1.00000000001,'r':-0.00000000001,'phi':-0.00000000001,'theta':-0.00000000001}
def simple_surface(origin,axes,grid):
matched=False
gk = set(grid.keys())
for c in range(3):
if gk==simple_surface_coords[c]:
matched=True
coord=c
#end if
#end for
if not matched:
print 'Error in simple_surface: invalid coordinate system provided'
print ' provided coordinates:',gk
print ' permitted coordinates:'
for c in range(3):
print ' ',simple_surface_coords[c]
#end for
sys.exit()
#end if
for k,v in grid.iteritems():
if min(v)<simple_surface_min[k]:
print 'Error in simple surface: '+k+' cannot be less than '+str(simple_surface_min[k])
print ' actual minimum: '+str(min(v))
sys.exit()
#end if
if max(v)>1.00000000001:
print 'Error in simple surface: '+k+' cannot be more than 1'
print ' actual maximum: '+str(max(v))
sys.exit()
#end if
#end if
u=empty((3,))
r=empty((3,))
if coord==0:
xl = grid['x']
yl = grid['y']
zl = grid['z']
dim = (len(xl),len(yl),len(zl))
npoints = prod(dim)
points = empty((npoints,3))
n=0
for i in xrange(dim[0]):
for j in xrange(dim[1]):
for k in xrange(dim[2]):
r[0] = xl[i]
r[1] = yl[j]
r[2] = zl[k]
points[n,:] = dot(axes,r) + origin
n+=1
#end for
#end for
#end for
elif coord==1:
rl = grid['r']
phil = 2.*pi*grid['phi']
zl = grid['z']
dim = (len(rl),len(phil),len(zl))
npoints = prod(dim)
points = empty((npoints,3))
n=0
for i in xrange(dim[0]):
for j in xrange(dim[1]):
for k in xrange(dim[2]):
u[0] = rl[i]
u[1] = phil[j]
u[2] = zl[k]
r[0] = u[0]*cos(u[1])
r[1] = u[0]*sin(u[1])
r[2] = u[2]
points[n,:] = dot(axes,r) + origin
n+=1
#end for
#end for
#end for
elif coord==2:
rl = grid['r']
phil = 2.*pi*grid['phi']
thetal = pi*grid['theta']
dim = (len(rl),len(phil),len(thetal))
if dim[0]==1:
sgn = -1. #this is to 'fix' surface normals
#sgn = 1. #this is to 'fix' surface normals
else:
sgn = 1.
#end if
npoints = prod(dim)
points = empty((npoints,3))
n=0
for i in xrange(dim[0]):
for j in xrange(dim[1]):
for k in xrange(dim[2]):
u[0] = rl[i]
u[1] = phil[j]
u[2] = thetal[k]
r[0] = sgn*u[0]*sin(u[2])*cos(u[1])
r[1] = sgn*u[0]*sin(u[2])*sin(u[1])
r[2] = sgn*u[0]*cos(u[2])
points[n,:] = dot(axes,r) + origin
n+=1
#end for
#end for
#end for
#end if
if min(dim)!=1:
print 'Error in simple_surface: minimum dimension must be 1'
print ' actual minimum dimension:',str(min(dim))
sys.exit()
#end if
dm = []
for d in dim:
if d>1:
dm.append(d)
#end if
#end for
dm=tuple(dm)
x = points[:,0].reshape(dm)
y = points[:,1].reshape(dm)
z = points[:,2].reshape(dm)
return x,y,z
#end def simple_surface
#least_squares = lambda p,x,y,f: ((f(p,x)-y)**2).sum()
def func_fit(x,y,fitting_function,p0,minimizer=least_squares):
f = fitting_function
p = fmin(minimizer,p0,args=(x,y,f),maxiter=10000,maxfun=10000)
return p
#end def func_fit
def distance_table(p1,p2,ordering=0):
n1 = len(p1)
n2 = len(p2)
if not isinstance(p1,ndarray):
p1=array(p1)
#end if
if not isinstance(p2,ndarray):
p2=array(p2)
#end if
dt = zeros((n1,n2))
for i1 in xrange(n1):
for i2 in xrange(n2):
dt[i1,i2] = norm(p1[i1]-p2[i2])
#end for
#end for
if ordering==0:
return dt
else:
if ordering==1:
n=n1
elif ordering==2:
n=n2
dt=dt.T
else:
print 'distance_table Error: ordering must be 1 or 2,\n you provided '+str(ordering)+'\nexiting.'
exit()
#end if
order = empty(dt.shape,dtype=int)
for i in xrange(n):
o = dt[i].argsort()
order[i] = o
dt[i,:] = dt[i,o]
#end for
return dt,order
#end if
#end def distance_table
def nearest_neighbors(n,points,qpoints=None,return_distances=False,slow=False):
extra = 0
if qpoints==None:
qpoints=points
if len(points)>1:
extra=1
elif return_distances:
return array([]),array([])
else:
return array([])
#end if
#end if
if n>len(qpoints)-extra:
print 'nearest_neighbors Error: requested more than the total number of neighbors\n maximum is: {0}\n you requested: {1}\nexiting.'.format(len(qpoints)-extra,n)
exit()
#end if
slow = slow or scipy_unavailable
if not slow:
kt = KDTree(points)
dist,ind = kt.query(qpoints,n+extra)
else:
dtable,order = distance_table(points,qpoints,ordering=2)
dist = dtable[:,0:n+extra]
ind = order[:,0:n+extra]
#end if
if extra==0 and n==1 and not slow:
nn = atleast_2d(ind).T
else:
nn = ind[:,extra:]
#end if
if not return_distances:
return nn
else:
return nn,dist
#end if
#end def nearest_neighbors
def voronoi_neighbors(points):
vor = Voronoi(points)
neighbor_pairs = vor.ridge_points
return neighbor_pairs
#end def voronoi_neighbors
def convex_hull(points,dimension=None,tol=None):
if dimension is None:
np,dimension = points.shape
#end if
d1 = dimension+1
tri = Delaunay(points)
all_inds = empty((d1,),dtype=bool)
all_inds[:] = True
verts = []
have_tol = tol!=None
for ni in range(len(tri.neighbors)):
n = tri.neighbors[ni]
ns = list(n)
if -1 in ns:
i = ns.index(-1)
inds = all_inds.copy()
inds[i] = False
v = tri.vertices[ni]
if have_tol:
iv = range(d1)
iv.pop(i)
c = points[v[iv[1]]]
a = points[v[i]]-c
b = points[v[iv[0]]]-c
bn = norm(b)
d = norm(a-dot(a,b)/(bn*bn)*b)
if d<tol:
inds[i]=True
#end if
#end if
verts.extend(v[inds])
#end if
#end for
verts = list(set(verts))
return verts
#end def convex_hull

14
nexus/library/opium.py
View File

@ -37,7 +37,7 @@ from simulation import Simulation,SimulationInput,SimulationAnalyzer
# in a form OpiumInput understands (this depends on your implementation)
#
# required functions to be implemented:
# OpiumInput: read_contents, write_contents
# OpiumInput: read_text, write_text
# Opium: app_command, check_sim_status
#
# 2) generated standalone simulation
@ -45,7 +45,7 @@ from simulation import Simulation,SimulationInput,SimulationAnalyzer
# generate functions provide a short-hand of minimal vars for input
#
# required functions to be implemented:
# OpiumInput: read_contents, write_contents
# OpiumInput: read_text, write_text
# Opium: app_command, check_sim_status
# generate_opium_input
#
@ -57,7 +57,7 @@ from simulation import Simulation,SimulationInput,SimulationAnalyzer
# this information is used by the others to populate input files
#
# required functions to be implemented:
# OpiumInput: read_contents, write_contents
# OpiumInput: read_text, write_text
# Opium: app_command,check_sim_status,
# check_result, get_result
#
@ -413,7 +413,7 @@ class OpiumInput(SimulationInput):
#end def __init__
def read_contents(self,contents):
def read_text(self,contents,filepath=None):
lines = contents.splitlines()
sections = obj()
sec=None
@ -446,10 +446,10 @@ class OpiumInput(SimulationInput):
section.read(sectext)
self[secname] = section
#end for
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self,filepath=None):
contents = ''
for secname in self.section_order:
secname = secname.lower()
@ -458,7 +458,7 @@ class OpiumInput(SimulationInput):
#end if
#end for
return contents
#end def write_contents
#end def write_text
#end class OpiumInput

10
nexus/library/periodic_table.py
View File

@ -1878,15 +1878,17 @@ def is_element(name,symbol=False):
nlen = len(name)
if name.find('_')!=-1:
s,n = name.split('_',1)
#iselem = n.isdigit() and s in periodic_table.elements
iselem = s in periodic_table.elements
elif nlen>1 and name[1:].isdigit():
s = name[0:1]
iselem = s in periodic_table.elements
elif nlen>2 and name[2:].isdigit():
s = name[0:2]
iselem = s in periodic_table.elements
#end if
if len(s)==1:
s = s.upper()
elif len(s)==2:
s = s[0].upper()+s[1].lower()
#end if
iselem = s in periodic_table.elements
#end if
#end if
if symbol:

88
nexus/library/plotter.py
View File

@ -1,88 +0,0 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# plotter.py #
# Interface with mayavi to plot isosurfaces and surfaces slices #
# of volume data. Used in the past to visualize energy densities.#
# See qmcpack_quantity_analyzers.py and spacegrid.py. #
# #
# Content summary: #
# Plotter #
# Wrapper class for mayavi visualization. #
# #
#====================================================================#
from abilities import AllAbilities
import numpy as np
class Plotter(AllAbilities):
def __init__(self):
self.initialized = False
return
#end def __init__
def ensure_init(self):
if not self.initialized:
from enthought.mayavi import mlab
from enthought.tvtk.api import tvtk
self.mlab = mlab
self.tvtk = tvtk
self.show = mlab.show
self.plot3d = mlab.plot3d
self.mesh = mlab.mesh
self.initialized = True
#end if
#end def ensure_init
def isosurface(self,points,scalars,contours,dimensions,name='val'):
self.ensure_init()
mlab = self.mlab
tvtk = self.tvtk
sg=tvtk.StructuredGrid(dimensions=dimensions,points=points)
sg.point_data.scalars = scalars
sg.point_data.scalars.name = name
d = mlab.pipeline.add_dataset(sg)
iso = mlab.pipeline.iso_surface(d)
if isinstance(contours,int):
iso.contour.number_of_contours = contours
elif isinstance(contours,list):
iso.contour.auto_contours = False
iso.contour.contours = contours
else:
self.error('isosurface contours must be an int or list\n a '+str(type(contours))+' was provided instead')
#end if
return
#end def isosurface
def surface_slice(self,x,y,z,scalars,options=None):
scale = 1.0
opacity= 1.0
if options!=None:
if 'norm_height' in options:
scale = options.norm_height/abs(scalars.max())
if 'scale' in options:
scale = options.scale
if 'opacity' in options:
opacity = options.opacity
#end if
self.ensure_init()
from extended_numpy import surface_normals
self.mesh(x,y,z,opacity=.2)
surfnorm = scale*surface_normals(x,y,z)
xs=x.copy()
ys=y.copy()
zs=z.copy()
xs[...] = x[...] + surfnorm[...,0]*scalars[...]
ys[...] = y[...] + surfnorm[...,1]*scalars[...]
zs[...] = z[...] + surfnorm[...,2]*scalars[...]
self.mesh(xs,ys,zs,scalars=scalars,opacity=opacity)
return
#end def surface_slice
#end class Plotter

271
nexus/library/project.py
View File

@ -1,271 +0,0 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# project.py #
# Gathering place for user-facing Nexus functions. Management #
# of user-provided settings. #
# #
# Content summary: #
# Settings #
# Class to set 'global' Nexus data. #
# #
# settings #
# A single Settings instance users interact with as a function. #
# #
# run_project #
# User interface to the ProjectManager. #
# Runs all simulations generated by the user. #
# #
#====================================================================#
import os
from generic import obj
from project_base import Pobj,modes
from project_manager import ProjectManager
from machines import Job,job,Machine,Supercomputer,get_machine
from structure import Structure,generate_structure,generate_cell,read_structure
from physical_system import PhysicalSystem,generate_physical_system
from pseudopotential import Pseudopotential,Pseudopotentials
from simulation import generate_simulation,input_template,multi_input_template,generate_template_input,generate_multi_template_input
from bundle import bundle
from opium import Opium , OpiumInput , OpiumAnalyzer
from sqd import Sqd , SqdInput , SqdAnalyzer , generate_sqd_input , generate_sqd, hunds_rule_filling
from pwscf import Pwscf , PwscfInput , PwscfAnalyzer , generate_pwscf_input , generate_pwscf
from pw2casino import Pw2casino , Pw2casinoInput , Pw2casinoAnalyzer , generate_pw2casino_input , generate_pw2casino
from pw2qmcpack import Pw2qmcpack , Pw2qmcpackInput , Pw2qmcpackAnalyzer , generate_pw2qmcpack_input , generate_pw2qmcpack
from wfconvert import Wfconvert , WfconvertInput , WfconvertAnalyzer , generate_wfconvert_input , generate_wfconvert
from gamess import Gamess , GamessInput , GamessAnalyzer , generate_gamess_input , generate_gamess, FormattedGroup
from convert4qmc import Convert4qmc, Convert4qmcInput, Convert4qmcAnalyzer, generate_convert4qmc_input, generate_convert4qmc
from qmcpack import Qmcpack , QmcpackInput , QmcpackAnalyzer , generate_qmcpack_input , generate_qmcpack
from vasp import Vasp , VaspInput , VaspAnalyzer , generate_vasp_input , generate_vasp
from qmcpack import loop,linear,cslinear,vmc,dmc
from qmcpack import generate_jastrows,generate_jastrow,generate_jastrow1,generate_jastrow2,generate_jastrow3,generate_opt,generate_opts
from qmcpack import generate_cusp_correction
from auxiliary import *
from debug import *
#set the machine if known, otherwise user will provide
hostmachine = Machine.get_hostname()
if Machine.exists(hostmachine):
Job.machine = hostmachine
ProjectManager.machine = Machine.get(hostmachine)
#end if
class Settings(Pobj):
singleton = None
project_vars = set([
'machine','account','generate_only','verbose','debug','mode',
'local_directory','remote_directory','runs','results','sleep',
'file_locations','load_images','trace','machine_mode','stages',
'pseudo_dir','skip_submit','interactive_cores','monitor',
'status_only','machine_info',
])
gamess_vars = set('ericfmt mcppath'.split())
all_vars = project_vars | gamess_vars
@staticmethod
def kw_set(vars,source=None):
kw = obj()
if source!=None:
for n in vars:
if n in source:
kw[n]=source[n]
del source[n]
#end if
#end for
#end if
return kw
#end def null_kw_set
def __init__(self):
if Settings.singleton is None:
Settings.singleton = self
else:
self.error('attempted to create a second Settings object\n please just use the original')
#end if
#end def __init__
def error(self,message,header='settings',exit=True,trace=True):
Pobj.error(self,message,header,exit,trace)
#end def error
def __call__(self,**kwargs):
# guard against invalid settings
not_allowed = set(kwargs.keys()) - self.all_vars
if len(not_allowed)>0:
self.error('unrecognized variables provided.\n You provided: '+str(list(not_allowed))+'\n Allowed variables are: '+str(list(vars)))
#end if
# extract settings based on keyword groups
kw = Settings.kw_set(self.project_vars,kwargs) # project settings
gamess_kw = Settings.kw_set(self.gamess_vars,kwargs) # gamess settings
if len(kwargs)>0:
self.error('some settings keywords have not been accounted for\nleftover keywords: {0}\nthis is a developer error'.format(sorted(kwargs.keys())))
#end if
# transfer project variables to project base class
for name,value in kw.iteritems():
Pobj.__dict__[name] = value
#end for
# process project manager settings
if 'debug' in kw and kw.debug:
Pobj.verbose = True
#end if
if 'mode' in kw:
Pobj.set_mode(kw.mode)
#end if
# process machine settings
if 'machine_info' in kw:
machine_info = Pobj.machine_info
del Pobj.machine_info
if isinstance(machine_info,dict) or isinstance(machine_info,obj):
for machine_name,minfo in machine_info.iteritems():
mname = machine_name.lower()
if Machine.exists(mname):
machine = Machine.get(mname)
machine.incorporate_user_info(minfo)
else:
self.error('machine {0} is unknown\n cannot set machine_info'.format(machine_name))
#end if
#end for
else:
self.error('machine_info must be a dict or obj\n you provided type '+machine_info.__class__.__name__)
#end if
#end if
if 'machine' in kw:
machine_name = kw.machine
if not Machine.exists(machine_name):
Pobj.class_error('machine {0} is unknown'.format(machine_name))
#end if
Job.machine = machine_name
ProjectManager.machine = Machine.get(machine_name)
#del Pobj.machine
if 'account' in kw:
account = Pobj.account
#del Pobj.account
if not isinstance(account,str):
self.error('account for {0} must be a string\n you provided: {1}'.format(machine_name,account))
#end if
ProjectManager.machine.account = account
#end if
if 'machine_mode' in kw:
machine_mode = kw.machine_mode
if machine_mode in Machine.modes:
machine_mode = Machine.modes[machine_mode]
#end if
if machine_mode==Machine.modes.interactive:
if ProjectManager.machine==None:
ProjectManager.class_error('no machine specified for interactive mode')
#end if
if not isinstance(ProjectManager.machine,Supercomputer):
self.error('interactive mode is not supported for machine type '+ProjectManager.machine.__class__.__name__)
#end if
if not 'interactive_cores' in kw:
self.error('interactive mode requested, but interactive_cores not set')
#end if
ProjectManager.machine = ProjectManager.machine.interactive_representation(Pobj.interactive_cores)
del Pobj.interactive_cores
#end if
del Pobj.machine_mode
#end if
#end if
# process simulation settings
if 'local_directory' in kw:
Pobj.file_locations.append(kw.local_directory)
#end if
if 'skip_submit' in kw:
Pobj.skip_submission = Pobj.skip_submit
del Pobj.skip_submit
#end if
if 'file_locations' in kw:
fl = kw.file_locations
if isinstance(fl,str):
Pobj.file_locations.extend([fl])
else:
Pobj.file_locations.extend(list(fl))
#end if
#end if
if not 'pseudo_dir' in kw:
Pobj.pseudopotentials = Pseudopotentials()
else:
pseudo_dir = kw.pseudo_dir
Pobj.file_locations.append(pseudo_dir)
if not os.path.exists(pseudo_dir):
self.error('pseudo_dir "{0}" does not exist'.format(pseudo_dir),trace=False)
#end if
files = os.listdir(pseudo_dir)
ppfiles = []
for f in files:
pf = os.path.join(pseudo_dir,f)
if os.path.isfile(pf):
ppfiles.append(pf)
#end if
#end for
Pobj.pseudopotentials = Pseudopotentials(ppfiles)
#end if
# more simulation/project manager settings processing
mode = Pobj.mode
modes = Pobj.modes
if mode==modes.stages:
stages = Pobj.stages
elif mode==modes.all:
stages = list(Pobj.primary_modes)
else:
stages = [kw.mode]
#end if
allowed_stages = set(Pobj.primary_modes)
if isinstance(stages,str):
stages = [stages]
#end if
if len(stages)==0:
stages = list(Pobj.primary_modes)
#self.error('variable stages must be a list of primary modes.\n Options are '+str(list(allowed_stages)))
elif 'all' in stages:
stages = list(Pobj.primary_modes)
else:
forbidden = set(Pobj.stages)-allowed_stages
if len(forbidden)>0:
self.error('some stages provided are not primary stages.\n You provided '+str(list(forbidden))+'\n Options are '+str(list(allowed_stages)))
#end if
#end if
Pobj.mode = modes.stages
Pobj.stages = stages
Pobj.stages_set = set(Pobj.stages)
# process gamess settings
Gamess.settings(**gamess_kw)
return
#end def __call__
#end class Settings
settings = Settings()
def run_project(*args,**kwargs):
pm = ProjectManager()
pm.add_simulations(*args,**kwargs)
pm.run_project()
#end def run_project

161
nexus/library/project_base.py
View File

@ -1,161 +0,0 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# project_base.py #
# Provides base class functionality and access to 'global' data #
# for core Nexus classes. #
# #
# Content summary: #
# Pobj #
# Base class for core Nexus classes. #
# 'Global' data is assigned by the Settings class. #
# #
#====================================================================#
import os
import traceback
import gc as garbage_collector
from memory import resident
from generic import obj
from developer import DevBase
modes = obj(
none = 0,
setup = 1,
send_files = 2,
submit = 3,
get_output = 4,
analyze = 5,
stages = 6,
all = 7
)
class Pobj(DevBase):
gc = garbage_collector
gc.enable()
mode = modes.stages
generate_only = False
status_only = False
monitor = True
skip_submission = False
emulate = False #override application and use application_emulator
verbose = True
debug = False
trace = False
load_images = True
sleep = 3
local_directory = './'
remote_directory = local_directory
file_locations = [local_directory]
runs = 'runs'
results = 'results'
#pseudo_dir = os.path.join(local_directory,'pseudopotentials')
pseudo_dir = None
pseudopotentials = None
modes = modes
primary_modes = ['setup','send_files','submit','get_output','analyze']
dependent_modes = set(['submit'])
stages = []
stages_set = set(stages)
wrote_something = False # for pretty printing
@staticmethod
def set_mode(mode):
if mode in Pobj.modes:
Pobj.mode = Pobj.modes[mode]
else:
print 'settings Error: invalid mode specified: '+mode+'\n valid modes are '+str(Pobj.modes.keys())
#end if
#end def set_mode
def mem_usage(self):
return int(resident()/1e6)
#end def mem_usage
indent = ' '
def log(self,*texts,**kwargs):
if self.verbose:
if len(kwargs)>0:
n = kwargs['n']
else:
n=0
#end if
text=''
for t in texts:
text+=str(t)+' '
#end for
pad = n*self.indent
self.logfile.write(pad+text.replace('\n','\n'+pad)+'\n')
#end if
Pobj.wrote_something = True
#end def log
#@classmethod
#def class_error(cls,msg,source=None,n=0,trace=True):
# if source==None:
# source = cls.__name__
# #end if
# pad = n*cls.indent
# text=pad+source+' Error: '+msg
# text = '\n'+text.replace('\n','\n'+pad)+'\n\n'
# cls.logfile.write(text)
# if trace:
# traceback.print_stack()
# #end if
# exit()
##end def class_error
#
#@classmethod
#def class_warn(cls,msg,source=None,n=0):
# if source==None:
# source = cls.__name__
# #end if
# pad = n*cls.indent
# text=pad+source+' Warning: '+msg
# cls.logfile.write(text.replace('\n','\n'+pad)+'\n')
##end def class_warn
def dlog(self,*texts,**kwargs):
if self.debug:
#self.log('mem_usage',self.mem_usage(),n=5)
self.log(*texts,**kwargs)
#end if
#end def dlog
def tlog(self,*texts,**kwargs):
if self.trace:
self.log(*texts,**kwargs)
w,s,j,f,g,a=int(self.setup),int(self.submitted),int(self.job.finished),int(self.finished),int(self.got_output),int(self.analyzed)
self.log('w,s,j,f,g,a',w,s,j,f,g,a,n=kwargs['n']+1)
#self.log('dependencies',self.dependencies.keys(),n=kwargs['n']+1)
#self.log('dependents ',self.dependents.keys(),n=kwargs['n']+1)
#end if
#end def tlog
working_directory = None
def enter(self,directory,changedir=True,msg=''):
self.working_directory = os.getcwd()
self.log(' Entering '+directory,msg)
if changedir:
os.chdir(directory)
#end if
pad = ' '
return pad
#end def enter
def leave(self):
os.chdir(self.working_directory)
#end def leave
#end class Pobj

48
nexus/library/project_manager.py
View File

@ -20,7 +20,7 @@ import os
import time
import memory
from generic import obj
from project_base import Pobj
from nexus_base import NexusCore,nexus_core
from simulation import Simulation
@ -29,10 +29,9 @@ def trivial(sim,*args,**kwargs):
#end def trivial
class ProjectManager(Pobj):
class ProjectManager(NexusCore):
def __init__(self):
#variables determined by self
modes = self.modes
modes = nexus_core.modes
self.persistent_modes = set([modes.submit,modes.all])
self.simulations = obj()
self.cascades = obj()
@ -63,7 +62,7 @@ class ProjectManager(Pobj):
self.resolve_file_collisions()
self.propagate_blockages()
self.log('loading cascade images',n=1)
if self.load_images:
if nexus_core.load_images:
self.load_cascades()
else:
self.log('cascades',n=1)
@ -80,32 +79,29 @@ class ProjectManager(Pobj):
def run_project(self,status=False,status_only=False):
self.log('\nProject starting',n=0)
self.init_cascades()
status_only = status_only or self.status_only
status_only = status_only or nexus_core.status_only
status = status or status_only
if status:
self.write_simulation_status()
if status_only:
NexusCore.write_end_splash()
return
#end if
#end if
self.log('\nstarting runs:\n'+30*'~',n=1)
if self.dependent_modes <= self.stages_set:
if self.monitor:
if nexus_core.dependent_modes <= nexus_core.stages_set:
if nexus_core.monitor:
ipoll = 0
while len(self.progressing_cascades)>0:
self.dlog('\n\n\n'+70*'=',n=1)
self.log('poll',ipoll,' memory %3.2f MB'%(memory.resident()/1e6),n=1)
Pobj.wrote_something = False
self.dlog('cascades',self.progressing_cascades.keys(),n=2)
NexusCore.wrote_something = False
ipoll+=1
self.machine.query_queue()
self.progress_cascades()
self.machine.submit_jobs()
self.update_process_ids()
self.dlog('sleeping',self.sleep,n=2)
time.sleep(self.sleep)
self.dlog('awake',n=2)
if Pobj.wrote_something:
time.sleep(nexus_core.sleep)
if NexusCore.wrote_something:
self.log()
#end if
#end while
@ -119,23 +115,20 @@ class ProjectManager(Pobj):
self.progress_cascades()
#end if
self.log('Project finished\n')
NexusCore.write_end_splash()
#end def run_project
def load_cascades(self):
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~load cascades',n=1)
cascades = obj()
progressing_cascades = obj()
for cid,cascade in self.cascades.iteritems():
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~reconstruct cascade',n=1)
rc = cascade.reconstruct_cascade()
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end reconstruct cascade',n=1)
cascades[rc.simid] = rc
progressing_cascades[rc.simid] = rc
#end for
self.cascades = cascades
self.progressing_cascades = progressing_cascades
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end load cascades',n=1)
#end def load_cascades
@ -204,7 +197,7 @@ class ProjectManager(Pobj):
def write_simulation_status(self):
self.log('cascade status',n=1)
self.log('\ncascade status',n=1)
self.log('setup, sent_files, submitted, finished, got_output, analyzed',n=2)
indicators = ('setup','sent_files','submitted','finished','got_output','analyzed')
for isim in self.simulations.keys():
@ -214,14 +207,14 @@ class ProjectManager(Pobj):
for stat in stats:
status+=str(int(stat))
#end for
self.log('{0} {1} {2}'.format(status,sim.identifier,sim.locdir),n=2)
self.log('{0} {1:<6} {2}'.format(status,sim.identifier,sim.locdir),n=2)
#end for
self.log('setup, sent_files, submitted, finished, got_output, analyzed',n=2)
#end def write_simulation_status
def write_cascade_status(self):
self.log('cascade status',n=1)
self.log('\ncascade status',n=1)
self.log('setup, sent_files, submitted, finished, got_output, analyzed',n=2)
def write_status(sim):
@ -319,31 +312,22 @@ class ProjectManager(Pobj):
#end def check_dependencies
def progress_cascades(self):
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~progress cascades',n=1)
self.gc.collect()
NexusCore.gc.collect()
finished = []
progressing_cascades = self.progressing_cascades
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~reset wait ids',n=1)
for cid,cascade in progressing_cascades.iteritems():
cascade.reset_wait_ids()
#end for
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end reset wait ids',n=1)
for cid,cascade in progressing_cascades.iteritems():
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~progress',n=1)
cascade.progress()
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end progress',n=1)
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~check subcascade',n=1)
cascade.check_subcascade()
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end check subcascade',n=1)
if cascade.subcascade_finished:
finished.append(cid)
#end if
#end for
for cid in finished:
self.dlog('removing cascade:',cid,n=1)
del progressing_cascades[cid]
#end for
self.dlog('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~end progress cascades',n=1)
#end def progress_cascades

894
nexus/library/pseudopotential.py
View File

File diff suppressed because it is too large Load Diff

9
nexus/library/pw2casino.py
View File

@ -30,6 +30,7 @@
import os
from generic import obj
from unit_converter import convert
from simulation import Simulation,SimulationInput,SimulationAnalyzer
@ -88,7 +89,7 @@ class Pw2casinoInput(SimulationInput):
allowed = set(ints+floats+strs+bools)
def read_contents(self,contents):
def read_text(self,contents,filepath=None):
lines = contents.split('\n')
inside = False
for l in lines:
@ -120,9 +121,9 @@ class Pw2casinoInput(SimulationInput):
inside=False
#end if
#end for
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self,filepath=None):
contents = ''
for sname,section in self.iteritems():
contents+='&'+sname+'\n'
@ -133,7 +134,7 @@ class Pw2casinoInput(SimulationInput):
contents+='/\n'
#end for
return contents
#end def write_contents
#end def write_text
def __init__(self,filepath=None,**vars):

360
nexus/library/pw2qmcpack.py
View File

@ -1,360 +0,0 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# pw2qmcpack.py #
# Nexus interface for the pw2qmcpack tool distributed with PWSCF. #
# pseudopotentials between file formats. #
# #
# Content summary: #
# generate_pw2qmcpack #
# User-facing function to create pw2qmcpack simulation objects. #
# #
# generate_pw2qmcpack_input #
# User-facing funcion to create input for pw2qmcpack. #
# #
# Pw2qmcpack #
# Simulation class for pw2qmcpack. #
# #
# Pw2qmcpackInput #
# SimulationInput class for pw2qmcpack. #
# #
# Pw2qmcpackAnalyzer #
# SimulationAnalyzer class for pw2qmcpack. #
# #
#====================================================================#
import os
from generic import obj
from simulation import Simulation,SimulationInput,SimulationAnalyzer
def read_str(sv):
return sv.strip('"').strip("'")
#end def read_str
def read_int(sv):
return int(sv)
#end def read_int
def read_float(sv):
return float(sv.replace('d','e').replace('D','e'))
#end def read_float
bconv = {'.true.':True,'.false.':False}
def read_bool(sv):
return bconv[sv]
#end def read_bool
def write_str(val):
return "'"+val+"'"
#end def write_str
def write_int(val):
return str(val)
#end def write_int
def write_float(val):
return str(val)
#end def write_float
def write_bool(val):
return '.'+str(val).lower()+'.'
#end def write_bool
readval={str:read_str,int:read_int,float:read_float,bool:read_bool}
writeval={str:write_str,int:write_int,float:write_float,bool:write_bool}
class Pw2qmcpackInput(SimulationInput):
ints = []
floats = []
strs = ['outdir','prefix']
bools = ['write_psir']
var_types = dict()
for v in ints:
var_types[v]=int
for v in floats:
var_types[v]=float
for v in strs:
var_types[v]=str
for v in bools:
var_types[v]=bool
allowed = set(ints+floats+strs+bools)
def read_contents(self,contents):
lines = contents.split('\n')
inside = False
for l in lines:
if inside:
tokens = l.split(',')
for t in tokens:
ts = t.strip()
if ts!='' and ts!='/':
name,value = ts.split('=')
name = name.strip()
value= value.strip()
if name in self.allowed:
vtype = self.var_types[name]
value = readval[vtype](value)
sobj[name]=value
else:
self.error('encountered unknown variable: '+name)
#end if
#end if
#end for
#end if
if '&' in l:
inside=True
section = l.lstrip('&').lower()
sobj = obj()
self[section]=sobj
elif l.strip()=='/':
inside=False
#end if
#end for
#end def read_contents
def write_contents(self):
contents = ''
for sname,section in self.iteritems():
contents+='&'+sname+'\n'
for name,value in section.iteritems():
vtype = type(value)
contents += ' '+name+' = '+writeval[vtype](value)+'\n'
#end for
contents+='/\n'
#end for
return contents
#end def write_contents
def __init__(self,filepath=None,**vars):
if filepath!=None:
self.read(filepath)
else:
inputpp = obj()
for name,value in vars.iteritems():
inputpp[name] = value
#end for
self.inputpp = inputpp
#end if
#end def __init__
#end class Pw2qmcpackInput
def generate_pw2qmcpack_input(prefix='pwscf',outdir='pwscf_output',write_psir=True):
pw = Pw2qmcpackInput(
prefix = prefix,
outdir = outdir,
write_psir = write_psir
)
return pw
#end def generate_pw2qmcpack_input
class Pw2qmcpackAnalyzer(SimulationAnalyzer):
def __init__(self,arg0):
if isinstance(arg0,Simulation):
sim = arg0
self.infile = sim.infile
prefix,outdir = sim.input.inputpp.tuple('prefix','outdir')
self.dir = sim.locdir
self.h5file = os.path.join(sim.locdir,outdir,prefix+'.pwscf.h5')
else:
self.infile = arg0
#end if
#end def __init__
def analyze(self):
if False:
import h5py
self.log('Fixing h5 file',n=5)
path = os.path.split(self.h5file)[0]
print os.getcwd()
print os.listdir('./')
if os.path.exists(path):
print os.listdir(path)
#end if
print self.h5file
h = h5py.File(self.h5file)
if 'electrons' in h:
elec = h['electrons']
nkpoints = 0
for name,val in elec.iteritems():
if name.startswith('kpoint'):
nkpoints+=1
#end for
#end if
nkold = elec['number_of_kpoints'][0]
self.log('Were',nkold,'kpoints, now',nkpoints,'kpoints',n=6)
elec['number_of_kpoints'][0] = nkpoints
#end for
#end if
#end def analyze
def get_result(self,result_name):
self.not_implemented()
#end def get_result
#end class Pw2qmcpackAnalyzer
class Pw2qmcpack(Simulation):
input_type = Pw2qmcpackInput
analyzer_type = Pw2qmcpackAnalyzer
generic_identifier = 'pw2qmcpack'
application = 'pw2qmcpack.x'
application_properties = set(['serial'])
application_results = set(['orbitals'])
def check_result(self,result_name,sim):
calculating_result = False
inputpp = self.input.inputpp
if result_name=='orbitals':
calculating_result = True
else:
calculating_result = False
self.error('ability to check for result '+result_name+' has not been implemented')
#end if
return calculating_result
#end def check_result
def get_result(self,result_name,sim):
result = obj()
inputpp = self.input.inputpp
prefix = 'pwscf'
outdir = './'
if 'prefix' in inputpp:
prefix = inputpp.prefix
#end if
if 'outdir' in inputpp:
outdir = inputpp.outdir
#end if
if outdir.startswith('./'):
outdir = outdir[2:]
#end if
if result_name=='orbitals':
result.h5file = os.path.join(self.locdir,outdir,prefix+'.pwscf.h5')
result.ptcl_xml = os.path.join(self.locdir,outdir,prefix+'.ptcl.xml')
result.wfs_xml = os.path.join(self.locdir,outdir,prefix+'.wfs.xml')
else:
self.error('ability to get result '+result_name+' has not been implemented')
#end if
return result
#end def get_result
def incorporate_result(self,result_name,result,sim):
if result_name=='orbitals':
pwin = sim.input.control
p2in = self.input.inputpp
pwprefix = 'pwscf'
p2prefix = 'pwscf'
pwoutdir = './'
p2outdir = './'
if 'prefix' in pwin:
pwprefix = pwin.prefix
if 'prefix' in p2in:
p2prefix = p2in.prefix
if 'outdir' in pwin:
pwoutdir = pwin.outdir
if 'outdir' in p2in:
p2outdir = p2in.outdir
if pwoutdir.startswith('./'):
pwoutdir = pwoutdir[2:]
if p2outdir.startswith('./'):
p2outdir = p2outdir[2:]
pwdir = os.path.abspath(os.path.join(sim.locdir ,pwoutdir))
p2dir = os.path.abspath(os.path.join(self.locdir,p2outdir))
errors = False
if pwdir!=p2dir:
self.error('to use orbitals, '+self.generic_identifier+' must have the same outdir as pwscf\n pwscf outdir: '+pwdir+'\n '+self.generic_identifier+' outdir: '+p2dir,exit=False)
errors = True
#end if
if pwprefix!=p2prefix:
self.error('to use orbitals, '+self.generic_identifier+' must have the same prefix as pwscf\n pwscf prefix: '+pwprefix+'\n '+self.generic_identifier+' prefix: '+p2prefix,exit=False)
errors = True
#end if
if errors:
self.error(self.generic_identifier+' cannot use orbitals from pwscf')
#end if
else:
self.error('ability to incorporate result '+result_name+' has not been implemented')
#end if
#end def incorporate_result
def check_sim_status(self):
outfile = os.path.join(self.locdir,self.outfile)
fobj = open(outfile,'r')
output = fobj.read()
fobj.close()
inputpp = self.input.inputpp
prefix = 'pwscf'
outdir = './'
if 'prefix' in inputpp:
prefix = inputpp.prefix
#end if
if 'outdir' in inputpp:
outdir = inputpp.outdir
#end if
if outdir.startswith('./'):
outdir = outdir[2:]
#end if
h5file = os.path.join(self.locdir,outdir,prefix+'.pwscf.h5')
ptcl_xml = os.path.join(self.locdir,outdir,prefix+'.ptcl.xml')
wfs_xml = os.path.join(self.locdir,outdir,prefix+'.wfs.xml')
must_exist = [h5file,ptcl_xml,wfs_xml]
files_exist = True
for file in must_exist:
files_exist = files_exist and os.path.exists(file)
#end for
outfin = True
#outfin = outfin and 'esh5 create' in output
#outfin = outfin and 'Creating electrons' in output
outfin = outfin and 'npw=' in output
outfin = outfin and 'ik=' in output
outfin = outfin or 'JOB DONE' in output
success = files_exist and outfin
#self.finished = success and self.job.finished
# pw2qmcpack has too many variants to assess completion based on log output
# assume (optimistically) that job completion indicates success
self.finished = files_exist and self.job.finished
#end def check_sim_status
def get_output_files(self):
output_files = []
return output_files
#end def get_output_files
def app_command(self):
return self.app_name+'<'+self.infile
#end def app_command
#end class Pw2qmcpack
def generate_pw2qmcpack(**kwargs):
sim_args,inp_args = Simulation.separate_inputs(kwargs)
if not 'input' in sim_args:
sim_args.input = generate_pw2qmcpack_input(**inp_args)
#end if
pw2qmcpack = Pw2qmcpack(**sim_args)
return pw2qmcpack
#end def generate_pw2qmcpack

2
nexus/library/pwscf.py
View File

@ -21,7 +21,7 @@ import os
from numpy import array
from generic import obj
from physical_system import PhysicalSystem
from converters import Pw2qmcpack
from qmcpack_converters import Pw2qmcpack
from simulation import Simulation
from pwscf_input import PwscfInput,generate_pwscf_input
from pwscf_analyzer import PwscfAnalyzer

146
nexus/library/pwscf_input.py
View File

@ -45,7 +45,7 @@
import os
import types as pytypes
import inspect
from copy import deepcopy
from superstring import string2val
from numpy import fromstring,empty,array,float64,ones,pi,dot
@ -54,7 +54,7 @@ from unit_converter import convert
from generic import obj
from structure import Structure,kmesh
from physical_system import PhysicalSystem
from project_base import Pobj
from developer import DevBase,warn
from simulation import SimulationInput
from debug import *
@ -153,7 +153,7 @@ def array_to_string(a,pad=' ',format=pwscf_array_format,converter=noconv,rowse
class PwscfInputBase(Pobj):
class PwscfInputBase(DevBase):
ints=['nstep','iprint','gdir','nppstr','nberrycyc','ibrav','nat','ntyp','nbnd','tot_charge','nr1','nr2','nr3','nr1s','nr2s','nr3s','nspin','multiplicity','tot_magnetization','edir','report','electron_maxstep','mixing_ndim','mixing_fixed_ns','ortho_para','diago_cg_maxiter','diago_david_ndim','nraise','bfgs_ndim','num_of_images','fe_nstep','sw_nstep','modenum','n_charge_compensation','nlev','lda_plus_u_kind']
floats=['dt','max_seconds','etot_conv_thr','forc_conv_thr','celldm','A','B','C','cosAB','cosAC','cosBC','nelec','ecutwfc','ecutrho','degauss','starting_magnetization','nelup','neldw','ecfixed','qcutz','q2sigma','Hubbard_alpha','Hubbard_U','Hubbard_J','starting_ns_eigenvalue','emaxpos','eopreg','eamp','angle1','angle2','fixed_magnetization','lambda','london_s6','london_rcut','conv_thr','mixing_beta','diago_thr_init','efield','tempw','tolp','delta_t','upscale','trust_radius_max','trust_radius_min','trust_radius_ini','w_1','w_2','temp_req','ds','k_max','k_min','path_thr','fe_step','g_amplitude','press','wmass','cell_factor','press_conv_thr','xqq','ecutcoarse','mixing_charge_compensation','comp_thr','exx_fraction','ecutfock']
strs=['calculation','title','verbosity','restart_mode','outdir','wfcdir','prefix','disk_io','pseudo_dir','occupations','smearing','input_dft','U_projection_type','constrained_magnetization','mixing_mode','diagonalization','startingpot','startingwfc','ion_dynamics','ion_positions','phase_space','pot_extrapolation','wfc_extrapolation','ion_temperature','opt_scheme','CI_scheme','cell_dynamics','cell_dofree','which_compensation','assume_isolated','exxdiv_treatment']
@ -285,22 +285,22 @@ class Card(Element):
def read(self,lines):
self.get_specifier(lines[0])
self.read_contents(lines[1:])
self.read_text(lines[1:])
#end def read
def write(self,parent):
c = self.name.upper()+' '+self.specifier+'\n'
c += self.write_contents()+'\n'
c += self.write_text()+'\n'
return c
#end def write
def read_contents(self,lines):
def read_text(self,lines):
self.not_implemented()
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self):
self.not_implemented()
#end def write_contents
#end def write_text
def change_specifier(self,new_specifier):
self.not_implemented()
@ -513,7 +513,7 @@ def check_section_classes(*sections):
class atomic_species(Card):
name = 'atomic_species'
def read_contents(self,lines):
def read_text(self,lines):
atoms = []
masses = obj()
pseudopotentials = obj()
@ -525,22 +525,22 @@ class atomic_species(Card):
pseudopotentials[atom] = tokens[2]
#end for
self.add(atoms=atoms,masses=masses,pseudopotentials=pseudopotentials)
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self):
c = ''
for at in self.atoms:
c += ' '+'{0:2}'.format(at)+' '+str(self.masses[at])+' '+self.pseudopotentials[at]+'\n'
#end for
return c
#end def write_contents
#end def write_text
#end class atomic_species
class atomic_positions(Card):
name = 'atomic_positions'
def read_contents(self,lines):
def read_text(self,lines):
npos = len(lines)
dim = 3
atoms = []
@ -562,10 +562,10 @@ class atomic_positions(Card):
if has_relax_directions:
self.add(relax_directions=relax_directions)
#end if
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self):
c = ''
has_relax_directions = 'relax_directions' in self
if has_relax_directions:
@ -581,7 +581,7 @@ class atomic_positions(Card):
#end if
#end for
return c
#end def write_contents
#end def write_text
def change_specifier(self,new_specifier,pwi):
@ -627,7 +627,7 @@ class atomic_positions(Card):
class k_points(Card):
name = 'k_points'
def read_contents(self,lines):
def read_text(self,lines):
if self.specifier in ['tpiba','crystal','tpiba_b','crystal_b','']:
self.nkpoints = int(lines[0])
a = array_from_lines(lines[1:])
@ -643,10 +643,10 @@ class k_points(Card):
else:
self.error('k_points specifier '+self.specifier+' is unrecognized')
#end if
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self):
c = ''
if self.specifier in ('tpiba','crystal','tpiba_b','crystal_b',''):
self.nkpoints = len(self.kpoints)
@ -665,7 +665,7 @@ class k_points(Card):
self.error('k_points specifier '+self.specifier+' is unrecognized')
#end if
return c
#end def write_contents
#end def write_text
def change_specifier(self,new_specifier,pwi):
@ -716,13 +716,13 @@ class k_points(Card):
class cell_parameters(Card):
name = 'cell_parameters'
def read_contents(self,lines):
def read_text(self,lines):
self.vectors = array_from_lines(lines)
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self):
return array_to_string(self.vectors)
#end def write_contents
#end def write_text
#end class cell_parameters
@ -731,17 +731,17 @@ class cell_parameters(Card):
class climbing_images(Card):
name = 'climbing_images'
def read_contents(self,lines):
def read_text(self,lines):
self.images = array_from_lines(lines)
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self):
c=' '
for n in self.images:
c+=str(int(n))+' '
#end for
return c
#end def write_contents
#end def write_text
#end class climbing_images
@ -750,7 +750,7 @@ class climbing_images(Card):
class constraints(Card):
name = 'constraints'
def read_contents(self,lines):
def read_text(self,lines):
tokens = lines[0].split()
self.ncontraints = int(tokens[0])
if len(tokens)>1:
@ -764,9 +764,9 @@ class constraints(Card):
cons.parameters = array(tokens[1:],dtype=float64)
self.constraints[i] = cons
#end for
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self):
c = ' '+str(self.nconstraints)
if 'tolerance' in self:
c+=' '+str(self.tolerance)
@ -775,7 +775,7 @@ class constraints(Card):
c+=' '+cons.type+' '+array_to_string(cons.parameters,pad='')
#end for
return c
#end def write_contents
#end def write_text
#end class constraints
@ -784,7 +784,7 @@ class constraints(Card):
class collective_vars(Card):
name = 'collective_vars'
def read_contents(self,lines):
def read_text(self,lines):
tokens = lines[0].split()
self.ncontraints = int(tokens[0])
if len(tokens)>1:
@ -798,9 +798,9 @@ class collective_vars(Card):
collv.parameters = array(tokens[1:],dtype=float64)
self.collective_vars[i] = collv
#end for
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self):
c= ' '+str(self.ncollective_vars)
if 'tolerance' in self:
c+=' '+str(self.tolerance)
@ -809,7 +809,7 @@ class collective_vars(Card):
c+=' '+collv.type+' '+array_to_string(collv.parameters,pad='')
#end for
return c
#end def write_contents
#end def write_text
#end class collective_vars
@ -818,13 +818,13 @@ class collective_vars(Card):
class occupations(Card):
name = 'occupations'
def read_contents(self,lines):
def read_text(self,lines):
self.occupations = array_from_lines(lines)
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self):
return array_to_string(self.occupations)
#end def write_contents
#end def write_text
#end class occupations
@ -893,7 +893,7 @@ class PwscfInput(SimulationInput):
#end def __init__
def read_contents(self,contents):
def read_text(self,contents,filepath=None):
lines = contents.splitlines()
in_element = False
elem_type = None
@ -948,10 +948,10 @@ class PwscfInput(SimulationInput):
for element in self:
element.post_process_read(self)
#end for
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self,filepath=None):
contents = ''
for s in self.sections:
if s in self:
@ -966,7 +966,7 @@ class PwscfInput(SimulationInput):
#end for
contents+='\n'
return contents
#end def write_contents
#end def write_text
def get_common_vars(self,*vars):
@ -1225,22 +1225,6 @@ class PwscfInput(SimulationInput):
ion_charge += atoms.count(atom)*valency[atom]
#end for
#ion_charge = 0
#valency = dict()
#atoms = list(self.atomic_positions.atoms)
#for atom in self.atomic_species.atoms:
# pseudo_file = self.atomic_species.pseudopotentials[atom]
# if self.pseudopotentials!=None and pseudo_file in self.pseudopotentials:
# pseudopot = self.pseudopotentials[pseudo_file]
# element = pseudopot.element
# valence = int(pseudopot.Z)
# ion_charge += atoms.count(atom)*valence
# valency[atom] = valence
# else:
# self.error('file '+pseudo_file+' was not listed in Pseudopotentials object\n please specify pseudopotentials with the settings function',trace=False)
# #end if
##end for
if 'nelup' in self.system:
nup = self.system.nelup
ndn = self.system.neldw
@ -1392,7 +1376,7 @@ def generate_any_pwscf_input(**kwargs):
for name,default in defaults.iteritems():
if not name in kwargs:
deftype = type(default)
if deftype==pytypes.ClassType or deftype==pytypes.FunctionType:
if inspect.isclass(default) or inspect.isfunction(default):
kwargs[name] = default()
else:
kwargs[name] = default
@ -1427,13 +1411,14 @@ def generate_any_pwscf_input(**kwargs):
#end for
#process other keywords
pseudos = kwargs.delete('pseudos')
system = kwargs.delete('system')
use_folded = kwargs.delete('use_folded')
hubbard_u = kwargs.delete('hubbard_u')
start_mag = kwargs.delete('start_mag')
kgrid = kwargs.delete('kgrid')
kshift = kwargs.delete('kshift')
pseudos = kwargs.delete_required('pseudos')
system = kwargs.delete_required('system')
use_folded = kwargs.delete_required('use_folded')
hubbard_u = kwargs.delete_required('hubbard_u')
start_mag = kwargs.delete_required('start_mag')
kgrid = kwargs.delete_required('kgrid')
kshift = kwargs.delete_required('kshift')
nogamma = kwargs.delete_optional('nogamma',False)
# pseudopotentials
pseudopotentials = obj()
@ -1507,17 +1492,32 @@ def generate_any_pwscf_input(**kwargs):
else:
zero_shift = tuple(kshift)==(0,0,0)
#end if
if zero_shift and (kgrid!=None and tuple(kgrid)==(1,1,1) or kgrid is None and system is None):
single_point = kgrid != None and tuple(kgrid)==(1,1,1)
no_info = kgrid is None and system is None
at_gamma = zero_shift and (single_point or no_info)
at_gamma |= 'specifier' in pw.k_points and pw.k_points.specifier=='gamma'
auto = kgrid!=None and kshift!=None
shifted = not zero_shift and kshift!=None
if at_gamma and not nogamma:
pw.k_points.clear()
pw.k_points.specifier = 'gamma'
elif kgrid!=None and kshift!=None:
elif at_gamma and nogamma:
pw.k_points.clear()
pw.k_points.set(
specifier = 'automatic',
grid = (1,1,1),
shift = (0,0,0)
)
elif auto:
pw.k_points.clear()
pw.k_points.set(
specifier = 'automatic',
grid = kgrid,
shift = kshift
)
elif not zero_shift and kshift!=None:
elif shifted:
pw.k_points.clear()
pw.k_points.set(
specifier = 'automatic',
@ -1897,7 +1897,7 @@ def generate_relax_input(prefix = 'pwscf',
if system is not None:
structure = system.structure
if group_atoms:
self.warn('requested grouping by atomic species, but pwscf does not group atoms anymore!')
warn('requested grouping by atomic species, but pwscf does not group atoms anymore!','generate_relax_input')
#end if
#if group_atoms: #don't group atoms for pwscf, any downstream consequences?
# structure.group_atoms()

149
nexus/library/qaobject.py
View File

@ -1,149 +0,0 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# qaobject.py #
# Provides base class functionality for all classes associated #
# with the QmcpackAnalyzer class. Maintains data global to #
# these classes. #
# #
# Content summary: #
# QAobject #
# Base class for all QMCPACK analyzer components. #
# Exposes settings options to the user. #
# #
# Checks #
# Class to assess overall validity based on stored results of #
# many checks (boolean values). Only use so far is to validate #
# the structure of Trace files. See qmcpack_method_analyzers.py.#
# #
#====================================================================#
from generic import obj
from developer import DevBase
from project_base import Pobj
from plotter import Plotter
class QAobj_base(Pobj):
None
#end class QAobj_base
class QAobject(QAobj_base):
_global = obj()
_global.dynamic_methods_objects=[]
plotter = Plotter()
opt_methods = set(['opt','linear','cslinear'])
def __init__(self):
return
#end def __init__
@staticmethod
def condense_name(name):
return name.strip().lower().replace(' ','_').replace('-','_').replace('__','_')
#end def condense_name
def _register_dynamic_methods(self):
QAobject._global.dynamic_methods_objects.append(self)
return
#end def _register_dynamic_methods
def _unlink_dynamic_methods(self):
for o in QAobject._global.dynamic_methods_objects:
o._unset_dynamic_methods()
#end for
return
#end def _unlink_dynamic_methods
def _relink_dynamic_methods(self):
for o in QAobject._global.dynamic_methods_objects:
o._reset_dynamic_methods()
#end for
return
#end def _relink_dynamic_methods
_allowed_settings = set(['optimize'])
_default_settings = obj(
#optimize = 'variance'
optimize = 'lastcost'
)
for var,val in _default_settings.iteritems():
QAobj_base.__dict__[var] = val
#end for
@classmethod
def settings(cls,**kwargs):
vars = set(kwargs.keys())
invalid = vars-cls._allowed_settings
if len(invalid)>0:
allowed = list(cls._allowed_settings)
allowed.sort()
invalid = list(invalid)
invalid.sort()
cls.class_error('attempted to set unknown variables\n unknown variables: {0}\n valid options are: {1}'.format(invalid,allowed))
#end if
for var,val in kwargs.iteritems():
QAobj_base.__dict__[var] = val
#end for
#end settings
#end class QAobject
class Checks(Pobj):
def __init__(self,label=''):
self._label = label
self._exclusions = set()
#end def __init__
def exclude(self,value):
self._exclusions.add(value)
#end def exclude
def valid(self):
valid = True
for name,value in self.iteritems():
if not (isinstance(name,str) and name.startswith('_')):
if not value in self._exclusions:
valid = valid and value
#end if
#end if
#end if
self._valid = valid
return valid
#end def valid
def write(self,pad=''):
pad2 = pad+' '
if not '_valid' in self:
self.valid()
#end if
valid = self._valid
if valid:
self.log(pad+self._label+' is valid')
else:
self.log(pad+self._label+' is invalid')
for name,value in self.iteritems():
if not (isinstance(name,str) and name.startswith('_')):
if value in self._exclusions:
self.log(pad2+name+' could not be checked')
elif value:
self.log(pad2+name+' is valid')
else:
self.log(pad2+name+' is invalid')
#end if
#end if
#end for
#end if
#end def write
#end class Checks

240
nexus/library/qmcpack.py
View File

@ -32,14 +32,12 @@ from periodic_table import periodic_table
from physical_system import PhysicalSystem
from simulation import Simulation
from qmcpack_input import QmcpackInput,generate_qmcpack_input
from qmcpack_input import BundledQmcpackInput,TracedQmcpackInput
from qmcpack_input import QmcpackInputTemplate
from qmcpack_input import TracedQmcpackInput
from qmcpack_input import loop,linear,cslinear,vmc,dmc,collection,determinantset,hamiltonian,init,pairpot,bspline_builder
from qmcpack_input import generate_jastrows,generate_jastrow,generate_jastrow1,generate_jastrow2,generate_jastrow3
from qmcpack_input import generate_opt,generate_opts
from qmcpack_analyzer import QmcpackAnalyzer
from converters import Pw2qmcpack,Wfconvert
from convert4qmc import Convert4qmc
from qmcpack_converters import Pw2qmcpack,Wfconvert,Convert4qmc
from sqd import Sqd
from debug import ci,ls,gs
from developer import unavailable
@ -60,7 +58,6 @@ class Qmcpack(Simulation):
application = 'qmcapp_complex' # always use complex version until kpoint handling is fixed
application_properties = set(['serial','omp','mpi'])
application_results = set(['jastrow','cuspcorr','wavefunction'])
preserve = Simulation.preserve | set(['should_twist_average'])
def post_init(self):
@ -71,12 +68,19 @@ class Qmcpack(Simulation):
#if self.system is None:
# self.error('system must be specified to determine type of run')
##end if
generic_input = self.has_generic_input()
if self.system is None:
if not isinstance(self.input,QmcpackInputTemplate):
if not generic_input:
self.warn('system must be specified to determine whether to twist average\n proceeding under the assumption of no twist averaging')
#end if
self.should_twist_average = False
else:
if generic_input:
cls = self.__class__
self.error('cannot twist average generic or templated input\nplease provide {0} instead of {1} for input'.format(cls.input_type.__class__.__name__,self.input.__class__.__name__))
#end if
self.system.group_atoms()
self.system.change_units('B')
twh = self.input.get_host('twist')
@ -91,7 +95,7 @@ class Qmcpack(Simulation):
def propagate_identifier(self):
if not isinstance(self.input,QmcpackInputTemplate):
if not self.has_generic_input():
self.input.simulation.project.id = self.identifier
#end if
#end def propagate_identifier
@ -356,45 +360,43 @@ class Qmcpack(Simulation):
def check_sim_status(self):
outfile = os.path.join(self.locdir,self.outfile)
errfile = os.path.join(self.locdir,self.errfile)
fobj = open(outfile,'r')
output = fobj.read()
fobj.close()
fobj = open(errfile,'r')
errors = fobj.read()
fobj.close()
cusp_run = self.input.cusp_correction()
if cusp_run:
sd = self.input.get('slaterdeterminant')
if sd!=None:
cuspfiles = []
for d in sd.determinants:
cuspfiles.append(d.id+'.cuspInfo.xml')
#end for
else: # assume multideterminant sposet names
cuspfiles = ['spo-up.cuspInfo.xml','spo-dn.cuspInfo.xml']
#end if
outfiles = cuspfiles
else:
outfiles = self.input.get_output_info('outfiles')
#end if
output = self.outfile_text()
errors = self.errfile_text()
ran_to_end = 'Total Execution' in output
aborted = 'Fatal Error' in errors
files_exist = True
for file in outfiles:
file_loc = os.path.join(self.locdir,file)
files_exist = files_exist and os.path.exists(file_loc)
#end for
if ran_to_end and not files_exist:
self.warn('run finished successfully, but output files do not seem to exist')
print outfiles
print os.listdir(self.locdir)
cusp_run = False
if not self.has_generic_input():
cusp_run = self.input.cusp_correction()
if cusp_run:
sd = self.input.get('slaterdeterminant')
if sd!=None:
cuspfiles = []
for d in sd.determinants:
cuspfiles.append(d.id+'.cuspInfo.xml')
#end for
else: # assume multideterminant sposet names
cuspfiles = ['spo-up.cuspInfo.xml','spo-dn.cuspInfo.xml']
#end if
outfiles = cuspfiles
else:
outfiles = self.input.get_output_info('outfiles')
#end if
for file in outfiles:
file_loc = os.path.join(self.locdir,file)
files_exist = files_exist and os.path.exists(file_loc)
#end for
if ran_to_end and not files_exist:
self.warn('run finished successfully, but output files do not seem to exist')
print outfiles
print os.listdir(self.locdir)
#end if
#end if
aborted = 'Fatal Error' in errors
self.succeeded = ran_to_end
self.failed = aborted
@ -408,30 +410,22 @@ class Qmcpack(Simulation):
#end for
#end if
#print
#print self.__class__.__name__
#print 'identifier ',self.identifier
#print 'ran_to_end ',ran_to_end
#print 'files_exist',files_exist
#print 'aborted ',aborted
#print 'job done ',self.job.finished
#print 'finished ',self.finished
#print
#end def check_sim_status
def get_output_files(self):
if self.should_twist_average and not isinstance(self.input,TracedQmcpackInput):
self.twist_average(range(len(self.system.structure.kpoints)))
br = self.bundle_request
input = self.input.trace(br.quantity,br.values)
input.generate_filenames(self.infile)
self.input = input
if self.has_generic_input():
output_files = []
else:
if self.should_twist_average and not isinstance(self.input,TracedQmcpackInput):
self.twist_average(range(len(self.system.structure.kpoints)))
br = self.bundle_request
input = self.input.trace(br.quantity,br.values)
input.generate_filenames(self.infile)
self.input = input
#end if
output_files = self.input.get_output_info('outfiles')
#end if
output_files = self.input.get_output_info('outfiles')
return output_files
#end def get_output_files
@ -446,16 +440,14 @@ class Qmcpack(Simulation):
br.quantity = 'twistnum'
br.values = list(twistnums)
self.bundle_request = br
#self.app_name = 'qmcapp_complex'
#print 'twist_average'
#print ' setting bundle request:'
#print self.bundle_request
#end def twist_average
def write_prep(self):
if self.got_dependencies:
if 'bundle_request' in self and not isinstance(self.input,TracedQmcpackInput):
traced_input = isinstance(self.input,TracedQmcpackInput)
generic_input = self.has_generic_input()
if 'bundle_request' in self and not traced_input and not generic_input:
br = self.bundle_request
input = self.input.trace(br.quantity,br.values)
input.generate_filenames(self.infile)
@ -475,124 +467,6 @@ class Qmcpack(Simulation):
class BundledQmcpack(Qmcpack):
infile_extension = '.in'
application_results = set([])
preserve = set(Simulation.preserve)
preserve.add('sims')
def __init__(self,**kwargs):
if not 'sims' in kwargs:
self.error('sims must be provided')
#end if
sims = kwargs['sims']
self.sims = sims
del kwargs['sims']
files = set()
for sim in sims:
files = files | sim.files
#end for
kwargs['files'] = files
inputs = []
filenames = []
for sim in sims:
inputs.append(sim.input)
filenames.append(sim.infile)
#end for
kwargs['input'] = BundledQmcpackInput(inputs=inputs,filenames=filenames)
Simulation.__init__(self,**kwargs)
deps = []
for sim in sims:
for dep in sim.dependencies:
deps.append((dep.sim,'other'))
#end for
#end for
self.depends(*deps)
#end def __init__
def propagate_identifier(self):
for sim in self.sims:
sim.propagate_identifier()
#end for
#end def propagate_identifier
def check_result(self,result_name,sim):
return False
#end def check_result
def get_result(self,result_name,sim):
self.error(result_name+' is not calculated by BundledQmcpack')
#end def get_result
def check_dependencies(self,result):
for sim in self.sims:
sim.check_dependencies(results)
#end for
Simulation.check_dependencies(self,result)
#end def check_dependencies
def get_dependencies(self):
for sim in self.sims:
sim.get_dependencies()
#end for
Simulation.get_dependencies(self)
#end def get_dependencies
def check_sim_status(self):
outfile = os.path.join(self.locdir,self.outfile)
errfile = os.path.join(self.locdir,self.errfile)
fobj = open(outfile,'r')
output = fobj.read()
fobj.close()
fobj = open(errfile,'r')
errors = fobj.read()
fobj.close()
ran_to_end = 'Total Execution' in output
files_exist = True
outfiles = self.input.get_output_info('outfiles')
for file in outfiles:
file_loc = os.path.join(self.locdir,file)
files_exist = files_exist and os.path.exists(file_loc)
#end for
if ran_to_end and not files_exist:
self.warn('run finished successfully, but output files do not seem to exist')
print outfiles
print os.listdir(self.locdir)
#end if
aborted = 'Fatal Error' in errors
self.failed = aborted
self.finished = files_exist and self.job.finished and not aborted
#print
#print self.__class__.__name__
#print 'identifier ',self.identifier
#print 'ran_to_end ',ran_to_end
#print 'files_exist',files_exist
#print 'aborted ',aborted
#print 'job done ',self.job.finished
#print 'finished ',self.finished
#print
#
#import code
#code.interact(local=dict(locals(),**globals()))
#end def check_sim_status
#end class BundledQmcpack
def generate_qmcpack(**kwargs):
sim_args,inp_args = Simulation.separate_inputs(kwargs)

15
nexus/library/qmcpack_analyzer.py
View File

@ -21,18 +21,6 @@
#====================================================================#
#
# bugs
# -equilibration length algorithm incorrectly picks first crossing of
# mean+sigma as starting place
# it should be the last crossing of mean+sigma preceeding the first
# crossing of mean-sigma
# -twist averaging energy density gets information from the first twist only
# -integral(Edens) /= E_tot for twist
# -generalized data averaging algorithm must not be finding all the data
# see /oahu/home/projects/project6_germanium/ed_paper_wrapup/final_analysis/data/pure/tile611_kg133/qmc/*
#
#
from time import time
@ -48,8 +36,7 @@ from developer import unavailable
from xmlreader import XMLreader
from plotting import *
#QmcpackAnalyzer classes imports
from qaobject import QAobject
from qmcpack_analyzer_base import QAanalyzer,QAanalyzerCollection
from qmcpack_analyzer_base import QAobject,QAanalyzer,QAanalyzerCollection
from qmcpack_property_analyzers \
import WavefunctionAnalyzer
from qmcpack_quantity_analyzers \

280
nexus/library/qmcpack_analyzer_base.py
View File

@ -6,8 +6,23 @@
#====================================================================#
# qmcpack_analyzer_base.py #
# Data object and analyzer base classes for QmcpackAnalyzer. #
# Maintains data global to these classes. #
# #
# Content summary: #
# QAobject #
# Base class for all QMCPACK analyzer components. #
# Exposes settings options to the user. #
# #
# Checks #
# Class to assess overall validity based on stored results of #
# many checks (boolean values). Only use so far is to validate #
# the structure of Trace files. See qmcpack_method_analyzers.py.#
# #
# Plotter #
# Wrapper class for mayavi visualization of isosurfaces and #
# surface slices. Previously used to visualize energy densities.#
# See qmcpack_quantity_analyzers.py and spacegrid.py. #
# #
# QAdata #
# Represents stored data from QMCPACK's output files. #
# Classification marks it as a target for potential merging, #
@ -17,22 +32,214 @@
# Specialization of QAdata for data from HDF files. #
# #
# QAanalyzer #
# Base class for analyzer classes. #
# Analyzers load and analyze data. #
# Base class functionality includes recursive traversal of #
# nested analyzer object structures for loading and analyzing #
# data. #
# Base class for analyzer classes. Analyzers load and analyze #
# data. Base class functionality includes recursive traversal #
# of nested analyzer object structures for loading and #
# analyzing data. #
# #
#====================================================================#
from numpy import minimum,resize
from generic import obj
from developer import DevBase
from hdfreader import HDFgroup
from qaobject import QAobject
from debug import *
import numpy as np
class Plotter(DevBase):
def __init__(self):
self.initialized = False
return
#end def __init__
def ensure_init(self):
if not self.initialized:
from enthought.mayavi import mlab
from enthought.tvtk.api import tvtk
self.mlab = mlab
self.tvtk = tvtk
self.show = mlab.show
self.plot3d = mlab.plot3d
self.mesh = mlab.mesh
self.initialized = True
#end if
#end def ensure_init
def isosurface(self,points,scalars,contours,dimensions,name='val'):
self.ensure_init()
mlab = self.mlab
tvtk = self.tvtk
sg=tvtk.StructuredGrid(dimensions=dimensions,points=points)
sg.point_data.scalars = scalars
sg.point_data.scalars.name = name
d = mlab.pipeline.add_dataset(sg)
iso = mlab.pipeline.iso_surface(d)
if isinstance(contours,int):
iso.contour.number_of_contours = contours
elif isinstance(contours,list):
iso.contour.auto_contours = False
iso.contour.contours = contours
else:
self.error('isosurface contours must be an int or list\n a '+str(type(contours))+' was provided instead')
#end if
return
#end def isosurface
def surface_slice(self,x,y,z,scalars,options=None):
scale = 1.0
opacity= 1.0
if options!=None:
if 'norm_height' in options:
scale = options.norm_height/abs(scalars.max())
if 'scale' in options:
scale = options.scale
if 'opacity' in options:
opacity = options.opacity
#end if
self.ensure_init()
from numerics import surface_normals
self.mesh(x,y,z,opacity=.2)
surfnorm = scale*surface_normals(x,y,z)
xs=x.copy()
ys=y.copy()
zs=z.copy()
xs[...] = x[...] + surfnorm[...,0]*scalars[...]
ys[...] = y[...] + surfnorm[...,1]*scalars[...]
zs[...] = z[...] + surfnorm[...,2]*scalars[...]
self.mesh(xs,ys,zs,scalars=scalars,opacity=opacity)
return
#end def surface_slice
#end class Plotter
class QAobj_base(DevBase):
None
#end class QAobj_base
class QAobject(QAobj_base):
_global = obj()
_global.dynamic_methods_objects=[]
plotter = Plotter()
opt_methods = set(['opt','linear','cslinear'])
def __init__(self):
return
#end def __init__
@staticmethod
def condense_name(name):
return name.strip().lower().replace(' ','_').replace('-','_').replace('__','_')
#end def condense_name
def _register_dynamic_methods(self):
QAobject._global.dynamic_methods_objects.append(self)
return
#end def _register_dynamic_methods
def _unlink_dynamic_methods(self):
for o in QAobject._global.dynamic_methods_objects:
o._unset_dynamic_methods()
#end for
return
#end def _unlink_dynamic_methods
def _relink_dynamic_methods(self):
for o in QAobject._global.dynamic_methods_objects:
o._reset_dynamic_methods()
#end for
return
#end def _relink_dynamic_methods
_allowed_settings = set(['optimize'])
_default_settings = obj(
#optimize = 'variance'
optimize = 'lastcost'
#optimize = 'energy_within_variance_tol' # also ewvt
)
QAobj_base.class_set(**_default_settings)
@classmethod
def settings(cls,**kwargs):
vars = set(kwargs.keys())
invalid = vars-cls._allowed_settings
if len(invalid)>0:
allowed = list(cls._allowed_settings)
allowed.sort()
invalid = list(invalid)
invalid.sort()
cls.class_error('attempted to set unknown variables\n unknown variables: {0}\n valid options are: {1}'.format(invalid,allowed))
#end if
QAobj_base.class_set(**kwargs)
#end settings
#end class QAobject
class Checks(DevBase):
def __init__(self,label=''):
self._label = label
self._exclusions = set()
#end def __init__
def exclude(self,value):
self._exclusions.add(value)
#end def exclude
def valid(self):
valid = True
for name,value in self.iteritems():
if not (isinstance(name,str) and name.startswith('_')):
if not value in self._exclusions:
valid = valid and value
#end if
#end if
#end if
self._valid = valid
return valid
#end def valid
def write(self,pad=''):
pad2 = pad+' '
if not '_valid' in self:
self.valid()
#end if
valid = self._valid
if valid:
self.log(pad+self._label+' is valid')
else:
self.log(pad+self._label+' is invalid')
for name,value in self.iteritems():
if not (isinstance(name,str) and name.startswith('_')):
if value in self._exclusions:
self.log(pad2+name+' could not be checked')
elif value:
self.log(pad2+name+' is valid')
else:
self.log(pad2+name+' is invalid')
#end if
#end if
#end for
#end if
#end def write
#end class Checks
class QAinformation(obj):
None
#end class QAinformation
@ -147,19 +354,6 @@ class QAanalyzer(QAobject):
dmc_methods = set(['dmc'])
allowed_settings = []
@classmethod
def settings(cls,**kwargs):
invalid = list(set(kwargs.keys())-set(allowed_settings))
if len(invalid)>0:
cls.class_error('invalid variable names encountered in settings\n invalid names: {0}\n valid options are: {1}'.format(invalid,allowed_settings))
#end if
for name,value in kwargs.iteritems():
cls.__dict__[name] = value
#end for
#end def settings
def __init__(self,nindent=0):
self.info = QAinformation(
initialized = False,
@ -218,30 +412,30 @@ class QAanalyzer(QAobject):
None
#end def unset_global_info
def traverse(self,function,block_name=None,callpost=True,**kwargs):
if not callpost:
cls.__dict__[func_name](self,**kwargs)
#end if
if block_name is None or not self.info[block_name]:
for name,value in self.iteritems():
if isinstance(value,QAanalyzer):
value.traverse(value,func_name,block_name,callpost,**kwargs)
elif isinstance(value,QAanalyzerCollection):
for n,v in value.iteritems():
if isinstance(v,QAanalyzer):
v.traverse(v,func_name,block_name,callpost,**kwargs)
#end if
#end for
#end if
#end for
#end if
if block_name!=None:
self.info[block_name] = True
#end if
if callpost:
cls.__dict__[func_name](self,**kwargs)
#end if
#end def traverse
#def traverse(self,function,block_name=None,callpost=True,**kwargs):
# if not callpost:
# cls.__dict__[func_name](self,**kwargs)
# #end if
# if block_name is None or not self.info[block_name]:
# for name,value in self.iteritems():
# if isinstance(value,QAanalyzer):
# value.traverse(value,func_name,block_name,callpost,**kwargs)
# elif isinstance(value,QAanalyzerCollection):
# for n,v in value.iteritems():
# if isinstance(v,QAanalyzer):
# v.traverse(v,func_name,block_name,callpost,**kwargs)
# #end if
# #end for
# #end if
# #end for
# #end if
# if block_name!=None:
# self.info[block_name] = True
# #end if
# if callpost:
# cls.__dict__[func_name](self,**kwargs)
# #end if
##end def traverse
def propagate_indicators(self,**kwargs):
self.reset_indicators(**kwargs)

786
nexus/library/qmcpack_converters.py Normal file
View File

@ -0,0 +1,786 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# qmcpack_converters.py #
# Nexus interfaces for orbital converter tools of QMCPACK: #
# pw2qmcpack, wfconvert, and convert4qmc #
# #
# Content summary: #
# generate_pw2qmcpack #
# User-facing function to create pw2qmcpack simulation objects. #
# #
# generate_pw2qmcpack_input #
# User-facing funcion to create input for pw2qmcpack. #
# #
# Pw2qmcpack #
# Simulation class for pw2qmcpack. #
# #
# Pw2qmcpackInput #
# SimulationInput class for pw2qmcpack. #
# #
# Pw2qmcpackAnalyzer #
# SimulationAnalyzer class for pw2qmcpack. #
# #
# #
# WfconvertInput #
# SimulationInput class for wfconvert. #
# #
# WfconvertAnalyzer #
# SimulationAnalyzer class for wfconvert. #
# #
# Wfconvert #
# Simulation class for wfconvert. #
# #
# generate_wfconvert #
# User-facing function to generate wfconvert simulation objects.#
# #
# #
# Convert4qmcInput #
# Class representing command line interface of convert4qmc. #
# #
# Convert4qmcAnalyzer #
# Placeholder class for output analysis. #
# #
# Convert4qmc #
# Class representing conver4qmc instance. #
# #
# generate_convert4qmc_input #
# Function to generate arbitrary convert4qmc input. #
# #
# generate_convert4qmc #
# Function to generate Convert4qmc simulation object. #
# #
#====================================================================#
import os
from generic import obj
from simulation import Simulation,SimulationInput,SimulationAnalyzer
# read/write functions associated with pw2qmcpack only
def read_str(sv):
return sv.strip('"').strip("'")
#end def read_str
def read_int(sv):
return int(sv)
#end def read_int
def read_float(sv):
return float(sv.replace('d','e').replace('D','e'))
#end def read_float
bconv = {'.true.':True,'.false.':False}
def read_bool(sv):
return bconv[sv]
#end def read_bool
def write_str(val):
return "'"+val+"'"
#end def write_str
def write_int(val):
return str(val)
#end def write_int
def write_float(val):
return str(val)
#end def write_float
def write_bool(val):
return '.'+str(val).lower()+'.'
#end def write_bool
readval={str:read_str,int:read_int,float:read_float,bool:read_bool}
writeval={str:write_str,int:write_int,float:write_float,bool:write_bool}
class Pw2qmcpackInput(SimulationInput):
ints = []
floats = []
strs = ['outdir','prefix']
bools = ['write_psir']
var_types = dict()
for v in ints:
var_types[v]=int
for v in floats:
var_types[v]=float
for v in strs:
var_types[v]=str
for v in bools:
var_types[v]=bool
allowed = set(ints+floats+strs+bools)
def read_text(self,contents,filepath=None):
lines = contents.split('\n')
inside = False
for l in lines:
if inside:
tokens = l.split(',')
for t in tokens:
ts = t.strip()
if ts!='' and ts!='/':
name,value = ts.split('=')
name = name.strip()
value= value.strip()
if name in self.allowed:
vtype = self.var_types[name]
value = readval[vtype](value)
sobj[name]=value
else:
self.error('encountered unknown variable: '+name)
#end if
#end if
#end for
#end if
if '&' in l:
inside=True
section = l.lstrip('&').lower()
sobj = obj()
self[section]=sobj
elif l.strip()=='/':
inside=False
#end if
#end for
#end def read_text
def write_text(self,filepath=None):
contents = ''
for sname,section in self.iteritems():
contents+='&'+sname+'\n'
for name,value in section.iteritems():
vtype = type(value)
contents += ' '+name+' = '+writeval[vtype](value)+'\n'
#end for
contents+='/\n'
#end for
return contents
#end def write_text
def __init__(self,filepath=None,**vars):
if filepath!=None:
self.read(filepath)
else:
inputpp = obj()
for name,value in vars.iteritems():
inputpp[name] = value
#end for
self.inputpp = inputpp
#end if
#end def __init__
#end class Pw2qmcpackInput
def generate_pw2qmcpack_input(prefix='pwscf',outdir='pwscf_output',write_psir=True):
pw = Pw2qmcpackInput(
prefix = prefix,
outdir = outdir,
write_psir = write_psir
)
return pw
#end def generate_pw2qmcpack_input
class Pw2qmcpackAnalyzer(SimulationAnalyzer):
def __init__(self,arg0):
if isinstance(arg0,Simulation):
sim = arg0
self.infile = sim.infile
prefix,outdir = sim.input.inputpp.tuple('prefix','outdir')
self.dir = sim.locdir
self.h5file = os.path.join(sim.locdir,outdir,prefix+'.pwscf.h5')
else:
self.infile = arg0
#end if
#end def __init__
def analyze(self):
if False:
import h5py
self.log('Fixing h5 file',n=5)
path = os.path.split(self.h5file)[0]
print os.getcwd()
print os.listdir('./')
if os.path.exists(path):
print os.listdir(path)
#end if
print self.h5file
h = h5py.File(self.h5file)
if 'electrons' in h:
elec = h['electrons']
nkpoints = 0
for name,val in elec.iteritems():
if name.startswith('kpoint'):
nkpoints+=1
#end for
#end if
nkold = elec['number_of_kpoints'][0]
self.log('Were',nkold,'kpoints, now',nkpoints,'kpoints',n=6)
elec['number_of_kpoints'][0] = nkpoints
#end for
#end if
#end def analyze
def get_result(self,result_name):
self.not_implemented()
#end def get_result
#end class Pw2qmcpackAnalyzer
class Pw2qmcpack(Simulation):
input_type = Pw2qmcpackInput
analyzer_type = Pw2qmcpackAnalyzer
generic_identifier = 'pw2qmcpack'
application = 'pw2qmcpack.x'
application_properties = set(['serial'])
application_results = set(['orbitals'])
def check_result(self,result_name,sim):
calculating_result = False
inputpp = self.input.inputpp
if result_name=='orbitals':
calculating_result = True
else:
calculating_result = False
self.error('ability to check for result '+result_name+' has not been implemented')
#end if
return calculating_result
#end def check_result
def get_result(self,result_name,sim):
result = obj()
inputpp = self.input.inputpp
prefix = 'pwscf'
outdir = './'
if 'prefix' in inputpp:
prefix = inputpp.prefix
#end if
if 'outdir' in inputpp:
outdir = inputpp.outdir
#end if
if outdir.startswith('./'):
outdir = outdir[2:]
#end if
if result_name=='orbitals':
result.h5file = os.path.join(self.locdir,outdir,prefix+'.pwscf.h5')
result.ptcl_xml = os.path.join(self.locdir,outdir,prefix+'.ptcl.xml')
result.wfs_xml = os.path.join(self.locdir,outdir,prefix+'.wfs.xml')
else:
self.error('ability to get result '+result_name+' has not been implemented')
#end if
return result
#end def get_result
def incorporate_result(self,result_name,result,sim):
if result_name=='orbitals':
pwin = sim.input.control
p2in = self.input.inputpp
pwprefix = 'pwscf'
p2prefix = 'pwscf'
pwoutdir = './'
p2outdir = './'
if 'prefix' in pwin:
pwprefix = pwin.prefix
if 'prefix' in p2in:
p2prefix = p2in.prefix
if 'outdir' in pwin:
pwoutdir = pwin.outdir
if 'outdir' in p2in:
p2outdir = p2in.outdir
if pwoutdir.startswith('./'):
pwoutdir = pwoutdir[2:]
if p2outdir.startswith('./'):
p2outdir = p2outdir[2:]
pwdir = os.path.abspath(os.path.join(sim.locdir ,pwoutdir))
p2dir = os.path.abspath(os.path.join(self.locdir,p2outdir))
errors = False
if pwdir!=p2dir:
self.error('to use orbitals, '+self.generic_identifier+' must have the same outdir as pwscf\n pwscf outdir: '+pwdir+'\n '+self.generic_identifier+' outdir: '+p2dir,exit=False)
errors = True
#end if
if pwprefix!=p2prefix:
self.error('to use orbitals, '+self.generic_identifier+' must have the same prefix as pwscf\n pwscf prefix: '+pwprefix+'\n '+self.generic_identifier+' prefix: '+p2prefix,exit=False)
errors = True
#end if
if errors:
self.error(self.generic_identifier+' cannot use orbitals from pwscf')
#end if
else:
self.error('ability to incorporate result '+result_name+' has not been implemented')
#end if
#end def incorporate_result
def check_sim_status(self):
outfile = os.path.join(self.locdir,self.outfile)
fobj = open(outfile,'r')
output = fobj.read()
fobj.close()
inputpp = self.input.inputpp
prefix = 'pwscf'
outdir = './'
if 'prefix' in inputpp:
prefix = inputpp.prefix
#end if
if 'outdir' in inputpp:
outdir = inputpp.outdir
#end if
if outdir.startswith('./'):
outdir = outdir[2:]
#end if
h5file = os.path.join(self.locdir,outdir,prefix+'.pwscf.h5')
ptcl_xml = os.path.join(self.locdir,outdir,prefix+'.ptcl.xml')
wfs_xml = os.path.join(self.locdir,outdir,prefix+'.wfs.xml')
must_exist = [h5file,ptcl_xml,wfs_xml]
files_exist = True
for file in must_exist:
files_exist = files_exist and os.path.exists(file)
#end for
outfin = True
#outfin = outfin and 'esh5 create' in output
#outfin = outfin and 'Creating electrons' in output
outfin = outfin and 'npw=' in output
outfin = outfin and 'ik=' in output
outfin = outfin or 'JOB DONE' in output
success = files_exist and outfin
#self.finished = success and self.job.finished
# pw2qmcpack has too many variants to assess completion based on log output
# assume (optimistically) that job completion indicates success
self.finished = files_exist and self.job.finished
#end def check_sim_status
def get_output_files(self):
output_files = []
return output_files
#end def get_output_files
def app_command(self):
return self.app_name+'<'+self.infile
#end def app_command
#end class Pw2qmcpack
def generate_pw2qmcpack(**kwargs):
sim_args,inp_args = Simulation.separate_inputs(kwargs)
if not 'input' in sim_args:
sim_args.input = generate_pw2qmcpack_input(**inp_args)
#end if
pw2qmcpack = Pw2qmcpack(**sim_args)
return pw2qmcpack
#end def generate_pw2qmcpack
class WfconvertInput(SimulationInput):
def __init__(self,app_name='wfconvert',h5in='MISSING.h5',h5out='wfconvert.h5',spline=False,format='eshdf',factor=None):
self.app_name = app_name
self.h5in = h5in
self.h5out= h5out
self.spline = spline
self.format = format
self.factor = factor
#end def __init__
#wfconvert --nospline --eshdf diamond.h5 out/diamond.pwscf.h5 >& diamond-wfconvert.out
def set_app_name(self,app_name):
self.app_name = app_name
#end def set_app_name
def app_command(self):
c = self.app_name+' '
if not self.spline:
c+= '--nospline '
#end if
c+='--'+self.format+' '+self.h5out+' '+self.h5in
return c
#end def app_command
def read(self,filepath):
None
#end def read
def write_text(self,filepath=None):
return self.app_command()
#end def write_text
#end class WfconvertInput
def generate_wfconvert_input(app_name='wfconvert',h5in='MISSING.h5',h5out='wfconvert.h5',spline=False,format='eshdf',factor=None):
wi = WfconvertInput(
app_name = app_name,
h5in = h5in,
h5out = h5out,
spline = spline,
format = format,
factor = factor
)
return wi
#end def generate_wfconvert_input
class WfconvertAnalyzer(SimulationAnalyzer):
def __init__(self,arg0):
if isinstance(arg0,Simulation):
sim = arg0
self.infile = sim.infile
self.dir = sim.locdir
self.h5file = os.path.join(sim.locdir,sim.input.h5out)
else:
self.infile = arg0
#end if
#end def __init__
def analyze(self):
if False:
import h5py
self.log('Fixing h5 file',n=5)
h = h5py.File(self.h5file)
if 'electrons' in h:
elec = h['electrons']
nkpoints = 0
for name,val in elec.iteritems():
if name.startswith('kpoint'):
nkpoints+=1
#end for
#end if
nkold = elec['number_of_kpoints'][0]
self.log('Were',nkold,'kpoints, now',nkpoints,'kpoints',n=6)
elec['number_of_kpoints'][0] = nkpoints
#end for
#end if
#end def analyze
#end class WfconvertAnalyzer
class Wfconvert(Simulation):
input_type = WfconvertInput
analyzer_type = WfconvertAnalyzer
generic_identifier = 'wfconvert'
application = 'wfconvert'
application_properties = set(['serial'])
application_results = set(['orbitals'])
def set_app_name(self,app_name):
self.app_name = app_name
self.input.set_app_name(app_name)
#end def set_app_name
def check_result(self,result_name,sim):
calculating_result = False
if result_name=='orbitals':
calculating_result = True
else:
calculating_result = False
self.error('ability to check for result '+result_name+' has not been implemented')
#end if
return calculating_result
#end def check_result
def get_result(self,result_name,sim):
result = obj()
if result_name=='orbitals':
result.h5file = os.path.join(self.locdir,self.input.h5out)
result.outfile = os.path.join(self.locdir,self.outfile)
else:
self.error('ability to get result '+result_name+' has not been implemented')
#end if
return result
#end def get_result
def incorporate_result(self,result_name,result,sim):
if result_name=='orbitals':
self.input.h5in = os.path.relpath(result.h5file,self.locdir)
self.job.app_command = self.input.app_command()
else:
self.error('ability to incorporate result '+result_name+' has not been implemented')
#end if
#end def incorporate_result
def check_sim_status(self):
outfile = os.path.join(self.locdir,self.outfile)
errfile = os.path.join(self.locdir,self.errfile)
fobj = open(outfile,'r')
output = fobj.read()
fobj.close()
fobj = open(errfile,'r')
errors = fobj.read()
fobj.close()
h5file = os.path.join(self.locdir,self.input.h5out)
file_exists = os.path.exists(h5file)
outfin = 'Successfully read' in errors and 'numSpins' in errors
outfin = outfin and 'Writing laplacians' in output
success = file_exists and outfin
self.finished = success
#end def check_sim_status
def get_output_files(self):
output_files = []
return output_files
#end def get_output_files
def app_command(self):
# app_name is passed along in post_init
return self.input.app_command()
#end def app_command
#end class Wfconvert
def generate_wfconvert(**kwargs):
sim_args,inp_args = Simulation.separate_inputs(kwargs)
if not 'input' in sim_args:
sim_args.input = generate_wfconvert_input(**inp_args)
#end if
wfconvert = Wfconvert(**sim_args)
return wfconvert
#end def generate_wfconvert
class Convert4qmcInput(SimulationInput):
def __init__(self,
app_name = 'convert4qmc',
prefix = None,
gamess_ascii = None,
ci = None,
read_initial_guess = None,
natural_orbitals = None,
threshold = None,
zero_ci = False,
add_3body_J = False
):
self.prefix = prefix
self.app_name = app_name
self.gamess_ascii = gamess_ascii
self.ci = ci
self.read_initial_guess = read_initial_guess
self.natural_orbitals = natural_orbitals
self.threshold = threshold
self.zero_ci = zero_ci
self.add_3body_J = add_3body_J
#end def __init__
def set_app_name(self,app_name):
self.app_name = app_name
#end def set_app_name
def app_command(self):
c = self.app_name
if self.prefix!=None:
c += ' -prefix '+self.prefix
#end if
if self.gamess_ascii!=None:
c += ' -gamessAscii '+self.gamess_ascii
#end if
if self.ci!=None:
c += ' -ci '+self.ci
#end if
if self.threshold!=None:
c += ' -threshold '+str(self.threshold)
#end if
if self.zero_ci:
c += ' -zeroCi'
#end if
if self.read_initial_guess!=None:
c += ' -readInitialGuess '+str(self.read_initial_guess)
#end if
if self.natural_orbitals!=None:
c += ' -NaturalOrbitals '+str(self.natural_orbitals)
#end if
if self.add_3body_J:
c += ' -add3BodyJ'
#end if
return c
#end def app_command
def read(self,filepath):
None
#end def read
def write_text(self,filepath=None):
return self.app_command()
#end def write_text
def output_files(self):
prefix = 'sample'
if self.prefix!=None:
prefix = self.prefix
#end if
wfn_file = prefix+'.Gaussian-G2.xml'
ptcl_file = prefix+'.Gaussian-G2.ptcl.xml'
return wfn_file,ptcl_file
#end def output_files
#end class Convert4qmcInput
def generate_convert4qmc_input(**kwargs):
return Convert4qmcInput(**kwargs)
#end def generate_convert4qmc_input
class Convert4qmcAnalyzer(SimulationAnalyzer):
def __init__(self,arg0):
if isinstance(arg0,Simulation):
self.infile = arg0.infile
else:
self.infile = arg0
#end if
#end def __init__
def analyze(self):
None
#end def analyze
#end class Convert4qmcAnalyzer
class Convert4qmc(Simulation):
input_type = Convert4qmcInput
analyzer_type = Convert4qmcAnalyzer
generic_identifier = 'convert4qmc'
application = 'convert4qmc'
application_properties = set(['serial'])
application_results = set(['orbitals'])
def set_app_name(self,app_name):
self.app_name = app_name
self.input.set_app_name(app_name)
#end def set_app_name
def propagate_identifier(self):
None
#self.input.prefix = self.identifier
#end def propagate_identifier
def check_result(self,result_name,sim):
calculating_result = False
if result_name=='orbitals':
calculating_result = True
else:
calculating_result = False
self.error('ability to check for result '+result_name+' has not been implemented')
#end if
return calculating_result
#end def check_result
def get_result(self,result_name,sim):
result = obj()
input = self.input
if result_name=='orbitals':
wfn_file,ptcl_file = input.output_files()
result.location = os.path.join(self.locdir,wfn_file)
else:
self.error('ability to get result '+result_name+' has not been implemented')
#end if
return result
#end def get_result
def incorporate_result(self,result_name,result,sim):
if result_name=='orbitals':
orbpath = os.path.relpath(result.location,self.locdir)
if result.scftyp=='mcscf':
self.input.gamess_ascii = orbpath
self.input.ci = orbpath
elif result.scftyp=='none': # cisd, etc
self.input.gamess_ascii = orbpath
self.input.ci = orbpath
if result.mos>0:
self.input.read_initial_guess = result.mos
elif result.norbitals>0:
self.input.read_initial_guess = result.norbitals
#end if
else:
self.input.gamess_ascii = orbpath
#end if
self.job.app_command = self.input.app_command()
else:
self.error('ability to incorporate result '+result_name+' has not been implemented')
#end if
#end def incorporate_result
def check_sim_status(self):
output = open(os.path.join(self.locdir,self.outfile),'r').read()
#errors = open(os.path.join(self.locdir,self.errfile),'r').read()
success = 'QMCGaussianParserBase::dump' in output
for filename in self.input.output_files():
success &= os.path.exists(os.path.join(self.locdir,filename))
#end for
self.failed = not success
self.finished = self.job.finished
#end def check_sim_status
def get_output_files(self):
output_files = []
return output_files
#end def get_output_files
def app_command(self):
return self.input.app_command()
#end def app_command
#end class Convert4qmc
def generate_convert4qmc(**kwargs):
sim_args,inp_args = Simulation.separate_inputs(kwargs)
if 'identifier' in sim_args and not 'prefix' in inp_args:
inp_args.prefix = sim_args.identifier
#end if
if not 'input' in sim_args:
sim_args.input = generate_convert4qmc_input(**inp_args)
#end if
convert4qmc = Convert4qmc(**sim_args)
return convert4qmc
#end def generate_convert4qmc

69
nexus/library/qmcpack_input.py
View File

@ -133,7 +133,7 @@
import os
import types as pytypes
import inspect
import keyword
from numpy import fromstring,empty,array,float64,\
loadtxt,ndarray,dtype,sqrt,pi,arange,exp,eye,\
@ -142,12 +142,11 @@ from StringIO import StringIO
from superstring import string2val
from generic import obj,hidden
from xmlreader import XMLreader,XMLelement
from project_base import Pobj
from developer import DevBase
from structure import Structure,Jellium
from physical_system import PhysicalSystem
from simulation import SimulationInput,SimulationInputTemplate
from pwscf_input import array_to_string as pwscf_array_string
from debug import *
from debug import ci as interact
@ -238,7 +237,7 @@ bool_write_types = set([yesno,onezero,truefalse])
class QIobj(Pobj):
class QIobj(DevBase):
None
#end class QIobj
@ -488,30 +487,25 @@ class QIxml(Names):
@classmethod
def init_class(cls):
vars = cls.__dict__.keys()
init_vars = dict(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
)
for var,val in init_vars.iteritems():
if not var in vars:
cls.__dict__[var] = val
#end if
#end for
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
)
for v in ['attributes','elements','parameters','attribs','costs','h5tags']:
names = cls.__dict__[v]
names = cls.class_get(v)
for i in range(len(names)):
if names[i] in cls.escape_names:
names[i]+='_'
@ -599,7 +593,7 @@ class QIxml(Names):
if e in self:
elem = self[e]
if isinstance(elem,QIxml):
c += self[e].write(indent_level+1)
c += elem.write(indent_level+1)
else:
begin = '<'+e+'>'
contents = param.write(elem)
@ -629,6 +623,7 @@ class QIxml(Names):
c+=indent+'</'+expanded_tag+'>\n'
#end if
param.reset_precision()
return c
#end def write
@ -797,7 +792,6 @@ class QIxml(Names):
def incorporate_defaults(self,elements=False,overwrite=False,propagate=True):
for name,value in self.defaults.iteritems():
valtype = type(value)
defval=None
if isinstance(value,classcollection):
if elements:
@ -809,11 +803,11 @@ class QIxml(Names):
#end for
defval = make_collection(coll)
#end if
elif valtype==pytypes.ClassType:
elif inspect.isclass(value):
if elements:
defval = value()
#end if
elif valtype==pytypes.FunctionType:
elif inspect.isfunction(value):
defval = value()
else:
defval = value
@ -1393,7 +1387,7 @@ class QIxml(Names):
#end def get_host
def get_precision(self):
return self.__class__.__dict__['precision']
return self.__class__.class_get('precision')
#end def get_precision
#end class QIxml
@ -2825,7 +2819,7 @@ class QmcpackInput(SimulationInput,Names):
#end def read
def write_contents(self):
def write_text(self,filepath=None):
c = ''
header = '''<?xml version="1.0"?>
'''
@ -2839,7 +2833,7 @@ class QmcpackInput(SimulationInput,Names):
c+=base.write(first=True)
Param.metadata = None
return c
#end def write_contents
#end def write_text
def unroll_calculations(self,modify=True):
@ -3757,6 +3751,8 @@ class QmcpackInput(SimulationInput,Names):
# base class for bundled qmcpack input
# not used on its own
class BundledQmcpackInput(SimulationInput):
def __init__(self,inputs,filenames):
@ -3799,8 +3795,8 @@ class BundledQmcpackInput(SimulationInput):
#end def get_output_info
def generate_filenames(self):
None
def generate_filenames(self,infile):
self.not_implemented()
#end def generate_filenames
@ -3845,6 +3841,7 @@ class TracedQmcpackInput(BundledQmcpackInput):
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

3
nexus/library/qmcpack_method_analyzers.py
View File

@ -21,8 +21,7 @@ import os
import re
from generic import obj
from hdfreader import HDFreader
from qaobject import Checks
from qmcpack_analyzer_base import QAanalyzer,QAdata,QAHDFdata
from qmcpack_analyzer_base import Checks,QAanalyzer,QAdata,QAHDFdata
from qmcpack_property_analyzers import WavefunctionAnalyzer
from qmcpack_quantity_analyzers import HDFAnalyzer
from debug import *

3
nexus/library/qmcpack_property_analyzers.py
View File

@ -38,8 +38,7 @@
import os
from numpy import loadtxt,array,ones,dot,floor,zeros
from qmcpack_input import QmcpackInput
from qaobject import QAobject
from qmcpack_analyzer_base import QAanalyzer
from qmcpack_analyzer_base import QAobject,QAanalyzer
from developer import unavailable
from debug import *
try:

1373
nexus/library/qmcpack_quantity_analyzers.py
View File

File diff suppressed because it is too large Load Diff

121
nexus/library/qmcpack_result_analyzers.py
View File

@ -26,12 +26,12 @@
#====================================================================#
from numpy import array,empty,zeros,sqrt
from numpy import array,empty,zeros,sqrt,arange
from generic import obj
from unit_converter import convert
from qmcpack_input import QmcpackInput
from qaobject import QAobject
from qmcpack_analyzer_base import QAanalyzer
from qmcpack_analyzer_base import QAobject,QAanalyzer
from debug import *
@ -94,43 +94,75 @@ class OptimizationAnalyzer(ResultAnalyzer):
vw = self.variance_weight
Efail = 1e6
EVratio_fail = 0.10
Vfail = 1e3
EVratio_fail = 0.30
EVratio_soft_fail = 0.15
#save the energies and variances of opt iterations
nopt = len(opts)
en = zeros((nopt,))+1e99
enerr = zeros((nopt,))+1e99
va = zeros((nopt,))+1e99
vaerr = zeros((nopt,))+1e99
res = obj()
variance_present = False
any_complete = False
all_complete = True
unstable = False
for i in range(nopt):
o = opts[i]
complete = o.info.complete
any_stable = False
for s,opt in opts.iteritems():
complete = opt.info.complete
any_complete |= complete
all_complete &= complete
if complete:
le = o.scalars.LocalEnergy
en[i] = le.mean
enerr[i] = le.error
if 'LocalEnergyVariance' in o.scalars:
fail = False
le = opt.scalars.LocalEnergy
en = le.mean
enerr = le.error
fail |= abs(en)>Efail
if 'LocalEnergyVariance' in opt.scalars:
variance_present = True
lev = o.scalars.LocalEnergyVariance
va[i] = lev.mean
vaerr[i] = lev.error
lev = opt.scalars.LocalEnergyVariance
va = lev.mean
vaerr = lev.error
fail |= abs(va)>Vfail or abs(va/en)>EVratio_fail
#end if
unstable |= abs(le.mean)>Efail
if not fail:
any_stable = True
sres = obj()
sres.en = en
sres.enerr = enerr
if variance_present:
sres.va = va
sres.vaerr = vaerr
#end if
res[s] = sres
#end if
unstable|=fail
#end if
#end for
unstable |= not any_complete
nseries = len(res)
en = zeros((nseries,),dtype=float)
enerr = zeros((nseries,),dtype=float)
va = zeros((nseries,),dtype=float)
vaerr = zeros((nseries,),dtype=float)
series = array(sorted(res.keys()),dtype=int)
i = 0
for s in series:
sres = res[s]
en[i] = sres.en
enerr[i] = sres.enerr
if variance_present:
va[i] = sres.va
vaerr[i] = sres.vaerr
#end if
i+=1
#end for
self.set(
any_complete = any_complete,
all_complete = all_complete,
unstable = unstable,
series = series,
energy = en,
energy_error = enerr,
variance = va,
@ -146,44 +178,37 @@ class OptimizationAnalyzer(ResultAnalyzer):
elif optimize=='energy':
ew = 1.0
vw = 0.0
elif optimize=='energy_within_variance_tol' or optimize=='ewvt':
None
elif isinstance(optimize,(tuple,list)) and len(optimize)==2:
ew,vw = optimize
else:
self.error('selection for optimization is invalid\n optimize setting: {0}\n valid options are: energy, variance, or a length 2 tuple containing the cost of each, e.g. (.5,.5)'.format(optimize))
self.error('selection for optimization is invalid\noptimize setting: {0}\nvalid options are: energy, variance, energy_within_variance_tol, or a length 2 tuple containing the cost of energy and variance, e.g. (.5,.5)'.format(optimize))
#end if
self.failed = True
self.optimal_series = None
self.optimal_file = None
self.optimal_wavefunction = None
if any_complete:
mincost = 1e99
index = -1
for i,opt in opts.iteritems():
if opt.info.complete:
s = opt.scalars
e = s.LocalEnergy.mean
if 'LocalEnergyVariance' in s:
v = s.LocalEnergyVariance.mean
else:
v = 0
#end if
cost = ew*e+vw*v
if cost<mincost:
mincost = cost
index = i
#end if
#end if
#end for
if index!=-1:
failed = abs(en[index])>Efail or (va[index]<1e98 and abs(va[index]/en[index])>EVratio_fail)
self.failed = failed
self.optimal_series = index
self.optimal_file = opts[index].info.files.opt
self.optimal_wavefunction = opts[index].wavefunction.info.wfn_xml.copy()
if any_stable:
if optimize=='energy_within_variance_tol' or optimize=='ewvt':
indices = arange(len(series),dtype=int)
vartol = 0.2
vmin = va.min()
vind = indices[abs(va-vmin)/vmin<vartol]
index = vind[en[vind].argmin()]
opt_series = series[index]
else:
cost = en*ew+va*vw
index = cost.argmin()
opt_series = series[index]
#end if
failed = abs(en[index])>Efail or abs(va[index])>Vfail or abs(va[index]/en[index])>EVratio_soft_fail
self.failed = failed
self.optimal_series = opt_series
self.optimal_file = opts[opt_series].info.files.opt
self.optimal_wavefunction = opts[opt_series].wavefunction.info.wfn_xml.copy()
#end if
#end def analyze_local

752
nexus/library/simulation.py
View File

File diff suppressed because it is too large Load Diff

1335
nexus/library/spacegrid.py
View File

File diff suppressed because it is too large Load Diff

55
nexus/library/sqd_input.py
View File

@ -23,7 +23,7 @@
import os
import types as pytypes
import inspect
import keyword
from numpy import fromstring,empty,array,float64,\
loadtxt,ndarray,dtype,sqrt,pi,arange,exp,eye,ceil,mod
@ -32,7 +32,7 @@ from superstring import string2val
from generic import obj
from xmlreader import XMLreader,XMLelement
from periodic_table import pt as periodic_table
from project_base import Pobj
from developer import DevBase
from structure import Structure
from physical_system import PhysicalSystem
from simulation import SimulationInput
@ -127,7 +127,7 @@ def string_to_val(s):
#end def string_to_val
class SQDobj(Pobj):
class SQDobj(DevBase):
None
#end class SQDobj
@ -312,24 +312,19 @@ class SQDxml(Names):
@classmethod
def init_class(cls):
vars = cls.__dict__.keys()
init_vars = dict(tag = cls.__name__,
attributes = [],
elements = [],
text = None,
parameters = [],
attribs = [],
costs = [],
h5tags = [],
defaults = obj()
)
for var,val in init_vars.iteritems():
if not var in vars:
cls.__dict__[var] = val
#end if
#end for
cls.class_set_optional(
tag = cls.__name__,
attributes = [],
elements = [],
text = None,
parameters = [],
attribs = [],
costs = [],
h5tags = [],
defaults = obj()
)
for v in ['attributes','elements','parameters','attribs','costs','h5tags']:
names = cls.__dict__[v]
names = cls.class_get(v)
for i in range(len(names)):
if names[i] in cls.escape_names:
names[i]+='_'
@ -600,7 +595,6 @@ class SQDxml(Names):
def incorporate_defaults(self,elements=False,overwrite=False,propagate=True):
for name,value in self.defaults.iteritems():
valtype = type(value)
defval=None
if isinstance(value,classcollection):
if elements:
@ -612,11 +606,11 @@ class SQDxml(Names):
#end for
defval = make_collection(coll)
#end if
elif valtype==pytypes.ClassType:
elif inspect.isclass(value):
if elements:
defval = value()
#end if
elif valtype==pytypes.FunctionType:
elif inspect.isfunction(value):
defval = value()
else:
defval = value
@ -1483,7 +1477,7 @@ class SqdInput(SimulationInput,Names):
#end def read
def write_contents(self):
def write_text(self,filepath=None):
c = ''
header = '''<?xml version="1.0"?>
'''
@ -1497,7 +1491,7 @@ class SqdInput(SimulationInput,Names):
c+=base.write(first=True)
Param.metadata = None
return c
#end def write_contents
#end def write_text
def unroll_calculations(self,modify=True):
qmc = []
@ -1705,7 +1699,7 @@ class SqdInput(SimulationInput,Names):
class Shells(Pobj):
class Shells(DevBase):
channel_names = tuple('spdfghik')
channels = obj()
for l in range(len(channel_names)):
@ -1898,7 +1892,7 @@ class Shells(Pobj):
self.error('unexpected key values for shells\n expected values: '+str(shells.keys())+'\n you provided '+str(self.all_shells.keys()))
#end if
else:
self.error('must provide a string, dict, or obj describing atomic shells\n you provided '+str(shell_str))
self.error('must provide a string, dict, or obj describing atomic shells\n you provided '+str(shells))
#end if
self.check_shells()
#end def __init__
@ -1958,7 +1952,7 @@ class Shells(Pobj):
m=v1
ma = abs(array(m))
if ma.max()>self.channel_indices[l]:
self.error('maximum |m| for {0} channel is {1}\n you requested {2}: {3}'.format(l,channel_indices[l],ma.max(),m))
self.error('maximum |m| for {0} channel is {1}\n you requested {2}: {3}'.format(l,self.channel_indices[l],ma.max(),m))
#end if
channel = m
m=None
@ -2082,11 +2076,6 @@ class Shells(Pobj):
#end for
return orbitals
#end def orbitals
def error(self,msg,exit=True):
Pobj.error(self,msg,self.location,exit=exit)
#end def error
#end class Shells
hunds_rule_filling = Shells.hunds_rule_filling

22
nexus/library/structure.py
View File

@ -32,7 +32,7 @@ from numpy import array,floor,empty,dot,diag,sqrt,pi,mgrid,exp,append,arange,cei
from numpy.linalg import inv,det,norm
from types import NoneType
from unit_converter import convert
from extended_numpy import nearest_neighbors,convex_hull,voronoi_neighbors
from numerics import nearest_neighbors,convex_hull,voronoi_neighbors
from periodic_table import pt,is_element
from generic import obj
from developer import DevBase,unavailable,error,warn
@ -603,9 +603,12 @@ class Structure(Sobj):
# while remaining periodically correct
# note that the unshearing procedure is not unique
# it depends on the order of unshearing operations
def unsheared_axes(self,distances=False):
def unsheared_axes(self,axes=None,distances=False):
if self.dim!=3:
self.error('rinscribe is currently only implemented for 3 dimensions')
self.error('unsheared_axes is currently only implemented for 3 dimensions')
#end if
if axes is None:
axes = self.axes
#end if
dim=3
axbar = identity(dim)
@ -654,11 +657,6 @@ class Structure(Sobj):
return self.face_vectors(distances=True)[1]
#end def face_distances
def set_orig(self):
self.orig_pos = pos.copy()
#end def set_orig
def rescale(self,scale):
self.scale *= scale
@ -2957,19 +2955,19 @@ class Structure(Sobj):
#end if
#end if
c = open(filepath,'r').read()
self.read_contents(c,format,elem=elem)
self.read_text(c,format,elem=elem)
return c
#end def read
def read_contents(self,contents,format,elem=None):
def read_text(self,contents,format,elem=None):
format = format.lower()
if format=='poscar':
self.read_poscar(contents,elem=elem,contents=True)
else:
self.error('unrecognized file format: {0}'.format(format))
#end if
#end def read_contents
#end def read_text
def read_poscar(self,filepath,elem=None,contents=False):
@ -4182,7 +4180,7 @@ def generate_crystal_structure(lattice=None,cell=None,centering=None,
basis_vectors=None,tiling=None,cscale=None,
axes=None,units=None,angular_units='degrees',
magnetization=None,magnetic_order=None,magnetic_prim=True,
kpoints=None,kgrid=None,kshift=(0,0,0),permute=None,
kpoints=None,kweights=None,kgrid=None,kshift=(0,0,0),permute=None,
structure=None,shape=None,element=None,scale=None, #legacy inputs
operations=None,
struct_type=Crystal,elem=None,pos=None,frozen=None,

32
nexus/library/template_simulation.py
View File

@ -54,7 +54,7 @@ from simulation import Simulation,SimulationInput,SimulationAnalyzer
# in a form TemplateSimulationInput understands (this depends on your implementation)
#
# required functions to be implemented:
# TemplateSimulationInput: read_contents, write_contents
# TemplateSimulationInput: read_text, write_text
# TemplateSimulation: app_command, check_sim_status
#
# 2) generated standalone simulation
@ -64,7 +64,7 @@ from simulation import Simulation,SimulationInput,SimulationAnalyzer
# a standard PhysicalSystem object
#
# required functions to be implemented:
# TemplateSimulationInput: read_contents, write_contents, incorporate_system
# TemplateSimulationInput: read_text, write_text, incorporate_system
# TemplateSimulation: app_command, check_sim_status
# generate_template_simulation_input
#
@ -76,7 +76,7 @@ from simulation import Simulation,SimulationInput,SimulationAnalyzer
# this information is used by the others to populate input files
#
# required functions to be implemented:
# TemplateSimulationInput: read_contents, write_contents
# TemplateSimulationInput: read_text, write_text
# TemplateSimulation: app_command,check_sim_status,
# check_result, get_result
#
@ -89,7 +89,7 @@ from simulation import Simulation,SimulationInput,SimulationAnalyzer
# info is used to populate own input file
#
# required functions to be implemented:
# TemplateSimulationInput: read_contents, write_contents
# TemplateSimulationInput: read_text, write_text
# TemplateSimulation: app_command,check_sim_status,
# check_result, get_result,
# incorporate_result
@ -111,10 +111,10 @@ class TemplateSimulationInput(SimulationInput):
#end def __init__
def read_contents(self,contents):
def read_text(self,text,filepath=None):
# required
# the string 'contents' contains the text of an input file
# translate contents into an internal representation of the input
# the string 'text' contains the text of an input file
# translate text into an internal representation of the input
# for example, an input file looking like:
# section_a
# var_a 1
@ -132,7 +132,7 @@ class TemplateSimulationInput(SimulationInput):
# sec['var_d'] = True
# self['section_b'] = sec
# with some line parsing:
# for line in contents.splitlines():
# for line in text.splitlines():
# # parse lines here
# #end for
# the resulting input object can then be worked with in a simple way:
@ -153,23 +153,23 @@ class TemplateSimulationInput(SimulationInput):
#
# >>> input.section_b.var_c = 25.0
None
#end def read_contents
#end def read_text
def write_contents(self):
def write_text(self,filepath=None):
# required
# translate the internal representation of input into a string
# for the above example, this might look like:
# contents = ''
# text = ''
# for secname,sec in self.iteritems():
# contents += secname+'\n'
# text += secname+'\n'
# for val,val in sec.iteritems():
# contents += ' '+var+' = '+val2str(val)
# text += ' '+var+' = '+val2str(val)
# #end for
# #end for
contents = ''
return contents
#end def write_contents
text = ''
return text
#end def write_text
def incorporate_system(self,system):

245
nexus/library/unit_converter.py
View File

@ -20,31 +20,44 @@
#====================================================================#
from abilities import Callable
from numpy import dot,zeros
from numpy.linalg import inv
from generic import obj
from developer import DevBase
class Unit:
def __init__(self,name,symbol,type,value,shift=0):
self.name = name
class Unit(DevBase):
unit_dicts = obj(all=obj())
def __init__(self,type,name,symbol,value,shift=0):
self.type = type
self.name = name
self.symbol = symbol
self.type = type
self.value = value
self.shift = shift
self.value = value
self.shift = shift
ud = Unit.unit_dicts
if type not in ud:
ud[type] = obj()
#end if
td = ud[type]
td[symbol] = self
ud.all[symbol] = self
#end def __init__
#end class Unit
class UnitConverter:
class UnitConverter(DevBase):
unassigned = None
kb = 1.3806503e-23 #J/K
count_set = set(['mol'])
# count
mol = 6.0221415e23
distance_set = set(['m','A','B','nm','pm','fm','a','b','c'])
# distance
m = 1.e0
A = 1.e-10*m
B = .52917720859e-10*m
@ -52,21 +65,21 @@ class UnitConverter:
pm = 1.e-12*m
fm = 1.e-15*m
time_set = set(['s','ms','ns','ps','fs'])
# time
s = 1.e0
ms = 1.e-3*s
ns = 1.e-9*s
ps = 1.e-12*s
fs = 1.e-15*s
mass_set = set(['kg','me','mp','amu','Da'])
# mass
kg = 1.e0
me = 9.10938291e-31*kg
mp = 1.672621777e-27*kg
amu = 1.660538921e-27*kg
Da = amu
energy_set = set(['J','eV','Ry','Ha','kJ_mol','K','degC','degF','kcal_mol'])
# energy
J = 1.e0
eV = 1.60217646e-19*J
Ry = 13.6056923*eV
@ -80,166 +93,85 @@ class UnitConverter:
degC_shift = -273.15
degF_shift = -459.67
charge_set = set(['C','e'])
# charge
C = 1.e0
e = 1.60217646e-19*C
pressure_set = set(['Pa','bar','Mbar','GPa','atm'])
# pressure
Pa = 1.e0
bar = 1.e5*Pa
Mbar = 1.e6*bar
GPa = 1.e9*Pa
atm = 1.01325e5*Pa
force_set = set(['N','pN'])
# force
N = 1.e0
pN = 1e-12*N
therm_cond_set = set(['W_mK'])
# thermal conductivity
W_mK = 1.0
alatt = unassigned
blatt = unassigned
clatt = unassigned
meter = Unit('meter' ,'m' ,'distance',m)
angstrom = Unit('Angstrom' ,'A' ,'distance',A)
bohr = Unit('Bohr' ,'B' ,'distance',B)
nanometer = Unit('nanometer' ,'nm','distance',nm)
picometer = Unit('picometer' ,'pm','distance',pm)
femtometer = Unit('femtometer','pm','distance',fm)
a = Unit('a' ,'a' ,'distance',alatt)
b = Unit('b' ,'b' ,'distance',blatt)
c = Unit('c' ,'c' ,'distance',clatt)
second = Unit('second' ,'s' ,'time',s )
millisecond = Unit('millisecond','ms','time',ms)
nanosecond = Unit('nanosecond' ,'ns','time',ns)
picosecond = Unit('picosecond' ,'ps','time',ps)
femtosecond = Unit('femtosecond','fs','time',fs)
kilogram = Unit('kilogram' ,'kg' ,'mass',kg)
electron_mass = Unit('electron mass' ,'me' ,'mass',me)
proton_mass = Unit('proton mass' ,'mp' ,'mass',mp)
atomic_mass_unit = Unit('atomic mass unit','amu','mass',amu)
dalton = Unit('Dalton' ,'Da' ,'mass',Da)
meter = Unit('distance' ,'meter' ,'m' ,m )
angstrom = Unit('distance' ,'Angstrom' ,'A' ,A )
bohr = Unit('distance' ,'Bohr' ,'B' ,B )
nanometer = Unit('distance' ,'nanometer' ,'nm' ,nm )
picometer = Unit('distance' ,'picometer' ,'pm' ,pm )
femtometer = Unit('distance' ,'femtometer' ,'pm' ,fm )
second = Unit('time' ,'second' ,'s' ,s )
millisecond = Unit('time' ,'millisecond' ,'ms' ,ms )
nanosecond = Unit('time' ,'nanosecond' ,'ns' ,ns )
picosecond = Unit('time' ,'picosecond' ,'ps' ,ps )
femtosecond = Unit('time' ,'femtosecond' ,'fs' ,fs )
kilogram = Unit('mass' ,'kilogram' ,'kg' ,kg )
electron_mass = Unit('mass' ,'electron mass' ,'me' ,me )
proton_mass = Unit('mass' ,'proton mass' ,'mp' ,mp )
atomic_mass_unit = Unit('mass' ,'atomic mass unit','amu' ,amu )
dalton = Unit('mass' ,'Dalton' ,'Da' ,Da )
joule = Unit('Joule' ,'J' ,'energy',J)
electron_volt = Unit('electron Volt','eV' ,'energy',eV)
rydberg = Unit('Rydberg' ,'Ry' ,'energy',Ry)
hartree = Unit('Hartree' ,'Ha' ,'energy',Ha)
kJ_mole = Unit('kJ_mole' ,'kJ_mol','energy',kJ_mol)
kcal_mole = Unit('kcal_mole' ,'kcal_mol','energy',kcal_mol)
kelvin = Unit('Kelvin' ,'K' ,'energy',K)
celcius = Unit('Celcius' ,'degC' ,'energy',degC,degC_shift)
fahrenheit = Unit('Fahrenheit' ,'degF' ,'energy',degF,degF_shift)
joule = Unit('energy' ,'Joule' ,'J' ,J )
electron_volt = Unit('energy' ,'electron Volt' ,'eV' ,eV )
rydberg = Unit('energy' ,'Rydberg' ,'Ry' ,Ry )
hartree = Unit('energy' ,'Hartree' ,'Ha' ,Ha )
kJ_mole = Unit('energy' ,'kJ_mole' ,'kJ_mol' ,kJ_mol )
kcal_mole = Unit('energy' ,'kcal_mole' ,'kcal_mol',kcal_mol )
kelvin = Unit('energy' ,'Kelvin' ,'K' ,K )
celcius = Unit('energy' ,'Celcius' ,'degC' ,degC,degC_shift)
fahrenheit = Unit('energy' ,'Fahrenheit' ,'degF' ,degF,degF_shift)
coulomb = Unit('charge' ,'Coulomb' ,'C' ,C )
electron_charge = Unit('charge' ,'electron charge' ,'e' ,e )
coulomb = Unit('Coulomb' ,'C','charge',C)
electron_charge = Unit('electron charge','e','charge',e)
pascal = Unit('Pascal' ,'Pa' ,'pressure',Pa)
bar = Unit('bar' ,'bar' ,'pressure',bar)
megabar = Unit('megabar' ,'Mbar','pressure',Mbar)
gigapascal = Unit('gigaPascal','Gpa' ,'pressure',GPa)
atmosphere = Unit('atmosphere','atm' ,'pressure',atm)
newton = Unit('Newton' ,'N' ,'force',N)
piconewton = Unit('picoNewton','pN','force',pN)
W_per_mK = Unit('W/(m*K)','W_mK','thermal_cond',W_mK)
pascal = Unit('pressure' ,'Pascal' ,'Pa' ,Pa )
bar = Unit('pressure' ,'bar' ,'bar' ,bar )
megabar = Unit('pressure' ,'megabar' ,'Mbar' ,Mbar )
gigapascal = Unit('pressure' ,'gigaPascal' ,'Gpa' ,GPa )
atmosphere = Unit('pressure' ,'atmosphere' ,'atm' ,atm )
newton = Unit('force' ,'Newton' ,'N' ,N )
piconewton = Unit('force' ,'picoNewton' ,'pN' ,pN )
W_per_mK = Unit('therm_cond','W/(m*K)' ,'W_mK' ,W_mK )
distance_dict = dict([('A',angstrom),\
('B',bohr),\
('a',a),\
('b',b),\
('c',c),\
('m',meter),\
('nm',nanometer),\
('pm',picometer),\
('fm',femtometer),\
])
distance_dict = Unit.unit_dicts.distance
mass_dict = Unit.unit_dicts.mass
energy_dict = Unit.unit_dicts.energy
charge_dict = Unit.unit_dicts.charge
pressure_dict = Unit.unit_dicts.pressure
force_dict = Unit.unit_dicts.force
therm_cond_dict = Unit.unit_dicts.therm_cond
mass_dict = dict([ ('kg',kilogram),\
('me',electron_mass),\
('mp',proton_mass),\
('amu',atomic_mass_unit),\
('Da',dalton),\
])
energy_dict = dict([\
('J',joule),\
('eV',electron_volt),\
('Ry',rydberg),\
('Ha',hartree),\
('kJ_mol',kJ_mole),\
('kcal_mol',kcal_mole),\
('K',kelvin),\
('degC',celcius),\
('degF',fahrenheit),\
])
charge_dict = dict([\
('C',coulomb),\
('e',electron_charge),\
])
pressure_dict = dict([\
('Pa',pascal),\
('bar',bar),\
('Mbar',megabar),\
('GPa',gigapascal),\
('atm',atmosphere),\
])
force_dict = dict([\
('N',newton),\
('pN',piconewton),\
])
therm_cond_dict = dict([\
('W_mK',W_per_mK),\
])
unit_dict = dict([ ('A',angstrom),\
('B',bohr),\
('a',a),\
('b',b),\
('c',c),\
('m',meter),\
('nm',nanometer),\
('pm',picometer),\
('fm',femtometer),\
('kg',kilogram),\
('me',electron_mass),\
('mp',proton_mass),\
('amu',atomic_mass_unit),\
('J',joule),\
('eV',electron_volt),\
('Ry',rydberg),\
('Ha',hartree),\
('kJ_mol',kJ_mole),\
('kcal_mol',kcal_mole),\
('K',kelvin),\
('degC',celcius),\
('degF',fahrenheit),\
('C',coulomb),\
('e',electron_charge),\
('Pa',pascal),\
('bar',bar),\
('Mbar',megabar),\
('GPa',gigapascal),\
('atm',atmosphere),\
('N',newton),\
('pN',piconewton),\
('W_mK',W_per_mK),\
])
unit_dict = Unit.unit_dicts.all
def __init(self):
None
self.error('UnitConverter should not be instantiated')
#def __init__
@staticmethod
def convert(value,source_unit,target_unit):
ui = UnitConverter.unit_dict[source_unit]
uo = UnitConverter.unit_dict[target_unit]
@ -253,11 +185,9 @@ class UnitConverter:
#end if
return (value_out,target_unit)
#end def convert
convert = Callable(convert)
@staticmethod
def convert_scalar_to_all(units,value_orig):
unit_type = UnitConverter.unit_dict[units].type
@ -272,11 +202,10 @@ class UnitConverter:
return value
#end def convert_scalar_to_all
convert_scalar_to_all = Callable(convert_scalar_to_all)
@staticmethod
def convert_array_to_all(units,arr,vectors=None):
A = dict()
A['orig'] = arr.copy()
@ -293,7 +222,7 @@ class UnitConverter:
arr_use = arr
ui = UnitConverter.unit_dict[units]
if(vectors!=None):
A['lattice'] = dot(arr,linalg.inv(vectors.A[units]))
A['lattice'] = dot(arr,inv(vectors.A[units]))
#end if
#end if
@ -311,12 +240,6 @@ class UnitConverter:
return A
#end def convert_array_to_all
convert_array_to_all = Callable(convert_array_to_all)
def submit_unit(uo):
UnitConverter.unit_dict[uo.symbol] = uo
#end def submit_unit
submit_unit = Callable(submit_unit)
#end class UnitConverter

8
nexus/library/vasp.py
View File

@ -28,7 +28,7 @@ import os
from generic import obj
from simulation import Simulation,SimulationInput,SimulationAnalyzer
from vasp_input import VaspInput,generate_vasp_input,generate_poscar,Poscar
from vasp_analyzer import VaspAnalyzer
from vasp_analyzer import VaspAnalyzer,read_vxml
from structure import Structure
from debug import *
@ -114,10 +114,10 @@ class Vasp(Simulation):
input = self.input
if result_name=='structure':
if input.performing_neb():
if 'neb_structures' not in self.temp:
self.temp.neb_structures = []
if 'neb_structures' not in self:
self.neb_structures = []
#end if
neb_structures = self.temp.neb_structures
neb_structures = self.neb_structures
if len(neb_structures)>1:
self.error('NEB simulation at {0} depends on more than two structures\n please check your inputs'.format(self.locdir))
#end if

17
nexus/library/vasp_input.py
View File

@ -41,6 +41,7 @@ import os
from numpy import array,abs,empty,ndarray
from generic import obj
from periodic_table import is_element
from nexus_base import nexus_noncore
from simulation import SimulationInput
from structure import interpolate_structures,Structure
from physical_system import PhysicalSystem
@ -339,18 +340,18 @@ class VKeywordFile(VFile):
@classmethod
def class_init(cls):
for kw_field in cls.kw_fields:
if not kw_field in cls.__dict__:
cls.__dict__[kw_field] = set()
if not cls.class_has(kw_field):
cls.class_set_single(kw_field,set())
#end if
#end for
#cls.check_consistency()
cls.scalar_keywords = set()
for scalar_field in cls.kw_scalars:
cls.scalar_keywords |= cls.__dict__[scalar_field]
cls.scalar_keywords |= cls.class_get(scalar_field)
#end for
cls.array_keywords = set()
for array_field in cls.kw_arrays:
cls.array_keywords |= cls.__dict__[array_field]
cls.array_keywords |= cls.class_get(array_field)
#end for
cls.keywords = cls.scalar_keywords | cls.array_keywords
cls.type = obj()
@ -358,7 +359,7 @@ class VKeywordFile(VFile):
cls.write_value = obj()
cls.assign_value = obj()
for type in cls.kw_scalars + cls.kw_arrays:
for name in cls.__dict__[type]:
for name in cls.class_get(type):
cls.type[name] = type
cls.read_value[name] = read_value_functions[type]
cls.write_value[name] = write_value_functions[type]
@ -374,14 +375,14 @@ class VKeywordFile(VFile):
types = cls.kw_scalars+cls.kw_arrays
untyped = cls.keywords
for type in types:
untyped -= cls.__dict__[type]
untyped -= cls.class_get(type)
#end for
if len(untyped)>0:
fail = True
msg += 'variables without a type: {0}\n'.format(sorted(untyped))
#end if
for type in types:
unknown = cls.__dict__[type]-cls.keywords
unknown = cls.class_get(type)-cls.keywords
if len(unknown)>0:
fail = True
msg += 'unknown {0}: {1}\n'.format(type,sorted(unknown))
@ -1348,7 +1349,7 @@ class VaspInput(SimulationInput,Vobj):
ordered_pseudos.append(pseudo_map[symbol])
#end for
#end if
self.potcar = Potcar(VaspInput.pseudo_dir,ordered_pseudos)
self.potcar = Potcar(nexus_noncore.pseudo_dir,ordered_pseudos)
#end def set_potcar

206
nexus/library/wfconvert.py
View File

@ -1,206 +0,0 @@
##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# wfconvert.py #
# Nexus interface to the wfconvert tool packaged with QMCPACK. #
# #
# Content summary: #
# WfconvertInput #
# SimulationInput class for wfconvert. #
# #
# WfconvertAnalyzer #
# SimulationAnalyzer class for wfconvert. #
# #
# Wfconvert #
# Simulation class for wfconvert. #
# #
# generate_wfconvert #
# User-facing function to generate wfconvert simulation objects.#
# #
#====================================================================#
import os
from generic import obj
from simulation import Simulation,SimulationInput,SimulationAnalyzer
class WfconvertInput(SimulationInput):
def __init__(self,app_name='wfconvert',h5in='MISSING.h5',h5out='wfconvert.h5',spline=False,format='eshdf',factor=None):
self.app_name = app_name
self.h5in = h5in
self.h5out= h5out
self.spline = spline
self.format = format
self.factor = factor
#end def __init__
#wfconvert --nospline --eshdf diamond.h5 out/diamond.pwscf.h5 >& diamond-wfconvert.out
def set_app_name(self,app_name):
self.app_name = app_name
#end def set_app_name
def app_command(self):
c = self.app_name+' '
if not self.spline:
c+= '--nospline '
#end if
c+='--'+self.format+' '+self.h5out+' '+self.h5in
return c
#end def app_command
def read(self,filepath):
None
#end def read
def write_contents(self):
return self.app_command()
#end def write_contents
#end class WfconvertInput
def generate_wfconvert_input(app_name='wfconvert',h5in='MISSING.h5',h5out='wfconvert.h5',spline=False,format='eshdf',factor=None):
wi = WfconvertInput(
app_name = app_name,
h5in = h5in,
h5out = h5out,
spline = spline,
format = format,
factor = factor
)
return wi
#end def generate_wfconvert_input
class WfconvertAnalyzer(SimulationAnalyzer):
def __init__(self,arg0):
if isinstance(arg0,Simulation):
sim = arg0
self.infile = sim.infile
self.dir = sim.locdir
self.h5file = os.path.join(sim.locdir,sim.input.h5out)
else:
self.infile = arg0
#end if
#end def __init__
def analyze(self):
if False:
import h5py
self.log('Fixing h5 file',n=5)
h = h5py.File(self.h5file)
if 'electrons' in h:
elec = h['electrons']
nkpoints = 0
for name,val in elec.iteritems():
if name.startswith('kpoint'):
nkpoints+=1
#end for
#end if
nkold = elec['number_of_kpoints'][0]
self.log('Were',nkold,'kpoints, now',nkpoints,'kpoints',n=6)
elec['number_of_kpoints'][0] = nkpoints
#end for
#end if
#end def analyze
#end class WfconvertAnalyzer
class Wfconvert(Simulation):
input_type = WfconvertInput
analyzer_type = WfconvertAnalyzer
generic_identifier = 'wfconvert'
application = 'wfconvert'
application_properties = set(['serial'])
application_results = set(['orbitals'])
def set_app_name(self,app_name):
self.app_name = app_name
self.input.set_app_name(app_name)
#end def set_app_name
def check_result(self,result_name,sim):
calculating_result = False
if result_name=='orbitals':
calculating_result = True
else:
calculating_result = False
self.error('ability to check for result '+result_name+' has not been implemented')
#end if
return calculating_result
#end def check_result
def get_result(self,result_name,sim):
result = obj()
if result_name=='orbitals':
result.h5file = os.path.join(self.locdir,self.input.h5out)
result.outfile = os.path.join(self.locdir,self.outfile)
else:
self.error('ability to get result '+result_name+' has not been implemented')
#end if
return result
#end def get_result
def incorporate_result(self,result_name,result,sim):
if result_name=='orbitals':
self.input.h5in = os.path.relpath(result.h5file,self.locdir)
self.job.app_command = self.input.app_command()
else:
self.error('ability to incorporate result '+result_name+' has not been implemented')
#end if
#end def incorporate_result
def check_sim_status(self):
outfile = os.path.join(self.locdir,self.outfile)
errfile = os.path.join(self.locdir,self.errfile)
fobj = open(outfile,'r')
output = fobj.read()
fobj.close()
fobj = open(errfile,'r')
errors = fobj.read()
fobj.close()
h5file = os.path.join(self.locdir,self.input.h5out)
file_exists = os.path.exists(h5file)
outfin = 'Successfully read' in errors and 'numSpins' in errors
outfin = outfin and 'Writing laplacians' in output
success = file_exists and outfin
self.finished = success
#end def check_sim_status
def get_output_files(self):
output_files = []
return output_files
#end def get_output_files
def app_command(self):
# app_name is passed along in post_init
return self.input.app_command()
#end def app_command
#end class Wfconvert
def generate_wfconvert(**kwargs):
sim_args,inp_args = Simulation.separate_inputs(kwargs)
if not 'input' in sim_args:
sim_args.input = generate_wfconvert_input(**inp_args)
#end if
wfconvert = Wfconvert(**sim_args)
return wfconvert
#end def generate_wfconvert

6
nexus/library/xmlreader.py
View File

@ -38,11 +38,11 @@ from superstring import \
valid_variable_name,\
string2val
#from abilities import AllAbilities
from generic import obj
from developer import DevBase
class XMLelement(obj):
class XMLelement(DevBase):
def _escape_name(self,name):
if name in self._escape_names:
name=name+'_'
@ -214,7 +214,7 @@ class XMLelement(obj):
class XMLReader
reads an xml file and creates a dynamic object out of its contents
'''
class XMLreader(obj):
class XMLreader(DevBase):
def __init__(self,fpath=None,element_joins=None,element_aliases=None,contract_names=False,strip_prefix=None,warn=True,xml=None):
if element_joins is None:
element_joins = []