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

update new_optimizer branch with recent trunk changes. 

git-svn-id: https://subversion.assembla.com/svn/qmcdev/branches/new_optimizer@7231 e5b18d87-469d-4833-9cc0-8cdfa06e9491
This commit is contained in:
Ye Luo 2016-11-16 15:16:01 +00:00
parent 43a54af4dc 4c9560cab8
commit 685b69d5ae
50 changed files with 4040 additions and 778 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
config/BGQ_Clang++11_ToolChain.cmake Normal file
View File

@ -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)

0
config/BGQToolChain.cmake → config/BGQ_XL_ToolChain.cmake
View File

10
manual/installation.tex
View File

@ -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}

BIN
manual/qmcpack_manual.pdf
View File

Binary file not shown.

2
manual/qmcpack_manual.tex
View File

@ -168,6 +168,8 @@
\input{lab_condensed_matter}
\input{unit_testing}
\input{contributing}
\newpage

132
manual/unit_testing.tex Normal file
View File

@ -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<double> grid;
OneDimGridFunctor<double> 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.

2
nexus/library/gamess.py
View File

@ -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)

99
nexus/library/generic.py
View File

@ -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]

51
nexus/library/machines.py
View File

@ -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')

13
nexus/library/nexus.py
View File

@ -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__

13
nexus/library/physical_system.py
View File

@ -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)

27
nexus/library/project_manager.py
View File

@ -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:

171
nexus/library/pseudopotential.py
View File

@ -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:

26
nexus/library/pwscf.py
View File

@ -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

472
nexus/library/pwscf_input.py
View File

@ -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

25
nexus/library/qmcpack.py
View File

@ -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

9
nexus/library/qmcpack_input.py
View File

@ -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

2685
nexus/library/qmcpack_workflows.py
View File

File diff suppressed because it is too large Load Diff

156
nexus/library/simulation.py
View File

@ -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,

176
nexus/library/structure.py
View File

@ -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

2
nexus/library/vasp.py
View File

@ -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)

10
src/CMakeLists.txt
View File

@ -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})

18
src/Estimators/PostProcessor.cpp
View File

@ -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)

13
src/Lattice/CrystalLattice.h
View File

@ -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<SingleParticlePos_t,Tensor_t,ORTHO>::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

6
src/Lattice/Uniform3DGridLayout.cpp
View File

@ -284,11 +284,11 @@ int Uniform3DGridLayout::connectGrid(value_type int_rad, value_type con_rad)
for(int ixyz=0; ixyz<RShell.size(); ixyz++)
{
SingleParticleIndex_t d0 = org+RShell[ixyz],d;
if(!BoxBConds[0] && d[0]<0 || d[0]>= 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];

4
src/LongRange/LRHandlerBase.h
View File

@ -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
{

5
src/LongRange/LRHandlerTemp.h
View File

@ -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 = "<<chisqr<<std::endl;
}
void fillXk(std::vector<TinyVector<mRealType,2> >& KList)

8
src/Message/CommOperatorsMPI.h
View File

@ -93,28 +93,28 @@ template<typename T> 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<typename T> 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<typename T> 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<typename T> 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<>

1
src/Message/Communicate.h
View File

@ -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) \

18
src/Numerics/OneDimGridBase.h
View File

@ -27,6 +27,7 @@
#include <cmath>
#include "OhmmsPETE/OhmmsVector.h"
#include "Numerics/GridTraits.h"
#include "einspline/bspline_base.h"
namespace qmcplusplus
{
@ -68,6 +69,8 @@ struct OneDimGridBase
virtual OneDimGridBase<T,CT>* 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<T,CT>
for(int i=0; i<n; i++)
X[i] = ri+Delta_DP*i;
}
inline Ugrid einspline_grid()
{
Ugrid grid;
grid.start = lower_bound;
grid.end = upper_bound;
grid.num = num_points;
grid.delta = Delta;
grid.delta_inv= DeltaInv;
return grid;
}
};
/** One-Dimensional logarithmic-grid.
@ -505,7 +519,7 @@ struct NumericalGrid: public OneDimGridBase<T,CT>
{
int k;
int klo=0;
int khi=num_points-1;
int khi=num_points;
//int khi=this->size()-1;
while(khi-klo > 1)
{

8
src/OOMPI/Group.cc
View File

@ -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;
}

8
src/OOMPI/Op.cc
View File

@ -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;
}

2
src/Optimize/NRCOptimization.h
View File

@ -56,7 +56,7 @@ struct sign2<float>
{
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);
}
};

28
src/Optimize/testDerivOptimization.h
View File

@ -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 ((xyit<xycleanup) & (std::abs(s-t)>GradTol) & (s*t<0.0) &(TargetFunc->IsValid))
if ((xyit<xycleanup) && (std::abs(s-t)>GradTol) && (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) & (xyit<xycleanup)) ;
while ((TargetFunc->IsValid) && (std::abs(s-t)>GradTol) && (xyit<xycleanup)) ;
}
else
if ((s*t>0.0) | (!TargetFunc->IsValid))
if ((s*t>0.0) || (!TargetFunc->IsValid))
{
return -1.0;
}

2
src/QMCDrivers/QMCCorrelatedSamplingLinearOptimize.h
View File

@ -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"

8
src/QMCDrivers/QMCCostFunctionOMP.cpp
View File

@ -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 "<<includeNonlocalH<< std::endl;
H_KE_Node[ip]->addOperator(a,includeNonlocalH.c_str());
H_KE_Node[ip]->addOperator(a,includeNonlocalH);
}
else
app_log()<<" Did not find non-local Hamiltonian element named "<<includeNonlocalH<< std::endl;
@ -365,7 +365,7 @@ void QMCCostFunctionOMP::checkConfigurations()
HDerivRecords[ip]->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();

1
src/QMCHamiltonians/CMakeLists.txt
View File

@ -30,6 +30,7 @@ SET(HAMSRCS
DensityEstimator.cpp
SkPot.cpp
SkEstimator.cpp
SkAllEstimator.cpp
MomentumEstimator.cpp
ForceBase.cpp
HamiltonianFactory.cpp

8
src/QMCHamiltonians/DensityMatrices1B.cpp
View File

@ -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

25
src/QMCHamiltonians/HamiltonianFactory.cpp
View File

@ -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"<<std::endl;
}
}
#endif
else if(potType == "Pressure")
{

2
src/QMCHamiltonians/QMCHamiltonianBase.cpp
View File

@ -80,7 +80,7 @@ void QMCHamiltonianBase::two_body_quantum_domain(const ParticleSet& P1,const Par
bool q2 = P2.is_quantum();
if(c1 && c2)
quantum_domain = classical_classical;
else if(q1 && c2 || c1 && q2)
else if((q1 && c2) || (c1 && q2))
quantum_domain = quantum_classical;
else if(q1 && q2)
quantum_domain = quantum_quantum;

324
src/QMCHamiltonians/SkAllEstimator.cpp Normal file
View File

@ -0,0 +1,324 @@
//////////////////////////////////////////////////////////////////
// (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++ -*-
#include <QMCHamiltonians/SkAllEstimator.h>
#include <LongRange/StructFact.h>
#include <Utilities/IteratorUtility.h>
#include <OhmmsData/AttributeSet.h>
#include "Configuration.h"
#include <fstream>
#include <sstream>
#include <string>
namespace qmcplusplus
{
SkAllEstimator::SkAllEstimator(ParticleSet& source, ParticleSet& target)
{
app_log()<<"SkAllEstimator Constructor\n";
elns= &target;
ions= &source;
NumeSpecies=elns->getSpeciesSet().getTotalNum();
NumIonSpecies=ions->getSpeciesSet().getTotalNum();
UpdateMode.set(COLLECTABLE,1);
NumK=source.SK->KLists.numk;
OneOverN=1.0/static_cast<RealType>(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; ks<MaxKshell; ks++)
{
Kmag[ks]=std::sqrt(source.SK->KLists.ksq[Kshell[ks]]);
OneOverDnk[ks]=1.0/static_cast<RealType>(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; i<NumIonSpecies; i++)
ss<<" rho_"<<i<<"_r"<<" rho_"<<i<<"_i";
filebuffer<<ss.str()<<std::endl;
for (int k=0; k<NumK; k++)
{
PosType kvec= ions->SK->KLists.kpts_cart[k];
filebuffer<<kvec;
for (int i=0; i<NumIonSpecies; i++)
{
double rho_i(0.0);
double rho_r(0.0);
#if defined(USE_REAL_STRUCT_FACTOR)
rho_r=ions->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<<" "<<rho_r<<" "<<rho_i;
}
filebuffer<<std::endl;
}
skfile<<filebuffer.str();
skfile.close();
}
SkAllEstimator::Return_t SkAllEstimator::evaluate(ParticleSet& P)
{
#if defined(USE_REAL_STRUCT_FACTOR)
//sum over species
std::copy(P.SK->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; i<NumeSpecies; ++i)
accumulate_elements(P.SK->rhok_r[i],P.SK->rhok_r[i]+NumK,RhokTot_r.begin());
for(int i=1; i<NumeSpecies; ++i)
accumulate_elements(P.SK->rhok_i[i],P.SK->rhok_i[i]+NumK,RhokTot_i.begin());
for(int k=0; k<NumK; k++)
values[k]=RhokTot_r[k]*RhokTot_r[k]+RhokTot_i[k]*RhokTot_i[k];
// for (int ionSpec=0; ionSpec<NumIonSpecies; ionSpec++)
// {
// for(int k=0; k<NumK; k++)
// {
// RealType rhok_A_r(ions->SK->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; k<NumK; k++)
{
values[NumK+2*k]=RhokTot_r[k];
values[NumK+2*k+1]=RhokTot_i[k];
}
#else
//sum over species
std::copy(P.SK->rhok[0],P.SK->rhok[0]+NumK,RhokTot.begin());
for(int i=1; i<NumeSpecies; ++i)
accumulate_elements(P.SK->rhok[i],P.SK->rhok[i]+NumK,RhokTot.begin());
Vector<ComplexType>::const_iterator iit(RhokTot.begin()),iit_end(RhokTot.end());
for(int k=0; k<NumK; k++)
values[k]=rhok[k].real()*rhok[k].real()+rhok[k].imag()*rhok[k].imag();
// for (int ionSpec=0; ionSpec<NumIonSpecies; ionSpec++)
// {
// for(int k=0; k<NumK; k++)
// {
// RealType rhok_A_r(ions->SK->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<NumK; k++)
{
value[NumK+count]=rhok[k].real();
count++;
value[NumK+count]=rhok[k].imag();
count++;
}
#endif
return 0.0;
}
void SkAllEstimator::addObservables(PropertySetType& plist, BufferType& collectables)
{
// if(hdf5_out)
// {
// myIndex=collectables.size();
// vector<RealType> tmp(NumK);
// collectables.add(tmp.begin(),tmp.end());
// }
// else
// {
myIndex=plist.size();
//First the electron structure factor
for (int i=0; i<NumK; i++)
{
std::stringstream sstr;
sstr << "rhok_e_e_" <<i;
int id=plist.add(sstr.str());
}
//Now the e-Ion structure factors. IonIon are dumped to file, and not evaluated.
// for (int ionSpec=0; ionSpec<NumIonSpecies; ionSpec++)
// {
// for (int i=0; i<NumK; i++)
// {
// std::stringstream sstr;
// sstr << "rhok_e_" <<ionSpec<<"_"<<i;
// int id=plist.add(sstr.str());
// }
// }
for(int k=0; k<NumK; k++)
{
std::stringstream sstr1,sstr2;
sstr1 << "rhok_e_r_"<<k;
sstr2 << "rhok_e_i_"<<k;
int id=plist.add(sstr1.str());
int id2=plist.add(sstr2.str());
}
// }
}
void SkAllEstimator::addObservables(PropertySetType& plist )
{
myIndex=plist.size();
//First the electron structure factor
for (int i=0; i<NumK; i++)
{
std::stringstream sstr;
sstr << "rhok_e_e_" <<i;
int id=plist.add(sstr.str());
}
//Now the e-Ion structure factors. IonIon are dumped to file, and not evaluated.
// for (int ionSpec=0; ionSpec<NumIonSpecies; ionSpec++)
// {
// for (int i=0; i<NumK; i++)
// {
// std::stringstream sstr;
// sstr << "rhok_e_" <<ionSpec<<"_"<<i;
// int id=plist.add(sstr.str());
// }
// }
for(int k=0; k<NumK; k++)
{
std::stringstream sstr1,sstr2;
sstr1 << "rhok_e_r_"<<k;
sstr2 << "rhok_e_i_"<<k;
int id=plist.add(sstr1.str());
int id2=plist.add(sstr2.str());
}
}
void SkAllEstimator::setObservables(PropertySetType& plist)
{
if (!hdf5_out)
std::copy(values.begin(),values.end(),plist.begin()+myIndex);
}
void SkAllEstimator::setParticlePropertyList(PropertySetType& plist
, int offset)
{
if (!hdf5_out)
std::copy(values.begin(),values.end(),plist.begin()+myIndex+offset);
}
void SkAllEstimator::registerCollectables(std::vector<observable_helper*>& h5desc
, hid_t gid) const
{
/* if (hdf5_out)
{
vector<int> 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 $
***************************************************************************/

99
src/QMCHamiltonians/SkAllEstimator.h Normal file
View File

@ -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 <QMCHamiltonians/QMCHamiltonianBase.h>
#include <vector>
namespace qmcplusplus
{
/** SkAllEstimator evaluate the structure factor of the target particleset
*
* <estimator name="sk" type="sk" debug="no"/>
*/
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<NonLocalData>& Txy)
{
return evaluate(P);
}
void evaluateIonIon();
void addObservables(PropertySetType& plist);
void addObservables(PropertySetType& plist, BufferType& collectables);
void registerCollectables(std::vector<observable_helper*>& 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<int> Kshell;
/** instantaneous structure factor */
std::vector<RealType> Kmag;
/** 1.0/degenracy for a ksell */
std::vector<RealType> OneOverDnk;
/** \f$rho_k = \sum_{\alpha} \rho_k^{\alpha} \f$ for species index \f$\alpha\f$ */
#if defined(USE_REAL_STRUCT_FACTOR)
Vector<RealType> RhokTot_r,RhokTot_i;
#else
Vector<ComplexType> RhokTot;
#endif
Vector<RealType> 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 $
***************************************************************************/

16
src/QMCHamiltonians/VHXC.cpp
View File

@ -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()

8
src/QMCTools/trace_density.cpp
View File

@ -465,11 +465,17 @@ struct Input
for(int p=0;p<tot_procs;++p)
procs.push_back(p);
if(groups.size()==0)
{
if(tot_groups>0)
{
for(int p=0;p<tot_groups;++p)
groups.push_back(p);
}
else
{
groups.push_back(-1);
}
}
Peln = get_particleset(elns);
Pion = get_particleset(ions);
if(tot_series<1)
@ -689,6 +695,8 @@ struct QuantityAnalyzer
block_size = 0;
}
virtual ~QuantityAnalyzer() {}
void init(xmlNodePtr xml,int s)
{
series_index = s;

2
src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp
View File

@ -541,6 +541,7 @@ void EinsplineSetBuilder::OccupyBands_ESHDF(int spin, int sortBands, int numOrbs
for(int ien=0; ien<SortBands.size(); ien++)
{
if((Occ[cnt] == SortBands[ien].TwistIndex)&&(Occ[cnt+1] == SortBands[ien].BandIndex))
{
if(cnt<particle_hole_pairs*2)
{
gsOcc[ien]-=1;
@ -553,6 +554,7 @@ void EinsplineSetBuilder::OccupyBands_ESHDF(int spin, int sortBands, int numOrbs
app_log()<<"adding orbital "<<ien<< std::endl;
cnt+=2;
}
}
}
std::vector<BandInfo> ReOrderedBands;
std::vector<BandInfo> RejectedBands;

6
src/QMCWaveFunctions/Jastrow/DiffOneBodySpinJastrowOrbital.h
View File

@ -80,14 +80,20 @@ public:
delete_iter(gradLogPsi.begin(),gradLogPsi.end());
delete_iter(lapLogPsi.begin(),lapLogPsi.end());
if(Spin)
{
for(int sg=0; sg<F.rows(); ++sg)
for(int tg=0; tg<F.cols(); ++tg)
if(F(sg,tg))
{
delete F(sg,tg);
}
else
{
for(int sg=0; sg<F.rows(); ++sg)
if(F(sg,0))
delete F(sg,0);
}
}
}
/** Add a radial functor for a group

6
src/QMCWaveFunctions/Jastrow/OneBodySpinJastrowOrbital.h
View File

@ -94,14 +94,20 @@ public:
~OneBodySpinJastrowOrbital()
{
if(Spin)
{
for(int sg=0; sg<F.rows(); ++sg)
for(int tg=0; tg<F.cols(); ++tg)
if(F(sg,tg))
{
delete F(sg,tg);
}
else
{
for(int sg=0; sg<F.rows(); ++sg)
if(F(sg,0))
delete F(sg,0);
}
}
}
//evaluate the distance table with P

11
src/QMCWaveFunctions/OptimizableSPOSet.cpp
View File

@ -230,20 +230,31 @@ OptimizableSPOSet::put (xmlNodePtr node, SPOPool_t &spo_pool)
params.resize(asize,0.0);
int ipm(0);
if (t_params.size()>0)
{
for (int ib=0; ib<M; ib++)
{
#if defined(QMC_COMPLEX)
if((allowedOrbs(state,ib).real()>0)&&(std::abs(t_params[ipm])>thr))
{
params[ib]=t_params[ipm++];
#else
if((allowedOrbs(state,ib)>0)&&(std::abs(t_params[ipm])>thr))
{
params[ib]=t_params[ipm++];
#endif
}
else
{
params[ib]=0;
}
}
}
}
else
{
if (params.size())
asize=params.size();
}
// for (int i=0; i< params.size(); i++) {
int indx=0;
for (int i=0; i< M; i++)

3
src/config/stdlib/math.h
View File

@ -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)

35
tests/test_automation/nightly_anl_cetus.sub
View File

@ -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"