diff --git a/config/BGQ_Clang++11_ToolChain.cmake b/config/BGQ_Clang++11_ToolChain.cmake new file mode 100644 index 000000000..f1581ef00 --- /dev/null +++ b/config/BGQ_Clang++11_ToolChain.cmake @@ -0,0 +1,62 @@ +set(CMAKE_C_COMPILER mpiclang) +set(CMAKE_CXX_COMPILER mpiclang++11) + +set(GNU_OPTS "-O3 -g -ffast-math -fopenmp -fstrict-aliasing -Wno-deprecated -Wswitch -Wunused-value") +set(GNU_FLAGS "-Drestrict=__restrict__ -DADD_ -DINLINE_ALL=inline -DDISABLE_TIMER=1 -DHAVE_MASS -DHAVE_MASSV -DUSE_REAL_STRUCT_FACTOR -DSPLINEFLOAT -DBGQPX") +set(CMAKE_CXX_FLAGS "${GNU_FLAGS} ${GNU_OPTS} -ftemplate-depth-60") +set(CMAKE_C_FLAGS "${GNU_FLAGS} ${GNU_OPTS} -std=c99" ) +SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Wl,--allow-multiple-definition") + +SET(QMC_CUDA 0) +SET(QMC_COMPLEX 0) +SET(ENABLE_OPENMP 1) +SET(HAVE_MPI 1) +set(HAVE_CUDA 0) +SET(QMC_BUILD_STATIC 1) +SET(HAVE_LIBESSL 1) +SET(HAVE_EINSPLINE 1) +SET(HAVE_EINSPLINE_EXT 0) +SET(HAVE_ADIOS 0) +SET(BUILD_QMCTOOLS 1) +#SET(BUILD_SANDBOX 1) + +SET(MPIEXEC "sh") +SET(MPIEXEC_NUMPROC_FLAG "${qmcpack_SOURCE_DIR}/utils/bgrunjobhelper.sh") +SET(QE_BIN /soft/applications/quantum_espresso/5.3.0-bgq-omp/bin) + +SET(BOOST_ROOT /home/projects/qmcpack/boost_1_45_0) + +SET(CMAKE_FIND_ROOT_PATH + /home/projects/qmcpack/libXML2-2.9.1 + /soft/libraries/hdf5/1.8.14/cnk-gcc/current + /soft/libraries/alcf/current/gcc/FFTW3 + /soft/libraries/alcf/current/gcc/ZLIB +) + +include_directories($ENV{IBM_MAIN_DIR}/xlmass/bg/7.3/include) + +SET(LAPACK_LIBRARY /soft/libraries/alcf/current/xl/LAPACK/lib/liblapack.a) +SET(BLAS_LIBRARY /soft/libraries/essl/current/essl/5.1/lib64/libesslbg.a) +SET(FORTRAN_LIBRARIES +$ENV{IBM_MAIN_DIR}/xlmass/bg/7.3/bglib64/libmass.a +$ENV{IBM_MAIN_DIR}/xlmass/bg/7.3/bglib64/libmassv.a +$ENV{IBM_FCMP_DIR}/bglib64/libxlf90_r.a +$ENV{IBM_FCMP_DIR}/bglib64/libxlopt.a +$ENV{IBM_FCMP_DIR}/bglib64/libxl.a +/soft/compilers/bgclang/xlsmp-nonconflicting/ibmcmp-feb2015/libxlsmp.a +) + +#LINK_LIBRARIES( +#/soft/perftools/hpctw/libmpihpm_smp.a +#/bgsys/drivers/ppcfloor/bgpm/lib/libbgpm.a +#/bgsys/drivers/ppcfloor/spi/lib/libSPI_upci_cnk.a +#-pg +#) + +FOREACH(type SHARED_LIBRARY SHARED_MODULE EXE) + SET(CMAKE_${type}_LINK_STATIC_C_FLAGS "-Wl,-Bstatic") + SET(CMAKE_${type}_LINK_DYNAMIC_C_FLAGS "-Wl,-Bstatic") + SET(CMAKE_${type}_LINK_STATIC_CXX_FLAGS "-Wl,-Bstatic") + SET(CMAKE_${type}_LINK_DYNAMIC_CXX_FLAGS "-Wl,-Bstatic") +ENDFOREACH(type) + diff --git a/config/BGQToolChain.cmake b/config/BGQ_XL_ToolChain.cmake similarity index 100% rename from config/BGQToolChain.cmake rename to config/BGQ_XL_ToolChain.cmake diff --git a/manual/installation.tex b/manual/installation.tex index def896f3d..274d8f6bb 100644 --- a/manual/installation.tex +++ b/manual/installation.tex @@ -502,8 +502,8 @@ Mira/Cetus is a Blue Gene/Q supercomputer at Argonne National Laboratory's Argon \begin{verbatim} cd build -cmake -DCMAKE_TOOLCHAIN_FILE=../config/BGQToolChain.cmake .. -cmake -DCMAKE_TOOLCHAIN_FILE=../config/BGQToolChain.cmake .. +cmake -DCMAKE_TOOLCHAIN_FILE=../config/BGQ_XL_ToolChain.cmake .. +cmake -DCMAKE_TOOLCHAIN_FILE=../config/BGQ_XL_ToolChain.cmake .. make -j 16 ls -l bin/qmcpack \end{verbatim} @@ -770,8 +770,10 @@ unit = 2.20 sec Total Test time (real) = 2.31 sec \end{verbatim} -Individual unit test executables can be found in build/tests/bin. -The source for the unit tests is located in the \texttt{tests} directory under each directory in \texttt{src} (e.g. src/QMCWavefunctions/tests). +Individual unit test executables can be found in \texttt{build/tests/bin}. +The source for the unit tests is located in the \texttt{tests} directory under each directory in \texttt{src} (e.g. \texttt{src/QMCWavefunctions/tests}). + +See Chapter \ref{chap:unit_testing} for more details about unit tests. \subsection{Automatic tests of QMCPACK} diff --git a/manual/qmcpack_manual.pdf b/manual/qmcpack_manual.pdf index ca1aedb4f..39a2a3609 100644 Binary files a/manual/qmcpack_manual.pdf and b/manual/qmcpack_manual.pdf differ diff --git a/manual/qmcpack_manual.tex b/manual/qmcpack_manual.tex index 41c8ef23b..3b1567aa8 100644 --- a/manual/qmcpack_manual.tex +++ b/manual/qmcpack_manual.tex @@ -168,6 +168,8 @@ \input{lab_condensed_matter} + +\input{unit_testing} \input{contributing} \newpage diff --git a/manual/unit_testing.tex b/manual/unit_testing.tex new file mode 100644 index 000000000..9c6107d87 --- /dev/null +++ b/manual/unit_testing.tex @@ -0,0 +1,132 @@ +\chapter{Unit Testing} +\label{chap:unit_testing} + +Unit testing is a standard software engineering practice to aid in ensuring a quality product. A good suite of unit tests provides confidence in refactoring and changing code, provides some documentation on how classes and functions are used, and can drive a more decoupled design. + +If unit tests do not already exist for a section of code, you are encouraged to add them when modifying that section of code. New code additions should also include unit tests. +When possible, fixes for specific bugs should also include a unit test that would have caught the bug. + +\section {Unit testing framework} The Catch framework is used for unit testing. +See the project site for a tutorial and documentation: \url{https://github.com/philsquared/Catch} + +Catch consists solely of header files. It is distributed as a single include file about 400KB in size. In QMCPACK, it is stored in \texttt{external\_codes/catch}. + +\section{Unit test organization} + +The source for the unit tests is located in the \texttt{tests} directory under each directory in \texttt{src} (e.g. \texttt{src/QMCWavefunctions/tests}). +All of the tests in each \texttt{tests} directory get compiled into an executable. +After building the project, the individual unit test executables can be found in \texttt{build/tests/bin}. +For example, the tests in \texttt{src/QMCWavefunctions/tests} are compiled into \texttt{build/tests/bin/test\_wavefunction}. + +All the unit test executables are collected under ctest with the \texttt{unit} label. +When checking the whole code, it's useful to run through cmake (\texttt{cmake -L unit}). +When working on an individual directory, it's useful to run the individual executable. + +Some of the tests reference input files. The unit test CMake setup places those input files in particular locations under the \texttt{tests} directory (e.g. \texttt{tests/xml\_test}). The individual test needs to be run from that directory to find the expected input files. + +Command line options are available on the unit test executables. Some of the more useful ones are +\begin{description} +\item[\texttt{-h}] List command line options. +\item [\texttt{--list-tests}] List all the tests in the executable. +\end{description} + +A test name can be given on the command line to execute just that test. This is useful when iterating +on a particular test, or when running in the debugger. Test names often contain spaces, so most command line environments require enclosing the test name in single or double quotes. + + + +\section{Example} + +The first example is one test from \texttt{src/Numerics/tests/test\_grid\_functor.cpp} + +\begin{minipage}{\linewidth} +\begin{lstlisting}[language=C++,caption={Unit test example using Catch},label=CatchExample,basicstyle=\ttfamily] +TEST_CASE("double_1d_grid_functor", "[numerics]") +{ + LinearGrid grid; + OneDimGridFunctor f(&grid); + + grid.set(0.0, 1.0, 3); + + REQUIRE(grid.size() == 3); + REQUIRE(grid.rmin() == 0.0); + REQUIRE(grid.rmax() == 1.0); + REQUIRE(grid.dh() == Approx(0.5)); + REQUIRE(grid.dr(1) == Approx(0.5)); +} +\end{lstlisting} +\end{minipage} + +The test function declaration is +\texttt{TEST\_CASE("double\_1d\_grid\_functor","[numerics]")}. +The first argument is the test name, and it must be unique in the test suite. +The second argument is an optional list of tags. Each tag is a name surrounded by brackets (\texttt{"[tag1][tag2]"}). It can also be the empty string. + +The \texttt{REQUIRE} macro accepts expressions with C++ comparison operators and records an error if the value of the expression is false. + +Floating point numbers may have small differences due to roundoff, etc. The \texttt{Approx} class adds some tolerance to the comparison. Place it on either side of the comparison (e.g. \texttt{Approx(a) == 0.3} or \texttt{a = Approx(0.3)}). To adjust the tolerance, use the \texttt{epsilon} and \texttt{scale} methods to \texttt{Approx} (\texttt{REQUIRE(Approx(a).epsilon(0.001) = 0.3);}. + +\subsection{Expected output} + +When running the test executables individually, the output of a run with no failures should look like +\begin{shade} +=============================================================================== +All tests passed (26 assertions in 4 test cases) +\end{shade} + +A test with failures will look like + +\begin{minipage}{\linewidth} +\begin{shade} +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +test_numerics is a Catch v1.4.0 host application. +Run with -? for options + +------------------------------------------------------------------------------- +double_1d_grid_functor +------------------------------------------------------------------------------- +/home/user/qmcpack/src/Numerics/tests/test_grid_functor.cpp:29 +............................................................................... + +/home/user/qmcpack/src/Numerics/tests/test_grid_functor.cpp:39: FAILED: + REQUIRE( grid.dh() == Approx(0.6) ) +with expansion: + 0.5 == Approx( 0.6 ) + +=============================================================================== +test cases: 4 | 3 passed | 1 failed +assertions: 25 | 24 passed | 1 failed +\end{shade} +\end{minipage} + + +\section{Adding tests} +There are three scenarios covered here: adding a new test in an existing file, adding a new test file, or adding a new \texttt{test} directory. + +\subsection{Adding a test to existing file} +Copy an existing test, or from the example shown here. Be sure to change the test name. + +\subsection{Adding a test file} +When adding a new test file, +create a file in the test directory, or copy from an existing file. Add the file name to the \texttt{ADD\_EXECUTABLE} in the \texttt{CMakeLists.txt} file in that directory. + +One (and only one) file must define the \texttt{main} function for the test executable by defining \texttt{CATCH\_CONFIG\_MAIN} before including the Catch header. If more than one file defines this value, there will be linking errors about multiply defined values. + +Some of the tests need to shut down MPI properly to avoid extraneous error messages. Those tests include \texttt{Message/catch\_mpi\_main.hpp} instead of defining \texttt{CATCH\_CONFIG\_MAIN}. + + +\subsection{Adding a test directory} +Copy the CMakeLists.txt file from an existing \texttt{tests} directory. +Change the \texttt{SRC\_DIR} name and the files in the \texttt{ADD\_EXECUTABLES} line. The libraries to link in \texttt{TARGET\_LINK\_LIBRARIES} may need to be updated. + +Add the new test directory to \texttt{src/CMakeLists.txt} in the \texttt{BUILD\_UNIT\_TESTS} section near the end. + + +\section{Testing with random numbers} +Many algorithms and parts of the code depend on random numbers, which makes validating the results difficult. +One solution is to verify that certain properties hold for any random number. +This approach is valuable at some levels of testing, but is unsatisfying at the unit test level. + +The \texttt{Utilities} directory contains a 'fake' random number generator that can be used for deterministic tests of these parts of the code. +Currently it outputs a single, fixed value every time it is called, but it could be expanded to produce more varied, but still deterministic, sequences. +See \texttt{src/QMCDrivers/test\_vmc.cpp} for an example of using the fake random number generator. diff --git a/nexus/library/gamess.py b/nexus/library/gamess.py index df432d5e7..82a577a53 100644 --- a/nexus/library/gamess.py +++ b/nexus/library/gamess.py @@ -183,7 +183,7 @@ class Gamess(Simulation): def generate_gamess(**kwargs): - sim_args,inp_args = Simulation.separate_inputs(kwargs,copy_pseudos=False) + sim_args,inp_args = Gamess.separate_inputs(kwargs,copy_pseudos=False) if not 'input' in sim_args: sim_args.input = generate_gamess_input(**inp_args) diff --git a/nexus/library/generic.py b/nexus/library/generic.py index cedec31a9..8a531cdd5 100644 --- a/nexus/library/generic.py +++ b/nexus/library/generic.py @@ -220,6 +220,72 @@ class object_interface(object): return eq #end def __eq__ + def tree(self,depth=None,all=False,types=False,nindent=1): + if depth==nindent-1: + return '' + #end if + pad = ' ' + npad = nindent*pad + s='' + normal = [] + qable = [] + for k,v in self._iteritems(): + if not isinstance(k,str) or k[0]!='_': + if isinstance(v,object_interface): + qable.append(k) + else: + normal.append(k) + #end if + #end if + #end for + normal.sort() + qable.sort() + indent = npad+18*' ' + if all: + for k in normal: + v = self[k] + if types: + s+=npad+'{0:<15} = '.format(k) + if hasattr(v,'__class__'): + s+='{0:<20}'.format(v.__class__.__name__) + else: + s+='{0:<20}'.format(type(v)) + #end if + else: + s+=npad+str(k) + #end if + s+='\n' + #end for + #end if + if all and depth!=nindent: + for k in qable: + v = self[k] + s+=npad+str(k)+'\n' + s+=v.tree(depth,all,types,nindent+1) + if isinstance(k,str): + s+=npad+'end '+k+'\n' + #end if + #end for + else: + for k in qable: + v = self[k] + if types: + s+=npad+'{0:<15} = '.format(k) + if hasattr(v,'__class__'): + s+='{0:<20}'.format(v.__class__.__name__) + else: + s+='{0:<20}'.format(type(v)) + #end if + else: + s+=npad+str(k) + #end if + s+='\n' + s+=v.tree(depth,all,types,nindent+1) + #end for + #end if + return s + #end def tree + # dict interface def keys(self): @@ -590,6 +656,39 @@ class obj(object_interface): o[path[-1]] = value #end def set_path + def get_path(self,path,value=None): + o = self + if isinstance(path,str): + path = path.split('/') + #end if + for p in path[0:-1]: + if not p in o: + return value + #end if + o = o[p] + #end for + lp = path[-1] + if lp not in o: + return value + else: + return o[lp] + #end if + #end def get_path + + def path_exists(self,path): + o = self + if isinstance(path,str): + path = path.split('/') + #end if + for p in path: + if not p in o: + return False + #end if + o = o[p] + #end for + return True + #end def path_exists + def get(self,key,value=None): # follow dict interface, no plural if key in self: value = self[key] diff --git a/nexus/library/machines.py b/nexus/library/machines.py index f002ed103..0547beabf 100644 --- a/nexus/library/machines.py +++ b/nexus/library/machines.py @@ -210,7 +210,8 @@ class Job(NexusCore): procs = None, processes = None, processes_per_proc = None, - processes_per_node = None + processes_per_node = None, + fake = False, ): self.directory = directory @@ -259,6 +260,7 @@ class Job(NexusCore): self.overtime = False self.successful = False self.finished = False + self.fake_job = fake if app != None: self.app_name = app @@ -1364,11 +1366,14 @@ class Supercomputer(Machine): def process_job(self,job): + if job.fake_job: + return + #end if job.subfile = job.name+'.'+self.sub_launcher+'.in' no_cores = job.cores==None no_nodes = job.nodes==None if no_cores and no_nodes: - self.error('job did not specify cores or nodes\n At least one must be provided') + self.error('job did not specify cores or nodes\nAt least one must be provided') elif no_cores: job.cores = self.cores_per_node*job.nodes elif no_nodes: @@ -1379,7 +1384,11 @@ class Supercomputer(Machine): else: job.cores = min(job.cores,job.nodes*self.cores_per_node) #end if - job.processes = max(1,int(float(job.cores)/job.threads)) + if job.processes_per_node!=None: + job.processes=job.nodes*job.processes_per_node + else: + job.processes = max(1,int(float(job.cores)/job.threads)) + #end if job.tot_cores = job.nodes*self.cores_per_node job.procs = job.nodes*self.procs_per_node @@ -2500,6 +2509,41 @@ class Mira(ALCF_Machine): #end class Mira +class Cooley(Supercomputer): + name = 'cooley' + requires_account = True + batch_capable = True + executable_subfile = True + + prefixed_output = True + outfile_extension = '.output' + errfile_extension = '.error' + + def process_job_extra(self,job): + if job.processes_per_node is None and job.threads!=1: + self.error('threads must be 1,2,3,4,6, or 12 on Cooley\nyou provided: {0}'.format(job.threads)) + #end if + + job.run_options.add( + f = '-f $COBALT_NODEFILE', + ppn = '-ppn {0}'.format(job.processes_per_node), + ) + #end def process_job_extra + + def write_job_header(self,job): + if job.queue is None: + job.queue = 'default' + #end if + c= '#!/bin/bash\n' + c+='#COBALT -q {0}\n'.format(job.queue) + c+='#COBALT -A {0}\n'.format(job.account) + c+='#COBALT -n {0}\n'.format(job.nodes) + c+='#COBALT -t {0}\n'.format(job.total_minutes()) + c+='#COBALT -O {0}\n'.format(job.identifier) + return c + #end def write_job_header +#end class Cooley + class Lonestar(Supercomputer): # Lonestar contribution from Paul Young @@ -2677,6 +2721,7 @@ EOS( 744, 2, 8, 64, 1000, 'aprun', 'qsub', 'qstat', 'q Vesta( 2048, 1, 16, 16, 10, 'runjob', 'qsub', 'qstata', 'qdel') Cetus( 1024, 1, 16, 16, 10, 'runjob', 'qsub', 'qstata', 'qdel') Mira( 49152, 1, 16, 16, 10, 'runjob', 'qsub', 'qstata', 'qdel') +Cooley( 126, 2, 6, 384, 10, 'mpirun', 'qsub', 'qstata', 'qdel') Lonestar( 22656, 2, 6, 12, 128, 'ibrun', 'qsub', 'qstat', 'qdel') Matisse( 20, 2, 8, 64, 2, 'mpirun', 'sbatch', 'sacct', 'scancel') Komodo( 24, 2, 6, 48, 2, 'mpirun', 'sbatch', 'sacct', 'scancel') diff --git a/nexus/library/nexus.py b/nexus/library/nexus.py index bd8308a32..4d95c719d 100755 --- a/nexus/library/nexus.py +++ b/nexus/library/nexus.py @@ -33,7 +33,7 @@ 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 pseudopotential import Pseudopotential,Pseudopotentials,ppset from basisset import BasisSets from bundle import bundle @@ -101,10 +101,14 @@ class Settings(NexusCore): ericfmt mcppath '''.split()) + pwscf_vars = set(''' + vdw_table + '''.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 + allowed_vars = nexus_vars | machine_vars | gamess_vars | pwscf_vars @staticmethod @@ -167,6 +171,7 @@ class Settings(NexusCore): 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) + pwscf_kw = Settings.kw_set(Settings.pwscf_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 @@ -175,6 +180,7 @@ class Settings(NexusCore): # copy input settings self.transfer_from(mach_kw.copy()) self.transfer_from(gamess_kw.copy()) + self.transfer_from(pwscf_kw.copy()) # process machine settings self.process_machine_settings(mach_kw) @@ -197,6 +203,9 @@ class Settings(NexusCore): # process gamess settings Gamess.settings(**gamess_kw) + # process pwscf settings + Pwscf.settings(**pwscf_kw) + return #end def __call__ diff --git a/nexus/library/physical_system.py b/nexus/library/physical_system.py index bf89d91b3..ad9c77d6e 100644 --- a/nexus/library/physical_system.py +++ b/nexus/library/physical_system.py @@ -208,6 +208,18 @@ class Particles(Matter): #end for return nelectrons #end def count_electrons + + def electron_counts(self): + counts = [] + for electron in ('up_electron','down_electron'): + if electron in self: + counts.append(self[electron].count) + else: + counts.append(0) + #end if + #end for + return counts + #end def electron_counts #end class Particles me_amu = convert(1.,'me','amu') @@ -613,7 +625,6 @@ def generate_physical_system(**kwargs): is_str = isinstance(s,str) if is_str and os.path.exists(s):# and '.' in os.path.split(s)[1]: if 'elem' in kwargs: - print 'system using elem' kwargs['structure'] = read_structure(s,elem=kwargs['elem']) else: kwargs['structure'] = read_structure(s) diff --git a/nexus/library/project_manager.py b/nexus/library/project_manager.py index 70dd660f4..7803b7737 100644 --- a/nexus/library/project_manager.py +++ b/nexus/library/project_manager.py @@ -44,10 +44,12 @@ class ProjectManager(NexusCore): simulations = simulations[0] #end if for sim in simulations: - if len(sim.dependencies)==0: - self.add_cascade(sim) + if not sim.fake(): + if len(sim.dependencies)==0: + self.add_cascade(sim) + #end if + self.simulations[sim.simid]=sim #end if - self.simulations[sim.simid]=sim #end for #end def add_simulations @@ -59,6 +61,7 @@ class ProjectManager(NexusCore): def init_cascades(self): + self.screen_fake_sims() self.resolve_file_collisions() self.propagate_blockages() self.log('loading cascade images',n=1) @@ -172,6 +175,24 @@ class ProjectManager(NexusCore): #end def save_cascades + def screen_fake_sims(self): + def collect_fake(sim,fake): + if sim.fake(): + fake.append(sim) + #end if + #end def collect_fake + fake = [] + self.traverse_cascades(collect_fake,fake) + if len(fake)>0: + msg = 'fake/temporary simulation objects detected in cascade\nthis is a developer error\nlist of fake sims and directories:\n' + for sim in fake: + msg +=' {0:>8} {1}\n'.format(sim.simid,sim.locdir) + #end for + self.error(msg) + #end if + #end def screen_fake_sims + + def propagate_blockages(self): def collect_blocked(sim,blocked): if sim.block or sim.block_subcascade: diff --git a/nexus/library/pseudopotential.py b/nexus/library/pseudopotential.py index 137792383..8d5847f82 100644 --- a/nexus/library/pseudopotential.py +++ b/nexus/library/pseudopotential.py @@ -51,6 +51,7 @@ from unit_converter import convert from generic import obj from developer import DevBase,unavailable,error from basisset import process_gaussian_text,GaussianBasisSet +from physical_system import PhysicalSystem from plotting import * from debug import * try: @@ -61,6 +62,26 @@ except: +def pp_elem_label(filename,guard=False): + el = '' + for c in filename: + if c=='.' or c=='_' or c=='-': + break + #end if + el+=c + #end for + elem_label = el + is_elem,symbol = is_element(el,symbol=True) + if guard: + if not is_elem: + error('cannot determine element for pseudopotential file: {0}\npseudopotential file names must be prefixed by an atomic symbol or label\n(e.g. Si, Si1, etc)'.format(filename)) + #end if + return elem_label,symbol + else: + return elem_label,symbol,is_elem + #end if +#end def pp_elem_label + # basic interface for nexus, only gamess really needs this for now class PseudoFile(DevBase): @@ -72,8 +93,9 @@ class PseudoFile(DevBase): if filepath!=None: self.filename = os.path.basename(filepath) self.location = os.path.abspath(filepath) - elem_label = self.filename.split('.')[0] - is_elem,symbol = is_element(elem_label,symbol=True) + elem_label,symbol,is_elem = pp_elem_label(self.filename) + #elem_label = self.filename.split('.')[0] + #is_elem,symbol = is_element(elem_label,symbol=True) if not is_elem: self.error('cannot determine element for pseudopotential file: {0}\npseudopotential file names must be prefixed by an atomic symbol or label\n(e.g. Si, Si1, etc)'.format(filepath)) #end if @@ -223,7 +245,104 @@ class Pseudopotentials(DevBase): +# user interface to group sets of pseudopotentials together and refer to them by labels +# labeling should eliminate the need to provide lists of pseudopotential files to each +# simulation object (e.g. via a generate_* call) separately +class PPset(DevBase): + instance_counter = 0 + known_codes = set('pwscf gamess vasp qmcpack'.split()) + + default_extensions = obj( + pwscf = ['ncpp','upf'], + gamess = ['gms'], + vasp = ['potcar'], + qmcpack = ['xml'], + ) + + def __init__(self): + if PPset.instance_counter!=0: + self.error('cannot instantiate more than one PPset object\nintended use follows a singleton pattern') + #end if + PPset.instance_counter+=1 + self.pseudos = obj() + #end def __init__ + + def supports_code(self,code): + return code in PPset.known_codes + #end def supports_code + + def __call__(self,label,**code_pps): + if not isinstance(label,str): + self.error('incorrect use of ppset\nlabel provided must be a string\nreceived type instead: {0}\nwith value: {1}'.format(label.__class__.__name__,label)) + #end if + if label in self.pseudos: + self.error('incorrect use of ppset\npseudopotentials with label "{0}" have already been added to ppset'.format(label)) + #end if + pseudos = obj() + self.pseudos[label]=pseudos + for code,pps in code_pps.iteritems(): + clow = code.lower() + if clow not in self.known_codes: + self.error('incorrect use of ppset\ninvalid simulation code "{0}" provided with set labeled "{1}"\nknown simulation codes are: {2}'.format(code,label,sorted(self.known_codes))) + #end if + if not isinstance(pps,(list,tuple)): + self.error('incorrect use of ppset\nmust provide a list of pseudopotentials for code "{0}" in set labeled "{1}"\ntype provided instead of list: {2}'.format(code,label,pps.__class__.__name__)) + #end if + ppcoll = obj() + for pp in pps: + if not isinstance(pp,str): + self.error('incorrect use of ppset\nnon-filename provided with set labeled "{0}" for simulation code "{1}"\neach pseudopential file name must be a string\nreceived type: {2}\nwith value: {3}'.format(label,code,pp.__class__.__name__,pp)) + #end if + elem_label,symbol,is_elem = pp_elem_label(pp) + if not is_elem: + self.error('invalid filename provided to ppset\ncannot determine element for pseudopotential file: {0}\npseudopotential file names must be prefixed by an atomic symbol or label\n(e.g. Si, Si1, etc)'.format(pp)) + elif symbol in ppcoll: + self.error('incorrect use of ppset\nmore than one pseudopotential file provided for element "{0}" for code "{1}" in set labeled "{2}"\nfirst file: {3}\nsecond file: {4}'.format(symbol,code,label,ppcoll[symbol],pp)) + #end if + ppcoll[symbol] = pp + #end for + pseudos[clow] = ppcoll + #end for + #end def __call__ + + def has_set(self,label): + return label in self.pseudos + #end def has_set + + def get(self,label,code,system): + if system is None or not system.pseudized: + return [] + #end if + if not isinstance(system,PhysicalSystem): + self.error('system object must be of type PhysicalSystem') + #end if + species_labels,species = system.structure.species(symbol=True) + if not isinstance(label,str): + self.error('incorrect use of ppset\nlabel provided must be a string\nreceived type instead: {0}\nwith value: {1}'.format(label.__class__.__name__,label)) + #end if + if not self.has_set(label): + self.error('incorrect use of ppset\npseudopotential set labeled "{0}" is not present in ppset\nset labels present: {1}\nplease either provide pseudopotentials with label "{0}" or correct the provided label'.format(label,sorted(self.pseudos.keys()))) + #end if + pseudos = self.pseudos[label] + clow = code.lower() + if clow not in self.known_codes: + self.error('simulation code "{0}" is not known to ppset\nknown codes are: {1}'.format(code,sorted(self.known_codes))) + elif clow not in pseudos: + self.error('incorrect use of ppset\npseudopotentials were not provided for simulation code "{0}" in set labeled "{1}"\npseudopotentials are required for physical system with pseudo-elements: {2}\nplease add these pseudopotentials for code "{0}" in set "{1}"'.format(code,label,sorted(species))) + #end if + ppcoll = pseudos[clow] + pps = [] + for symbol in species: + if symbol not in ppcoll: + self.error('incorrect use of ppset\npseudopotentials were not provided for element "{0}" code "{1}" in set labeled "{2}"\nphysical system encountered with pseudo-elements: {3}\nplease ensure that pseudopotentials are provided for these elements in set "{2}" for code "{1}"'.format(symbol,code,label,sorted(species))) + #end if + pps.append(ppcoll[symbol]) + #end for + return pps + #end def get +#end class PPset +ppset = PPset() @@ -326,6 +445,11 @@ class SemilocalPP(Pseudopotential): formats = ['qmcpack'] + channel_indices = obj() + for i,c in enumerate(l_channels): + channel_indices[c] = i + #end for + def __init__(self,filepath=None,format=None,name=None,src=None): self.name = name self.rcut = None @@ -509,7 +633,7 @@ class SemilocalPP(Pseudopotential): if self.name!=None: lab = self.name+' '+lab #end if - v = self.evaluate(r,c,with_local) + v = self.evaluate(r,c,with_local=with_local) if metric=='r2': v = r**2*v if c==self.local: @@ -518,6 +642,14 @@ class SemilocalPP(Pseudopotential): elif metric!=None: self.error('invalid metric for plotting: {0}\nvalid options are: r2'.format(metric)) #end if + #l = self.channel_indices[c] + #if l==0: + # lmult = 1 + #else: + # lmult = (l*(l+1)) + #vs = self.evaluate(r,'s',with_local=with_local) + #v-=vs + #v/=lmult plot(r,v,color+linestyle,label=lab) #end for #end for @@ -535,6 +667,39 @@ class SemilocalPP(Pseudopotential): #end if #end def plot + + def plot_L2(self,show=True,fig=True,r=None,rmin=0.01,rmax=5.0,linestyle='-',title=None): + if r is None and self.numeric: + r = self.r + elif r is None: + r = linspace(rmin,rmax,1000) + #end if + vs = self.evaluate(r,'s',with_local=True) + channels = self.all_channels + for c in channels[2:]: + if c in self.channels: + color = self.channel_colors[c] + v = self.evaluate(r,c,with_local=True) + l = self.channel_indices[c] + vL2 = (v-vs)/(l*(l+1)) + plot(r,vL2,color+linestyle,label='(v{0}-vs)/(l(l+1))'.format(c)) + #end if + #end for + if fig: + xlim([0,rmax]) + if title is None: + title = 'Semilocal {0} PP ({1} core)'.format(self.element,self.core) + #end if + set_title(title) + ylabel('vL2 for channels above s') + xlabel('r') + legend() + if show: + show_plots() + #end if + #end if + #end def plot_L2 + def gaussian_fit(self,r=None,pmax=3,maxfev=100000,verbose=False,filepath=None,format=None,offset=0): if verbose: diff --git a/nexus/library/pwscf.py b/nexus/library/pwscf.py index b0090b13e..23d77afff 100644 --- a/nexus/library/pwscf.py +++ b/nexus/library/pwscf.py @@ -36,6 +36,16 @@ class Pwscf(Simulation): application_properties = set(['serial','mpi']) application_results = set(['charge_density','orbitals','structure']) + + vdw_table = None + + @staticmethod + def settings(vdw_table=None): + # van der Waals family of functional require the vdW table generated by + # generate_vdW_kernel_table.x: specify 'vdw_table' in settings + Pwscf.vdw_table = vdw_table + #end def settings + #def propagate_identifier(self): # self.input.control.prefix = self.identifier ##end def propagate_identifier @@ -53,7 +63,7 @@ class Pwscf(Simulation): self.warn('requested grouping by atomic species, but pwscf does not group atoms anymore!') #self.system.structure.group_atoms() #end if - #end def post_init + #end def __init__ def write_prep(self): @@ -62,6 +72,18 @@ class Pwscf(Simulation): if not os.path.exists(outdir): os.makedirs(outdir) #end if + #copy over vdw_table for vdW-DF functional + if self.path_exists('input/system/input_dft'): + if self.input.system.input_dft.lower().startswith('vdw-df'): + if self.vdw_table is None: + self.error('attempting to run vdW-DF functional, but vdw_table is missing\nplease provide path to table file via "vdw_table" parameter in settings') + #end if + cd_rel = os.path.relpath(self.vdw_table,self.locdir) + # copy instead of link to vdw_table to avoid file-lock from multiple pw.x instances + cp_cmd = 'cd '+self.locdir+';cp '+cd_rel+' .' + os.system(cp_cmd) + #end if + #end if #end def write_prep @@ -196,7 +218,7 @@ class Pwscf(Simulation): def generate_pwscf(**kwargs): - sim_args,inp_args = Simulation.separate_inputs(kwargs) + sim_args,inp_args = Pwscf.separate_inputs(kwargs) if not 'input' in sim_args: input_type = inp_args.input_type diff --git a/nexus/library/pwscf_input.py b/nexus/library/pwscf_input.py index 966cd5998..f0213fd77 100644 --- a/nexus/library/pwscf_input.py +++ b/nexus/library/pwscf_input.py @@ -48,7 +48,7 @@ import os import inspect from copy import deepcopy from superstring import string2val -from numpy import fromstring,empty,array,float64,ones,pi,dot +from numpy import fromstring,empty,array,float64,ones,pi,dot,ceil from numpy.linalg import inv from unit_converter import convert from generic import obj @@ -76,7 +76,7 @@ def read_float(sv): bconv = {'.true.':True,'.false.':False} def read_bool(sv): - return bconv[sv] + return bconv[sv.lower()] #end def read_bool @@ -154,10 +154,69 @@ def array_to_string(a,pad=' ',format=pwscf_array_format,converter=noconv,rowse 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'] - bools=['wf_collect','tstress','tprnfor','lkpoint_dir','tefield','dipfield','lelfield','lberry','nosym','nosym_evc','noinv','force_symmorphic','noncolin','lda_plus_u','lspinorb','do_ee','london','diago_full_acc','tqr','remove_rigid_rot','refold_pos','first_last_opt','use_masses','use_freezing','la2F'] + ints=[ + # pre 5.4 + 'nstep','iprint','gdir','nppstr','nberrycyc','ibrav','nat','ntyp', + 'nbnd','nr1','nr2','nr3','nr1s','nr2s','nr3s','nspin', + 'multiplicity','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', + # 5.4 additions + 'nqx1','nqx2','nqx3','esm_nfit','space_group','origin_choice', + ] + floats=[ + # pre 5.4 + 'dt','max_seconds','etot_conv_thr','forc_conv_thr','celldm','A','B','C', + 'cosAB','cosAC','cosBC','nelec','ecutwfc','ecutrho','degauss', + 'tot_charge','tot_magnetization','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', + # 5.4 additions + 'conv_thr_init','conv_thr_multi','efield_cart','screening_parameter', + 'ecutvcut','Hubbard_J0','Hubbard_beta','Hubbard_J','esm_w', + 'esm_efield','fcp_mu','london_c6','london_rvdw','xdm_a1','xdm_a2', + ] + strs=[ + # pre 5.4 + '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', + # 5.4 additions + 'esm_bc','vdw_corr', + ] + bools=[ + # pre 5.4 + 'wf_collect','tstress','tprnfor','lkpoint_dir','tefield','dipfield', + 'lelfield','lberry','nosym','nosym_evc','noinv','force_symmorphic', + 'noncolin','lda_plus_u','lspinorb','do_ee','london','diago_full_acc', + 'tqr','remove_rigid_rot','refold_pos','first_last_opt','use_masses', + 'use_freezing','la2F', + # 5.4 additions + 'lorbm','lfcpopt','scf_must_converge','adaptive_thr','no_t_rev', + 'use_all_frac','one_atom_occupations','starting_spin_angle', + 'x_gamma_extrapolation','xdm','uniqueb','rhombohedral', + ] + + # real arrays: celldm,starting_magnetization, hubbard_alpha, hubbard_u, + # hubbard_j0, hubbard_beta, hubbard_j, + # starting_ns_eigenvalue, angle1, angle2, fixed_magnetization + # fe_step,efield_cart + # london_c6 london_rvdw + + # species arrays: starting_magnetization, hubbard_alpha, hubbard_u, hubbard_j0, hubbard_beta, hubbard_j, angle1, angle2, london_c6, london_rvdw + # multidimensional arrays: starting_ns_eigenvalue(3), hubbard_j(2) + # all_variables = set(ints+floats+strs+bools) @@ -165,16 +224,16 @@ class PwscfInputBase(DevBase): var_types = dict() for v in ints: - var_types[v]=int + var_types[v.lower()]=int #end for for v in floats: - var_types[v]=float + var_types[v.lower()]=float #end for for v in strs: - var_types[v]=str + var_types[v.lower()]=str #end for for v in bools: - var_types[v]=bool + var_types[v.lower()]=bool #end for #end class PwscfInputBase @@ -204,29 +263,74 @@ class Element(PwscfInputBase): class Section(Element): + @classmethod + def class_init(cls): + cls.varlist = list(cls.variables) + cls.variables = set([v.lower() for v in cls.varlist]) + cls.case_map = obj() + for vname in cls.varlist: + cls.case_map[vname.lower()] = vname + #end for + #end if + def assign(self,**variables): self.transfer_from(variables) #end def assign def read(self,lines): for l in lines: - tokens = l.split(',') + # exclude comments + cloc = l.find('!') + if cloc!=-1: + l = l[:cloc] + #end if + # parse tokens accounting for significant whitespace + # replace commas outside array brackets with spaces + if '(' not in l: + lout = l + else: + lout = '' + inparen=False + for c in l: + if c=='(': + inparen=True + lout+=c + elif c==')': + inparen=False + lout+=c + elif inparen and c==',': + lout+='|' # new delimiter + else: + lout += c + #end if + #end for + #end if + tokens = lout.split(',') for t in tokens: if len(t)>0: - var,val = t.split('=') - var = var.strip() - val = val.strip() - varname = var.split('(')[0] + tsplt = t.split('=') + if len(tsplt)!=2: + self.error('attempted to read misformatted line\nmisformatted line: {0}\ntokens: {1}'.format(l,tsplt)) + #end if + var,val = tsplt + var = var.strip().lower() + val = val.strip() + if '(' not in var: + varname = var + else: + var = var.replace('|',',') + varname = var.split('(')[0] + #end if if not varname in self.variables: - self.error('pwscf input section {0} does not have a variable named {1}, please check your input\nif correct, please add a new variable ({1}) to the {0} PwscfInput class'.format(self.__class__.__name__,varname),trace=False) + self.error('pwscf input section {0} does not have a variable named "{1}", please check your input\nif correct, please add a new variable ({1}) to the {0} PwscfInput class'.format(self.__class__.__name__,varname),trace=False) #end if if not varname in self.var_types: - self.error('a type has not been specified for variable {0}\nplease add it to PwscfInputBase'.format(varname),trace=False) + self.error('a type has not been specified for variable "{0}"\nplease add it to PwscfInputBase'.format(varname),trace=False) #end if vtype = self.var_types[varname] val = readval[vtype](val) self[var] = val - #end if + #end for #end for #end for #end def read @@ -279,7 +383,7 @@ class Card(Element): def get_specifier(self,line): tokens = line.split() if len(tokens)>1: - self.specifier = tokens[1].strip('{}').lower() + self.specifier = tokens[1].strip('{}()').lower() #end if #end def get_specifier @@ -311,32 +415,96 @@ class Card(Element): class control(Section): name = 'control' - variables = set( - ['calculation','title','verbosity','restart_mode','wf_collect','nstep', - 'iprint','tstress','tprnfor','dt','outdir','wfcdir','prefix', - 'lkpoint_dir','max_seconds','etot_conv_thr','forc_conv_thr','disk_io', - 'pseudo_dir','tefield','dipfield','lelfield','lberry','gdir','nppstr', - 'nberrycyc']) -#end class control + # all known keywords + variables = [ + 'calculation','title','verbosity','restart_mode','wf_collect','nstep', + 'iprint','tstress','tprnfor','dt','outdir','wfcdir','prefix', + 'lkpoint_dir','max_seconds','etot_conv_thr','forc_conv_thr','disk_io', + 'pseudo_dir','tefield','dipfield','lelfield','nberrycyc','lorbm','lberry', + 'gdir','nppstr','lfcpopt' + ] + + # 5.4 keyword spec + #variables = [ + # 'calculation','title','verbosity','restart_mode','wf_collect','nstep', + # 'iprint','tstress','tprnfor','dt','outdir','wfcdir','prefix', + # 'lkpoint_dir','max_seconds','etot_conv_thr','forc_conv_thr','disk_io', + # 'pseudo_dir','tefield','dipfield','lelfield','nberrycyc','lorbm','lberry', + # 'gdir','nppstr','lfcpopt' + # ] + + # sometime prior to 5.4 + #variables = [ + # 'calculation','title','verbosity','restart_mode','wf_collect','nstep', + # 'iprint','tstress','tprnfor','dt','outdir','wfcdir','prefix', + # 'lkpoint_dir','max_seconds','etot_conv_thr','forc_conv_thr','disk_io', + # 'pseudo_dir','tefield','dipfield','lelfield','lberry','gdir','nppstr', + # 'nberrycyc' + # ] +#end class control class system(Section): name = 'system' - variables = set( - ['ibrav','celldm','A','B','C','cosAB','cosAC','cosBC','nat','ntyp', - 'nbnd','nelec','tot_charge','ecutwfc','ecutrho','nr1','nr2','nr3', - 'nr1s','nr2s','nr3s','nosym','nosym_evc','noinv','force_symmorphic', - 'occupations','degauss','smearing','nspin','noncolin', - 'starting_magnetization','nelup','neldw','multiplicity', - 'tot_magnetization','ecfixed','qcutz','q2sigma','input_dft', - 'lda_plus_u','Hubbard_alpha','Hubbard_U','starting_ns_eigenvalue', - 'U_projection_type','edir','emaxpos','eopreg','eamp','angle1', - 'angle2','constrained_magnetization','fixed_magnetization','lambda', - 'report','lspinorb','assume_isolated','do_ee','london','london_s6', - 'london_rcut','exx_fraction','ecutfock', - 'lda_plus_u_kind','Hubbard_J','exxdiv_treatment','la2F']) + + # all known keywords + variables = [ + 'ibrav','celldm','A','B','C','cosAB','cosAC','cosBC','nat','ntyp', + 'nbnd','tot_charge','tot_magnetization','starting_magnetization', + 'ecutwfc','ecutrho','ecutfock','nr1','nr2','nr3','nr1s','nr2s','nr3s', + 'nosym','nosym_evc','noinv','no_t_rev','force_symmorphic','use_all_frac', + 'occupations','one_atom_occupations','starting_spin_angle','degauss', + 'smearing','nspin','noncolin','ecfixed','qcutz','q2sigma','input_dft', + 'exx_fraction','screening_parameter','exxdiv_treatment', + 'x_gamma_extrapolation','ecutvcut','nqx1','nqx2','nqx3','lda_plus_u', + 'lda_plus_u_kind','Hubbard_U','Hubbard_J0','Hubbard_alpha', + 'Hubbard_beta','Hubbard_J','starting_ns_eigenvalue','U_projection_type', + 'edir','emaxpos','eopreg','eamp','angle1','angle2', + 'constrained_magnetization','fixed_magnetization','lambda','report', + 'lspinorb','assume_isolated','esm_bc','esm_w','esm_efield','esm_nfit', + 'fcp_mu','vdw_corr','london','london_s6','london_c6','london_rvdw', + 'london_rcut','xdm','xdm_a1','xdm_a2','space_group','uniqueb', + 'origin_choice','rhombohedral', + 'nelec','nelup','neldw','multiplicity','do_ee','la2F', + ] + + # 5.4 keyword spec + #variables = [ + # 'ibrav','celldm','A','B','C','cosAB','cosAC','cosBC','nat','ntyp', + # 'nbnd','tot_charge','tot_magnetization','starting_magnetization', + # 'ecutwfc','ecutrho','ecutfock','nr1','nr2','nr3','nr1s','nr2s','nr3s', + # 'nosym','nosym_evc','noinv','no_t_rev','force_symmorphic','use_all_frac', + # 'occupations','one_atom_occupations','starting_spin_angle','degauss', + # 'smearing','nspin','noncolin','ecfixed','qcutz','q2sigma','input_dft', + # 'exx_fraction','screening_parameter','exxdiv_treatment', + # 'x_gamma_extrapolation','ecutvcut','nqx1','nqx2','nqx3','lda_plus_u', + # 'lda_plus_u_kind','Hubbard_U','Hubbard_J0','Hubbard_alpha', + # 'Hubbard_beta','Hubbard_J','starting_ns_eigenvalue','U_projection_type', + # 'edir','emaxpos','eopreg','eamp','angle1','angle2', + # 'constrained_magnetization','fixed_magnetization','lambda','report', + # 'lspinorb','assume_isolated','esm_bc','esm_w','esm_efield','esm_nfit', + # 'fcp_mu','vdw_corr','london','london_s6','london_c6','london_rvdw', + # 'london_rcut','xdm','xdm_a1','xdm_a2','space_group','uniqueb', + # 'origin_choice','rhombohedral' + # ] + + # sometime prior to 5.4 + #variables = [ + # 'ibrav','celldm','A','B','C','cosAB','cosAC','cosBC','nat','ntyp', + # 'nbnd','nelec','tot_charge','ecutwfc','ecutrho','nr1','nr2','nr3', + # 'nr1s','nr2s','nr3s','nosym','nosym_evc','noinv','force_symmorphic', + # 'occupations','degauss','smearing','nspin','noncolin', + # 'starting_magnetization','nelup','neldw','multiplicity', + # 'tot_magnetization','ecfixed','qcutz','q2sigma','input_dft', + # 'lda_plus_u','Hubbard_alpha','Hubbard_U','starting_ns_eigenvalue', + # 'U_projection_type','edir','emaxpos','eopreg','eamp','angle1', + # 'angle2','constrained_magnetization','fixed_magnetization','lambda', + # 'report','lspinorb','assume_isolated','do_ee','london','london_s6', + # 'london_rcut','exx_fraction','ecutfock', + # 'lda_plus_u_kind','Hubbard_J','exxdiv_treatment','la2F' + # ] atomic_variables = obj( hubbard_u = 'Hubbard_U', @@ -354,9 +522,12 @@ class system(Section): akeys = [] for key in keys: if key.startswith(name): - has_hubbard_u = True akeys.append(key) - index = int(key.replace(name,'').strip('()')) + if '(' not in key: + index=1 + else: + index = int(key.replace(name,'').strip('()')) + #end if avals[index] = self[key] #end if #end for @@ -387,6 +558,9 @@ class system(Section): for var in vars: val = self[var] if var in self.atomic_variables: + if val is None: # jtk mark: patch fix for generate_pwscf_input + continue # i.e. hubbard_u = None + #end if if 'atomic_species' in parent: atoms = parent.atomic_species.atoms avar = self.atomic_variables[var] @@ -438,60 +612,139 @@ class system(Section): class electrons(Section): name = 'electrons' - variables = set( - ['electron_maxstep','conv_thr','mixing_mode','mixing_beta','mixing_ndim', - 'mixing_fixed_ns','diagonalization','ortho_para','diago_thr_init', - 'diago_cg_maxiter','diago_david_ndim','diago_full_acc','efield', - 'startingpot','startingwfc','tqr']) + + # all known keywords + variables = [ + 'electron_maxstep','scf_must_converge','conv_thr','adaptive_thr', + 'conv_thr_init','conv_thr_multi','mixing_mode','mixing_beta', + 'mixing_ndim','mixing_fixed_ns','diagonalization','ortho_para', + 'diago_thr_init','diago_cg_maxiter','diago_david_ndim','diago_full_acc', + 'efield','efield_cart','startingpot','startingwfc','tqr' + ] + + # 5.4 keyword spec + #variables = [ + # 'electron_maxstep','scf_must_converge','conv_thr','adaptive_thr', + # 'conv_thr_init','conv_thr_multi','mixing_mode','mixing_beta', + # 'mixing_ndim','mixing_fixed_ns','diagonalization','ortho_para', + # 'diago_thr_init','diago_cg_maxiter','diago_david_ndim','diago_full_acc', + # 'efield','efield_cart','startingpot','startingwfc','tqr' + # ] + + # sometime prior to 5.4 + #variables = [ + # 'electron_maxstep','conv_thr','mixing_mode','mixing_beta','mixing_ndim', + # 'mixing_fixed_ns','diagonalization','ortho_para','diago_thr_init', + # 'diago_cg_maxiter','diago_david_ndim','diago_full_acc','efield', + # 'startingpot','startingwfc','tqr' + # ] #end class electrons class ions(Section): name = 'ions' - variables = set( - ['ion_dynamics','ion_positions','phase_space','pot_extrapolation', - 'wfc_extrapolation','remove_rigid_rot','ion_temperature','tempw', - 'tolp','delta_t','nraise','refold_pos','upscale','bfgs_ndim', - 'trust_radius_max','trust_radius_min','trust_radius_ini','w_1','w_2', - 'num_of_images','opt_scheme','CI_scheme','first_last_opt','temp_req', - 'ds','k_max','k_min','path_thr','use_masses','use_freezing','fe_step', - 'g_amplitude','fe_nstep','sw_nstep']) + + # all known keywords + variables = [ + 'ion_dynamics','ion_positions','pot_extrapolation','wfc_extrapolation', + 'remove_rigid_rot','ion_temperature','tempw','tolp','delta_t','nraise', + 'refold_pos','upscale','bfgs_ndim','trust_radius_max','trust_radius_min', + 'trust_radius_ini','w_1','w_2', + 'num_of_images','opt_scheme','CI_scheme','first_last_opt','temp_req', + 'ds','k_max','k_min','path_thr','use_masses','use_freezing','fe_step', + 'g_amplitude','fe_nstep','sw_nstep','phase_space', + ] + + # 5.4 keyword spec + #variables = [ + # 'ion_dynamics','ion_positions','pot_extrapolation','wfc_extrapolation', + # 'remove_rigid_rot','ion_temperature','tempw','tolp','delta_t','nraise', + # 'refold_pos','upscale','bfgs_ndim','trust_radius_max','trust_radius_min', + # 'trust_radius_ini','w_1','w_2' + # ] + + # sometime prior to 5.4 + #variables = [ + # 'ion_dynamics','ion_positions','phase_space','pot_extrapolation', + # 'wfc_extrapolation','remove_rigid_rot','ion_temperature','tempw', + # 'tolp','delta_t','nraise','refold_pos','upscale','bfgs_ndim', + # 'trust_radius_max','trust_radius_min','trust_radius_ini','w_1','w_2', + # 'num_of_images','opt_scheme','CI_scheme','first_last_opt','temp_req', + # 'ds','k_max','k_min','path_thr','use_masses','use_freezing','fe_step', + # 'g_amplitude','fe_nstep','sw_nstep' + # ] #end class ions class cell(Section): name = 'cell' - variables = set( - ['cell_dynamics','press','wmass','cell_factor','press_conv_thr', - 'cell_dofree']) + + # all known keywords + variables = [ + 'cell_dynamics','press','wmass','cell_factor','press_conv_thr', + 'cell_dofree' + ] + + # 5.4 keyword spec + #variables = [ + # 'cell_dynamics','press','wmass','cell_factor','press_conv_thr', + # 'cell_dofree' + # ] + + # sometime prior to 5.4 + #variables = [ + # 'cell_dynamics','press','wmass','cell_factor','press_conv_thr', + # 'cell_dofree' + # ] #end class cell class phonon(Section): name = 'phonon' - variables = set(['modenum','xqq']) + # all known keywords + variables = ['modenum','xqq'] + + # sometime prior to 5.4 + #variables = ['modenum','xqq'] #end class phonon class ee(Section): name = 'ee' - variables = set( - ['which_compensation','ecutcoarse','mixing_charge_compensation', - 'n_charge_compensation','comp_thr','nlev']) + # all known keywords + variables = [ + 'which_compensation','ecutcoarse','mixing_charge_compensation', + 'n_charge_compensation','comp_thr','nlev' + ] + + # sometime prior to 5.4 + #variables = [ + # 'which_compensation','ecutcoarse','mixing_charge_compensation', + # 'n_charge_compensation','comp_thr','nlev' + # ] #end class ee +section_classes = [ + control,system,electrons,ions,cell,phonon,ee + ] +for sec in section_classes: + sec.class_init() +#end for + def check_section_classes(*sections): all_variables = PwscfInputBase.all_variables global_missing = set(all_variables) local_missing = obj() locs_missing = False + secs = obj() for section in sections: variables = section.variables global_missing -= variables loc_missing = variables - all_variables local_missing[section.name] = loc_missing locs_missing |= len(loc_missing)>0 + secs[section.name] = section #end for if len(global_missing)>0 or locs_missing: msg = 'PwscfInput: variable information is not consistent for section classes\n' @@ -499,15 +752,27 @@ def check_section_classes(*sections): msg+=' some typed variables have not been assigned to a section:\n {0}\n'.format(sorted(global_missing)) #end if if locs_missing: - for name in sorted(locs_missing.keys()): - msg+=' some variables in section {0} have not been assigned a type:\n {1}\n'.format(name,sorted(locs_missing[name])) + for name in sorted(local_missing.keys()): + lmiss = local_missing[name] + if len(lmiss)>0: + vmiss = [] + for vname in secs[name].varlist: + if vname in lmiss: + vmiss.append(vname) + #end if + #end for + msg+=' some variables in section {0} have not been assigned a type\n missing variable counts: {1} {2}\n missing variables: {3}\n'.format(name,len(lmiss),len(vmiss),vmiss) + #end if #end for #end if - PwscfInput.class_error(msg) + print msg + else: + print 'pwscf input checks passed' #end if - ci(ls(),gs()) + exit() #end def check_section_classes -#check_section_classes(control,system,electrons,ions,cell,phonon,ee) +#check_section_classes(*section_classes) + class atomic_species(Card): @@ -521,7 +786,7 @@ class atomic_species(Card): tokens = l.split() atom = tokens[0] atoms.append(tokens[0]) - masses[atom] = float(tokens[1]) + masses[atom] = read_float(tokens[1]) pseudopotentials[atom] = tokens[2] #end for self.add(atoms=atoms,masses=masses,pseudopotentials=pseudopotentials) @@ -537,6 +802,7 @@ class atomic_species(Card): #end class atomic_species + class atomic_positions(Card): name = 'atomic_positions' @@ -623,6 +889,38 @@ class atomic_positions(Card): +class atomic_forces(Card): + name = 'atomic_forces' + + def read_text(self,lines): + npos = len(lines) + dim = 3 + atoms = [] + forces = empty((npos,dim)) + i=0 + for l in lines: + tokens = l.split() + atoms.append(tokens[0]) + forces[i,:] = array(tokens[1:4],dtype=float64) + i+=1 + #end for + self.add(atoms=atoms,forces=forces) + #end def read_text + + + def write_text(self): + c = '' + rowsep = '\n' + for i in range(len(self.atoms)): + c +=' '+'{0:2}'.format(self.atoms[i])+' ' + c += array_to_string(self.forces[i],pad='',rowsep=rowsep) + #end for + return c + #end def write_text +#end class atomic_forces + + + class k_points(Card): name = 'k_points' @@ -839,8 +1137,8 @@ class occupations(Card): class PwscfInput(SimulationInput): sections = ['control','system','electrons','ions','cell','phonon','ee'] - cards = ['atomic_species','atomic_positions','k_points', - 'cell_parameters','climbing_images','constraints', + cards = ['atomic_species','atomic_positions','atomic_forces', + 'k_points','cell_parameters','climbing_images','constraints', 'collective_vars','occupations'] section_types = obj( @@ -855,6 +1153,7 @@ class PwscfInput(SimulationInput): card_types = obj( atomic_species = atomic_species , atomic_positions = atomic_positions, + atomic_forces = atomic_forces , k_points = k_points , cell_parameters = cell_parameters , climbing_images = climbing_images , @@ -910,8 +1209,7 @@ class PwscfInput(SimulationInput): elem_type = 'section' c=[] else: - print 'PwscfInput Error: '+l[1:]+' is not a recognized pwscf section, exiting.' - exit() + self.error('encountered unrecognized input section during read\n{0} is not a recognized pwscf section\nfile read failed'.format(l[1:])) #end if elif tokens[0].lower() in self.cards and '=' not in l: if elem_type == 'card': @@ -933,9 +1231,7 @@ class PwscfInput(SimulationInput): elif in_element: c.append(l) else: - print 'PwscfInput Error: the following line is invalid, exiting.' - print l - exit() + self.error('invalid line encountered during read\ninvalid line: {0}\nfile read failed'.format(l)) #end if #end if #end for @@ -1024,9 +1320,7 @@ class PwscfInput(SimulationInput): ndn = p.down_electron.count self.system.ibrav = 0 - if 'celldm(1)' not in self.system: - self.system['celldm(1)'] = 1.0e0 - #end if + self.system['celldm(1)'] = 1.0e0 nions,nspecies = p.count_ions(species=True) self.system.nat = nions self.system.ntyp = nspecies @@ -1258,7 +1552,6 @@ class PwscfInput(SimulationInput): - def generate_pwscf_input(selector,**kwargs): if 'system' in kwargs: system = kwargs['system'] @@ -1390,14 +1683,10 @@ def generate_any_pwscf_input(**kwargs): #end for #copy certain keywords - tot_magnetization = None - if 'tot_magnetization' in kwargs: - tot_magnetization = kwargs.tot_magnetization - #end if - nspin = None - if 'nspin' in kwargs: - nspin = kwargs.nspin - #end if + tot_magnetization = kwargs.get_optional('tot_magnetization',None) + nspin = kwargs.get_optional('nspin',None) + nbnd = kwargs.get_optional('nbnd',None) + hubbard_u = kwargs.get_optional('hubbard_u',None) #make an empty input file pw = PwscfInput() @@ -1419,12 +1708,12 @@ def generate_any_pwscf_input(**kwargs): 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) totmag_sys = kwargs.delete_optional('totmag_sys',False) + bandfac = kwargs.delete_optional('bandfac',None) # pseudopotentials pseudopotentials = obj() @@ -1483,6 +1772,12 @@ def generate_any_pwscf_input(**kwargs): pw.system.nspin = 2 #end if + # set nbnd using bandfac, if provided + if nbnd is None and bandfac is not None: + nocc = max(system.particles.electron_counts()) + pw.system.nbnd = int(ceil(nocc*bandfac)) + #end if + # Hubbard U if hubbard_u!=None: if not isinstance(hubbard_u,(dict,obj)): @@ -1974,6 +2269,7 @@ def generate_relax_input(prefix = 'pwscf', return pw #end def generate_relax_input + def generate_vcrelax_input( press = None, # None = use pw.x default cell_factor = None, @@ -2075,3 +2371,5 @@ def generate_vcrelax_input( # # return pw ##end def generate_nscf_input + + diff --git a/nexus/library/qmcpack.py b/nexus/library/qmcpack.py index 901b7343b..8ff10c197 100644 --- a/nexus/library/qmcpack.py +++ b/nexus/library/qmcpack.py @@ -41,6 +41,7 @@ from qmcpack_converters import Pw2qmcpack,Wfconvert,Convert4qmc from sqd import Sqd from debug import ci,ls,gs from developer import unavailable +from nexus_base import nexus_core try: import h5py except ImportError: @@ -81,6 +82,17 @@ class Qmcpack(Simulation): user_twist_given |= isinstance(tnh,htypes) and tnh.twistnum!=None many_kpoints = len(self.system.structure.kpoints)>1 self.should_twist_average = many_kpoints and not user_twist_given + if self.should_twist_average: + # correct the job app command to account for the change in input file name + # this is necessary for twist averaged runs in bundles + app_comm = self.app_command() + prefix,ext = self.infile.split('.',1) + self.infile = prefix+'.in' + app_comm_new = self.app_command() + if self.job.app_command==app_comm: + self.job.app_command=app_comm_new + #end if + #end if #end if #end def post_init @@ -92,6 +104,17 @@ class Qmcpack(Simulation): #end def propagate_identifier + def pre_write_inputs(self,save_image): + # fix to make twist averaged input file under generate_only + if nexus_core.generate_only: + twistnums = range(len(self.system.structure.kpoints)) + if self.should_twist_average: + self.twist_average(twistnums) + #end if + #end if + #end def pre_write_inputs + + def check_result(self,result_name,sim): calculating_result = False if result_name=='jastrow' or result_name=='wavefunction': @@ -493,7 +516,7 @@ class Qmcpack(Simulation): def generate_qmcpack(**kwargs): - sim_args,inp_args = Simulation.separate_inputs(kwargs) + sim_args,inp_args = Qmcpack.separate_inputs(kwargs) if not 'input' in sim_args: input_type = inp_args.input_type diff --git a/nexus/library/qmcpack_input.py b/nexus/library/qmcpack_input.py index 4c1f64789..8b7abeeb8 100644 --- a/nexus/library/qmcpack_input.py +++ b/nexus/library/qmcpack_input.py @@ -2574,7 +2574,12 @@ localenergy.defaults.set( dm1b.defaults.set( type = 'dm1b',name='DensityMatrices' ) - +density.defaults.set( + type='density',name='Density' + ) +spindensity.defaults.set( + type='spindensity',name='SpinDensity' + ) @@ -5098,6 +5103,8 @@ def generate_basic_input(id = 'qmc', #end if if corrections=='default' and tuple(bconds)==tuple('ppp'): corrections = ['mpc','chiesa'] + elif isinstance(corrections,(list,tuple)): + None else: corrections = [] #end if diff --git a/nexus/library/qmcpack_workflows.py b/nexus/library/qmcpack_workflows.py index d815eba91..8fdf36d27 100644 --- a/nexus/library/qmcpack_workflows.py +++ b/nexus/library/qmcpack_workflows.py @@ -1,7 +1,12 @@ import os -from numpy import ndarray +import itertools +from time import clock +from numpy import ndarray,ceil from developer import obj,ci,error as dev_error,devlog,DevBase +from physical_system import generate_physical_system +from simulation import Simulation,GenericSimulation,graph_sims +from machines import job from vasp import generate_vasp from pwscf import generate_pwscf from pwscf_postprocessors import generate_projwfc @@ -10,6 +15,51 @@ from qmcpack_input import generate_jastrow,loop,linear,cslinear,vmc,dmc from qmcpack import generate_qmcpack + +# remaining dev tasks +# enable scans over shared parameters (e.g. pseudopotential sets) +# introduce relax sims (pwscf and vasp) into qmcpack_chain +# make ppset and add inputs +# eliminate detected fake sims in project manager +# probably do traversals to collect fake ones, then iterate to eliminate +# +# track association between input section and generated sims in qmcpack_chain +# this is needed in SimRep.__init__ +# set first input of each scan into fake chain prior to generating +# make "i" the default pattern if scanned variables are not str/int/float +# enable user provided labels for scan variables +# enables non-integer directory labeling for non-str/int/float variables +# consider how to replace deep copies in qmcpack_workflow (marked) +# test robustness of shallow copy, may be ok +# enable text print out of workflow tree +# useful on machines that do not allow plotting +# +# longer term +# make way to limit downstream branching action, not just "merge" +# perhaps at the level of segments based on key matching +# determine how best to extract output data from scanned simulations +# + +# +# enabling scan over system params +# add system_inputs to qmcpack chain +# process_qmcpack_chain_kwargs should essentially leave it alone +# if system_inputs is present +# generate a system object from inputs (via gen_phys_sys or user function) +# set kw.system to this +# if kw.fake then also create a SystemHolder(Simulation) object +# is fake sim by definition +# will need to have simlabel='system' for current secrep association +# make all sims in chain depend on it via 'other' +# this way, will appear in segmented representation and all variations +# should it go in the chain sims object or not? it has to +# SystemHolder's should be eliminated when qmcpack_workflow completes +# will also need to remove them from workflow sims object +# perhaps via implementing a 'recurse' function in obj class +# should probably make temp simlist at qmcpack_workflow start + + + def error(msg,loc=None,exit=True,trace=True,indent=' ',logfile=devlog): header = 'qmcpack_workflows' if loc!=None: @@ -81,7 +131,7 @@ default = Defaults() class Requirement(CallBase): def __call__(self,name,value): if missing(value): - error('a value has not been provided for required variable named {0}'.format(name),self.location) + error('an input value has not been provided for required variable named {0}'.format(name),self.location) #end if return value #end def __call__ @@ -89,6 +139,20 @@ class Requirement(CallBase): require = Requirement() +class RequireAny(CallBase): + def __call__(self,*name_values): + any_present = False + for name,value in name_values: + any_present |= not missing(value) + #end for + if not any_present: + error('an input value must be provided for at least one of these variables: {0}'.format([name for (name,value) in name_values]),self.location) + #end if + #end def __call__ +#end class RequireAny +require_any = RequireAny() + + class Assignment(CallBase): def __call__(self,o,name,value): if not missing(value): @@ -102,7 +166,7 @@ assign = Assignment() class RequireAssignment(CallBase): def __call__(self,o,name,value): if missing(value): - error('a value has not been provided for required variable named {0}'.format(name),self.location) + error('a value has not been provided for required input variable named "{0}"'.format(name),self.location) #end if o[name] = value #end def __call__ @@ -117,7 +181,7 @@ class DefaultAssignment(CallBase): elif name in self.defs: o[name] = self.defs[name] elif self.strict: - error('default value is missing for variable named {0}'.format(name),self.location) + error('default value is missing for variable named "{0}"'.format(name),self.location) #end if #end def __call__ #end class DefaultAssignment @@ -168,12 +232,12 @@ def extract_keywords(o,names,optional=False): #end def extract_keywords -def prevent_invalid_input(invalid,loc): +def prevent_invalid_input(invalid,valid,loc): if len(invalid)>0: if isinstance(invalid,(dict,obj)): invalid = invalid.keys() #end if - error('invalid input keywords encountered\ninvalid keywords: {0}\nfunction location: {1}'.format(sorted(invalid),loc)) + error('invalid input keywords encountered\ninvalid keywords: {0}\nvalid options are: {1}\nfunction location: {2}'.format(sorted(invalid),valid,loc)) #end if #end def prevent_invalid_input @@ -215,6 +279,80 @@ def capture_inputs(inputs): #end def capture_inputs +class SimSet(DevBase): + def __setitem__(self,simname,sim): + if simname in self: + self.error('simulation {0} is already set\nthis is a developer error'.format(simname)) + elif not isinstance(sim,(Simulation,SimSet)): + self.error('only simulation objects are allowed in simset\nreceived type: {0}\nwith name: {1}\nthis is a developer error'.format(sim.__class__.__name__,simname)) + else: + self.__dict__[simname] = sim + #end if + #end def + + def remove_fake(self): + fake_sims = [] + for k,v in self.iteritems(): + if isinstance(v,Simulation): + if v.fake(): + fake_sims.append(k) + #end if + elif isinstance(v,SimSet): + v.remove_fake() + #end if + #end for + for k in fake_sims: + del self[k] + #end for + #end def remove_fake +#end class SimSet + + +class SectionPlaceHolder(DevBase): + None +#end class SectionPlaceHolder +placeholder = SectionPlaceHolder + + +class SimHolder(GenericSimulation): + cls_simlabel = None + def __init__(self,*args,**kwargs): + cls = self.__class__ + Simulation.__init__(self,path='',job=job(app_command='',fake=True),fake_sim=True) + if cls.cls_simlabel!=None: + self.simlabel = cls.cls_simlabel + #end if + #end def __init__ +#end class SimHolder + + +class SystemHolder(SimHolder): + cls_simlabel = 'system' +#end class SystemHolder + + + +from memory import resident +tprev = clock() +mprev = resident() +def tpoint(loc,n=0,indent=' '): + global tprev + global mprev + tnow = clock() + mnow = resident() + #print '{0} {1} elapsed={2}'.format(n*indent,loc,tnow-tprev) + print '{0} {1} elapsed={2} {3:3.2f}'.format(n*indent,loc,tnow-tprev,(mnow-mprev)/1e6) + tprev = tnow + mprev = mnow +#end def tpoint + + + + + + + + jastrow_factor_keys = ['J1','J2','J3', 'J1_size','J1_rcut', @@ -226,9 +364,9 @@ jastrow_factor_defaults = obj( J1 = True, J2 = True, J3 = False, - J1_size = 10, + J1_size = None, J1_rcut = None, - J2_size = 10, + J2_size = None, J2_rcut = None, J2_init = 'zero', J3_isize = 3, @@ -437,6 +575,13 @@ dmc_sections_defaults = obj( +system_workflow_keys = [] +system_input_defaults = obj( + v1 = obj(), + ) + + + vasp_workflow_keys = [] vasp_input_defaults = obj( none = obj( @@ -451,9 +596,15 @@ vasp_input_defaults = obj( -scf_workflow_keys = [] +scf_workflow_keys = [ + 'struct_src', + ] +fixed_defaults = obj( + struct_src = None, + ) scf_input_defaults = obj( none = obj( + **fixed_defaults ), minimal = obj( identifier = 'scf', @@ -463,6 +614,7 @@ scf_input_defaults = obj( wf_collect = True, use_folded = True, nogamma = True, + **fixed_defaults ), v1 = obj( identifier = 'scf', @@ -481,12 +633,20 @@ scf_input_defaults = obj( wf_collect = True, # should be False if nscf is next use_folded = True, nogamma = True, + **fixed_defaults ), ) -nscf_workflow_keys = [] +nscf_workflow_keys = [ + 'struct_src','dens_src', + ] +fixed_defaults = obj( + struct_src = None, + dens_src = None, + ) nscf_input_defaults = obj( none = obj( + **fixed_defaults ), minimal = obj( identifier = 'nscf', @@ -494,6 +654,7 @@ nscf_input_defaults = obj( calculation = 'nscf', use_folded = True, nogamma = True, + **fixed_defaults ), v1 = obj( identifier = 'nscf', @@ -510,35 +671,57 @@ nscf_input_defaults = obj( degauss = 0.0001, use_folded = True, nogamma = True, + **fixed_defaults ), ) -p2q_workflow_keys = [] +p2q_workflow_keys = [ + 'orb_src', + ] +fixed_defaults = obj( + orb_src = None, + ) p2q_input_defaults = obj( minimal = obj( identifier = 'p2q', write_psir = False, + **fixed_defaults ), v1 = obj( identifier = 'p2q', write_psir = False, + **fixed_defaults ), ) -pwf_workflow_keys = [] +pwf_workflow_keys = [ + 'src', + ] +fixed_defaults = obj( + src = None, + ) pwf_input_defaults = obj( minimal = obj( identifier = 'pwf', + **fixed_defaults ), v1 = obj( identifier = 'pwf', + **fixed_defaults ), ) -opt_workflow_keys = ['J2_prod','J3_prod','J_defaults'] +opt_workflow_keys = [ + 'J2_run','J3_run','J_defaults', + 'struct_src','orb_src','J2_src','J3_src', + ] fixed_defaults = obj( - J2_prod = False, - J3_prod = False, + J2_run = False, + J3_run = False, + struct_src = None, + orb_src = None, + J2_src = None, + J3_src = None, J_defaults = defaults_version, ) opt_input_defaults = obj( @@ -556,19 +739,24 @@ opt_input_defaults = obj( ) vmc_workflow_keys = [ - 'J0_prod','J2_prod','J3_prod', + 'J0_run','J2_run','J3_run', 'J0_test','J2_test','J3_test', + 'struct_src','orb_src','J2_src','J3_src', 'job','test_job', ] fixed_defaults = obj( - J0_prod = False, - J2_prod = False, - J3_prod = False, - J0_test = False, - J2_test = False, - J3_test = False, - job = None, - test_job = None, + J0_run = False, + J2_run = False, + J3_run = False, + J0_test = False, + J2_test = False, + J3_test = False, + struct_src = None, + orb_src = None, + J2_src = None, + J3_src = None, + job = None, + test_job = None, ) vmc_input_defaults = obj( minimal = obj( @@ -598,19 +786,27 @@ vmc_input_defaults = obj( ) dmc_workflow_keys = [ - 'J0_prod','J2_prod','J3_prod', + 'J0_run','J2_run','J3_run', 'J0_test','J2_test','J3_test', - 'tmoves','locality' + 'struct_src','orb_src','J2_src','J3_src', + 'job','test_job', + 'tmoves','locality', ] fixed_defaults = obj( - J0_prod = False, - J2_prod = False, - J3_prod = False, - J0_test = False, - J2_test = False, - J3_test = False, - tmoves = False, - locality = False, + J0_run = False, + J2_run = False, + J3_run = False, + J0_test = False, + J2_test = False, + J3_test = False, + struct_src = None, + orb_src = None, + J2_src = None, + J3_src = None, + job = None, + test_job = None, + tmoves = False, + locality = False, ) dmc_input_defaults = obj( minimal = obj( @@ -628,39 +824,138 @@ dmc_input_defaults = obj( +def resolve_label(name,ch,loc='resolve_label',inp=False,quant=None,index=None,ind=False): + ind_req = ind + del ind + kw = ch.kw + task = ch.task + wf = task.workflow + if index is None: + index = 1 + if name in kw.def_srcs: + index = kw.def_srcs[name] + #end if + if quant is not None: + if quant=='structure' and 'struct_src' in wf: + ind = wf.struct_src + elif quant=='charge_density' and 'dens_src' in wf: + ind = wf.dens_src + elif quant=='orbitals' and 'orb_src' in wf: + ind = wf.orb_src + elif quant=='jastrow': + if 'J2' in name and 'J2_src' in wf: + ind = wf.J2_src + elif 'J3' in name and 'J3_src' in wf: + ind = wf.J3_src + #end if + elif quant=='other' and 'src' in wf: + ind = wf.src + #end if + if ind is not None: + index = ind + #end if + #end if + #end if + + if isinstance(index,Simulation): + sim = index + if not inp: + ret = sim + else: + ret = sim,name+'_inputs' + #end if + else: + if index==1: + sind = '' + else: + sind = str(index) + #end if + if '_' not in name: + label = name+sind + else: + name,postfix = name.split('_',1) + label = name+sind+'_'+postfix + #end if + inputs = name+'_inputs'+sind + if not inp: + ret = label + else: + ret = label,inputs + #end if + #end if + if ind_req: + if isinstance(ret,tuple): + ret = tuple(list(ret)+[index]) + else: + ret = ret,index + #end if + #end if + return ret +#end def resolve_label +def resolve_path(name,ch,loc='resolve_path'): + kw = ch.kw + task = ch.task + wf = task.workflow + index = None + if 'scf' in name: + if task.name=='nscf': + index = wf.dens_src + elif task.name=='p2q': + index = wf.orb_src + elif task.name=='pwf': + index = wf.src + #end if + #end if + if isinstance(index,Simulation): + sim = index + path = sim.locdir + else: + label = resolve_label(name,ch,loc,index=index) + path = os.path.join(wf.basepath,label) + #end if + return path +#end def resolve_path - - - - - -def resolve_deps(name,sims,deps,loc='resolve_deps'): +def resolve_deps(name,deps,ch,loc='resolve_deps'): + kw = ch.kw + task = ch.task + sims = ch.sims deplist = [] missing = [] + missing_inp = [] for depname,depquant in deps: - if depname in sims: - deplist.append((sims[depname],depquant)) + label,inp = resolve_label(depname,ch,loc,inp=True,quant=depquant) + if isinstance(label,Simulation): + sim = label + deplist.append((sim,depquant)) + elif label in sims: + deplist.append((sims[label],depquant)) else: - missing.append(depname) + missing.append(label) + missing_inp.append(inp) #end if #end for if len(missing)>0: - keywords = [] - for m in missing: - if len(m)>=3: - keywords.append(m[0:3]+'_inputs') - #end if - #end for - error('workflow cannot be run\nsimulation "{0}" depends on other simulations that have not been requested\nmissing simulations: {1}\nthe user needs to provide more detailed input\nthis issue can likely be fixed by providing the following keywords: {2}'.format(name,sorted(missing),sorted(set(keywords)))) + error('workflow cannot be run\nsimulation "{0}" depends on other simulations that have not been requested\nmissing simulations: {1}\nthe user needs to provide more detailed input\nthis issue can likely be fixed by providing the following keywords: {2}'.format(name,sorted(missing),sorted(set(missing_inp)))) #end if return deplist #end def resolve_deps - +def get_pseudos(type,shared): + if not missing(shared.pseudos): + return shared.pseudos + elif type=='pwscf': + return shared.dft_pseudos + elif type=='qmcpack': + return shared.qmc_pseudos + else: + error('pseudopotential type "{0}" is unrecognized\nthis is a developer error'.format(type),'get_pseudos') + #end if +#end def get_pseudos @@ -720,14 +1015,28 @@ def jastrow_factor( J3 = process_jastrow(J3,system) if J1==True: - if openbc and J1_rcut is None: - J1_rcut = J1_rcut_open + if J1_rcut is None: + if openbc: + J1_rcut = J1_rcut_open + else: + J1_rcut = system.structure.rwigner(1) + #end if + #end if + if J1_size is None: + J1_size = int(ceil(J1_rcut/0.5)) #end if J1 = generate_jastrow('J1','bspline',J1_size,J1_rcut,system=system) #end if if J2==True: - if openbc and J2_rcut is None: - J2_rcut = J2_rcut_open + if J2_rcut is None: + if openbc: + J2_rcut = J2_rcut_open + else: + J2_rcut = system.structure.rwigner(1) + #end if + #end if + if J2_size is None: + J2_size = int(ceil(J2_rcut/0.5)) #end if J2 = generate_jastrow('J2','bspline',J2_size,J2_rcut,init=J2_init,system=system) #end if @@ -856,7 +1165,7 @@ def opt_sections( error('invalid optimization cost function encountered\ninvalid cost fuction: {0}\nvalid options are: variance, energy, (0.95,0.05), etc'.format(cost),loc) #end if opt_calcs = [] - if abs(cost[0])>1e-6: + if abs(cost[0])>1e-6 and var_cycles>0: vmin_opt = opt( energy = 0.0, unreweightedvariance = 1.0, @@ -1082,564 +1391,1645 @@ def dmc_sections( +def process_system_inputs(inputs,shared,task,loc): + return inputs +#end def process_system_inputs +def process_vasp_inputs(inputs,shared,task,loc): + inputs.set( + system = shared.system, + ) + return inputs +#end def process_vasp_inputs -qmcpack_chain_required = ['system','sim_list','dft_pseudos','qmc_pseudos'] -qmcpack_chain_defaults = obj( - basepath = '', - orb_source = None, - J2_source = None, - J3_source = None, - vasp = False, - vasp_inputs = None, - vasp_defaults = defaults_version, - scf = False, - scf_inputs = None, - scf_defaults = defaults_version, - nscf = False, - nscf_inputs = None, - nscf_defaults = defaults_version, - nscf_transfer = True, - p2q = False, - p2q_inputs = None, - p2q_defaults = defaults_version, - pwf = False, - pwf_inputs = None, - pwf_defaults = defaults_version, - opt = False, - opt_inputs = None, - opt_defaults = defaults_version, - vmc = False, - vmc_inputs = None, - vmc_defaults = defaults_version, - dmc = False, - dmc_inputs = None, - dmc_defaults = defaults_version, - processed = False, + +def process_scf_inputs(inputs,shared,task,loc): + if shared.run.nscf: + if 'nosym' not in task.inputs_in: + inputs.nosym = False + #end if + if 'wf_collect' not in task.inputs_in: + inputs.wf_collect = False + #end if + #end if + inputs.set( + system = shared.system, + pseudos = get_pseudos('pwscf',shared), + ) + return inputs +#end def process_scf_inputs + + +def process_nscf_inputs(inputs,shared,task,loc): + inputs.set( + nosym = True, + wf_collect = True, + system = shared.system, + pseudos = get_pseudos('pwscf',shared), + ) + return inputs +#end def process_nscf_inputs + + +def process_p2q_inputs(inputs,shared,task,loc): + return inputs +#end def process_p2q_inputs + + +def process_pwf_inputs(inputs,shared,task,loc): + return inputs +#end def process_pwf_inputs + + +def process_opt_inputs(inputs,shared,task,loc): + inputs.set( + system = shared.system, + pseudos = get_pseudos('qmcpack',shared), + ) + + set_loc(loc+'opt_inputs jastrows') + jkw = extract_keywords(inputs,jastrow_factor_keys,optional=True) + jkw.system = shared.system + j2kw = jkw.copy() + j2kw.set(J1=1,J2=1,J3=0) + j3kw = jkw.copy() + j3kw.set(J1=1,J2=1,J3=1) + task.J2_inputs = j2kw + task.J3_inputs = j3kw + + set_loc(loc+'opt_inputs opt_methods') + task.sec_inputs = extract_keywords(inputs,opt_sections_keys,optional=True) + + return inputs +#end def process_opt_inputs + + +def process_vmc_inputs(inputs,shared,task,loc): + inputs.set( + system = shared.system, + pseudos = get_pseudos('qmcpack',shared), + ) + task.sec_inputs = extract_keywords(inputs,vmc_sections_keys,optional=True) + return inputs +#end def process_vmc_inputs + + +def process_dmc_inputs(inputs,shared,task,loc): + inputs.set( + system = shared.system, + pseudos = get_pseudos('qmcpack',shared), + ) + task.sec_inputs = extract_keywords(inputs,dmc_sections_keys,optional=True) + return inputs +#end def process_dmc_inputs + + +sim_input_defaults = obj( + system = system_input_defaults, + vasp = vasp_input_defaults, + scf = scf_input_defaults, + nscf = nscf_input_defaults, + p2q = p2q_input_defaults, + pwf = pwf_input_defaults, + opt = opt_input_defaults, + vmc = vmc_input_defaults, + dmc = dmc_input_defaults, ) +sim_workflow_keys = obj( + system = system_workflow_keys, + vasp = vasp_workflow_keys, + scf = scf_workflow_keys, + nscf = nscf_workflow_keys, + p2q = p2q_workflow_keys, + pwf = pwf_workflow_keys, + opt = opt_workflow_keys, + vmc = vmc_workflow_keys, + dmc = dmc_workflow_keys, + ) +process_sim_inp = obj( + system = process_system_inputs, + vasp = process_vasp_inputs, + scf = process_scf_inputs, + nscf = process_nscf_inputs, + p2q = process_p2q_inputs, + pwf = process_pwf_inputs, + opt = process_opt_inputs, + vmc = process_vmc_inputs, + dmc = process_dmc_inputs, + ) + +def process_sim_inputs(name,inputs_in,defaults,shared,task,loc): + set_loc(loc) + inputs = obj(**inputs_in) + assign_defaults(inputs,sim_input_defaults[name][defaults]) + inputs = process_sim_inp[name](inputs,shared,task,loc) + workflow = extract_keywords(inputs,sim_workflow_keys[name]) + task.inputs = inputs + task.workflow = workflow +#end def process_sim_inputs + + +qmcpack_chain_sim_names = ['system','vasp','scf','nscf','p2q','pwf','opt','vmc','dmc'] +def capture_qmcpack_chain_inputs(**kwargs): + kwcap = obj(**kwargs) + for name in qmcpack_chain_sim_names: + section = name+'_inputs' + if section in kwcap: + inputs = kwcap[section] + if inputs is not None: + kwcap[section] = obj(**inputs) + #end if + #end if + #end for + return kwcap +#end def capture_qmcpack_chain_inputs + def process_qmcpack_chain_kwargs( - basepath = missing, - system = missing, - sim_list = missing, - dft_pseudos = missing, - qmc_pseudos = missing, - orb_source = missing, - J2_source = missing, - J3_source = missing, - vasp = missing, - vasp_inputs = missing, - vasp_defaults = missing, - scf = missing, - scf_inputs = missing, - scf_defaults = missing, - nscf = missing, - nscf_inputs = missing, - nscf_defaults = missing, - nscf_transfer = missing, - p2q = missing, - p2q_inputs = missing, - p2q_defaults = missing, - pwf = missing, - pwf_inputs = missing, - pwf_defaults = missing, - opt = missing, - opt_inputs = missing, - opt_defaults = missing, - vmc = missing, - vmc_inputs = missing, - vmc_defaults = missing, - dmc = missing, - dmc_inputs = missing, - dmc_defaults = missing, - defaults = missing, - loc = 'process_qmcpack_chain_kwargs', - processed = False, - **invalid_kwargs + basepath = '', + orb_source = None, + J2_source = None, + J3_source = None, + system = missing, + system_function = missing, + sim_list = missing, + pseudos = missing, + dft_pseudos = missing, + qmc_pseudos = missing, + loc = 'process_qmcpack_chain_kwargs', + valid_inputs = None, + fake = False, + sims = None, + **other_kwargs ): - prevent_invalid_input(invalid_kwargs,loc) + if valid_inputs is None: + valid_inputs = [] + #end if + valid_inputs.extend([ + 'basepath','orb_source','J2_source','J3_source','system', + 'system_function','sim_list','pseudos','dft_pseudos', + 'qmc_pseudos','sims' + ]) + for sname in qmcpack_chain_sim_names: + valid_inputs.append(sname+'_inputs') + #end for - if missing(defaults): - defaults = qmcpack_chain_defaults + if 'system_inputs' in other_kwargs: + system_inputs = other_kwargs['system_inputs'] else: - defaults.set_optional(**qmcpack_chain_defaults) + system_inputs = missing #end if - set_def_loc(defaults,loc) + sim_names = set(qmcpack_chain_sim_names) - # capture scf_inputs for nscf, if present - scf_inputs_capture = capture_inputs(scf_inputs) + # collect other inputs into tasks + # organize all sim_inputs and sim_defaults input data + invalid = [] + tasks = obj() + task_defaults = obj() + for name in sim_names: + task_defaults[name] = defaults_version + #end for + for k,v in other_kwargs.iteritems(): + if '_inputs' in k: + name,tmp = k.split('_',1) + if name in sim_names: + index = tmp.replace('inputs','') + if len(index)==0: + index = 1 + else: + try: + index = int(index) + except: + error('inputs index must be an integer\ninput section: {0}\ninvalid index: {1}'.format(k,index),loc) + #end try + #end if + if not isinstance(v,(dict,obj)): + error('{0} must be a dict or obj\ntype received: {1}'.format(k,v.__class__.__name__),loc) + #end if + if name not in tasks: + tasks[name] = obj() + #end if + if index in tasks[name]: + error('repeated inputs index encountered\ninput section: {0}\ninvalid index: {1}'.format(k,index),loc) + #end if + if isinstance(v,SectionPlaceHolder): + tasks[name][index] = placeholder(inputs_in=obj()) + else: + tasks[name][index] = obj(inputs_in = obj(**v)) # capture the inputs + #end if + else: + invalid.append(k) + #end if + elif k.endswith('_defaults') and k.split('_')[0] in sim_names: + name = k.split('_')[0] + if name in sim_names: + task_defaults[name] = v + else: + invalid.append(k) + #end if + else: + invalid.append(k) + #end if + #end for + for name in sim_names: + if name not in tasks: + tasks[name] = None + #end if + #end for + prevent_invalid_input(invalid,valid_inputs,loc) + + set_loc(loc) kw = obj() - assign_require(kw,'system' ,system ) - assign_require(kw,'sim_list' ,sim_list ) - assign_default(kw,'basepath' ,basepath ) + run = obj() + for name in sim_names: + run[name] = tasks[name]!=None + #end for + kw.run = run - assign_default(kw,'orb_source' ,orb_source ) - assign_default(kw,'J2_source' ,J2_source ) - assign_default(kw,'J3_source' ,J3_source ) + kw.basepath = basepath + kw.orb_source = orb_source + kw.J2_source = J2_source + kw.J3_source = J3_source + kw.pseudos = pseudos + kw.fake = fake + kw.sims = sims - assign_default(kw,'vasp' ,vasp ) - assign_default(kw,'vasp_inputs' ,vasp_inputs ) - assign_default(kw,'vasp_defaults',vasp_defaults) - - assign_default(kw,'scf' ,scf ) - assign_default(kw,'scf_inputs' ,scf_inputs ) - assign_default(kw,'scf_defaults',scf_defaults) - - assign_default(kw,'nscf' ,nscf ) - assign_default(kw,'nscf_inputs' ,nscf_inputs ) - assign_default(kw,'nscf_defaults',nscf_defaults) - assign_default(kw,'nscf_transfer',nscf_transfer) - - assign_default(kw,'p2q' ,p2q ) - assign_default(kw,'p2q_inputs' ,p2q_inputs ) - assign_default(kw,'p2q_defaults',p2q_defaults) - - assign_default(kw,'pwf' ,pwf ) - assign_default(kw,'pwf_inputs' ,pwf_inputs ) - assign_default(kw,'pwf_defaults',pwf_defaults) - - assign_default(kw,'opt' ,opt ) - assign_default(kw,'opt_inputs' ,opt_inputs ) - assign_default(kw,'opt_defaults',opt_defaults) - - assign_default(kw,'vmc' ,vmc ) - assign_default(kw,'vmc_inputs' ,vmc_inputs ) - assign_default(kw,'vmc_defaults',vmc_defaults) - - assign_default(kw,'dmc' ,dmc ) - assign_default(kw,'dmc_inputs' ,dmc_inputs ) - assign_default(kw,'dmc_defaults',dmc_defaults) - - kw.vasp |= kw.vasp_inputs !=None - kw.scf |= kw.scf_inputs !=None - kw.nscf |= kw.nscf_inputs!=None - kw.p2q |= kw.p2q_inputs !=None - kw.pwf |= kw.pwf_inputs !=None - kw.opt |= kw.opt_inputs !=None - kw.vmc |= kw.vmc_inputs !=None - kw.dmc |= kw.dmc_inputs !=None - - if kw.scf or kw.nscf: - assign_require(kw,'dft_pseudos',dft_pseudos) - #end if - if kw.opt or kw.vmc or kw.dmc: - assign_require(kw,'qmc_pseudos',qmc_pseudos) - #end if - - - if kw.vasp: - # kw.vasp_inputs contains inputs to generate_vasp after this - set_loc(loc+' vasp_inputs') - kw.vasp_inputs = capture_inputs(kw.vasp_inputs) - assign_defaults(kw.vasp_inputs,vasp_input_defaults[kw.vasp_defaults]) - kw.vasp_inputs.set( - system = kw.system, - ) - kw.vasp_workflow = extract_keywords(kw.vasp_inputs,vasp_workflow_keys) - #end if - - if kw.scf: - # kw.scf_inputs contains inputs to generate_pwscf after this - set_loc(loc+' scf_inputs') - kw.scf_inputs = capture_inputs(kw.scf_inputs) - assign_defaults(kw.scf_inputs,scf_input_defaults[kw.scf_defaults]) - kw.scf_inputs.set( - system = kw.system, - pseudos = kw.dft_pseudos, - ) - kw.scf_workflow = extract_keywords(kw.scf_inputs,scf_workflow_keys) - #end if - - if kw.nscf: - # kw.nscf_inputs contains inputs to generate_pwscf after this - set_loc(loc+' nscf_inputs') - kw.nscf_inputs = capture_inputs(kw.nscf_inputs) - if kw.nscf_transfer: - scf_inputs_capture.delete_optional('kgrid') - scf_inputs_capture.delete_optional('kshift') - assign_defaults(kw.nscf_inputs,scf_inputs_capture) + kw.system_inputs_provided = not missing(system_inputs) + if not missing(system): + kw.system = system + elif not missing(system_inputs): + if missing(system_function): + system_function = generate_physical_system #end if - assign_defaults(kw.nscf_inputs,nscf_input_defaults[kw.nscf_defaults]) - kw.nscf_inputs.set( - nosym = True, - wf_collect = True, - system = kw.system, - pseudos = kw.dft_pseudos, - ) - kw.nscf_workflow = extract_keywords(kw.nscf_inputs,nscf_workflow_keys) + system = system_function(**system_inputs) + kw.system = system + else: + assign_require(kw,'system' ,system ) #end if - if kw.p2q: - # kw.p2q_inputs contains inputs to generate_pw2qmcpack after this - set_loc(loc+' p2q_inputs') - kw.p2q_inputs = capture_inputs(kw.p2q_inputs) - assign_defaults(kw.p2q_inputs,p2q_input_defaults[kw.p2q_defaults]) - kw.p2q_workflow = extract_keywords(kw.p2q_inputs,p2q_workflow_keys) + assign_require(kw,'sim_list' ,sim_list ) + + if missing(pseudos): + if run.scf or run.nscf: + assign_require(kw,'dft_pseudos',dft_pseudos) + #end if + if run.opt or run.vmc or run.dmc: + assign_require(kw,'qmc_pseudos',qmc_pseudos) + #end if #end if - if kw.pwf: - # kw.pwf_inputs contains inputs to generate_projwfc after this - set_loc(loc+' pwf_inputs') - kw.pwf_inputs = capture_inputs(kw.pwf_inputs) - assign_defaults(kw.pwf_inputs,pwf_input_defaults[kw.pwf_defaults]) - kw.pwf_workflow = extract_keywords(kw.pwf_inputs,pwf_workflow_keys) - #end if + for name in qmcpack_chain_sim_names: # order matters + if run[name]: + sim_tasks = tasks[name] + prev_index = None + for index in sorted(sim_tasks.keys()): + sim_task = sim_tasks[index] + if 'vary' in sim_task.inputs_in: + vary_index = sim_task.inputs_in.delete('vary') + else: + vary_index = prev_index + #end if + if vary_index!=None: + inputs_in = obj(**sim_tasks[vary_index].inputs_in) + inputs_in.set(**sim_task.inputs_in) + sim_task.inputs_in = inputs_in + #end if + sim_task.vary = vary_index + process_sim_inputs( + name = name, + inputs_in = sim_task.inputs_in, + defaults = task_defaults[name], + shared = kw, + task = sim_task, + loc = '{0} {1}_inputs'.format(loc,name), + ) + if index==1: + sind = '' + else: + sind = str(index) + #end if + label = name+sind + sim_task.name = name + sim_task.label = label + sim_task.index = index + sim_task.workflow.set( + basepath = basepath, + path = os.path.join(basepath,label), + vary_index = vary_index, + ) + prev_index = index + #end for + #end if + #end for - if kw.opt: - set_loc(loc+' opt_inputs') - kw.opt_inputs = capture_inputs(kw.opt_inputs) - assign_defaults(kw.opt_inputs,opt_input_defaults[kw.opt_defaults]) - kw.opt_inputs.set( - system = kw.system, - pseudos = kw.qmc_pseudos, - ) - kw.opt_workflow = extract_keywords(kw.opt_inputs,opt_workflow_keys) - - set_loc(loc+'opt_inputs jastrows') - #assign_defaults(kw.opt_inputs,jastrow_factor_defaults[kw.opt_workflow.J_defaults]) - jkw = extract_keywords(kw.opt_inputs,jastrow_factor_keys,optional=True) - jkw.system = kw.system - jkw.defaults = kw.opt_workflow.J_defaults - j2kw = jkw.copy() - j2kw.set(J1=1,J2=1,J3=0) - j3kw = jkw.copy() - j3kw.set(J1=1,J2=1,J3=1) - kw.J2_inputs = j2kw - kw.J3_inputs = j3kw - - set_loc(loc+'opt_inputs opt_methods') - kw.opt_sec_inputs = extract_keywords(kw.opt_inputs,opt_sections_keys,optional=True) - #end if - - if kw.vmc: - set_loc(loc+' vmc_inputs') - kw.vmc_inputs = capture_inputs(kw.vmc_inputs) - assign_defaults(kw.vmc_inputs,vmc_input_defaults[kw.vmc_defaults]) - kw.vmc_inputs.set( - system = kw.system, - pseudos = kw.qmc_pseudos, - ) - kw.vmc_workflow = extract_keywords(kw.vmc_inputs,vmc_workflow_keys) - - kw.vmc_sec_inputs = extract_keywords(kw.vmc_inputs,vmc_sections_keys,optional=True) - #end if - - if kw.dmc: - set_loc(loc+' dmc_inputs') - kw.dmc_inputs = capture_inputs(kw.dmc_inputs) - assign_defaults(kw.dmc_inputs,dmc_input_defaults[kw.dmc_defaults]) - kw.dmc_inputs.set( - system = kw.system, - pseudos = kw.qmc_pseudos, - ) - kw.dmc_workflow = extract_keywords(kw.dmc_inputs,dmc_workflow_keys) - - kw.dmc_sec_inputs = extract_keywords(kw.dmc_inputs,dmc_sections_keys,optional=True) - #end if - - del kw.vasp_defaults - del kw.scf_defaults - del kw.nscf_defaults - del kw.p2q_defaults - del kw.pwf_defaults - del kw.opt_defaults - del kw.vmc_defaults - del kw.dmc_defaults - - kw.processed = True + kw.tasks = tasks return kw #end def process_qmcpack_chain_kwargs -def qmcpack_chain(**kwargs): - loc = kwargs.pop('loc','qmcpack_chain') - processed = kwargs.pop('processed',False) - if processed: - kw = obj(**kwargs) - else: - kw = process_qmcpack_chain_kwargs(loc=loc,**kwargs) - #end if - basepath = kw.basepath - sim_list = kw.sim_list - sims = obj() - if kw.vasp: - vasp = generate_vasp( - path = os.path.join(basepath,'vasp'), - **kw.vasp_inputs +def gen_vasp_chain(ch,loc): + task = ch.task + wf = task.workflow + vasp = generate_vasp( + path = wf.path, + **task.inputs + ) + return vasp +#end def gen_vasp_chain + + +def gen_scf_chain(ch,loc): + task = ch.task + wf = task.workflow + scf = generate_pwscf( + path = wf.path, + **task.inputs + ) + return scf +#end def gen_scf_chain + + +def gen_nscf_chain(ch,loc): + task = ch.task + wf = task.workflow + deps = resolve_deps('nscf',[('scf','charge_density')],ch,loc) + # obtain user inputs to scf + scf_label,scf_index = resolve_label('scf',ch,loc,ind=True) + scf_inputs = obj(ch.tasks.scf[scf_index].inputs) + # pass inputs on to nscf after removing kpoint information + scf_inputs.delete_optional('kgrid') + scf_inputs.delete_optional('kshift') + task.inputs.set_optional(**scf_inputs) + # generate the nscf sim + path = resolve_path('scf',ch,loc) # added to always run in scf dir + nscf = generate_pwscf( + path = path, + dependencies = deps, + **task.inputs + ) + return nscf +#end def gen_nscf_chain + + +def gen_p2q_chain(ch,loc): + kw = ch.kw + task = ch.task + wf = task.workflow + run = kw.run + if run.nscf: + source = 'nscf' + else: + source = 'scf' + #end if + #path = resolve_path(source,ch,loc) + path = resolve_path('scf',ch,loc) # added to always run in scf dir + deps = resolve_deps('p2q',[(source,'orbitals')],ch,loc) + p2q = generate_pw2qmcpack( + path = path, + dependencies = deps, + **task.inputs + ) + return p2q +#end def gen_p2q_chain + + +def gen_pwf_chain(ch,loc): + kw = ch.kw + task = ch.task + wf = task.workflow + run = kw.run + if run.nscf: + source = 'nscf' + else: + source = 'scf' + #end if + path = resolve_path(source,ch,loc) + deps = resolve_deps('pwf',[(source,'other')],ch,loc) + pwf = generate_projwfc( + path = path, + dependencies = deps, + **task.inputs + ) + return pwf +#end def gen_pwf_chain + + +def gen_opt_chain(ch,loc): + kw = ch.kw + task = ch.task + sims = ch.sims + wf = task.workflow + orbdep = [('p2q','orbitals')] + J2dep = orbdep + [(task.label+'_J2','jastrow')] + if wf.J2_src is None: + optJ2_dep = orbdep + else: + optJ2_dep = orbdep + [('opt_J2','jastrow')] + #end if + if wf.J3_src is None: + if wf.J2_run: + optJ3_dep = J2dep + else: + optJ3_dep = orbdep + #end if + else: + optJ3_dep = orbdep + [('opt_J3','jastrow')] + #end if + if wf.J2_run: + label = task.label+'_J2' + deps = resolve_deps(label,optJ2_dep,ch,loc) + optJ2 = generate_qmcpack( + path = os.path.join(wf.basepath,label), + jastrows = jastrow_factor(**task.J2_inputs), + calculations = opt_sections(**task.sec_inputs), + dependencies = deps, + **task.inputs ) - sims.vasp = vasp + sims[label] = optJ2 + kw.def_srcs[task.name+'_J2'] = task.index + #end if + if wf.J3_run: + label = task.label+'_J3' + deps = resolve_deps(label,optJ3_dep,ch,loc) + optJ3 = generate_qmcpack( + path = os.path.join(wf.basepath,label), + jastrows = jastrow_factor(**task.J3_inputs), + calculations = opt_sections(**task.sec_inputs), + dependencies = deps, + **task.inputs + ) + sims[label] = optJ3 + kw.def_srcs[task.name+'_J3'] = task.index + #end if + return None +#end def gen_opt_chain + + +def gen_vmc(simlabel,ch,depset,J,test=0,loc=''): + task = ch.task + sims = ch.sims + wf = task.workflow + deps = resolve_deps(simlabel,depset[J],ch,loc) + if test: + qmcjob = wf.test_job + else: + qmcjob = wf.job + #end if + qmc = generate_qmcpack( + path = os.path.join(wf.basepath,simlabel), + job = qmcjob, + jastrows = [], + calculations = vmc_sections(test=test,J0=J=='J0',**task.sec_inputs), + dependencies = deps, + **task.inputs + ) + sims[simlabel] = qmc +#end def gen_vmc + + +def gen_vmc_chain(ch,loc): + task = ch.task + wf = task.workflow + orbdep = [('p2q','orbitals')] + J2dep = orbdep + [('opt_J2','jastrow')] + J3dep = orbdep + [('opt_J3','jastrow')] + depset = obj( + J0 = orbdep, + J2 = J2dep, + J3 = J3dep, + ) + lab = task.label + if wf.J0_test: + gen_vmc(lab+'_J0_test',ch,depset,J='J0',test=1,loc=loc) + #end if + if wf.J0_run: + gen_vmc(lab+'_J0',ch,depset,J='J0',loc=loc) + #end if + if wf.J2_test: + gen_vmc(lab+'_J2_test',ch,depset,J='J2',test=1,loc=loc) + #end if + if wf.J2_run: + gen_vmc(lab+'_J2',ch,depset,J='J2',loc=loc) + #end if + if wf.J3_test: + gen_vmc(lab+'_J3_test',ch,depset,J='J3',test=1,loc=loc) + #end if + if wf.J3_run: + gen_vmc(lab+'_J3',ch,depset,J='J3',loc=loc) + #end if + return None +#end def gen_vmc_chain + + +def gen_dmc(simlabel,ch,depset,J,nlmove=None,test=0,loc=''): + task = ch.task + sims = ch.sims + wf = task.workflow + deps = resolve_deps(simlabel,depset[J],ch,loc) + if test: + qmcjob = wf.test_job + else: + qmcjob = wf.job + #end if + other_inputs = obj(task.inputs) + if 'calculations' not in other_inputs: + other_inputs.calculations = dmc_sections(nlmove=nlmove,test=test,J0=J=='J0',**task.sec_inputs) + #end if + qmc = generate_qmcpack( + path = os.path.join(wf.basepath,simlabel), + job = qmcjob, + jastrows = [], + dependencies = deps, + **other_inputs + ) + sims[simlabel] = qmc +#end def gen_dmc + + +def gen_dmc_chain(ch,loc): + kw = ch.kw + task = ch.task + wf = task.workflow + orbdep = [('p2q','orbitals')] + J2dep = orbdep + [('opt_J2','jastrow')] + J3dep = orbdep + [('opt_J3','jastrow')] + depset = obj( + J0 = orbdep, + J2 = J2dep, + J3 = J3dep, + ) + nonloc_labels = {None:'',True:'_tm',False:'_la'} + nonlocalmoves_default = None + if kw.system.pseudized: + nonlocalmoves_default = True + #end if + nloc_moves = [] + if wf.tmoves: + nloc_moves.append(True) + #end if + if wf.locality: + nloc_moves.append(False) + #end if + if not wf.tmoves and not wf.locality: + nloc_moves.append(nonlocalmoves_default) + #end if + lab = task.label + for nlmove in nloc_moves: + nll = nonloc_labels[nlmove] + if wf.J0_test: + gen_dmc(lab+'_J0'+nll+'_test',ch,depset,J='J0',nlmove=nlmove,test=1,loc=loc) + #end if + if wf.J0_run: + gen_dmc(lab+'_J0'+nll,ch,depset,J='J0',nlmove=nlmove,loc=loc) + #end if + if wf.J2_test: + gen_dmc(lab+'_J2'+nll+'_test',ch,depset,J='J2',nlmove=nlmove,test=1,loc=loc) + #end if + if wf.J2_run: + gen_dmc(lab+'_J2'+nll,ch,depset,J='J2',nlmove=nlmove,loc=loc) + #end if + if wf.J3_test: + gen_dmc(lab+'_J3'+nll+'_test',ch,depset,J='J3',nlmove=nlmove,test=1,loc=loc) + #end if + if wf.J3_run: + gen_dmc(lab+'_J3'+nll,ch,depset,J='J3',nlmove=nlmove,loc=loc) + #end if + #end for + return None +#end def gen_dmc_chain + + +gen_sim_ch = obj( + vasp = gen_vasp_chain, + scf = gen_scf_chain, + nscf = gen_nscf_chain, + p2q = gen_p2q_chain, + pwf = gen_pwf_chain, + opt = gen_opt_chain, + vmc = gen_vmc_chain, + dmc = gen_dmc_chain, + ) + +def gen_sim_chain(name,ch,loc): + kw = ch.kw + tasks = ch.tasks + sims = ch.sims + sim_tasks = tasks[name] + for index in sorted(sim_tasks.keys()): + task = sim_tasks[index] + if not isinstance(task,SectionPlaceHolder): + ch.set( + name = name, + index = index, + task = task, + ) + sim = gen_sim_ch[name](ch,loc) + if sim is not None: + sims[task.label] = sim + #end if + #end if + kw.def_srcs[task.name] = task.index # only the last index will matter + #end for +#end def gen_sim_chain + + + +def qmcpack_chain(**kwargs): + loc = kwargs.pop('loc','qmcpack_chain') + + kw = process_qmcpack_chain_kwargs(loc=loc,**kwargs) + + sim_list = kw.sim_list + run = kw.run + tasks = kw.tasks + + if kw.fake: + Simulation.creating_fake_sims = True + #end if + + if kw.sims is None: + sims = SimSet() + else: + sims = kw.sims + #end if + kw.def_srcs = obj() + + ch = obj( + kw = kw, + tasks = tasks, + sims = sims, + ) + + if run.vasp: + gen_sim_chain('vasp',ch,loc) #end if if kw.orb_source!=None: sims.p2q = kw.orb_source else: - if kw.scf: - scf = generate_pwscf( - path = os.path.join(basepath,'scf'), - **kw.scf_inputs - ) - sims.scf = scf + if run.scf: + gen_sim_chain('scf',ch,loc) #end if - - if kw.nscf: - deps = resolve_deps('nscf',sims,[('scf','charge_density')],loc) - nscf = generate_pwscf( - path = os.path.join(basepath,'scf'), - dependencies = deps, - **kw.nscf_inputs - ) - sims.nscf = nscf + if run.nscf: + gen_sim_chain('nscf',ch,loc) #end if - - if kw.p2q: - if kw.nscf: - p2q_path = os.path.join(basepath,'scf') - deps = resolve_deps('p2q',sims,[('nscf','orbitals')],loc) - else: - p2q_path = os.path.join(basepath,'scf') - deps = resolve_deps('p2q',sims,[('scf','orbitals')],loc) - #end if - p2q = generate_pw2qmcpack( - path = p2q_path, - dependencies = deps, - **kw.p2q_inputs - ) - sims.p2q = p2q + if run.p2q: + gen_sim_chain('p2q',ch,loc) #end if - - if kw.pwf: - deps = resolve_deps('pwf',sims,[('scf','other')],loc) - pwf = generate_projwfc( - path = os.path.join(basepath,'scf'), - dependencies = deps, - **kw.pwf_inputs - ) - sims.pwf = pwf + if run.pwf: + gen_sim_chain('pwf',ch,loc) #end if #end if - - orbdep = [('p2q','orbitals')] - J2dep = orbdep + [('optJ2','jastrow')] - J3dep = orbdep + [('optJ3','jastrow')] - if kw.opt: - if kw.opt_workflow.J2_prod and kw.J2_source is None: - deps = resolve_deps('optJ2',sims,orbdep,loc) - optJ2 = generate_qmcpack( - path = os.path.join(basepath,'optJ2'), - jastrows = jastrow_factor(**kw.J2_inputs), - calculations = opt_sections(**kw.opt_sec_inputs), - dependencies = deps, - **kw.opt_inputs - ) - sims.optJ2 = optJ2 - #end if - if kw.opt_workflow.J3_prod and kw.J3_source is None: - deps = resolve_deps('optJ3',sims,J2dep,loc) - optJ3 = generate_qmcpack( - path = os.path.join(basepath,'optJ3'), - jastrows = jastrow_factor(**kw.J3_inputs), - calculations = opt_sections(**kw.opt_sec_inputs), - dependencies = deps, - **kw.opt_inputs - ) - sims.optJ3 = optJ3 - #end if + if run.opt: + gen_sim_chain('opt',ch,loc) #end if if kw.J2_source!=None: - sims.optJ2 = kw.J2_source + sims.opt_J2 = kw.J2_source #end if if kw.J3_source!=None: - sims.optJ3 = kw.J3_source + sims.opt_J3 = kw.J3_source #end if - if kw.vmc: - if kw.vmc_workflow.J0_test: - deps = resolve_deps('vmcJ0_test',sims,orbdep,loc) - vmcJ0_test = generate_qmcpack( - path = os.path.join(basepath,'vmcJ0_test'), - job = kw.vmc_workflow.test_job, - jastrows = [], - calculations = vmc_sections(test=1,J0=1,**kw.vmc_sec_inputs), - dependencies = deps, - **kw.vmc_inputs - ) - sims.vmcJ0_test = vmcJ0_test - #end if - if kw.vmc_workflow.J0_prod: - deps = resolve_deps('vmcJ0',sims,orbdep,loc) - vmcJ0 = generate_qmcpack( - path = os.path.join(basepath,'vmcJ0'), - job = kw.vmc_workflow.job, - jastrows = [], - calculations = vmc_sections(J0=1,**kw.vmc_sec_inputs), - dependencies = deps, - **kw.vmc_inputs - ) - sims.vmcJ0 = vmcJ0 - #end if - if kw.vmc_workflow.J2_test: - deps = resolve_deps('vmcJ2_test',sims,J2dep,loc) - vmcJ2_test = generate_qmcpack( - path = os.path.join(basepath,'vmcJ2_test'), - job = kw.vmc_workflow.test_job, - jastrows = [], - calculations = vmc_sections(test=1,**kw.vmc_sec_inputs), - dependencies = deps, - **kw.vmc_inputs - ) - sims.vmcJ2_test = vmcJ2_test - #end if - if kw.vmc_workflow.J2_prod: - deps = resolve_deps('vmcJ2',sims,J2dep,loc) - vmcJ2 = generate_qmcpack( - path = os.path.join(basepath,'vmcJ2'), - job = kw.vmc_workflow.job, - jastrows = [], - calculations = vmc_sections(**kw.vmc_sec_inputs), - dependencies = deps, - **kw.vmc_inputs - ) - sims.vmcJ2 = vmcJ2 - #end if - if kw.vmc_workflow.J3_test: - deps = resolve_deps('vmcJ3_test',sims,J3dep,loc) - vmcJ3_test = generate_qmcpack( - path = os.path.join(basepath,'vmcJ3_test'), - job = kw.vmc_workflow.test_job, - jastrows = [], - calculations = vmc_sections(test=1,**kw.vmc_sec_inputs), - dependencies = deps, - **kw.vmc_inputs - ) - sims.vmcJ3_test = vmcJ3_test - #end if - if kw.vmc_workflow.J3_prod: - deps = resolve_deps('vmcJ3',sims,J3dep,loc) - vmcJ3 = generate_qmcpack( - path = os.path.join(basepath,'vmcJ3'), - job = kw.vmc_workflow.job, - jastrows = [], - calculations = vmc_sections(**kw.vmc_sec_inputs), - dependencies = deps, - **kw.vmc_inputs - ) - sims.vmcJ3 = vmcJ3 - #end if + if run.vmc: + gen_sim_chain('vmc',ch,loc) #end if - if kw.dmc: - nonloc_labels = {None:'',True:'_tm',False:'_la'} - nonlocalmoves_default = None - if kw.system.pseudized: - nonlocalmoves_default = True - #end if - nloc_moves = [] - if kw.dmc_workflow.tmoves: - nloc_moves.append(True) - #end if - if kw.dmc_workflow.locality: - nloc_moves.append(False) - #end if - if not kw.dmc_workflow.tmoves and not kw.dmc_workflow.locality: - nloc_moves.append(nonlocalmoves_default) - #end if - for nlmove in nloc_moves: - nll = nonloc_labels[nlmove] - if kw.dmc_workflow.J0_test: - label = 'dmcJ0'+nll+'_test' - deps = resolve_deps(label,sims,orbdep,loc) - dmcJ0_test = generate_qmcpack( - path = os.path.join(basepath,label), - jastrows = [], - calculations = dmc_sections(nlmove=nlmove,test=1,J0=1,**kw.dmc_sec_inputs), - dependencies = deps, - **kw.dmc_inputs - ) - sims[label] = dmcJ0_test - #end if - if kw.dmc_workflow.J0_prod: - label = 'dmcJ0'+nll - deps = resolve_deps(label,sims,orbdep,loc) - dmcJ0 = generate_qmcpack( - path = os.path.join(basepath,label), - jastrows = [], - calculations = dmc_sections(nlmove=nlmove,J0=1,**kw.dmc_sec_inputs), - dependencies = deps, - **kw.dmc_inputs - ) - sims[label] = dmcJ0 - #end if - if kw.dmc_workflow.J2_test: - label = 'dmcJ2'+nll+'_test' - deps = resolve_deps(label,sims,J2dep,loc) - dmcJ2_test = generate_qmcpack( - path = os.path.join(basepath,label), - jastrows = [], - calculations = dmc_sections(nlmove=nlmove,test=1,**kw.dmc_sec_inputs), - dependencies = deps, - **kw.dmc_inputs - ) - sims[label] = dmcJ2_test - #end if - if kw.dmc_workflow.J2_prod: - label = 'dmcJ2'+nll - deps = resolve_deps(label,sims,J2dep,loc) - dmcJ2 = generate_qmcpack( - path = os.path.join(basepath,label), - jastrows = [], - calculations = dmc_sections(nlmove=nlmove,**kw.dmc_sec_inputs), - dependencies = deps, - **kw.dmc_inputs - ) - sims[label] = dmcJ2 - #end if - if kw.dmc_workflow.J3_test: - label = 'dmcJ3'+nll+'_test' - deps = resolve_deps(label,sims,J3dep,loc) - dmcJ3_test = generate_qmcpack( - path = os.path.join(basepath,label), - jastrows = [], - calculations = dmc_sections(nlmove=nlmove,test=1,**kw.dmc_sec_inputs), - dependencies = deps, - **kw.dmc_inputs - ) - sims[label] = dmcJ3_test - #end if - if kw.dmc_workflow.J3_prod: - label = 'dmcJ3'+nll - deps = resolve_deps(label,sims,J3dep,loc) - dmcJ3 = generate_qmcpack( - path = os.path.join(basepath,label), - jastrows = [], - calculations = dmc_sections(nlmove=nlmove,**kw.dmc_sec_inputs), - dependencies = deps, - **kw.dmc_inputs - ) - sims[label] = dmcJ3 + if run.dmc: + gen_sim_chain('dmc',ch,loc) + #end if + + # apply labeling to simulations + for simlabel,sim in sims.iteritems(): + sim.simlabel = simlabel + #end for + + if not kw.fake: + sim_list.extend(sims.list()) + else: + Simulation.creating_fake_sims = False + #end if + + # for scanned workflows only + # add fake system sim + # make all other sims depend on it + if kw.fake and kw.system_inputs_provided: + sys = SystemHolder(system=kw.system) + for sim in sims: + if len(sim.dependencies)==0: + sim.depends(sys,'other') #end if #end for + sims.system = sys #end if - sim_list.extend(sims.list()) + #print kw.def_srcs + #exit() return sims #end def qmcpack_chain -# alias to qmcpack_workflow for user interaction -qmcpack_workflow = qmcpack_chain + +def unpack_scan_inputs(scan,loc='unpack_scan_inputs'): + if not isinstance(scan,(tuple,list)): + error('parameter "scan" must be a list or tuple\ntype received: {0}'.format(scan.__class__.__name__),loc=loc) + #end if + scans = obj() + for scan_list in scan: + if not isinstance(scan_list,(tuple,list)): + error('scan list for parameter "scan" must be a tuple or list\ntype received: {0}'.format(scan_list.__class__.__name__),loc=loc) + #end if + misformatted = len(scan_list)<3 + misformatted|= (len(scan_list)-1)%2!=0 or not isinstance(scan_list[0],str) + if not misformatted: + section = scan_list[0] + vars = scan_list[1::2] + vals = scan_list[2::2] + values_in = obj() + for var,val_list in zip(vars,vals): + values_in[var] = val_list + misformatted |= not isinstance(var,str) or not isinstance(val_list,(tuple,list,ndarray)) + #end for + if not misformatted: + if len(vars)==1: + vals_in = vals[0] + vals = [] + inds = [] + i = 1 + for v in vals_in: + vals.append((v,)) + inds.append((i,)) + i+=1 + #end for + else: + vars = tuple(vars) + n=1 + for vals_in in vals[1:]: + if len(vals_in)!=len(vals[0]): + error('problem with "scan" input\nall value_lists must have the same length (they are covarying)\nvar_list "{0}" has length {1}\nvar_list "{2}" has length {3}'.format(vars[0],len(vals[0]),vars[n],len(vals[n])),loc=loc) + #end if + n+=1 + #end for + vals = zip(*vals) + r = range(1,len(vals)+1) + inds = zip(*[r for n in range(len(vars))]) + #end if + if section not in scans: + sec = obj(parameters=list(vars),values=vals,indices=inds,values_in=values_in) + scans[section] = sec + else: + sec = scans[section] + sec.values_in.transfer_from(values_in) + sec.parameters.extend(vars) + values = [] + for v in sec.values: + for v2 in vals: + val = tuple(list(v)+list(v2)) + values.append(val) + #end for + #end for + sec.values = values + indices = [] + for i in sec.indices: + for i2 in inds: + ind = tuple(list(i)+list(i2)) + indices.append(ind) + #end for + #end for + sec.indices = indices + #end if + #end if + #end if + if misformatted: + error('scan list for parameter "scan" is formatted improperly\nlist must have input section name (string) followed by variable-value_list pairs\nnumber of entries present (should be odd): {0}\nvalue_lists must be of type tuple/list/array\nscan list contents: {1}'.format(len(scan_list),scan_list,loc)) + #end if + #end for + return scans +#end def unpack_scan_inputs + + +def unpack_fix_inputs(fix,loc='unpack_fix_inputs'): + if not isinstance(fix,(tuple,list)): + error('parameter "fix" must be a list or tuple\ntype received: {0}'.format(fix.__class__.__name__),loc=loc) + #end if + constraints = obj() + for clist in fix: + if not isinstance(clist,(tuple,list)): + error('constraint list for parameter "fix" must be a tuple or list\ntype received: {0}'.format(clist.__class__.__name__),loc=loc) + #end if + misformatted = len(clist)<3 + misformatted|= (len(clist)-1)%2!=0 or not isinstance(clist[0],str) + if not misformatted: + sim_name = clist[0] + variables = list(clist[1::2]) + selectors = list(clist[2::2]) + for var in variables: + misformatted |= not isinstance(var,str) + #end for + if not misformatted: + if variables[-1]=='pattern': + pattern = list(selectors[-1]) + variables.pop() + selectors.pop() + if len(pattern)!=len(variables): + error('must have a single pattern element for each variable in contraint list for parameter "fix"\nnumber of pattern elements: {0}\nnumber of variables: {1}\npattern list: {2}\nvariable list: {3}'.format(len(pattern),len(variables),pattern,variables),loc=loc) + #end if + pattern_misformatted = False + pstr = '' + prob = '' + for v,s,p in zip(variables,selectors,pattern): + bad_int = p=='i' and not isinstance(s,int) + bad_val = p=='v' and not isinstance(s,(str,int,float)) + bad_pat = p not in ('v','i') + if bad_int or bad_val or bad_pat: + pattern_misformatted = True + pstr+='({0})'.format(p) + prob+='{0} (invalid pattern {1})'.format(v,p) + else: + pstr+=p + #end if + #end for + if pattern_misformatted: + error('constraint values do not match pattern provided\nthe problem variables are: {0}\npattern provided: {1}\nvariables provided: {2}\nfor pattern "v" only an integer,float, or string can be provided\nfor pattern "i" only an integer can be provided'.format(prob,pstr,variables)) + #end if + else: + pattern = list(len(variables)*'v') + #end if + #end if + #end if + if misformatted: + error('constraint list for parameter "fix" is formatted improperly\nlist must have sim name (string) followed by variable-value/index (string/int etc.) pairs\nnumber of entries present (should be odd): {0}\nlist contents: {1}'.format(len(clist),clist,loc)) + #end if + constraints[sim_name] = obj( + sim_name = sim_name, + variables = variables, + selectors = selectors, + pattern = pattern, + ) + #end for + return constraints +#end def unpack_fix_inputs +class WFRep(DevBase): + @classmethod + def get_id(cls): + if not cls.class_has('rep_count'): + cls.class_set(rep_count=0) + #end if + id = cls.rep_count + cls.rep_count+=1 + return id + #end def get_id + + def __init__(self,label=None): + self.identifier = label + self.id = self.__class__.get_id() + self.label = label + self.graphid = self.id + self.depth = 0 + self.dependents = obj() + self.dependencies = obj() + #end def __init__ + + def add_dependent(self,rep): + self.dependents[rep.id] = rep + rep.dependencies[self.id] = self + #end def add_dependent + + + def assign_depth(self,depth=None): + if depth is None: + depth = self.depth + #end if + self.depth = max(self.depth,depth) + for rep in self.dependents: + rep.assign_depth(depth+1) + #end for + #end def assign_depth + + + def assign_graphid_down(self,graphid=None): + if graphid is None: + graphid = self.graphid + else: + self.graphid = graphid + #end if + for rep in self.dependents: + rep.assign_graphid_down(graphid) + #end for + #end def assign_graphid_down + + + def assign_graphid_up(self,graphid=None): + if graphid is None: + graphid = self.graphid + else: + self.graphid = graphid + #end if + for rep in self.dependencies: + rep.assign_graphid_up(graphid) + #end for + #end def assign_graphid_up + + + def downstream_labels(self,collection=None): + if collection is None: + collection = set() + #end if + for rep in self.dependents: + collection.add(rep.label) + #end for + for rep in self.dependents: + rep.downstream_labels(collection) + #end for + return collection + #end def downstream_labels + + + def upstream_labels(self,collection=None): + if collection is None: + collection = set() + #end if + for rep in self.dependencies: + collection.add(rep.label) + #end for + for rep in self.dependencies: + rep.upstream_labels(collection) + #end for + return collection + #end def upstream_labels + + @property + def simid(self): + return self.id + #end def simid + + @property + def simlabel(self): + return self.label + #end def simlabel +#end class WFRep + + +class SimRep(WFRep): + def __init__(self,sim): + WFRep.__init__(self) + self.transfer_from(sim,['identifier']) + self.id = sim.simid + self.label = sim.simlabel + self.graphid = self.simid + self.secrep = None + label = self.simlabel + if '_' in label: + label = label.split('_',1)[0] + #end if + n=len(label) + while label[n-1].isdigit(): + n=-1 + #end while + # get the input section + # the association between label and section will not always hold + # jtk mark: need to replace this with something better + # probably should track this during qmcpack_chain + self.section = label[:n]+'_inputs'+label[n:] + #print self.simlabel,self.section + #end def __init__ +#end class SimRep + + +class SecRep(WFRep): + def __init__(self,section,simreps_in): + WFRep.__init__(self,section) + self.segrep = None + sims = obj() + for sim in simreps_in: + if sim.section==section: + sims[sim.simid] = sim + sim.secrep = self + #end if + #end for + self.simreps = sims + #end def __init__ + + @property + def section(self): + return self.label + #end def section +#end class SecRep + + +class SegRep(WFRep): + def __init__(self,label,secreps,invariant=False): + WFRep.__init__(self,label) + self.invariant = invariant + self.secreps = obj() + for sec in secreps: + self.secreps[sec.label] = sec + sec.segrep = self + #end for + self.scan_data = None # non-invariant only + self.fixed_inputs = None + self.reference_inputs = None + #end def __init__ + + + def set_reference(self,active_sections,scans,kwargs,placeholders): + if not self.invariant: + self.scan_data = scans[self.label] + #end if + fixed = obj() + for k,v in kwargs.iteritems(): + if k not in active_sections: + fixed[k] = v + #end if + #end for + self.fixed_inputs = fixed + ref = obj() + for section in self.secreps.keys(): + ref[section] = obj(**kwargs[section]) # shallow copy inputs data + #ref[section].delete_optional('job') + #end for + self.placeholders = obj() + for section,ph in placeholders.iteritems(): + if section not in ref: + self.placeholders[section] = ph + #end if + #end for + self.reference_inputs = ref + #end def set_reference + + + def check_constraint(self,constraint,tol,first=True,loc='check_constraint'): + if first: + constraint.indices = obj() + constraint.check = obj( + variables_found = set(), + matches_found = set(), + all_upstream_variables = set(), + ) + #end if + sd = self.scan_data + if sd is not None and len(sd.parameters)>0: + c = constraint + for v,s,p in zip(c.variables,c.selectors,c.pattern): + for k in sd.values_in.keys(): + c.check.all_upstream_variables.add(k) + #end for + if v in sd.values_in: + values = sd.values_in[v] + c.check.variables_found.add(v) + if p=='i': + if s>0 and s1: + fix_vars = sorted(set(constraint.variables)) + fix_inds = constraint.indices.tuple(fix_vars) + fixed_sims = obj() + fixed_simids = set() + grouped_sims = obj() + inds0 = scoll[0][1] + other_vars = sorted(set(inds0.keys())-set(fix_vars)) + for sim,inds,vals in scoll: + ifix = inds.tuple(fix_vars) + iother = inds.tuple(other_vars) + if ifix==fix_inds: + fixed_sims[iother] = sim + else: + if iother not in grouped_sims: + grouped_sims[iother] = [] + #end if + grouped_sims[iother].append(sim) + #end if + #end for + for iother,fixed_sim in fixed_sims.iteritems(): + if iother in grouped_sims: + for sim in grouped_sims[iother]: + fixed_sim.acquire_dependents(sim) + #end for + #end if + #end for + #end if + #end for + #end if + + # get rid of any fake/temporary simulations + sims.remove_fake() + + + if graph_workflow: + graph_sims(loc_sims,exit=False) + #end if + if write_workflow: + write_sims(loc_sims,exit=False) + #end if + + if graph_simreps: + graph_sims(simreps,quants=False,exit=False) + #end if + if write_simreps: + write_sims(simreps,quants=False,exit=False) + #end if + if graph_secreps: + graph_sims(secreps,quants=False,exit=False) + #end if + if write_secreps: + write_sims(secreps,quants=False,exit=False) + #end if + if graph_segreps: + graph_sims(segreps,quants=False,exit=False) + #end if + if write_segreps: + write_sims(segreps,quants=False,exit=False) + #end if + + + sim_list.extend(loc_sims) + + #tpoint('finish workflows') + + + #ci() + #exit() + + return sims +#end def qmcpack_workflow + + + +# deprecated ecut_scan_chain_defaults = obj( scf = True, p2q = True, @@ -1717,6 +3107,7 @@ def ecut_scan( +# deprecated def system_scan( basepath = missing, dirname = 'system_scan', @@ -1794,7 +3185,7 @@ def system_scan( #end def system_scan - +# deprecated def system_parameter_scan( basepath = missing, dirname = 'system_param_scan', @@ -1852,7 +3243,7 @@ def system_parameter_scan( - +# deprecated def input_parameter_scan( basepath = missing, dirname = 'input_param_scan', @@ -1916,7 +3307,7 @@ def input_parameter_scan( same_jastrow = not missing(jastrow_key) if same_jastrow: - jastrow_key = render_key(jastrow_key,loc) + jastrow_key = render_parameter_key(jastrow_key,loc) if not jastrow_key in set(paramkeys): error('input parameter value for fixed Jastrow not found\njastrow key provided: {0}\nparameter keys present: {1}'.format(jastrow_key,paramkeys),loc) #end if @@ -1978,6 +3369,192 @@ def input_parameter_scan( +# deprecated +def input_parameter_scan2( + basepath = missing, + dirname = 'input_param_scan', + section = missing, + parameter = missing, + parameterset = missing, + variable = missing, # same as parameter, to be deprecated + values = missing, + tags = missing, + jastrow_key = missing, + J2_source = None, + J3_source = None, + loc = 'input_parameter_scan', + **kwargs + ): + + set_loc(loc) + + if missing(parameter): + parameter = variable + #end if + + require('basepath' ,basepath ) + require('section' ,section ) + if missing(parameterset): + require('parameter',parameter) + require('values' ,values ) + if not missing(tags) and len(tags)!=len(values): + error('must provide one tag (directory label) per parameter value\nnumber of values: {0}\nnumber of tags: {1}\nvalues provided: {2}\ntags provided: {3}'.format(len(values),len(tags),values,tags),loc) + #end if + else: + require('tags',tags) + for pname,pvalues in parameterset.iteritems(): + if len(tags)!=len(pvalues): + error('must provide one tag (directory label) per parameter value\nparameter name: {0}\nnumber of values: {1}\nnumber of tags: {2}\nvalues provided: {3}\ntags provided: {4}'.format(pname,len(pvalues),len(tags),pvalues,tags),loc) + #end if + #end for + #end if + + + paramdirs = [] + paramkeys = [] + if missing(parameterset): + n = 0 + for v in values: + if not missing(tags): + vkey = tags[n] + else: + vkey = render_parameter_key(v,loc) + #end if + vlabel = render_parameter_label(vkey,loc) + paramdir = '{0}_{1}'.format(parameter,vlabel) + paramdirs.append(paramdir) + paramkeys.append(vkey) + n+=1 + #end for + parameterset = obj() + parameterset[parameter] = values + else: + paramdirs.extend(tags) + paramkeys.extend(tags) + #end if + + same_jastrow = not missing(jastrow_key) + if same_jastrow: + jastrow_key = render_key(jastrow_key,loc) + if not jastrow_key in set(paramkeys): + error('input parameter value for fixed Jastrow not found\njastrow key provided: {0}\nparameter keys present: {1}'.format(jastrow_key,paramkeys),loc) + #end if + index = paramkeys.index(jastrow_key) + paramdirs.insert(0,paramdirs.pop(index)) + paramkeys.insert(0,paramkeys.pop(index)) + for pname,pvalues in parameterset.iteritems(): + pvalues = list(pvalues) + pvalues.insert(0,pvalues.pop(index)) + parameterset[pname] = pvalues + #end for + #end if + + varying = obj() + for name in qmcpack_chain_sim_names: + varying[name] = section==name+'_inputs' + #end for + varying.nscf |= varying.scf + varying.p2q |= varying.nscf + varying.opt |= varying.p2q + varying.vmc |= varying.opt + varying.dmc |= varying.opt + + suppress_J0 = False + suppress_J2 = False + for pname in parameterset.keys(): + if pname.startswith('J3'): + suppress_J0 = True + suppress_J2 = True + elif pname.startswith('J2'): + suppress_J0 = True + #end if + #end for + + ip_sims = obj() + for n in xrange(len(paramkeys)): + + kwcap = capture_qmcpack_chain_inputs(**kwargs) + + if 'inputs' not in section: + error('section must be a name of the form *_inputs\nsection provided: {0}\ninput sections present: {1}'.format(section,[s for s in sorted(kwcap.keys()) if 'inputs' in s]),loc) + elif section not in kwcap: + error('input section not found\ninput section provided: {0}\ninput sections present: {1}'.format(section,[s for s in sorted(kwcap.keys()) if 'inputs' in s]),loc) + #end if + sec = kwcap[section] + for pname,pvalues in parameterset.iteritems(): + pvalue = pvalues[n] + if not novalue(pvalue): + sec[pname] = pvalue + #end if + #end for + + qckw = process_qmcpack_chain_kwargs( + defaults = qmcpack_chain_defaults, + loc = loc, + **kwcap + ) + + if n>0: + for name,vary in varying.iteritems(): + if not vary: + qckw[name] = False + qckw[name+'_inputs'] = None + #end if + #end for + if 'vmc_inputs' in qckw: + if suppress_J0: + qckw.vmc_inputs.J0_run = 0 + qckw.vmc_inputs.J0_test = 0 + #end if + if suppress_J2: + qckw.vmc_inputs.J2_run = 0 + qckw.vmc_inputs.J2_test = 0 + #end if + #end if + if 'opt_inputs' in qckw: + if suppress_J2: + qckw.opt_inputs.J2_run = 0 + qckw.opt_inputs.J2_test = 0 + #end if + #end if + if 'dmc_inputs' in qckw: + if suppress_J0: + qckw.dmc_inputs.J0_run = 0 + qckw.dmc_inputs.J0_test = 0 + #end if + if suppress_J2: + qckw.dmc_inputs.J2_run = 0 + qckw.dmc_inputs.J2_test = 0 + #end if + #end if + #end if + + qckw.basepath = os.path.join(basepath,dirname,paramdirs[n]) + qckw.J2_source = J2_source + qckw.J3_source = J3_source + qcsims = qmcpack_chain(**qckw) + if n==0: + if not varying.p2q: + orb_source = qcsims.get_optional('p2q',None) + elif not varying.opt: + J2_source = qcsims.get_optional('optJ2',None) + J3_source = qcsims.get_optional('optJ3',None) + elif varying.opt and suppress_J2: + J3_source = qcsims.get_optional('optJ3',None) + elif same_jastrow: + J2_source = qcsims.get_optional('optJ2',None) + J3_source = qcsims.get_optional('optJ3',None) + #end if + #end if + ip_sims[paramkeys[n]] = qcsims + #end for + + return ip_sims +#end def input_parameter_scan2 + + + + if __name__=='__main__': print 'simple driver for qmcpack_workflows functions' diff --git a/nexus/library/simulation.py b/nexus/library/simulation.py index 7b9449e85..c5d3043c9 100644 --- a/nexus/library/simulation.py +++ b/nexus/library/simulation.py @@ -74,6 +74,7 @@ from generic import obj from periodic_table import is_element from physical_system import PhysicalSystem from machines import Job +from pseudopotential import ppset from nexus_base import NexusCore,nexus_core @@ -257,7 +258,7 @@ class Simulation(NexusCore): allowed_inputs = set(['identifier','path','infile','outfile','errfile','imagefile', 'input','job','files','dependencies','analysis_request', 'block','block_subcascade','app_name','app_props','system', - 'skip_submit','force_write']) + 'skip_submit','force_write','simlabel','fake_sim']) sim_imagefile = 'sim.p' input_imagefile = 'input.p' analyzer_imagefile = 'analyzer.p' @@ -266,12 +267,17 @@ class Simulation(NexusCore): is_bundle = False sim_count = 0 + creating_fake_sims = False sim_directories = dict() + @classmethod + def code_name(cls): + return cls.generic_identifier + #end def code_name - @staticmethod - def separate_inputs(kwargs,overlapping_kw=-1,copy_pseudos=True): + @classmethod + def separate_inputs(cls,kwargs,overlapping_kw=-1,copy_pseudos=True): if overlapping_kw==-1: overlapping_kw = set(['system']) elif overlapping_kw is None: @@ -286,6 +292,22 @@ class Simulation(NexusCore): inp_args.transfer_from(kwargs,inp_kw) if 'pseudos' in inp_args and inp_args.pseudos!=None: pseudos = inp_args.pseudos + # support ppset labels + if isinstance(pseudos,str): + code = cls.code_name() + if not ppset.supports_code(code): + cls.class_error('ppset labeled pseudopotential groups are not supported for code "{0}"'.format(code)) + #end if + if 'system' not in inp_args: + cls.class_error('system must be provided when using a ppset label') + #end if + system = inp_args.system + pseudos = ppset.get(pseudos,code,system) + if 'pseudos' in sim_args: + sim_args.pseudos = pseudos + #end if + inp_args.pseudos = pseudos + #end if if copy_pseudos: if 'files' in sim_args: sim_args.files = list(sim_args.files) @@ -297,17 +319,17 @@ class Simulation(NexusCore): if 'system' in inp_args: system = inp_args.system if not isinstance(system,PhysicalSystem): - Simulation.class_error('system object must be of type PhysicalSystem') + cls.class_error('system object must be of type PhysicalSystem') #end if species_labels,species = system.structure.species(symbol=True) pseudopotentials = nexus_core.pseudopotentials for ppfile in pseudos: if not ppfile in pseudopotentials: - Simulation.class_error('pseudopotential file {0} cannot be found'.format(ppfile)) + cls.class_error('pseudopotential file {0} cannot be found'.format(ppfile)) #end if pp = pseudopotentials[ppfile] if not pp.element_label in species_labels and not pp.element in species: - Simulation.class_error('the element {0} for pseudopotential file {1} is not in the physical system provided'.format(pp.element,ppfile)) + cls.class_error('the element {0} for pseudopotential file {1} is not in the physical system provided'.format(pp.element,ppfile)) #end if #end for #end if @@ -329,6 +351,7 @@ class Simulation(NexusCore): #variables determined by self self.identifier = self.generic_identifier self.simid = Simulation.sim_count + self.simlabel = None Simulation.sim_count+=1 self.files = set() self.app_name = self.application @@ -364,6 +387,7 @@ class Simulation(NexusCore): self.errfile = None self.bundled = False self.bundler = None + self.fake_sim = Simulation.creating_fake_sims #variables determined by derived classes self.outputs = None #object representing output data @@ -377,12 +401,15 @@ class Simulation(NexusCore): self.init_job() #end if self.post_init() - if self.system!=None: - self.system.check_folded_system() - #end if + #end def __init__ + def fake(self): + return self.fake_sim + #end def fake + + def init_job(self): if self.job==None: self.error('job not provided. Input field job must be set to a Job object.') @@ -402,7 +429,12 @@ class Simulation(NexusCore): self.depends(*kw['dependencies']) del kw['dependencies'] #end if - allowed = set(kw.keys()) & self.allowed_inputs + kwset = set(kw.keys()) + invalid = kwset - self.allowed_inputs + if len(invalid)>0: + self.error('received invalid inputs\ninvalid inputs: {0}\nallowed inputs are: {1}'.format(sorted(invalid),sorted(self.allowed_inputs))) + #end if + allowed = kwset & self.allowed_inputs for name in allowed: self[name] = kw[name] #end for @@ -439,14 +471,18 @@ class Simulation(NexusCore): self.remdir = os.path.join(nexus_core.remote_directory,nexus_core.runs,self.path) self.resdir = os.path.join(nexus_core.local_directory,nexus_core.results,nexus_core.runs,self.path) - if not self.locdir in self.sim_directories: - self.sim_directories[self.locdir] = set([self.identifier]) - else: - idset = self.sim_directories[self.locdir] - if not self.identifier in idset: - idset.add(self.identifier) + if not self.fake(): + print ' creating sim {0} in {1}'.format(self.simid,self.locdir) + + if not self.locdir in self.sim_directories: + self.sim_directories[self.locdir] = set([self.identifier]) else: - self.error('multiple simulations in a single directory have the same identifier\nplease assign unique identifiers to each simulation\nsimulation directory: {0}\nrepeated identifier: {1}\nother identifiers: {2}\nbetween the directory shown and the identifiers listed, it should be clear which simulations are involved\nmost likely, you described two simulations with identifier {3}'.format(self.locdir,self.identifier,sorted(idset),self.identifier)) + idset = self.sim_directories[self.locdir] + if not self.identifier in idset: + idset.add(self.identifier) + else: + self.error('multiple simulations in a single directory have the same identifier\nplease assign unique identifiers to each simulation\nsimulation directory: {0}\nrepeated identifier: {1}\nother identifiers: {2}\nbetween the directory shown and the identifiers listed, it should be clear which simulations are involved\nmost likely, you described two simulations with identifier {3}'.format(self.locdir,self.identifier,sorted(idset),self.identifier)) + #end if #end if #end if @@ -660,6 +696,53 @@ class Simulation(NexusCore): #end def depends + def undo_depends(self,sim): + i=0 + for dep in self.ordered_dependencies: + if dep.sim.simid==sim.simid: + break + #end if + i+=1 + #end for + self.ordered_dependencies.pop(i) + del self.dependencies[sim.simid] + del sim.dependents[self.simid] + self.dependency_ids.remove(sim.simid) + if sim.simid in self.wait_ids: + self.wait_ids.remove(sim.simid) + #end if + #end def undo_depends + + + def acquire_dependents(self,sim): + # acquire the dependents from the other simulation + dsims = obj(sim.dependents) + for dsim in dsims: + dep = dsim.dependencies[sim.simid] + dsim.depends(self,*dep.result_names) + #end for + # eliminate the other simulation + # this renders it void (fake) and removes all dependency relationships + sim.eliminate() + #end def acquire_dependents + + + def eliminate(self): + # reverse relationship of dependents (downstream) + dsims = obj(self.dependents) + for dsim in dsims: + dsim.undo_depends(self) + #end for + # reverse relationship of dependencies (upstream) + deps = obj(self.dependencies) + for dep in deps: + self.undo_depends(dep.sim) + #end for + # mark sim to be ignored in all future interactions + self.fake_sim = True + #end def eliminate + + def check_dependencies(self,result): dep_satisfied = result.dependencies_satisfied for dep in self.dependencies: @@ -978,14 +1061,13 @@ class Simulation(NexusCore): if self.finished: self.enter(self.locdir,False,self.simid) self.log('analyzing'+self.idstr(),n=3) - analyzer = self.analyzer_type(self) if not nexus_core.generate_only: + analyzer = self.analyzer_type(self) analyzer.analyze() self.post_analyze(analyzer) + analyzer.save(os.path.join(self.imresdir,self.analyzer_image)) + del analyzer #end if - analyzer.save(os.path.join(self.imresdir,self.analyzer_image)) - - del analyzer self.analyzed = True self.save_image() #end if @@ -1541,16 +1623,25 @@ except: Image = unavailable('Image') #end try import tempfile -def graph_sims(sims): +exit_call = exit +def graph_sims(sims,useid=False,exit=True,quants=True): graph = Dot(graph_type='digraph') graph.set_label('simulation workflows') graph.set_labelloc('t') nodes = obj() for sim in sims: + if 'fake_sim' in sim and sim.fake_sim: + continue + #end if + if sim.simlabel is not None and not useid: + nlabel = sim.simlabel+' '+str(sim.simid) + else: + nlabel = sim.identifier+' '+str(sim.simid) + #end if node = obj( id = sim.simid, sim = sim, - node = Node(sim.identifier+' '+str(sim.simid),style='filled',shape='Mrecord'), + node = Node(nlabel,style='filled',shape='Mrecord'), edges = obj(), ) nodes[node.id] = node @@ -1559,10 +1650,15 @@ def graph_sims(sims): for node in nodes: for simid,dep in node.sim.dependencies.iteritems(): other = nodes[simid].node - for quantity in dep.result_names: - edge = Edge(other,node.node,label=quantity,fontsize='10.0') + if quants: + for quantity in dep.result_names: + edge = Edge(other,node.node,label=quantity,fontsize='10.0') + graph.add_edge(edge) + #end for + else: + edge = Edge(other,node.node) graph.add_edge(edge) - #end for + #end if #end for #end for @@ -1573,5 +1669,11 @@ def graph_sims(sims): image = Image.open(savefile) image.show() - exit() + if exit: + exit_call() + #end if #end def graph_sims + + + +#def write_sims(sims, diff --git a/nexus/library/structure.py b/nexus/library/structure.py index cdb703802..c52ca025f 100755 --- a/nexus/library/structure.py +++ b/nexus/library/structure.py @@ -460,6 +460,7 @@ class Structure(Sobj): magnetization=None,magnetic_order=None,magnetic_prim=True, operations=None,background_charge=0,frozen=None,bconds=None, posu=None): + if center is None: if axes is not None: center = array(axes,dtype=float).sum(0)/2 @@ -489,16 +490,16 @@ class Structure(Sobj): if mag is None: mag = len(elem)*[None] #end if - self.scale = 1. - self.units = units - self.dim = dim - self.center = array(center,dtype=float) - self.axes = array(axes,dtype=float) + self.scale = 1. + self.units = units + self.dim = dim + self.center = array(center,dtype=float) + self.axes = array(axes,dtype=float) self.set_bconds(bconds) self.set_elem(elem) - self.pos = array(pos,dtype=float) - self.frozen = None - self.mag = array(mag,dtype=object) + self.pos = array(pos,dtype=float) + self.frozen = None + self.mag = array(mag,dtype=object) self.kpoints = empty((0,dim)) self.kweights = empty((0,)) self.background_charge = background_charge @@ -674,6 +675,60 @@ class Structure(Sobj): #end if #end def reshape_axes + + def corners(self): + a = self.axes + c = array([(0,0,0), + a[0], + a[1], + a[2], + a[0]+a[1], + a[1]+a[2], + a[2]+a[0], + a[0]+a[1]+a[2], + ]) + return c + #end def corners + + + def miller_direction(self,h,k,l,normalize=False): + d = dot((h,k,l),self.axes) + if normalize: + d/=norm(d) + #end if + return d + #end def miller_direction + + + def miller_normal(self,h,k,l,normalize=False): + d = dot((h,k,l),self.kaxes) + if normalize: + d/=norm(d) + #end if + return d + #end def miller_normal + + + def project_plane(self,a1,a2,points=None): + # a1/a2: in plane vectors + if points is None: + points = self.pos + #end if + a1n = norm(a1) + a2n = norm(a2) + a1/=a1n + a2/=a2n + n = cross(a1,a2) + plane_coords = [] + for p in points: + p -= dot(n,p)*n # project point into plane + c1 = dot(a1,p)/a1n + c2 = dot(a2,p)/a2n + plane_coords.append((c1,c2)) + #end for + return array(plane_coords,dtype=float) + #end def project_plane + def bounding_box(self,scale=1.0,box='tight',recenter=False): pmin = self.pos.min(0) @@ -736,7 +791,7 @@ class Structure(Sobj): P[:,i] = pv[:] #end for self.center = dot(self.center,P) - if len(self.axes)>0: + if self.has_axes(): self.axes = dot(self.axes,P) #end if if len(self.pos)>0: @@ -798,19 +853,18 @@ class Structure(Sobj): def any_periodic(self): - has_cell = len(self.axes)>0 - has_kpoints = len(self.kpoints)>0 + has_cell = self.has_axes() pbc = False for bc in self.bconds: pbc |= bc=='p' #end if - periodic = has_cell and (pbc or has_kpoints) + periodic = has_cell and pbc return periodic #end def any_periodic def all_periodic(self): - has_cell = len(self.axes)>0 + has_cell = self.has_axes() pbc = True for bc in self.bconds: pbc &= bc=='p' @@ -839,7 +893,7 @@ class Structure(Sobj): def volume(self): - if len(self.axes)==0: + if not self.has_axes(): return None else: return abs(det(self.axes)) @@ -885,6 +939,22 @@ class Structure(Sobj): #end def rinscribe + def rwigner_cube(self,*args,**kwargs): + cube = Structure() + a = self.volume()**(1./3) + cube.set_axes([[a,0,0],[0,a,0],[0,0,a]]) + return cube.rwigner(*args,**kwargs) + #end def rwigner_cube + + + def rinscribe_cube(self,*args,**kwargs): + cube = Structure() + a = self.volume()**(1./3) + cube.set_axes([[a,0,0],[0,a,0],[0,0,a]]) + return cube.rinscribe(*args,**kwargs) + #end def rinscribe_cube + + def rmin(self): return self.rwigner() #end def rmin @@ -1038,7 +1108,7 @@ class Structure(Sobj): v = -v # preserve the normal direction for atom identification, but reverse the shift direction #end if self.recorner() # want box contents to be static - if len(self.axes)>0: + if self.has_axes(): components = 0 dim = self.dim axes = self.axes @@ -2627,7 +2697,7 @@ class Structure(Sobj): def tile_opt(self,volfac,dn=1,tol=1e-3): - Topt = self.opt_tilematrix(volfac,dn,tol) + Topt,ropt = self.opt_tilematrix(volfac,dn,tol) return self.tile(Topt) #end def tile_opt @@ -3332,9 +3402,14 @@ class Structure(Sobj): self.read_poscar(filepath,elem=elem) elif format=='cif': self.read_cif(filepath,block=block,grammar=grammar,cell=cell) + elif format=='fhi-aims': + self.read_fhi_aims(filepath) else: self.error('cannot read structure from file\nunsupported file format: {0}'.format(format)) #end if + if self.has_axes(): + self.set_bconds('ppp') + #end if #end def read @@ -3489,6 +3564,7 @@ class Structure(Sobj): def read_cif(self,filepath,block=None,grammar='1.1',cell='prim'): axes,elem,pos,units = read_cif(filepath,block,grammar,cell,args_only=True) + self.dim = 3 self.set_axes(axes) self.set_elem(elem) self.pos = pos @@ -3496,6 +3572,40 @@ class Structure(Sobj): #end def read_cif + def read_fhi_aims(self,filepath): + if os.path.exists(filepath): + lines = open(filepath,'r').read().splitlines() + else: + lines = filepath.splitlines() # "filepath" is contents + #end if + axes = [] + posu = [] + elem = [] + for line in lines: + ls = line.strip() + if len(ls)>0 and ls[0]!='#': + tokens = ls.split() + if ls.startswith('lattice_vector'): + axes.append(tokens[1:]) + elif ls.startswith('atom_frac'): + posu.append(tokens[1:4]) + elem.append(tokens[4]) + else: + None + #self.error('unrecogonized or not yet supported token in fhi-aims geometry file: {0}'.format(tokens[0])) + #end if + #end if + #end for + axes = array(axes,dtype=float) + pos = dot(array(posu,dtype=float),axes) + self.dim = 3 + self.set_axes(axes) + self.set_elem(elem) + self.pos = pos + self.units = 'A' + #end def read_fhi_aims + + def write(self,filepath=None,format=None): if filepath is None and format is None: self.error('please specify either the filepath or format arguments to write()') @@ -3511,6 +3621,8 @@ class Structure(Sobj): c = self.write_xyz(filepath) elif format=='xsf': c = self.write_xsf(filepath) + elif format=='fhi-aims': + c = self.write_fhi_aims(filepath) else: self.error('file format {0} is unrecognized'.format(format)) #end if @@ -3578,6 +3690,25 @@ class Structure(Sobj): #end def write_xsf + def write_fhi_aims(self,filepath=None): + s = self.copy() + s.change_units('A') + c = '' + c+='\n' + for a in s.axes: + c += 'lattice_vector {0: 12.8f} {1: 12.8f} {2: 12.8f}\n'.format(*a) + #end for + c+='\n' + for p,e in zip(self.pos,self.elem): + c += 'atom_frac {0: 12.8f} {1: 12.8f} {2: 12.8f} {3}\n'.format(p[0],p[1],p[2],e) + #end for + if filepath!=None: + open(filepath,'w').write(c) + #end if + return c + #end def write_fhi_aims + + def plot2d_ax(self,ix,iy,*args,**kwargs): if self.dim!=3: self.error('plot2d_ax is currently only implemented for 3 dimensions') @@ -4600,7 +4731,7 @@ def generate_structure(type='crystal',*args,**kwargs): -def generate_atom_structure(atom=None,units='A',Lbox=None,skew=0,axes=None,kgrid=(1,1,1),kshift=(0,0,0),struct_type=Structure): +def generate_atom_structure(atom=None,units='A',Lbox=None,skew=0,axes=None,kgrid=(1,1,1),kshift=(0,0,0),bconds=tuple('nnn'),struct_type=Structure): if atom is None: Structure.class_error('atom must be provided','generate_atom_structure') #end if @@ -4608,16 +4739,17 @@ def generate_atom_structure(atom=None,units='A',Lbox=None,skew=0,axes=None,kgrid axes = [[Lbox*(1-skew),0,0],[0,Lbox,0],[0,0,Lbox*(1+skew)]] #end if if axes is None: - s = Structure(elem=[atom],pos=[[0,0,0]],units=units) + s = Structure(elem=[atom],pos=[[0,0,0]],units=units,bconds=bconds) else: - s = Structure(elem=[atom],pos=[[0,0,0]],axes=axes,kgrid=kgrid,kshift=kshift,units=units) + s = Structure(elem=[atom],pos=[[0,0,0]],axes=axes,kgrid=kgrid,kshift=kshift,bconds=bconds,units=units) s.center_molecule() #end if + return s #end def generate_atom_structure -def generate_dimer_structure(dimer=None,units='A',separation=None,Lbox=None,skew=0,axes=None,kgrid=(1,1,1),kshift=(0,0,0),struct_type=Structure,axis='x'): +def generate_dimer_structure(dimer=None,units='A',separation=None,Lbox=None,skew=0,axes=None,kgrid=(1,1,1),kshift=(0,0,0),bconds=tuple('nnn'),struct_type=Structure,axis='x'): if dimer is None: Structure.class_error('dimer atoms must be provided to construct dimer','generate_dimer_structure') #end if @@ -4637,9 +4769,9 @@ def generate_dimer_structure(dimer=None,units='A',separation=None,Lbox=None,skew Structure.class_error('dimer orientation axis must be x,y,z\n you provided: {0}'.format(axis),'generate_dimer_structure') #end if if axes is None: - s = Structure(elem=dimer,pos=[[0,0,0],p2],units=units) + s = Structure(elem=dimer,pos=[[0,0,0],p2],units=units,bconds=bconds) else: - s = Structure(elem=dimer,pos=[[0,0,0],p2],axes=axes,kgrid=kgrid,kshift=kshift,units=units) + s = Structure(elem=dimer,pos=[[0,0,0],p2],axes=axes,kgrid=kgrid,kshift=kshift,units=units,bconds=bconds) s.center_molecule() #end if return s diff --git a/nexus/library/vasp.py b/nexus/library/vasp.py index 0d9edc26e..5b7c67a3e 100644 --- a/nexus/library/vasp.py +++ b/nexus/library/vasp.py @@ -186,7 +186,7 @@ class Vasp(Simulation): def generate_vasp(**kwargs): - sim_args,inp_args = Simulation.separate_inputs(kwargs,copy_pseudos=False) + sim_args,inp_args = Vasp.separate_inputs(kwargs,copy_pseudos=False) sim_args.input = generate_vasp_input(**inp_args) vasp = Vasp(**sim_args) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index dbe57c5d1..511efdd0e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -156,8 +156,12 @@ else() Particle/HDFWalkerOutput.cpp Particle/HDFWalkerInput_0_0.cpp Particle/HDFWalkerInput_0_4.cpp - ADIOS/ADIOS_verify.cpp ) + IF(HAVE_ADIOS) + SET(PARTICLEIO ${PARTICLEIO} + ADIOS/ADIOS_verify.cpp + ) + ENDIF(HAVE_ADIOS) IF(OHMMS_DIM MATCHES 3) SET(PARTICLEIO ${PARTICLEIO} ParticleIO/ESHDFParticleParser.cpp) @@ -181,14 +185,14 @@ else() SET(ADIOS ADIOS/ADIOS_config.cpp) - #if(HAVE_ADIOS) + IF(HAVE_ADIOS) SET(ADIOS ${ADIOS} Particle/AdiosWalkerInput.cpp ADIOS/ADIOS_profile.cpp ADIOS/ADIOS_verify.cpp ADIOS/ADIOSTrace.cpp ) - #endif(HAVE_ADIOS) + ENDIF(HAVE_ADIOS) ADD_LIBRARY(adios_config ${ADIOS}) diff --git a/src/Estimators/PostProcessor.cpp b/src/Estimators/PostProcessor.cpp index 625900d1b..7f15aeeba 100644 --- a/src/Estimators/PostProcessor.cpp +++ b/src/Estimators/PostProcessor.cpp @@ -67,15 +67,15 @@ namespace qmcplusplus APP_ABORT("PostProcessor::put xml input is incomplete, see messages above."); } - ParticleSet& Pq = *ppool.getParticleSet(Pqname); - ParticleSet& Pc = *ppool.getParticleSet(Pcname); - TrialWaveFunction& Psi = *wpool.getWaveFunction(Psiname); - QMCHamiltonian& H = *hpool.getHamiltonian(Hname); + ParticleSet* Pq = ppool.getParticleSet(Pqname); + ParticleSet* Pc = ppool.getParticleSet(Pcname); + TrialWaveFunction* Psi = wpool.getWaveFunction(Psiname); + QMCHamiltonian* H = hpool.getHamiltonian(Hname); - notPq = &Pq==0 || Pq.getName()!=Pqname; - notPc = &Pc==0 || Pc.getName()!=Pcname; - notPsi = &Psi==0 || Psi.getName()!=Psiname; - notH = &H==0 || H.getName()!=Hname; + notPq = Pq==0 || Pq->getName()!=Pqname; + notPc = Pc==0 || Pc->getName()!=Pcname; + notPsi = Psi==0 || Psi->getName()!=Psiname; + notH = H==0 || H->getName()!=Hname; if(notPq||notPc||notPsi||notH) { if(notPq) @@ -96,7 +96,7 @@ namespace qmcplusplus bool not_text = cname!="text"; PostProcessorBase* pp = 0; if(cname=="spindensity") - pp = new SpinDensityPostProcessor(Pq,Pc,H); + pp = new SpinDensityPostProcessor(*Pq,*Pc,*H); else if(not_text) APP_ABORT("PostProcessor::put "+cname+" is not a valid element of postprocess"); if(not_text) diff --git a/src/Lattice/CrystalLattice.h b/src/Lattice/CrystalLattice.h index 02a8e61c1..3d68b4580 100644 --- a/src/Lattice/CrystalLattice.h +++ b/src/Lattice/CrystalLattice.h @@ -303,6 +303,19 @@ struct CrystalLattice #endif } + /** conversion of a caresian reciprocal-vector to unit k-vector + *@param kin an input reciprocal vector in cartesian form + *@return k(reciprocal vector) as unit vector + */ + inline SingleParticlePos_t k_unit(const SingleParticlePos_t& kin) const + { +#ifdef OHMMS_LATTICEOPERATORS_H + return DotProduct::apply(R,kin)/TWOPI; +#else + return dot(R,kin)/TWOPI; +#endif + } + /** evaluate \f$k^2\f$ * *@param kin an input reciprocal vector in reciprocal-vector unit diff --git a/src/Lattice/Uniform3DGridLayout.cpp b/src/Lattice/Uniform3DGridLayout.cpp index 769bf7c65..17f55947d 100644 --- a/src/Lattice/Uniform3DGridLayout.cpp +++ b/src/Lattice/Uniform3DGridLayout.cpp @@ -284,11 +284,11 @@ int Uniform3DGridLayout::connectGrid(value_type int_rad, value_type con_rad) for(int ixyz=0; ixyz= NP[0]) + if((!BoxBConds[0] && d[0]<0) || d[0]>= NP[0]) continue; - if(!BoxBConds[1] && d[1]<0 || d[1]>= NP[1]) + if((!BoxBConds[1] && d[1]<0) || d[1]>= NP[1]) continue; - if(!BoxBConds[2] && d[2]<0 || d[2]>= NP[2]) + if((!BoxBConds[2] && d[2]<0) || d[2]>= NP[2]) continue; d[0]=d0[0]%NP[0]; d[0]= (d[0]<0)?d[0]+NP[0] : d[0]; diff --git a/src/LongRange/LRHandlerBase.h b/src/LongRange/LRHandlerBase.h index fa1644977..578a360d1 100644 --- a/src/LongRange/LRHandlerBase.h +++ b/src/LongRange/LRHandlerBase.h @@ -63,6 +63,10 @@ struct LRHandlerBase //constructor explicit LRHandlerBase(mRealType kc):LR_kc(kc) {} + + // virtual destructor + virtual ~LRHandlerBase() {} + //return r cutoff inline mRealType get_rc() const { diff --git a/src/LongRange/LRHandlerTemp.h b/src/LongRange/LRHandlerTemp.h index 8e9c33015..5d49b58d6 100644 --- a/src/LongRange/LRHandlerTemp.h +++ b/src/LongRange/LRHandlerTemp.h @@ -224,7 +224,10 @@ private: fillXk(breakuphandler.KList); //Allocate the space for the coefficients. coefs.resize(Basis.NumBasisElem()); //This must be after SetupKVecs. - breakuphandler.DoBreakup(Fk.data(),coefs.data()); //Fill array of coefficients. + + mRealType chisqr(0.0); + chisqr=breakuphandler.DoBreakup(Fk.data(),coefs.data()); //Fill array of coefficients. + app_log()<<"\n LR Breakup chi^2 = "< >& KList) diff --git a/src/Message/CommOperatorsMPI.h b/src/Message/CommOperatorsMPI.h index 0bc6ba8d4..1a23be9a9 100644 --- a/src/Message/CommOperatorsMPI.h +++ b/src/Message/CommOperatorsMPI.h @@ -93,28 +93,28 @@ template inline Communicate::request Communicate::irecv(int source, int tag, T& ) { APP_ABORT("Need specialization for irecv(int source, int tag, T& )"); - return 1; + return MPI_REQUEST_NULL; } template inline Communicate::request Communicate::isend(int dest, int tag, T&) { APP_ABORT("Need specialization for isend(int source, int tag, T& )"); - return 1; + return MPI_REQUEST_NULL; } template inline Communicate::request Communicate::irecv(int source, int tag, T* , int n) { APP_ABORT("Need specialization for irecv(int source, int tag, T*, int )"); - return 1; + return MPI_REQUEST_NULL; } template inline Communicate::request Communicate::isend(int dest, int tag, T*, int n) { APP_ABORT("Need specialization for isend(int source, int tag, T*, int )"); - return 1; + return MPI_REQUEST_NULL; } template<> diff --git a/src/Message/Communicate.h b/src/Message/Communicate.h index b8aa2451e..4b929c6f4 100644 --- a/src/Message/Communicate.h +++ b/src/Message/Communicate.h @@ -46,6 +46,7 @@ struct CommunicatorTraits typedef int status; typedef int request; typedef int intra_comm_type; + static const int MPI_REQUEST_NULL = 1; }; #define APP_ABORT(msg) \ diff --git a/src/Numerics/OneDimGridBase.h b/src/Numerics/OneDimGridBase.h index e5c34bd01..b0e9a6fc1 100644 --- a/src/Numerics/OneDimGridBase.h +++ b/src/Numerics/OneDimGridBase.h @@ -27,6 +27,7 @@ #include #include "OhmmsPETE/OhmmsVector.h" #include "Numerics/GridTraits.h" +#include "einspline/bspline_base.h" namespace qmcplusplus { @@ -68,6 +69,8 @@ struct OneDimGridBase virtual OneDimGridBase* makeClone() const =0; + virtual ~OneDimGridBase() {} + inline int getGridTag() const { return GridTag; @@ -83,7 +86,7 @@ struct OneDimGridBase { int k; int klo=0; - int khi=this->size()-1; + int khi=this->size(); while(khi-klo > 1) { k=(khi+klo) >> 1; @@ -337,7 +340,18 @@ struct LinearGrid: public OneDimGridBase for(int i=0; i { int k; int klo=0; - int khi=num_points-1; + int khi=num_points; //int khi=this->size()-1; while(khi-klo > 1) { diff --git a/src/OOMPI/Group.cc b/src/OOMPI/Group.cc index 9d0a2a194..52d6f9bcd 100644 --- a/src/OOMPI/Group.cc +++ b/src/OOMPI/Group.cc @@ -33,12 +33,14 @@ const OOMPI_Group OOMPI_GROUP_NULL(MPI_GROUP_NULL); static int Free_mpi_group(MPI_Group *ptr) { - if (ptr != 0 && *ptr != MPI_GROUP_NULL && *ptr != MPI_GROUP_EMPTY ) - if (OOMPI_COMM_WORLD.Finalized()) + if (ptr != 0 && *ptr != MPI_GROUP_NULL && *ptr != MPI_GROUP_EMPTY ) { + if (OOMPI_COMM_WORLD.Finalized()) { std::cerr << "Attempt to free group after finalize (ignored, probably resulting in memory leak)." << std::endl; - else + } else { MPI_Group_free(ptr); + } + } return MPI_SUCCESS; } diff --git a/src/OOMPI/Op.cc b/src/OOMPI/Op.cc index 1aa01c1d5..ac78de342 100644 --- a/src/OOMPI/Op.cc +++ b/src/OOMPI/Op.cc @@ -27,12 +27,14 @@ static int Free_mpi_op(MPI_Op *op_wrapper) { - if (op_wrapper != 0 && *op_wrapper != MPI_OP_NULL ) - if (OOMPI_COMM_WORLD.Finalized()) + if (op_wrapper != 0 && *op_wrapper != MPI_OP_NULL ) { + if (OOMPI_COMM_WORLD.Finalized()) { std::cerr << "Attempt to free operator after finalize (ignored, probably resulting in memory leak)." << std::endl; - else + } else { MPI_Op_free(op_wrapper); + } + } return MPI_SUCCESS; } diff --git a/src/Optimize/NRCOptimization.h b/src/Optimize/NRCOptimization.h index 9f5c9f448..71c8a80b8 100644 --- a/src/Optimize/NRCOptimization.h +++ b/src/Optimize/NRCOptimization.h @@ -56,7 +56,7 @@ struct sign2 { inline static float apply(float a, float b) { - return (b > 0.0)? abs(a): -abs(a); + return (b > 0.0)? std::abs(a): -std::abs(a); } }; diff --git a/src/Optimize/testDerivOptimization.h b/src/Optimize/testDerivOptimization.h index 91368fc1c..c2c8c0aef 100644 --- a/src/Optimize/testDerivOptimization.h +++ b/src/Optimize/testDerivOptimization.h @@ -222,7 +222,7 @@ public: for (int j = 0 ; j < NumParams ; j ++) gg += a_xi[j]*a_xi[j]; a_gtyp = sqrt(gg / (1.0*NumParams)); - for (int a_its = 0 ; ((a_its <= a_itmax)&(TargetFunc->IsValid)&(gg>CostTol)) ; ++a_its) + for (int a_its = 0 ; ((a_its <= a_itmax)&&(TargetFunc->IsValid)&(gg>CostTol)) ; ++a_its) { if (a_verbose > 0) printf("mac_it %d of %d : gg = %6.3g tol = %6.3g: \n", a_its , a_itmax , gg , CostTol) ; @@ -351,7 +351,7 @@ public: } if (a_verbose > 0) printf("FINAL : gg = %6.3g tol = %6.3g: \n", gg , CostTol) ; - return ((TargetFunc->IsValid) & (gg <= CostTol)); + return ((TargetFunc->IsValid) && (gg <= CostTol)); } @@ -369,13 +369,13 @@ public: int its(0); Return_t x = a_lastx / a_gtyp ; Return_t s = lineProduct(Parms , a_gx , x) ; - if (s < 0 & TargetFunc->IsValid) + if (s < 0 && TargetFunc->IsValid) { do { y = x * a_linmin_g1 ; t = lineProduct(Parms , a_gy , y) ; - if ((t >= 0.0) | (!TargetFunc->IsValid)) + if ((t >= 0.0) || (!TargetFunc->IsValid)) { if (a_verbose > 1) printf("s < 0: s = %6.3g: t = %6.3g; x = %6.3g y = %6.3g\n",s, t , x , y); @@ -399,16 +399,16 @@ public: } its++ ; } - while ((its <= a_linmin_maxits)&(TargetFunc->IsValid)) ; + while ((its <= a_linmin_maxits)&&(TargetFunc->IsValid)) ; } else - if (s > 0 & TargetFunc->IsValid) + if (s > 0 && TargetFunc->IsValid) { do { y = x * a_linmin_g3 ; t = lineProduct(Parms , a_gy , y) ; - if ((t <= 0.0) | (!TargetFunc->IsValid)) + if ((t <= 0.0) || (!TargetFunc->IsValid)) { if (a_verbose > 1) printf("s > 0: s = %6.3g: t = %6.3g; x = %6.3g y = %6.3g\n",s, t , x , y); @@ -432,7 +432,7 @@ public: } its ++ ; } - while ((its <= a_linmin_maxits)&(TargetFunc->IsValid)) ; + while ((its <= a_linmin_maxits)&&(TargetFunc->IsValid)) ; } else { @@ -444,7 +444,7 @@ public: Return_t u(0.0); Return_t m(0.5*(x+y)); int xyit(0); - if ((xyitGradTol) & (s*t<0.0) &(TargetFunc->IsValid)) + if ((xyitGradTol) && (s*t<0.0) &&(TargetFunc->IsValid)) { int XYBisectCounter=xybisect; do @@ -464,9 +464,9 @@ public: u=lineProduct(Parms , a_gunused , m) ; if (a_verbose > 1) printf("xyit %d of %d : u = %6.3g m = %6.3g: \n", xyit ,xycleanup, u , m) ; - if ((u==u)) + if (u==u) { - if (u*s >= 0.0 & TargetFunc->IsValid) + if (u*s >= 0.0 && TargetFunc->IsValid) { s=u; x=m; @@ -475,7 +475,7 @@ public: printf("s = %6.3g: t = %6.3g; x = %6.3g y = %6.3g\n",s, t , x , y); } else - if (u*t >= 0.0 & TargetFunc->IsValid) + if (u*t >= 0.0 && TargetFunc->IsValid) { t=u; y=m; @@ -499,10 +499,10 @@ public: } xyit++; } - while ((TargetFunc->IsValid) & (std::abs(s-t)>GradTol) & (xyitIsValid) && (std::abs(s-t)>GradTol) && (xyit0.0) | (!TargetFunc->IsValid)) + if ((s*t>0.0) || (!TargetFunc->IsValid)) { return -1.0; } diff --git a/src/QMCDrivers/QMCCorrelatedSamplingLinearOptimize.h b/src/QMCDrivers/QMCCorrelatedSamplingLinearOptimize.h index 95e12945d..3c80b5235 100644 --- a/src/QMCDrivers/QMCCorrelatedSamplingLinearOptimize.h +++ b/src/QMCDrivers/QMCCorrelatedSamplingLinearOptimize.h @@ -16,7 +16,7 @@ /** @file QMCCorrelatedSamplingLinearOptimize.h * @brief Definition of QMCDriver which performs VMC and optimization. */ -#ifndef QMCPLUSPLUS_QMCCSINEAROPTIMIZATION_VMCSINGLE_H +#ifndef QMCPLUSPLUS_QMCCSLINEAROPTIMIZATION_VMCSINGLE_H #define QMCPLUSPLUS_QMCCSLINEAROPTIMIZATION_VMCSINGLE_H #include "QMCDrivers/QMCLinearOptimize.h" diff --git a/src/QMCDrivers/QMCCostFunctionOMP.cpp b/src/QMCDrivers/QMCCostFunctionOMP.cpp index 69952fbdb..e32198c19 100644 --- a/src/QMCDrivers/QMCCostFunctionOMP.cpp +++ b/src/QMCDrivers/QMCCostFunctionOMP.cpp @@ -215,13 +215,13 @@ void QMCCostFunctionOMP::getConfigurations(const std::string& aroot) { H_KE_Node[ip]= new QMCHamiltonian; H_KE_Node[ip]->addOperator(hClones[ip]->getHamiltonian("Kinetic"),"Kinetic"); - if (includeNonlocalH.c_str()!="no") + if (includeNonlocalH!="no") { - QMCHamiltonianBase* a=hClones[ip]->getHamiltonian(includeNonlocalH.c_str()); + QMCHamiltonianBase* a=hClones[ip]->getHamiltonian(includeNonlocalH); if(a) { app_log()<<" Found non-local Hamiltonian element named "<addOperator(a,includeNonlocalH.c_str()); + H_KE_Node[ip]->addOperator(a,includeNonlocalH); } else app_log()<<" Did not find non-local Hamiltonian element named "<resize(wRef.getActiveWalkers(),NumOptimizables); } } - QMCHamiltonianBase* nlpp = (includeNonlocalH =="no")? 0: hClones[ip]->getHamiltonian(includeNonlocalH.c_str()); + QMCHamiltonianBase* nlpp = (includeNonlocalH =="no")? 0: hClones[ip]->getHamiltonian(includeNonlocalH); bool compute_nlpp=useNLPPDeriv && nlpp; //set the optimization mode for the trial wavefunction psiClones[ip]->startOptimization(); diff --git a/src/QMCHamiltonians/CMakeLists.txt b/src/QMCHamiltonians/CMakeLists.txt index 51d1f6e6a..e8c63ac22 100644 --- a/src/QMCHamiltonians/CMakeLists.txt +++ b/src/QMCHamiltonians/CMakeLists.txt @@ -30,6 +30,7 @@ SET(HAMSRCS DensityEstimator.cpp SkPot.cpp SkEstimator.cpp + SkAllEstimator.cpp MomentumEstimator.cpp ForceBase.cpp HamiltonianFactory.cpp diff --git a/src/QMCHamiltonians/DensityMatrices1B.cpp b/src/QMCHamiltonians/DensityMatrices1B.cpp index ed1084ef4..0ba6931a6 100644 --- a/src/QMCHamiltonians/DensityMatrices1B.cpp +++ b/src/QMCHamiltonians/DensityMatrices1B.cpp @@ -471,8 +471,6 @@ namespace qmcplusplus void DensityMatrices1B::get_required_traces(TraceManager& tm) { app_log()<<"dm1b get_required_traces"<< std::endl; - if(&Pq==NULL) - APP_ABORT("DensityMatrices1B::get_required_traces quantum particleset reference is NULL"); w_trace = tm.get_real_trace("weight"); if(energy_mat) { @@ -894,10 +892,16 @@ namespace qmcplusplus generate_density_samples(save,steps,rng); if(save) + { if(sampling==metropolis) + { sample_weights*=weight; + } else + { fill(sample_weights.begin(),sample_weights.end(),weight); + } + } // temporary check diff --git a/src/QMCHamiltonians/HamiltonianFactory.cpp b/src/QMCHamiltonians/HamiltonianFactory.cpp index 8322e29bd..c3ba3528b 100644 --- a/src/QMCHamiltonians/HamiltonianFactory.cpp +++ b/src/QMCHamiltonians/HamiltonianFactory.cpp @@ -46,6 +46,7 @@ #endif #if OHMMS_DIM == 3 #include "QMCHamiltonians/ChiesaCorrection.h" +#include "QMCHamiltonians/SkAllEstimator.h" #if defined(HAVE_LIBFFTW_LS) #include "QMCHamiltonians/ModInsKineticEnergy.h" #include "QMCHamiltonians/MomentumDistribution.h" @@ -468,6 +469,30 @@ bool HamiltonianFactory::build(xmlNodePtr cur, bool buildtree) ChiesaCorrection *chiesa = new ChiesaCorrection (source, psi); targetH->addOperator(chiesa,"KEcorr",false); } + else if(potType == "skall") + { + std::string SourceName = ""; + OhmmsAttributeSet attrib; + attrib.add(SourceName,"source"); + attrib.put(cur); + + PtclPoolType::iterator pit(ptclPool.find(SourceName)); + if(pit == ptclPool.end()) + { + APP_ABORT("Unknown source \"" + SourceName + "\" for LocalMoment."); + } + ParticleSet* source = (*pit).second; + + if(PBCType) + { + + SkAllEstimator* apot=new SkAllEstimator(*source, *targetPtcl); + apot->put(cur); + targetH->addOperator(apot,potName,false); + app_log()<<"Adding S(k) ALL estimator"< +#include +#include +#include +#include "Configuration.h" +#include +#include +#include +namespace qmcplusplus +{ + +SkAllEstimator::SkAllEstimator(ParticleSet& source, ParticleSet& target) +{ + app_log()<<"SkAllEstimator Constructor\n"; + elns= ⌖ + ions= &source; + NumeSpecies=elns->getSpeciesSet().getTotalNum(); + NumIonSpecies=ions->getSpeciesSet().getTotalNum(); + UpdateMode.set(COLLECTABLE,1); + + NumK=source.SK->KLists.numk; + OneOverN=1.0/static_cast(source.getTotalNum()); + Kshell=source.SK->KLists.kshell; + MaxKshell=Kshell.size()-1; + +#if defined(USE_REAL_STRUCT_FACTOR) + RhokTot_r.resize(NumK); + RhokTot_i.resize(NumK); + +#else + RhokTot.resize(NumK); +#endif + //for values, we are including e-e structure factor, and e-Ion. So a total of NumIonSpecies+1 structure factors. + //+2 for the real and imaginary parts of rho_k^e + values.resize(3*NumK); + Kmag.resize(MaxKshell); + OneOverDnk.resize(MaxKshell); + for(int ks=0, k=0; ksKLists.ksq[Kshell[ks]]); + OneOverDnk[ks]=1.0/static_cast(Kshell[ks+1]-Kshell[ks]); + } + hdf5_out=false; + +} + +void SkAllEstimator::resetTargetParticleSet(ParticleSet& P) +{ + elns = &P; +} + + +void SkAllEstimator::evaluateIonIon() +{ + std::stringstream ss; + ss << " kx ky kz"; + + std::ofstream skfile; + std::stringstream filebuffer; + skfile.open("rhok_IonIon.dat"); + for(int i=0; iSK->KLists.kpts_cart[k]; + + filebuffer<SK->rhok_r[i][k]; + rho_i=ions->SK->rhok_i[i][k]; +#else + rho_r=ions->SK->rhok[i][k].real(); + rho_i=ions->SK->rhok[i][k].imag(); +#endif + filebuffer<<" "<rhok_r[0],P.SK->rhok_r[0]+NumK,RhokTot_r.begin()); + std::copy(P.SK->rhok_i[0],P.SK->rhok_i[0]+NumK,RhokTot_i.begin()); + for(int i=1; irhok_r[i],P.SK->rhok_r[i]+NumK,RhokTot_r.begin()); + for(int i=1; irhok_i[i],P.SK->rhok_i[i]+NumK,RhokTot_i.begin()); + + for(int k=0; kSK->rhok_r[ionSpec][k]), rhok_A_i(ions->SK->rhok_i[ionSpec][k]); +// values[(ionSpec+1)*NumK+k]=RhokTot_r[k]*rhok_A_r+RhokTot_i[k]*rhok_A_i; +// } +// } + + for(int k=0,count=0; krhok[0],P.SK->rhok[0]+NumK,RhokTot.begin()); + for(int i=1; irhok[i],P.SK->rhok[i]+NumK,RhokTot.begin()); + Vector::const_iterator iit(RhokTot.begin()),iit_end(RhokTot.end()); + for(int k=0; kSK->rhok[ionSpec][k].real()), rhok_A_i(ions->SK->rhok[ionSpec][k].imag()); +// values[(ionSpec+1)*NumK+k]=rhok[k].real()*rhok_A_r+rho_k[k].imag()*rhok_A_i; +// } +// } + for(int k=0,count=0; k tmp(NumK); +// collectables.add(tmp.begin(),tmp.end()); +// } +// else +// { + + myIndex=plist.size(); +//First the electron structure factor + for (int i=0; i& h5desc + , hid_t gid) const +{ +/* if (hdf5_out) + { + vector ndim(1,NumK); + observable_helper* h5o=new observable_helper(myName); + h5o->set_dimensions(ndim,myIndex); + h5o->open(gid); + h5desc.push_back(h5o); + hsize_t kdims[2]; + kdims[0] = NumK; + kdims[1] = OHMMS_DIM; + string kpath = myName + "/kpoints"; + hid_t k_space = H5Screate_simple(2,kdims, NULL); + hid_t k_set = H5Dcreate (gid, kpath.c_str(), H5T_NATIVE_DOUBLE, k_space, H5P_DEFAULT); + hid_t mem_space = H5Screate_simple (2, kdims, NULL); + double *ptr = &(sourcePtcl->SK->KLists.kpts_cart[0][0]); + herr_t ret = H5Dwrite(k_set, H5T_NATIVE_DOUBLE, mem_space, k_space, H5P_DEFAULT, ptr); + H5Dclose (k_set); + H5Sclose (mem_space); + H5Sclose (k_space); + H5Fflush(gid, H5F_SCOPE_GLOBAL); + } */ +} + +bool SkAllEstimator::put(xmlNodePtr cur) +{ + OhmmsAttributeSet pAttrib; + std::string hdf5_flag="no"; + std::string write_ionion_flag="no"; + + pAttrib.add(hdf5_flag,"hdf5"); + pAttrib.add(hdf5_flag,"hdf5"); + pAttrib.add(write_ionion_flag,"writeionion"); + + pAttrib.put(cur); + if (hdf5_flag=="yes") + hdf5_out=true; + else + hdf5_out=false; + + if (write_ionion_flag=="Yes" || write_ionion_flag=="yes") + { + app_log()<<"SkAll: evaluateIonIon()\n"; + evaluateIonIon(); + } + + return true; +} + +bool SkAllEstimator::get(std::ostream& os) const +{ + return true; +} + +QMCHamiltonianBase* SkAllEstimator::makeClone(ParticleSet& qp + , TrialWaveFunction& psi) +{ + SkAllEstimator* myclone = new SkAllEstimator(*this); + myclone->hdf5_out=hdf5_out; + myclone->myIndex=myIndex; + return myclone; +} +} + +/*************************************************************************** + * $RCSfile$ $Author: jnkim $ + * $Revision: 2945 $ $Date: 2008-08-05 10:21:33 -0500 (Tue, 05 Aug 2008) $ + * $Id: ForceBase.h 2945 2008-08-05 15:21:33Z jnkim $ + ***************************************************************************/ diff --git a/src/QMCHamiltonians/SkAllEstimator.h b/src/QMCHamiltonians/SkAllEstimator.h new file mode 100644 index 000000000..c2d282f64 --- /dev/null +++ b/src/QMCHamiltonians/SkAllEstimator.h @@ -0,0 +1,99 @@ +////////////////////////////////////////////////////////////////// +// (c) Copyright 2008- by Jeongnim Kim +////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////// +// National Center for Supercomputing Applications & +// Materials Computation Center +// University of Illinois, Urbana-Champaign +// Urbana, IL 61801 +// e-mail: jnkim@ncsa.uiuc.edu +// +// Supported by +// National Center for Supercomputing Applications, UIUC +// Materials Computation Center, UIUC +////////////////////////////////////////////////////////////////// +// -*- C++ -*- +/** @file SkAllEstimator.h + * @brief Declare SkAllEstimator + */ +#ifndef QMCPLUSPLUS_SK_ALL_ESTIMATOR_H +#define QMCPLUSPLUS_SK_ALL_ESTIMATOR_H +#include +#include +namespace qmcplusplus +{ + +/** SkAllEstimator evaluate the structure factor of the target particleset + * + * + */ +class SkAllEstimator: public QMCHamiltonianBase +{ +public: + + SkAllEstimator(ParticleSet& ions, ParticleSet& elns); + + void resetTargetParticleSet(ParticleSet& P); + + Return_t evaluate(ParticleSet& P); + + inline Return_t evaluate(ParticleSet& P, std::vector& Txy) + { + return evaluate(P); + } + + void evaluateIonIon(); + + void addObservables(PropertySetType& plist); + void addObservables(PropertySetType& plist, BufferType& collectables); + void registerCollectables(std::vector& h5desc, hid_t gid) const ; + void setObservables(PropertySetType& plist); + void setParticlePropertyList(PropertySetType& plist, int offset); + bool put(xmlNodePtr cur); + bool get(std::ostream& os) const; + QMCHamiltonianBase* makeClone(ParticleSet& qp, TrialWaveFunction& psi); + +protected: +// ParticleSet *sourcePtcl; + ParticleSet *elns; + ParticleSet *ions; + /** number of species */ + int NumSpecies; + int NumeSpecies; + int NumIonSpecies; + /** number of kpoints */ + unsigned int NumK; + /** number of kshells */ + int MaxKshell; + /** normalization factor */ + RealType OneOverN; + /** kshell counters */ + std::vector Kshell; + /** instantaneous structure factor */ + std::vector Kmag; + /** 1.0/degenracy for a ksell */ + std::vector OneOverDnk; + /** \f$rho_k = \sum_{\alpha} \rho_k^{\alpha} \f$ for species index \f$\alpha\f$ */ +#if defined(USE_REAL_STRUCT_FACTOR) + Vector RhokTot_r,RhokTot_i; +#else + Vector RhokTot; +#endif + Vector values; + /** resize the internal data + * + * The argument list is not completed + */ + void resize(); + + bool hdf5_out; +}; + +} +#endif + +/*************************************************************************** + * $RCSfile$ $Author: jnkim $ + * $Revision: 2945 $ $Date: 2008-08-05 10:21:33 -0500 (Tue, 05 Aug 2008) $ + * $Id: ForceBase.h 2945 2008-08-05 15:21:33Z jnkim $ + ***************************************************************************/ diff --git a/src/QMCHamiltonians/VHXC.cpp b/src/QMCHamiltonians/VHXC.cpp index 38575bf68..2d5f1a541 100644 --- a/src/QMCHamiltonians/VHXC.cpp +++ b/src/QMCHamiltonians/VHXC.cpp @@ -42,8 +42,10 @@ VHXC::VHXC(ParticleSet& ptcl) : VHXC::~VHXC() { - if (VSpline) - destroy_Bspline(VSpline); + if (VSpline[0]) + destroy_Bspline(VSpline[0]); + if (VSpline[1]) + destroy_Bspline(VSpline[1]); } @@ -79,8 +81,12 @@ VHXC::init_spline() (grid0, grid1, grid2, bc0, bc1, bc2, VHXC_r_DP.data()); } - if (!VSpline[1]) - VSpline[1] = VSpline[1]; + // The following code does nothing, but I don't know what it is + // supposed to do. Maybe set VSpline[1] = VSpline[0]? If so, the + // destructor would need to be updated to hande the case where both + // VSpline elements points to the same spline structure. + //if (!VSpline[1]) + // VSpline[1] = VSpline[1]; } @@ -116,7 +122,7 @@ VHXC::evaluate (ParticleSet& P) VHXC::VHXC(ParticleSet& ptcl) : PtclRef(&ptcl), FirstTime(true) { - APP_ABORT("VHXC::VHXC cannot be used with einspline/fftw "); + APP_ABORT("VHXC::VHXC cannot be used without einspline/fftw "); } VHXC::~VHXC() diff --git a/src/QMCTools/trace_density.cpp b/src/QMCTools/trace_density.cpp index 2bcca7a43..c3ee06d95 100644 --- a/src/QMCTools/trace_density.cpp +++ b/src/QMCTools/trace_density.cpp @@ -465,11 +465,17 @@ struct Input for(int p=0;p0) + { for(int p=0;p ReOrderedBands; std::vector RejectedBands; diff --git a/src/QMCWaveFunctions/Jastrow/DiffOneBodySpinJastrowOrbital.h b/src/QMCWaveFunctions/Jastrow/DiffOneBodySpinJastrowOrbital.h index 735c9ebbe..a5068401d 100644 --- a/src/QMCWaveFunctions/Jastrow/DiffOneBodySpinJastrowOrbital.h +++ b/src/QMCWaveFunctions/Jastrow/DiffOneBodySpinJastrowOrbital.h @@ -80,14 +80,20 @@ public: delete_iter(gradLogPsi.begin(),gradLogPsi.end()); delete_iter(lapLogPsi.begin(),lapLogPsi.end()); if(Spin) + { for(int sg=0; sg0) + { for (int ib=0; ib0)&&(std::abs(t_params[ipm])>thr)) + { params[ib]=t_params[ipm++]; #else if((allowedOrbs(state,ib)>0)&&(std::abs(t_params[ipm])>thr)) + { params[ib]=t_params[ipm++]; #endif + } else + { params[ib]=0; + } + } + } } else + { if (params.size()) asize=params.size(); + } // for (int i=0; i< params.size(); i++) { int indx=0; for (int i=0; i< M; i++) diff --git a/src/config/stdlib/math.h b/src/config/stdlib/math.h index 3bf80bb0b..4a0c66762 100644 --- a/src/config/stdlib/math.h +++ b/src/config/stdlib/math.h @@ -23,11 +23,12 @@ #endif /* M_PI */ #endif /* TWOPI */ - +#if __cplusplus < 201103L inline float round(float x) { return roundf(x); } +#endif #if defined(HAVE_SINCOS) inline void sincos(float a, float* s, float* c) diff --git a/tests/test_automation/nightly_anl_cetus.sub b/tests/test_automation/nightly_anl_cetus.sub index 495ae491e..07d1fedc5 100755 --- a/tests/test_automation/nightly_anl_cetus.sub +++ b/tests/test_automation/nightly_anl_cetus.sub @@ -20,6 +20,13 @@ #Must be an absolute path place=/projects/QMCSim/QMCPACK_CI_BUILDS_DO_NOT_REMOVE +#define compiler +compiler=XL +#compiler=Clang++11 + +. /soft/environment/softenv/etc/softenv-aliases.sh +resoft ~/.soft.clang + if [ ! -e $place ]; then mkdir $place fi @@ -38,7 +45,7 @@ entry=trunk if [ -e $entry/CMakeLists.txt ]; then cd $entry -for sys in build_BGQ_XL_real #build_BGQ_XL_complex +for sys in build_BGQ_${compiler}_real build_BGQ_${compiler}_complex build_BGQ_${compiler}_MP_real build_BGQ_${compiler}_MP_complex do mkdir $sys @@ -53,15 +60,25 @@ then fi case $sys in -"build_BGQ_XL_real") -# Build the real version with XL compiler on BGQ - export QMCPACK_TEST_SUBMIT_NAME=BGQ-XL-Real-Release - ctest -DCMAKE_TOOLCHAIN_FILE=../config/BGQToolChain.cmake -DQMC_COMPLEX=0 -S $PWD/../CMake/ctest_script.cmake,release -R short -VV &> $place/log/$entry/$mydate/${QMCPACK_TEST_SUBMIT_NAME}.log +"build_BGQ_${compiler}_real") +# Build the real version with ${compiler} compiler on BGQ + export QMCPACK_TEST_SUBMIT_NAME=BGQ-${compiler}-Real-Release + ctest -DCMAKE_TOOLCHAIN_FILE=../config/BGQ_${compiler}_ToolChain.cmake -DQMC_COMPLEX=0 -S $PWD/../CMake/ctest_script.cmake,release -R short -E "dmc" -VV &> $place/log/$entry/$mydate/${QMCPACK_TEST_SUBMIT_NAME}.log ;; -"build_BGQ_XL_complex") -# Build the complex version with XL compiler on BGQ - export QMCPACK_TEST_SUBMIT_NAME=BGQ-XL-Complex-Release - ctest -DCMAKE_TOOLCHAIN_FILE=../config/BGQToolChain.cmake -DQMC_COMPLEX=1 -S $PWD/../CMake/ctest_script.cmake,release -R "short|unit|qe-LiH-unpolarized-no-collect" -VV &> $place/log/$entry/$mydate/${QMCPACK_TEST_SUBMIT_NAME}.log +"build_BGQ_${compiler}_complex") +# Build the complex version with ${compiler} compiler on BGQ + export QMCPACK_TEST_SUBMIT_NAME=BGQ-${compiler}-Complex-Release + ctest -DCMAKE_TOOLCHAIN_FILE=../config/BGQ_${compiler}_ToolChain.cmake -DQMC_COMPLEX=1 -S $PWD/../CMake/ctest_script.cmake,release -R "unit|short|qe-LiH-unpolarized-no-collect-np-16-12-12-nk-16-2-2" -E "dmc" -VV &> $place/log/$entry/$mydate/${QMCPACK_TEST_SUBMIT_NAME}.log + ;; +"build_BGQ_${compiler}_MP_real") +# Build the real version with ${compiler} compiler on BGQ + export QMCPACK_TEST_SUBMIT_NAME=BGQ-${compiler}-MixedPrec-Real-Release + ctest -DCMAKE_TOOLCHAIN_FILE=../config/BGQ_${compiler}_ToolChain.cmake -DQMC_COMPLEX=0 -DQMC_MIXED_PRECISION=1 -S $PWD/../CMake/ctest_script.cmake,release -R short -E "dmc" -VV &> $place/log/$entry/$mydate/${QMCPACK_TEST_SUBMIT_NAME}.log + ;; +"build_BGQ_${compiler}_MP_complex") +# Build the complex version with ${compiler} compiler on BGQ + export QMCPACK_TEST_SUBMIT_NAME=BGQ-${compiler}-MixedPrec-Complex-Release + ctest -DCMAKE_TOOLCHAIN_FILE=../config/BGQ_${compiler}_ToolChain.cmake -DQMC_COMPLEX=1 -DQMC_MIXED_PRECISION=1 -S $PWD/../CMake/ctest_script.cmake,release -R "unit|short|qe-LiH-unpolarized-no-collect-np-16-12-12-nk-16-2-2" -E "dmc" -VV &> $place/log/$entry/$mydate/${QMCPACK_TEST_SUBMIT_NAME}.log ;; *) echo "ERROR: Unknown build type $sys"