Ignore:
Timestamp:
Mar 9, 2017, 9:43:23 AM (8 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_StructOpt_integration_tests, AutomationFragmentation_failures, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_ChronosMutex, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, GeometryObjects, Gui_displays_atomic_force_velocity, IndependentFragmentGrids_IntegrationTest, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, PartialCharges_OrthogonalSummation, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, RotateToPrincipalAxisSystem_UndoRedo, StoppableMakroAction, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, Ubuntu_1604_changes, stable
Children:
7672284
Parents:
ba1152
git-author:
Frederik Heber <heber@…> (03/08/17 08:11:59)
git-committer:
Frederik Heber <heber@…> (03/09/17 09:43:23)
Message:

HUGE: Extracted libmolecuilder_mpqc that is now linked to poolworker.

  • This stops the problems with MemDebug and mpqc when linking against libLinearAlgebra in debug mode: static global variables in libLinAlg are allocated using MemDebug (prefixed with a checksum) but are deallocated using normal libc's free() on exit. This causes an invalid free() as the ptr given to free points inside the block and not at its start. The problem comes from mpqc's use of own new and delete implementation. I'm guessing its reference counting is the culprit. Hence, we need to separate these two compilations from another in different units/libraries. Therefore, we have split off libmolecuilder_mpqc, .._mpqc_extract and additionally place the MPQCJob::Work() implementation inside libMolecuilderJobs_Work.
  • libmolecuilder_mpqc contains all MPQC's code in mpqc.cc (and linked libraries) that is not the main() function.
  • libmolecuilder_mpqc_extract contains functions that extract results such as energies, forces, charge grids from the obtained mpqc solution. These were added by myself to the mpqc code before.
  • molecuilder_mpqc is then linked against a NoOp .._extract library version. Thereby, it does not use any of the Molecuilder or related libraries and does not come in contact with MemDebug.
  • molecuilder_poolworker however is linked with the full .._extract library and hence performs these extractions.
  • poolworker now executes MPQCJob, MPQCCommandJob, and VMGJob and therefore needs to enforce binding to all of them.
  • TESTS: renamed molecuilder_mpqc.in to molecuilder_poolworker.in.
Location:
ThirdParty/mpqc_open/src/bin/mpqc
Files:
4 added
2 edited

Legend:

Unmodified
Added
Removed

  • ThirdParty/mpqc_open/src/bin/mpqc/Makefile.am

    rba1152 rfbf005  
    1 AM_CPPFLAGS = -I$(top_srcdir)/include -I$(top_srcdir)/src/lib
     1AM_CPPFLAGS = \
     2        -I$(top_srcdir)/include -I$(top_srcdir)/src/lib \
     3        -DHAVE_JOBMARKET -I$(top_srcdir)/../JobMarket/src -I$(top_builddir)/../JobMarket -I$(top_srcdir)/../../src -I$(top_srcdir)/../LinearAlgebra/src -I$(top_srcdir)/../CodePatterns/src $(BOOST_SYSTEM_CFLAGS) \
     4        -DHAVE_MPQCDATA -I$(top_builddir)/../../src
    25
    36bin_PROGRAMS = molecuilder_mpqc
    47
    5 noinst_LTLIBRARIES = libmpqc.la
     8lib_LTLIBRARIES = \
     9        libmolecuilder_mpqc.la
    610
    7 libmpqc_la_SOURCES = \
     11noinst_LTLIBRARIES = \
     12        libmolecuilder_mpqc_extract_dummy.la
     13
     14LIBMPQCSOURCES = \
     15        mpqc.cc \
    816        mpqcin.cc \
     17        parse.cc \
     18        scan.cc
     19
     20LIBMPQCHEADERS = \
     21        mpqc.h \
    922        mpqcin.h \
    10         parse.cc \
    1123        parse.h \
    12         scan.cc \
    1324        scdirlist.h
    1425
     26libmolecuilder_mpqc_extract_dummy_la_SOURCES = \
     27        mpqc_extract_dummy.cc \
     28        mpqc_extract.h
     29
     30libmolecuilder_mpqc_la_includedir = $(includedir)/molecuilder_mpqc/
     31libmolecuilder_mpqc_la_CPPFLAGS = $(AM_CPPFLAGS)
     32libmolecuilder_mpqc_la_LDFLAGS = $(AM_LDFLAGS)
     33       
     34libmolecuilder_mpqc_la_LIBADD = \
     35        ../../lib/libSCmbpt.la ../../lib/libSCscf.la ../../lib/libSCdft.la ../../lib/libSCwfn.la ../../lib/libSCpsi.la ../../lib/libSCsolvent.la ../../lib/libSCintv3.la ../../lib/libSCbasis.la ../../lib/libSCoint3.la ../../lib/libSCmolecule.la ../../lib/libSCisosurf.la ../../lib/libSCoptimize.la ../../lib/libSCsymmetry.la ../../lib/libSCscmat.la ../../lib/libSCrender.la ../../lib/libSCgroup.la ../../lib/libSCmisc.la ../../lib/libSCstate.la ../../lib/libSCkeyval.la ../../lib/libSCclass.la ../../lib/libSCcontainer.la ../../lib/libSCref.la ../../lib/libSCoptions.la
     36
     37nobase_libmolecuilder_mpqc_la_include_HEADERS = ${LIBMPQCHEADERS}
     38
     39## Define the source file list for the "libexample-@MOLECUILDER_API_VERSION@.la"
     40## target.  Note that @MOLECUILDER_API_VERSION@ is not interpreted by Automake and
     41## will therefore be treated as if it were literally part of the target name,
     42## and the variable name derived from that.
     43## The file extension .cc is recognized by Automake, and makes it produce
     44## rules which invoke the C++ compiler to produce a libtool object file (.lo)
     45## from each source file.  Note that it is not necessary to list header files
     46## which are already listed elsewhere in a _HEADERS variable assignment.
     47libmolecuilder_mpqc_la_SOURCES = ${LIBMPQCSOURCES}
     48
     49## Instruct libtool to include ABI version information in the generated shared
     50## library file (.so).  The library ABI version is defined in configure.ac, so
     51## that all version information is kept in one place.
     52libmolecuilder_mpqc_la_LDFLAGS += -version-info $(MPQC_SO_VERSION)
     53
     54## The generated configuration header is installed in its own subdirectory of
     55## $(libdir).  The reason for this is that the configuration information put
     56## into this header file describes the target platform the installed library
     57## has been built for.  Thus the file must not be installed into a location
     58## intended for architecture-independent files, as defined by the Filesystem
     59## Hierarchy Standard (FHS).
     60## The nodist_ prefix instructs Automake to not generate rules for including
     61## the listed files in the distribution on 'make dist'.  Files that are listed
     62## in _HEADERS variables are normally included in the distribution, but the
     63## configuration header file is generated at configure time and should not be
     64## shipped with the source tarball.
     65libmolecuilder_mpqc_la_libincludedir = $(libdir)/molecuilder_mpqc/include
     66nodist_libmolecuilder_mpqc_la_libinclude_HEADERS = $(top_builddir)/src/lib/scconfig.h
     67
     68## Install the generated pkg-config file (.pc) into the expected location for
     69## architecture-dependent package configuration information.  Occasionally,
     70## pkg-config files are also used for architecture-independent data packages,
     71## in which case the correct install location would be $(datadir)/pkgconfig.
     72#pkgconfigdir = $(libdir)/pkgconfig
     73#pkgconfig_DATA = $(top_builddir)/molecuilder_mpqc.pc
     74
    1575molecuilder_mpqc_SOURCES = \
    16         mpqc.cc
     76        mpqc_main.cc
    1777
    18 molecuilder_mpqc_CPPFLAGS = \
    19         $(AM_CPPFLAGS) \
    20         -DHAVE_JOBMARKET -I$(top_srcdir)/../JobMarket/src -I$(top_builddir)/../JobMarket -I$(top_srcdir)/../../src -I$(top_srcdir)/../LinearAlgebra/src $(BOOST_SYSTEM_CFLAGS) \
    21         -DHAVE_MPQCDATA -I$(top_builddir)/../../src
    22 
    23 molecuilder_mpqc_LDFLAGS = \
    24         $(AM_LDFLAGS) \
    25         -L$(top_builddir)/../JobMarket/src/JobMarket/.libs -Wl,-rpath,$(top_builddir)/../JobMarket/src/JobMarket/.libs \
    26         -L$(top_builddir)/../../src/.libs -Wl,-rpath,$(top_builddir)/../../src/.libs \
    27         $(BOOST_SYSTEM_LDFLAGS)
    2878molecuilder_mpqc_LDADD = \
    29         libmpqc.la ../../lib/libSCmbpt.la ../../lib/libSCscf.la ../../lib/libSCdft.la ../../lib/libSCwfn.la ../../lib/libSCpsi.la ../../lib/libSCsolvent.la ../../lib/libSCintv3.la ../../lib/libSCbasis.la ../../lib/libSCoint3.la ../../lib/libSCmolecule.la ../../lib/libSCisosurf.la ../../lib/libSCoptimize.la ../../lib/libSCsymmetry.la ../../lib/libSCscmat.la ../../lib/libSCrender.la ../../lib/libSCgroup.la ../../lib/libSCmisc.la ../../lib/libSCstate.la ../../lib/libSCkeyval.la ../../lib/libSCclass.la ../../lib/libSCcontainer.la ../../lib/libSCref.la ../../lib/libSCoptions.la \
    30         -lJobMarketPoolWorker -lJobMarket  \
    31         -lMolecuilderJobs -lMolecuilderFragmentationSummation \
    32         -lboost_serialization
     79        libmolecuilder_mpqc.la libmolecuilder_mpqc_extract_dummy.la
    3380if COND_LIBINT
    3481molecuilder_mpqc_LDADD += -lint

  • ThirdParty/mpqc_open/src/bin/mpqc/mpqc.cc

    rba1152 rfbf005  
    123123#endif
    124124
    125 #ifdef HAVE_JOBMARKET
    126 // include headers that implement a archive in simple text format
    127 // otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect
    128 #include <boost/archive/text_oarchive.hpp>
    129 #include <boost/archive/text_iarchive.hpp>
    130 #include <boost/bind.hpp>
    131 #include <boost/tokenizer.hpp>
    132 
    133 #include "JobMarket/Results/FragmentResult.hpp"
    134 #include "JobMarket/poolworker_main.hpp"
    135 
    136 #include "chemistry/qc/scf/scfops.h"
    137 
    138 #ifdef HAVE_MPQCDATA
    139 #include "Jobs/MPQCJob.hpp"
    140 #include "Fragmentation/Summation/Containers/MPQCData.hpp"
    141 
    142 #include <chemistry/qc/basis/obint.h>
    143 #include <chemistry/qc/basis/symmint.h>
    144 #endif
    145 
    146 #include <algorithm>
    147 #include <stdlib.h>
    148 #endif
    149 
    150125using namespace std;
    151126using namespace sc;
    152127
    153128#include "mpqcin.h"
     129#include "mpqc.h"
     130#include "mpqc_extract.h"
    154131
    155132//////////////////////////////////////////////////////////////////////////
    … …  
    178155}
    179156
    180 static void
    181 clean_up(void)
    182 {
    183   MemoryGrp::set_default_memorygrp(0);
    184   ThreadGrp::set_default_threadgrp(0);
    185   SCMatrixKit::set_default_matrixkit(0);
    186   Integral::set_default_integral(0);
    187   RegionTimer::set_default_regiontimer(0);
    188 }
     157namespace detail {
     158
     159  void
     160  clean_up(void)
     161  {
     162    ::MemoryGrp::set_default_memorygrp(0);
     163    ::ThreadGrp::set_default_threadgrp(0);
     164    ::SCMatrixKit::set_default_matrixkit(0);
     165    ::Integral::set_default_integral(0);
     166    ::RegionTimer::set_default_regiontimer(0);
     167  }
     168
     169} /* namespace detail */
    189170
    190171static void
    … …  
    335316}
    336317
    337 /** Temporary structure for storing information from command-line
    338  *
    339  * This structure has been introduced to gather the various calls to GetLongOpts
    340  * at one (initial) place and to abstract it from the source of command-lines.
    341  * This temporary object can be set by other means, too. I.e. we become
    342  * independent of usage in command-line programs.
    343  */
    344 struct OptionValues {
    345   const char *keyvalue; // option "k"
    346   const char *debug; // option ""
    347   int limit; // option "l"
    348   const char *check; // option "c"
    349   const char *simple_input; // option "i"
    350   string executablename;
    351 
    352 #ifdef HAVE_CHEMISTRY_CCA
    353   string cca_load;  // option "cca-path"
    354   string cc_path; // option "cca-load"
    355 #endif
    356 };
    357 
    358318/** Parse remainder options not treated by ComputeOptions() into temporary storage.
    359319 *
    … …  
    364324void parseRemainderOptions(
    365325    GetLongOpt &options,
    366     struct OptionValues &values,
     326    detail::OptionValues &values,
    367327    int argc,
    368328    char **argv)
    … …  
    655615/** Parse the input configuration from char array into keyvalue container.
    656616 *
    657  * \param parsedkv key value container to foll
     617 * \param parsedkv key value container to fill
    658618 * \param values temporary options value structure
    659619 * \param in_char_array char array with input file
    … …  
    663623void parseIntoKeyValue(
    664624    Ref<ParsedKeyVal> &parsedkv,
    665     struct OptionValues &values,
     625    detail::OptionValues &values,
    666626    char *&in_char_array,
    667627    int use_simple_input)
    … …  
    699659    Ref<MessageGrp> &grp,
    700660    Ref<ParsedKeyVal> &parsedkv,
    701     struct OptionValues &values,
     661    detail::OptionValues &values,
    702662    const char *&input,
    703663    const char *&generic_input,
    … …  
    804764void prepareCCA(
    805765    Ref<KeyVal> &keyval,
    806     struct OptionValues &values
     766    detail::OptionValues &values
    807767    )
    808768{
    … …  
    851811    Ref<MessageGrp> &grp,
    852812    Ref<Debugger> &debugger,
    853     struct OptionValues &values)
     813    detail::OptionValues &values)
    854814{
    855815  debugger << keyval->describedclassvalue("debug");
    … …  
    971931int checkBasisSetLimit(
    972932    Ref<MolecularEnergy> &mole,
    973     struct OptionValues &values
     933    detail::OptionValues &values
    974934    )
    975935{
    … …  
    11821142}
    11831143
    1184 static int getCoreElectrons(const int z)
    1185 {
    1186   int n=0;
    1187   if (z > 2) n += 2;
    1188   if (z > 10) n += 8;
    1189   if (z > 18) n += 8;
    1190   if (z > 30) n += 10;
    1191   if (z > 36) n += 8;
    1192   if (z > 48) n += 10;
    1193   if (z > 54) n += 8;
    1194   return n;
    1195 }
    1196 
    11971144void init()
    11981145{
    … …  
    12001147
    12011148  int i;
    1202   atexit(clean_up);
     1149  atexit(detail::clean_up);
    12031150
    12041151#ifdef HAVE_FEENABLEEXCEPT
    … …  
    12331180}
    12341181
    1235 /** Finds the region index to a given timer region name.
    1236  *
    1237  * @param nregion number of regions
    1238  * @param region_names array with name of each region
    1239  * @param name name of desired region
    1240  * @return index of desired region in array
    1241  */
    1242 int findTimerRegion(const int &nregion, const char **&region_names, const char *name)
    1243 {
    1244   size_t region=0;
    1245   for (;region<nregion;++region) {
    1246     //std::cout << "Comparing " << region_names[region] << " and " << name << "." << std::endl;
    1247     if (strcmp(region_names[region], name) == 0)
    1248       break;
    1249   }
    1250   if (region == nregion)
    1251     region = 0;
    1252   return region;
    1253 }
    1254 
    12551182/** Performs the main work to calculate the ground state energies, gradients, etc.
    12561183 *
    1257  * @param grp message group
    1258  * @param values temporary value storage for parsed command-line
    1259  * @param input input file name
    1260  * @param generic_input generic input file name
    1261  * @param in_char_array either NULL or contains char array with read input
     1184 * @param _initvalues struct with all state variables
    12621185 * @param argc argument count
    12631186 * @param argv argument array
    1264  */
    1265 void mainFunction(
    1266     Ref<MessageGrp> grp,
    1267     struct OptionValues &values,
    1268     const char *&output,
    1269     const char *&input,
    1270     const char *&generic_input,
    1271     char *&in_char_array,
     1187 * @param data void ptr to extract structure
     1188 */
     1189void mpqc::mainFunction(
     1190    mpqc::InitValues &_initvalues,
    12721191    int argc,
    1273     char **argv
    1274 #ifdef HAVE_MPQCDATA
    1275     , MPQCData &data
    1276 #endif
     1192    char **argv,
     1193    void *data
    12771194    )
    12781195{
    12791196  // get the basename for output files
    1280   setOutputBaseName(input, output);
     1197  setOutputBaseName(_initvalues.input, _initvalues.output);
    12811198
    12821199  // parse input into keyvalue container
    12831200  Ref<ParsedKeyVal> parsedkv;
    12841201  int use_simple_input = 0; // default is object-oriented if in_char_array is given
    1285   if (!in_char_array) // obtain from file
    1286     parseInputfile(grp, parsedkv, values, input, generic_input, in_char_array, use_simple_input);
    1287   parseIntoKeyValue(parsedkv, values, in_char_array, use_simple_input);
    1288   delete[] in_char_array;
     1202  if (!_initvalues.in_char_array) // obtain from file
     1203    parseInputfile(
     1204        _initvalues.grp,
     1205        parsedkv,
     1206        _initvalues.values,
     1207        _initvalues.input,
     1208        _initvalues.generic_input,
     1209        _initvalues.in_char_array,
     1210        use_simple_input);
     1211  parseIntoKeyValue(
     1212      parsedkv,
     1213      _initvalues.values,
     1214      _initvalues.in_char_array,
     1215      use_simple_input);
     1216//  delete[] in_char_array;
    12891217
    12901218  // prefix parsed values wit "mpqc"
    1291   if (values.keyvalue) parsedkv->verbose(1);
     1219  if (_initvalues.values.keyvalue) parsedkv->verbose(1);
    12921220  Ref<KeyVal> keyval = new PrefixKeyVal(parsedkv.pointer(),"mpqc");
    12931221
    12941222  // set up output classes
    1295   setupSCFormIO(grp);
     1223  setupSCFormIO(_initvalues.grp);
    12961224
    12971225  // initialize timing for mpqc
    12981226  Ref<RegionTimer> tim;
    1299   initTimings(grp, keyval, tim);
     1227  initTimings(_initvalues.grp, keyval, tim);
    13001228 
    13011229  // announce ourselves
    … …  
    13111239
    13121240  ExEnv::out0() << indent
    1313        << "Using " << grp->class_name()
    1314        << " for message passing (number of nodes = " << grp->n() << ")." << endl
     1241       << "Using " << _initvalues.grp->class_name()
     1242       << " for message passing (number of nodes = " << _initvalues.grp->n() << ")." << endl
    13151243       << indent
    13161244       << "Using " << thread->class_name()
    … …  
    13201248       << " for distributed shared memory." << endl
    13211249       << indent
    1322        << "Total number of processors = " << grp->n() * thread->nthread() << endl;
     1250       << "Total number of processors = " << _initvalues.grp->n() * thread->nthread() << endl;
    13231251
    13241252  // prepare CCA if available
    1325   prepareCCA(keyval, values);
     1253  prepareCCA(keyval, _initvalues.values);
    13261254
    13271255  // now set up the debugger
    13281256  Ref<Debugger> debugger;
    1329   setupDebugger(keyval, grp, debugger, values);
     1257  setupDebugger(keyval, _initvalues.grp, debugger, _initvalues.values);
    13301258
    13311259  // now check to see what matrix kit to use
    … …  
    13701298
    13711299  // read in restart file if we do restart
    1372   performRestart(keyval, grp, opt, mole, restartfile);
     1300  performRestart(keyval, _initvalues.grp, opt, mole, restartfile);
    13731301
    13741302  // setup molecule checkpoint file
    1375   setMolecularCheckpointFile(keyval, grp, mole, mole_ckpt_file);
     1303  setMolecularCheckpointFile(keyval, _initvalues.grp, mole, mole_ckpt_file);
    13761304  delete[] mole_ckpt_file;
    13771305
    … …  
    13791307  if (checkpoint && opt.nonnull()) {
    13801308    opt->set_checkpoint();
    1381     if (grp->me() == 0) opt->set_checkpoint_file(ckptfile);
     1309    if (_initvalues.grp->me() == 0) opt->set_checkpoint_file(ckptfile);
    13821310    else opt->set_checkpoint_file(devnull);
    13831311  }
    … …  
    13901318
    13911319  // check basis set limit
    1392   const int check = checkBasisSetLimit(mole, values);
     1320  const int check = checkBasisSetLimit(mole, _initvalues.values);
    13931321  if (check) {
    13941322    ExEnv::out0() << endl << indent
    … …  
    14261354  // requiring optimized geometries, etc.
    14271355  if (renderer.null()) {
    1428     if (parsedkv.nonnull()) print_unseen(parsedkv, input);
     1356    if (parsedkv.nonnull()) print_unseen(parsedkv, _initvalues.input);
    14291357    keyval = 0;
    14301358    parsedkv = 0;
    … …  
    14651393
    14661394  // save this before doing the frequency stuff since that obsoletes the
    1467   saveState(wfn_file, savestate, opt, grp, mole, molname, ckptfile);
     1395  saveState(wfn_file, savestate, opt, _initvalues.grp, mole, molname, ckptfile);
    14681396 
    14691397  // Frequency calculation.
    … …  
    14731401
    14741402  if (renderer.nonnull()) {
    1475     renderObjects(renderer, keyval, tim, grp);
     1403    renderObjects(renderer, keyval, tim, _initvalues.grp);
    14761404
    14771405    Ref<MolFreqAnimate> molfreqanim;
    … …  
    14891417      mole->print(ExEnv::out0());
    14901418
    1491     saveToPdb(do_pdb, grp, mole, molname);
     1419    saveToPdb(do_pdb, _initvalues.grp, mole, molname);
    14921420
    14931421  }
    … …  
    14971425         << endl;
    14981426  }
    1499   if (parsedkv.nonnull()) print_unseen(parsedkv, input);
     1427  if (parsedkv.nonnull()) print_unseen(parsedkv, _initvalues.input);
    15001428
    15011429  // here, we may gather the results
    15021430  // we start to fill the MPQC_Data object
    1503   if (tim.nonnull()) tim->enter("gather");
    1504   {
    1505      Ref<Wavefunction> wfn;
    1506      wfn << mole;
    1507 //     ExEnv::out0() << "The number of atomic orbitals: " << wfn->ao_dimension()->n() << endl;
    1508 //     ExEnv::out0() << "The AO density matrix is ";
    1509 //     wfn->ao_density()->print(ExEnv::out0());
    1510 //     ExEnv::out0() << "The natural density matrix is ";
    1511 //     wfn->natural_density()->print(ExEnv::out0());
    1512 //     ExEnv::out0() << "The Gaussian basis is " << wfn->basis()->name() << endl;
    1513 //     ExEnv::out0() << "The Gaussians sit at the following centers: " << endl;
    1514 //     for (int nr = 0; nr< wfn->basis()->ncenter(); ++nr) {
    1515 //       ExEnv::out0() << nr << " basis function has its center at ";
    1516 //       for (int i=0; i < 3; ++i)
    1517 //           ExEnv::out0() << wfn->basis()->r(nr,i) << "\t";
    1518 //       ExEnv::out0() << endl;
    1519 //     }
    1520      // store accuracies
    1521      data.accuracy = mole->value_result().actual_accuracy();
    1522      data.desired_accuracy = mole->value_result().desired_accuracy();
    1523      // print the energy
    1524      data.energies.total = wfn->energy();
    1525      data.energies.nuclear_repulsion = wfn->nuclear_repulsion_energy();
    1526      {
    1527        CLHF *clhf = dynamic_cast<CLHF*>(wfn.pointer());
    1528        if (clhf != NULL) {
    1529          double ex, ec;
    1530          clhf->two_body_energy(ec, ex);
    1531          data.energies.electron_coulomb = ec;
    1532          data.energies.electron_exchange = ex;
    1533          clhf = NULL;
    1534        } else {
    1535          ExEnv::out0() << "INFO: There is no direct CLHF information available." << endl;
    1536          data.energies.electron_coulomb = 0.;
    1537          data.energies.electron_exchange = 0.;
    1538        }
    1539      }
    1540      SCF *scf = NULL;
    1541      {
    1542        MBPT2 *mbpt2 = dynamic_cast<MBPT2*>(wfn.pointer());
    1543        if (mbpt2 != NULL) {
    1544          data.energies.correlation = mbpt2->corr_energy();
    1545          scf = mbpt2->ref().pointer();
    1546          CLHF *clhf = dynamic_cast<CLHF*>(scf);
    1547          if (clhf != NULL) {
    1548            double ex, ec;
    1549            clhf->two_body_energy(ec, ex);
    1550            data.energies.electron_coulomb = ec;
    1551            data.energies.electron_exchange = ex;
    1552            clhf = NULL;
    1553          } else {
    1554            ExEnv::out0() << "INFO: There is no reference CLHF information available either." << endl;
    1555            data.energies.electron_coulomb = 0.;
    1556            data.energies.electron_exchange = 0.;
    1557          }
    1558          mbpt2 = 0;
    1559        } else {
    1560          ExEnv::out0() << "INFO: There is no MBPT2 information available." << endl;
    1561          data.energies.correlation = 0.;
    1562          scf = dynamic_cast<SCF*>(wfn.pointer());
    1563          if (scf == NULL)
    1564            abort();
    1565        }
    1566      }
    1567      {
    1568        // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)
    1569 
    1570        RefSymmSCMatrix t = scf->overlap();
    1571        RefSymmSCMatrix cl_dens_ = scf->ao_density();
    1572 
    1573        SCFEnergy *eop = new SCFEnergy;
    1574        eop->reference();
    1575        if (t.dim()->equiv(cl_dens_.dim())) {
    1576          Ref<SCElementOp2> op = eop;
    1577          t.element_op(op,cl_dens_);
    1578          op=0;
    1579        }
    1580        eop->dereference();
    1581 
    1582        data.energies.overlap = eop->result();
    1583 
    1584        delete eop;
    1585        t = 0;
    1586        cl_dens_ = 0;
    1587      }
    1588      {
    1589        // taken from Wavefunction::core_hamiltonian()
    1590        RefSymmSCMatrix hao(scf->basis()->basisdim(), scf->basis()->matrixkit());
    1591        hao.assign(0.0);
    1592        Ref<PetiteList> pl = scf->integral()->petite_list();
    1593        Ref<SCElementOp> hc =
    1594          new OneBodyIntOp(new SymmOneBodyIntIter(scf->integral()->kinetic(), pl));
    1595        hao.element_op(hc);
    1596        hc=0;
    1597 
    1598        RefSymmSCMatrix h(scf->so_dimension(), scf->basis_matrixkit());
    1599        pl->symmetrize(hao,h);
    1600 
    1601        // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)
    1602        RefSymmSCMatrix cl_dens_ = scf->ao_density();
    1603 
    1604        SCFEnergy *eop = new SCFEnergy;
    1605        eop->reference();
    1606        if (h.dim()->equiv(cl_dens_.dim())) {
    1607          Ref<SCElementOp2> op = eop;
    1608          h.element_op(op,cl_dens_);
    1609          op=0;
    1610        }
    1611        eop->dereference();
    1612 
    1613        data.energies.kinetic = 2.*eop->result();
    1614 
    1615        delete eop;
    1616        hao = 0;
    1617        h = 0;
    1618        cl_dens_ = 0;
    1619      }
    1620      {
    1621        // set to potential energy between nuclei and electron charge distribution
    1622        RefSymmSCMatrix hao(scf->basis()->basisdim(), scf->basis()->matrixkit());
    1623        hao.assign(0.0);
    1624        Ref<PetiteList> pl = scf->integral()->petite_list();
    1625        Ref<SCElementOp> hc =
    1626          new OneBodyIntOp(new SymmOneBodyIntIter(scf->integral()->nuclear(), pl));
    1627        hao.element_op(hc);
    1628        hc=0;
    1629 
    1630        RefSymmSCMatrix h(scf->so_dimension(), scf->basis_matrixkit());
    1631        pl->symmetrize(hao,h);
    1632 
    1633        // taken from clscf.cc: CLSCF::scf_energy() (but see also Szabo/Ostlund)
    1634        RefSymmSCMatrix cl_dens_ = scf->ao_density();
    1635 
    1636        SCFEnergy *eop = new SCFEnergy;
    1637        eop->reference();
    1638        if (h.dim()->equiv(cl_dens_.dim())) {
    1639          Ref<SCElementOp2> op = eop;
    1640          h.element_op(op,cl_dens_);
    1641          op=0;
    1642        }
    1643        eop->dereference();
    1644 
    1645        data.energies.hcore = 2.*eop->result();
    1646 
    1647        delete eop;
    1648        hao = 0;
    1649        h = 0;
    1650        cl_dens_ = 0;
    1651      }
    1652      ExEnv::out0() << "total is " << data.energies.total << endl;
    1653      ExEnv::out0() << "nuclear_repulsion is " << data.energies.nuclear_repulsion << endl;
    1654      ExEnv::out0() << "electron_coulomb is " << data.energies.electron_coulomb << endl;
    1655      ExEnv::out0() << "electron_exchange is " << data.energies.electron_exchange << endl;
    1656      ExEnv::out0() << "correlation is " << data.energies.correlation << endl;
    1657      ExEnv::out0() << "overlap is " << data.energies.overlap << endl;
    1658      ExEnv::out0() << "kinetic is " << data.energies.kinetic << endl;
    1659      ExEnv::out0() << "hcore is " << data.energies.hcore << endl;
    1660      ExEnv::out0() << "sum is " <<
    1661          data.energies.nuclear_repulsion
    1662          + data.energies.electron_coulomb
    1663          + data.energies.electron_exchange
    1664          + data.energies.correlation
    1665          + data.energies.kinetic
    1666          + data.energies.hcore
    1667          << endl;
    1668 
    1669      ExEnv::out0() << endl << indent
    1670           << scprintf("Value of the MolecularEnergy: %15.10f",
    1671                       mole->energy())
    1672           << endl;
    1673      // print the gradient
    1674      RefSCVector grad;
    1675      if (mole->gradient_result().computed()) {
    1676        grad = mole->gradient_result().result_noupdate();
    1677      }
    1678      // gradient calculation needs to be activated in the configuration
    1679      // some methods such as open shell MBPT2 do not allow for gradient calc.
    1680 //     else {
    1681 //       grad = mole->gradient();
    1682 //     }
    1683      if (grad.nonnull()) {
    1684        data.forces.resize(grad.dim()/3);
    1685        for (int j=0;j<grad.dim()/3; ++j) {
    1686          data.forces[j].resize(3, 0.);
    1687        }
    1688        ExEnv::out0() << "Gradient of the MolecularEnergy:" << std::endl;
    1689        for (int j=0;j<grad.dim()/3; ++j) {
    1690          ExEnv::out0() << "\t";
    1691          for (int i=0; i< 3; ++i) {
    1692            data.forces[j][i] = grad[3*j+i];
    1693            ExEnv::out0() << grad[3*j+i] << "\t";
    1694          }
    1695          ExEnv::out0() << endl;
    1696        }
    1697      } else {
    1698        ExEnv::out0() << "INFO: There is no gradient information available." << endl;
    1699      }
    1700      grad = NULL;
    1701 
    1702      {
    1703        // eigenvalues (this only works if we have a OneBodyWavefunction, i.e. SCF procedure)
    1704        // eigenvalues seem to be invalid for unrestricted SCF calculation
    1705        // (see UnrestrictedSCF::eigenvalues() implementation)
    1706        UnrestrictedSCF *uscf = dynamic_cast<UnrestrictedSCF*>(wfn.pointer());
    1707        if ((scf != NULL) && (uscf == NULL)) {
    1708          const double scfernergy = scf->energy();
    1709          RefDiagSCMatrix evals = scf->eigenvalues();
    1710 
    1711          ExEnv::out0() << "Eigenvalues:" << endl;
    1712          for(int i=0;i<wfn->oso_dimension(); ++i) {
    1713            data.energies.eigenvalues.push_back(evals(i));
    1714            ExEnv::out0() << i << "th eigenvalue is " << evals(i) << endl;
    1715          }
    1716        } else {
    1717            ExEnv::out0() << "INFO: There is no eigenvalue information available." << endl;
    1718        }
    1719      }
    1720      // we do sample the density only on request
    1721        {
    1722          // fill positions and charges (NO LONGER converting from bohr radii to angstroem)
    1723          const double AtomicLengthToAngstroem = 1.;//0.52917721;
    1724          data.positions.reserve(wfn->molecule()->natom());
    1725          data.atomicnumbers.reserve(wfn->molecule()->natom());
    1726          data.charges.reserve(wfn->molecule()->natom());
    1727          for (int iatom=0;iatom < wfn->molecule()->natom(); ++iatom) {
    1728            data.atomicnumbers.push_back(wfn->molecule()->Z(iatom));
    1729            double charge = wfn->molecule()->Z(iatom);
    1730            if (data.DoValenceOnly == MPQCData::DoSampleValenceOnly)
    1731              charge -= getCoreElectrons((int)charge);
    1732            data.charges.push_back(charge);
    1733            std::vector<double> pos(3, 0.);
    1734            for (int j=0;j<3;++j)
    1735              pos[j] = wfn->molecule()->r(iatom, j)*AtomicLengthToAngstroem;
    1736            data.positions.push_back(pos);
    1737          }
    1738          ExEnv::out0() << "We have "
    1739              << data.positions.size() << " positions and "
    1740              << data.charges.size() << " charges." << endl;
    1741        }
    1742        if (data.DoLongrange) {
    1743        if (data.sampled_grid.level != 0)
    1744        {
    1745          // we now need to sample the density on the grid
    1746          // 1. get max and min over all basis function positions
    1747          assert( scf->basis()->ncenter() > 0 );
    1748          SCVector3 bmin( scf->basis()->r(0,0), scf->basis()->r(0,1), scf->basis()->r(0,2) );
    1749          SCVector3 bmax( scf->basis()->r(0,0), scf->basis()->r(0,1), scf->basis()->r(0,2) );
    1750          for (int nr = 1; nr< scf->basis()->ncenter(); ++nr) {
    1751            for (int i=0; i < 3; ++i) {
    1752              if (scf->basis()->r(nr,i) < bmin(i))
    1753                bmin(i) = scf->basis()->r(nr,i);
    1754              if (scf->basis()->r(nr,i) > bmax(i))
    1755                bmax(i) = scf->basis()->r(nr,i);
    1756            }
    1757          }
    1758          ExEnv::out0() << "Basis min is at " << bmin << " and max is at " << bmax << endl;
    1759 
    1760          // 2. choose an appropriately large grid
    1761          // we have to pay attention to capture the right amount of the exponential decay
    1762          // and also to have a power of two size of the grid at best
    1763          SCVector3 boundaryV(5.);  // boundary extent around compact domain containing basis functions
    1764          bmin -= boundaryV;
    1765          bmax += boundaryV;
    1766          for (size_t i=0;i<3;++i) {
    1767            if (bmin(i) < data.sampled_grid.begin[i])
    1768              bmin(i) = data.sampled_grid.begin[i];
    1769            if (bmax(i) > data.sampled_grid.end[i])
    1770              bmax(i) = data.sampled_grid.end[i];
    1771          }
    1772          // set the non-zero window of the sampled_grid
    1773          data.sampled_grid.setWindow(bmin.data(), bmax.data());
    1774 
    1775          // for the moment we always generate a grid of full size
    1776          // (NO LONGER converting grid dimensions from angstroem to bohr radii)
    1777          const double AtomicLengthToAngstroem = 1.;//0.52917721;
    1778          SCVector3 min;
    1779          SCVector3 max;
    1780          SCVector3 delta;
    1781          size_t samplepoints[3];
    1782          // due to periodic boundary conditions, we don't need gridpoints-1 here
    1783          // TODO: in case of open boundary conditions, we need data on the right
    1784          // hand side boundary as well
    1785          const int gridpoints = data.sampled_grid.getGridPointsPerAxis();
    1786          for (size_t i=0;i<3;++i) {
    1787            min(i) = data.sampled_grid.begin_window[i]/AtomicLengthToAngstroem;
    1788            max(i) = data.sampled_grid.end_window[i]/AtomicLengthToAngstroem;
    1789            delta(i) = data.sampled_grid.getDeltaPerAxis(i)/AtomicLengthToAngstroem;
    1790            samplepoints[i] = data.sampled_grid.getWindowGridPointsPerAxis(i);
    1791          }
    1792          ExEnv::out0() << "Grid starts at " << min
    1793              << " and ends at " << max
    1794              << " with a delta of " << delta
    1795              << " to get "
    1796              << samplepoints[0] << ","
    1797              << samplepoints[1] << ","
    1798              << samplepoints[2] << " samplepoints."
    1799              << endl;
    1800          assert( data.sampled_grid.sampled_grid.size() == samplepoints[0]*samplepoints[1]*samplepoints[2]);
    1801 
    1802          // 3. sample the atomic density
    1803          const double element_volume_conversion =
    1804              1./AtomicLengthToAngstroem/AtomicLengthToAngstroem/AtomicLengthToAngstroem;
    1805          SCVector3 r = min;
    1806 
    1807          std::set<int> valence_indices;
    1808          RefDiagSCMatrix evals = scf->eigenvalues();
    1809          if (data.DoValenceOnly == MPQCData::DoSampleValenceOnly) {
    1810            // find valence orbitals
    1811 //           std::cout << "All Eigenvalues:" << std::endl;
    1812 //           for(int i=0;i<wfn->oso_dimension(); ++i)
    1813 //             std::cout << i << "th eigenvalue is " << evals(i) << std::endl;
    1814            int n_electrons = scf->nelectron();
    1815            int n_core_electrons =  wfn->molecule()->n_core_electrons();
    1816            std::set<double> evals_sorted;
    1817            {
    1818              int i=0;
    1819              double first_positive_ev = std::numeric_limits<double>::max();
    1820              for(i=0;i<wfn->oso_dimension(); ++i) {
    1821                if (evals(i) < 0.)
    1822                  evals_sorted.insert(evals(i));
    1823                else
    1824                  first_positive_ev = std::min(first_positive_ev, (double)evals(i));
    1825              }
    1826              // add the first positive for the distance
    1827              evals_sorted.insert(first_positive_ev);
    1828            }
    1829            std::set<double> evals_distances;
    1830            std::set<double>::const_iterator advancer = evals_sorted.begin();
    1831            std::set<double>::const_iterator iter = advancer++;
    1832            for(;advancer != evals_sorted.end(); ++advancer,++iter)
    1833              evals_distances.insert((*advancer)-(*iter));
    1834            const double largest_distance = *(evals_distances.rbegin());
    1835            ExEnv::out0()  << "Largest distance between EV is " << largest_distance << std::endl;
    1836            advancer = evals_sorted.begin();
    1837            iter = advancer++;
    1838            for(;advancer != evals_sorted.begin(); ++advancer,++iter)
    1839              if (fabs(fabs((*advancer)-(*iter)) - largest_distance) < 1e-10)
    1840                break;
    1841            assert( advancer != evals_sorted.begin() );
    1842            const double last_core_ev = (*iter);
    1843            ExEnv::out0()  << "Last core EV might be " << last_core_ev << std::endl;
    1844            ExEnv::out0()  << "First valence index is " << n_core_electrons/2 << std::endl;
    1845            for(int i=n_core_electrons/2;i<wfn->oso_dimension(); ++i)
    1846              if (evals(i) > last_core_ev)
    1847                valence_indices.insert(i);
    1848 //           {
    1849 //             int i=0;
    1850 //             std::cout << "Valence eigenvalues:" << std::endl;
    1851 //             for (std::set<int>::const_iterator iter = valence_indices.begin();
    1852 //                 iter != valence_indices.end(); ++iter)
    1853 //               std::cout << i++ << "th eigenvalue is " << (*iter) << std::endl;
    1854 //           }
    1855          } else {
    1856            // just insert all indices
    1857            for(int i=0;i<wfn->oso_dimension(); ++i)
    1858              valence_indices.insert(i);
    1859          }
    1860 
    1861          // testing alternative routine from SCF::so_density()
    1862          RefSCMatrix oso_vector = scf->oso_eigenvectors();
    1863          RefSCMatrix vector = scf->so_to_orthog_so().t() * oso_vector;
    1864          oso_vector = 0;
    1865          RefSymmSCMatrix occ(scf->oso_dimension(), scf->basis_matrixkit());
    1866          occ.assign(0.0);
    1867          for (std::set<int>::const_iterator iter = valence_indices.begin();
    1868              iter != valence_indices.end(); ++iter) {
    1869            const int i = *iter;
    1870            occ(i,i) = scf->occupation(i);
    1871            ExEnv::out0() << "# " << i << " has ev of " << evals(i) << ", occupied by " << scf->occupation(i) << std::endl;
    1872          }
    1873          RefSymmSCMatrix d2(scf->so_dimension(), scf->basis_matrixkit());
    1874          d2.assign(0.0);
    1875          d2.accumulate_transform(vector, occ);
    1876 
    1877          // taken from scf::density()
    1878          RefSCMatrix nos
    1879              = scf->integral()->petite_list()->evecs_to_AO_basis(scf->natural_orbitals());
    1880          RefDiagSCMatrix nd = scf->natural_density();
    1881          GaussianBasisSet::ValueData *valdat
    1882            = new GaussianBasisSet::ValueData(scf->basis(), scf->integral());
    1883          std::vector<double>::iterator griditer = data.sampled_grid.sampled_grid.begin();
    1884          const int nbasis = scf->basis()->nbasis();
    1885          double *bs_values = new double[nbasis];
    1886 
    1887          // TODO: need to take care when we have periodic boundary conditions.
    1888          for (size_t x = 0; x < samplepoints[0]; ++x, r.x() += delta(0)) {
    1889            std::cout << "Sampling now for x=" << r.x() << std::endl;
    1890            for (size_t y = 0; y < samplepoints[1]; ++y, r.y() += delta(1)) {
    1891              for (size_t z = 0; z < samplepoints[2]; ++z, r.z() += delta(2)) {
    1892                scf->basis()->values(r,valdat,bs_values);
    1893 
    1894                // loop over natural orbitals adding contributions to elec_density
    1895                double elec_density=0.0;
    1896                for (int i=0; i<nbasis; ++i) {
    1897                    double tmp = 0.0;
    1898                    for (int j=0; j<nbasis; ++j) {
    1899                        tmp += d2(j,i)*bs_values[j]*bs_values[i];
    1900                      }
    1901                    elec_density += tmp;
    1902                  }
    1903                const double dens_at_r = elec_density * element_volume_conversion;
    1904   //             const double dens_at_r = scf->density(r) * element_volume_conversion;
    1905 
    1906   //             if (fabs(dens_at_r) > 1e-4)
    1907   //               std::cout << "Electron density at " << r << " is " << dens_at_r << std::endl;
    1908                if (griditer != data.sampled_grid.sampled_grid.end())
    1909                  *griditer++ = dens_at_r;
    1910                else
    1911                  std::cerr << "PAST RANGE!" << std::endl;
    1912              }
    1913              r.z() = min.z();
    1914            }
    1915            r.y() = min.y();
    1916          }
    1917          delete[] bs_values;
    1918          delete valdat;
    1919          assert( griditer == data.sampled_grid.sampled_grid.end());
    1920          // normalization of electron charge to equal electron number
    1921          {
    1922            double integral_value = 0.;
    1923            const double volume_element = pow(AtomicLengthToAngstroem,3)*delta(0)*delta(1)*delta(2);
    1924            for (std::vector<double>::const_iterator diter = data.sampled_grid.sampled_grid.begin();
    1925                diter != data.sampled_grid.sampled_grid.end(); ++diter)
    1926                integral_value += *diter;
    1927            integral_value *= volume_element;
    1928            int n_electrons = scf->nelectron();
    1929            if (data.DoValenceOnly == MPQCData::DoSampleValenceOnly)
    1930              n_electrons -= wfn->molecule()->n_core_electrons();
    1931            const double normalization =
    1932                ((integral_value == 0) || (n_electrons == 0)) ?
    1933                    1. : n_electrons/integral_value;
    1934            std::cout << "Created " << data.sampled_grid.sampled_grid.size() << " grid points"
    1935                << " with integral value of " << integral_value
    1936                << " against " << ((data.DoValenceOnly == MPQCData::DoSampleValenceOnly) ? "n_valence_electrons" : "n_electrons")
    1937                << " of " << n_electrons << "." << std::endl;
    1938            // with normalization we also get the charge right : -1.
    1939            for (std::vector<double>::iterator diter = data.sampled_grid.sampled_grid.begin();
    1940                diter != data.sampled_grid.sampled_grid.end(); ++diter)
    1941                *diter *= -1.*normalization;
    1942          }
    1943        }
    1944      }
    1945      scf = 0;
    1946   }
    1947   if (tim.nonnull()) tim->exit("gather");
     1431  if (tim.nonnull()) tim->enter("extract");
     1432  extractResults(mole, data);
     1433  if (tim.nonnull()) tim->exit("extract");
    19481434
    19491435  if (print_timings)
    19501436    if (tim.nonnull()) tim->print(ExEnv::out0());
    19511437
    1952 
    1953   {
    1954     // times obtain from key "mpqc" which should be the first
    1955     const int nregion = tim->nregion();
    1956     //std::cout << "There are " << nregion << " timed regions." << std::endl;
    1957     const char **region_names = new const char*[nregion];
    1958     tim->get_region_names(region_names);
    1959     // find "gather"
    1960     size_t gather_region = findTimerRegion(nregion, region_names, "gather");
    1961     size_t mpqc_region = findTimerRegion(nregion, region_names, "mpqc");
    1962     delete[] region_names;
    1963 
    1964     // get timings
    1965     double *cpu_time = new double[nregion];
    1966     double *wall_time = new double[nregion];
    1967     double *flops = new double[nregion];
    1968     tim->get_cpu_times(cpu_time);
    1969     tim->get_wall_times(wall_time);
    1970     tim->get_flops(flops);
    1971     if (cpu_time != NULL) {
    1972       data.times.total_cputime = cpu_time[mpqc_region];
    1973       data.times.gather_cputime = cpu_time[gather_region];
    1974     }
    1975     if (wall_time != NULL) {
    1976       data.times.total_walltime = wall_time[mpqc_region];
    1977       data.times.gather_walltime = wall_time[gather_region];
    1978     }
    1979     if (flops != NULL) {
    1980       data.times.total_flops = flops[mpqc_region];
    1981       data.times.gather_flops = flops[gather_region];
    1982     }
    1983     delete[] cpu_time;
    1984     delete[] wall_time;
    1985     delete[] flops;
    1986   }
     1438  extractTimings(tim, data);
    19871439
    19881440  delete[] molname;
    … …  
    20011453  parsedkv = 0;
    20021454  memory = 0;
    2003   clean_up();
     1455  detail::clean_up();
    20041456
    20051457#if defined(HAVE_TIME) && defined(HAVE_CTIME)
    … …  
    20161468
    20171469// static values object
    2018 OptionValues values;
    2019 
    2020 #ifdef HAVE_JOBMARKET
    2021 FragmentResult::ptr MPQCJob::Work()
    2022 {
    2023   char mpqc[] = "mpqc" ;
    2024   char **argv = new char*[1];
    2025   argv[0] = &mpqc[0];
    2026   int argc = 1;
    2027 //  init();
    2028 //
    2029 //  ExEnv::init(argc, argv);
    2030 //
    2031 //  // parse commandline options
    2032 //  GetLongOpt options;
    2033 //  int optind = ParseOptions(options, argc, argv);
    2034 //  const char *output = 0;
    2035 //  ostream *outstream = 0;
    2036 //  ComputeOptions(options, output, outstream);
    2037 //  OptionValues values;
    2038 //  parseRemainderOptions(options, values, argc, argv);
    2039 //
    2040 //  // get the basename for output files
    2041 //  char filename_template[] = "mpqc_temp_XXXXXX";
    2042 //  output = mktemp(filename_template);
    2043 //  setOutputBaseName(NULL, output);
    2044 
    2045   // now comes the actual work
    2046   int nfilebase = (int) inputfile.length();
    2047   char *in_char_array = new char[nfilebase + 1];
    2048   strncpy(in_char_array, inputfile.c_str(), nfilebase);
    2049   in_char_array[nfilebase] = '\0';
    2050   const char *input = 0;
    2051   const char *generic_input = 0;
    2052   Ref<MessageGrp> grp = MessageGrp::get_default_messagegrp();
    2053   // create unique, temporary name and check whether it exists
    2054   const char *output = NULL;
    2055   char *tempfilename = NULL;
     1470detail::OptionValues values;
     1471
     1472
     1473namespace mpqc {
     1474
     1475  InitValues::InitValues() :
     1476      input(NULL),
     1477      object_input(NULL),
     1478      generic_input(NULL),
     1479      output(NULL),
     1480      outstream(NULL),
     1481      in_char_array(NULL),
     1482      values(::values)
     1483  {};
     1484
     1485  void init(InitValues &_initvalues, int argc, char *argv[]) {
     1486
     1487#if !defined(DEFAULT_MPI) && !defined(ALWAYS_USE_MPI) && defined(HAVE_MPI)
     1488    MPI_Init(&argc, &argv);
     1489#endif
     1490
     1491    ::init();
     1492
     1493    ExEnv::init(argc, argv);
     1494
     1495  }
     1496
     1497  void parseOptions(InitValues &_initvalues, int argc, char *argv[])
    20561498  {
    2057     std::ifstream test;
    2058     do {
    2059       if (output != NULL) // free buffer from last round
    2060         delete output;
    2061       char filename_template[] = "mpqc_temp_XXXXXX\0";
    2062       char filename_suffix[] = ".in\0";
    2063       tempfilename = (char *) malloc ( (strlen(filename_template)+strlen(filename_suffix)+2)*(sizeof(char)));
    2064       strncpy(tempfilename, mktemp(filename_template), strlen(filename_template));
    2065       tempfilename[strlen(filename_template)] = '\0';
    2066       strncat(tempfilename, filename_suffix, strlen(filename_suffix));
    2067       output = tempfilename;
    2068       //free(tempfilename); // don't free! output takes over pointer!
    2069       test.open(output);
    2070     } while (test.good()); // test whether file does not(!) exist
    2071     test.close();
    2072   }
    2073   // put info how to sample the density into MPQCData
    2074   MPQCData data(grid);
    2075   data.DoLongrange = DoLongrange; // set whether we sample the density
    2076   data.DoValenceOnly = DoValenceOnly; // set whether we sample just valence electron and nuclei densities
    2077 // now call work horse
    2078   try {
    2079     mainFunction(grp, values, output, input, generic_input, in_char_array, argc, argv, data);
    2080   }
    2081   catch (SCException &e) {
    2082     cout << argv[0] << ": ERROR: SC EXCEPTION RAISED:" << endl
    2083          << e.what()
    2084          << endl;
    2085     clean_up();
    2086   }
    2087 
    2088   //delete[] in_char_array; // is deleted in mainFunction()
    2089   if (output != 0) {
    2090     free(tempfilename);
    2091   }
    2092   delete[] argv;
    2093   grp = NULL;
    2094 
    2095   // place into returnstream
    2096   std::stringstream returnstream;
    2097   boost::archive::text_oarchive oa(returnstream);
    2098   oa << data;
    2099 
    2100   FragmentResult::ptr s( new FragmentResult(getId(), returnstream.str()) );
    2101   if (s->exitflag != 0)
    2102     cerr << "Job #" << s->getId() << " failed to reach desired accuracy." << endl;
    2103 
    2104   return s;
    2105 }
    2106 #endif
    2107 
    2108 // we need to explicitly instantiate the serialization functions as
    2109 // its is only serialized through its base class FragmentJob
    2110 BOOST_CLASS_EXPORT_IMPLEMENT(MPQCJob)
     1499    // parse commandline options
     1500    int optind = ParseOptions(_initvalues.options, argc, argv);
     1501    ComputeOptions(_initvalues.options, _initvalues.output, _initvalues.outstream);
     1502    parseRemainderOptions(_initvalues.options, _initvalues.values, argc, argv);
     1503
     1504    // get input file names, either object-oriented or generic
     1505    getInputFileNames(_initvalues.object_input, _initvalues.generic_input, _initvalues.options, optind, argc, argv);
     1506    if (_initvalues.object_input)
     1507      _initvalues.input = _initvalues.object_input;
     1508    if (_initvalues.generic_input)
     1509      _initvalues.input = _initvalues.generic_input;
     1510
     1511    // get the message group.  first try the commandline and environment
     1512    getMessageGroup(_initvalues.grp, argc, argv);
     1513  }
     1514
     1515  void cleanup(InitValues &_initvalues) {
     1516    if (_initvalues.output != 0) {
     1517      ExEnv::set_out(&cout);
     1518      delete _initvalues.outstream;
     1519    }
     1520
     1521    _initvalues.grp = 0;
     1522    final_clean_up();
     1523
     1524#if !defined(DEFAULT_MPI) && !defined(ALWAYS_USE_MPI) && defined(HAVE_MPI)
     1525    // there is an MPI_Finalize() hidden somewhere else
     1526//    MPI_Finalize();
     1527#endif
     1528  }
     1529
     1530} /* namespace mpqc */
     1531
     1532namespace detail {
    21111533
    21121534int
    21131535try_main(int argc, char *argv[])
    21141536{
    2115   init();
    2116 
    2117   ExEnv::init(argc, argv);
    2118 
    2119   // parse commandline options
    2120   GetLongOpt options;
    2121   int optind = ParseOptions(options, argc, argv);
    2122   const char *output = 0;
    2123   ostream *outstream = 0;
    2124   ComputeOptions(options, output, outstream);
    2125   parseRemainderOptions(options, values, argc, argv);
    2126 
    2127   // get input file names, either object-oriented or generic
    2128   const char *object_input = 0;
    2129   const char *generic_input = 0;
    2130   getInputFileNames(object_input, generic_input, options, optind, argc, argv);
    2131   const char *input = 0;
    2132   if (object_input) input = object_input;
    2133   if (generic_input) input = generic_input;
    2134 
    2135   // get the message group.  first try the commandline and environment
    2136   Ref<MessageGrp> grp;
    2137   getMessageGroup(grp, argc, argv);
     1537  mpqc::InitValues initvalues;
     1538
     1539  mpqc::init(initvalues, argc, argv);
     1540
     1541  mpqc::parseOptions(initvalues, argc, argv);
    21381542
    21391543  // check if we got option "-n"
    21401544  int exitflag = 0;
    2141 #ifdef HAVE_JOBMARKET
    2142   if (options.retrieve("n")) {
    2143     /// create new argc, argv and call by splitting with tokenizer
    2144     std::string networkstring(options.retrieve("n"));
    2145     typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
    2146     boost::char_separator<char> sep(" ");
    2147     tokenizer tok(networkstring, sep);
    2148     // simple count because tokenizer has now size()
    2149     int argc_network = 0;
    2150     for(tokenizer::iterator beg=tok.begin(); beg!=tok.end();++beg)
    2151       ++argc_network;
    2152     // argv[0] is program name
    2153     char **argv_network = new char*[argc_network+1];
    2154     argv_network[0] = new char[5];
    2155     strcpy(argv_network[0], "mpqc");
    2156     // then we place each token as a new argument
    2157     argc_network = 1;
    2158     for(tokenizer::iterator beg=tok.begin(); beg!=tok.end();++beg){
    2159       const size_t strlength = (*beg).length();
    2160       const char *strarray = beg->c_str();
    2161       //std::cout << "Token " << argc_network << " is " << strarray << ", length " << strlength << endl;
    2162       argv_network[argc_network] = new char[strlength+1];
    2163       strcpy(argv_network[argc_network], strarray);
    2164       argv_network[argc_network][strlength] = '\0';
    2165       for (size_t index = 0; index < strlength; ++index)
    2166         if (argv_network[argc_network][index] == '+')
    2167           argv_network[argc_network][index] = '-';
    2168       ++argc_network;
    2169     }
    2170     /// and start listening for MPQCJobs
    2171     exitflag = poolworker_main(argc_network, argv_network);
    2172     /// remove the artifical [argv,argv] again
    2173     for (int i=0;i<argc_network;++i)
    2174       delete[] argv_network[i];
    2175     delete[] argv_network;
    2176   } else
    2177 #endif
    2178   {
    2179     // if not, work on the command line input
    2180     char *in_char_array = 0;
    2181 #ifdef HAVE_MPQCDATA
    2182     MPQCData data;
    2183     mainFunction(grp, values, output, input, generic_input, in_char_array, argc, argv, data);
    2184 #else
    2185     mainFunction(grp, values, output, input, generic_input, in_char_array, argc, argv);
    2186 #endif
    2187   }
    2188 
    2189   if (output != 0) {
    2190     ExEnv::set_out(&cout);
    2191     delete outstream;
    2192   }
    2193 
    2194   grp = 0;
    2195   final_clean_up();
     1545  void *data = NULL;
     1546    mpqc::mainFunction(
     1547        initvalues,
     1548        argc, argv,
     1549        data);
     1550
     1551  mpqc::cleanup(initvalues);
    21961552
    21971553  return exitflag;
    21981554}
    21991555
     1556} /* namespace detail */
    22001557
    22011558double EvaluateDensity(SCVector3 &r, Ref<Integral> &intgrl, GaussianBasisSet::ValueData &vdat, Ref<Wavefunction> &wfn)
    … …  
    22121569}
    22131570
    2214 int
    2215 main(int argc, char *argv[])
    2216 {
    2217 #if !defined(DEFAULT_MPI) && !defined(ALWAYS_USE_MPI) && defined(HAVE_MPI)
    2218   MPI_Init(&argc, &argv);
    2219 #endif
    2220   size_t exitflag = 0;
    2221   try {
    2222     exitflag = try_main(argc, argv);
    2223   }
    2224   catch (SCException &e) {
    2225     cout << argv[0] << ": ERROR: SC EXCEPTION RAISED:" << endl
    2226          << e.what()
    2227          << endl;
    2228     clean_up();
    2229     throw;
    2230   }
    2231   catch (bad_alloc &e) {
    2232     cout << argv[0] << ": ERROR: MEMORY ALLOCATION FAILED:" << endl
    2233          << e.what()
    2234          << endl;
    2235     clean_up();
    2236     throw;
    2237   }
    2238   catch (exception &e) {
    2239     cout << argv[0] << ": ERROR: EXCEPTION RAISED:" << endl
    2240          << e.what()
    2241          << endl;
    2242     clean_up();
    2243     throw;
    2244   }
    2245   catch (...) {
    2246     cout << argv[0] << ": ERROR: UNKNOWN EXCEPTION RAISED" << endl;
    2247     clean_up();
    2248     throw;
    2249   }
    2250 #if !defined(DEFAULT_MPI) && !defined(ALWAYS_USE_MPI) && defined(HAVE_MPI)
    2251   // there is an MPI_Finalize() hidden somewhere else
    2252 //  MPI_Finalize();
    2253 #endif
    2254   return exitflag;
    2255 }
    2256 
    22571571/////////////////////////////////////////////////////////////////////////////
    22581572
Note: See TracChangeset for help on using the changeset viewer.