mpqc.cc

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

    debugger->set_exec(argv[0]);
    debugger->set_prefix(grp->me());
    if (options.retrieve("d"))
      debugger->debug("Starting debugger because -d given on command line.");
  }

  // now check to see what matrix kit to use
  if (keyval->exists("matrixkit"))
    SCMatrixKit::set_default_matrixkit(
      dynamic_cast<SCMatrixKit*>(
        keyval->describedclassvalue("matrixkit").pointer()));

  // get the integral factory. first try commandline and environment
  Ref<Integral> integral = Integral::initial_integral(argc, argv);
  if (integral.nonnull())
    Integral::set_default_integral(integral);
  else
    integral = Integral::get_default_integral();
  ExEnv::out0() << endl << indent
       << "Using " << integral->class_name()
       << " by default for molecular integrals evaluation" << endl << endl;
  
  // check for a molecular energy and optimizer
  KeyValValueString molnamedef(basename);
  char * molname = keyval->pcharvalue("filename", molnamedef);
  if (strcmp(molname, basename))
    SCFormIO::set_default_basename(molname);

  char * ckptfile = new char[strlen(molname)+6];
  sprintf(ckptfile,"%s.ckpt",molname);
  
  KeyValValueString restartfiledef(ckptfile);
  char * restartfile = keyval->pcharvalue("restart_file", restartfiledef);
  
  char * wfn_file = keyval->pcharvalue("wfn_file");
  if (wfn_file == 0) {
    wfn_file = new char[strlen(molname)+6];
    sprintf(wfn_file,"%s.wfn",molname);
  }
  char *mole_ckpt_file = new char[strlen(wfn_file)+1];
  sprintf(mole_ckpt_file,"%s",wfn_file);

  int restart = keyval->booleanvalue("restart",truevalue);

  int checkpoint = keyval->booleanvalue("checkpoint",truevalue);
  int checkpoint_freq = keyval->intvalue("checkpoint_freq",KeyValValueint(1));
  
  int savestate = keyval->booleanvalue("savestate",truevalue);

  struct stat sb;
  Ref<MolecularEnergy> mole;
  Ref<Optimize> opt;

  int statresult, statsize;
  if (restart) {
    if (grp->me() == 0) {
      statresult = stat(restartfile,&sb);
      statsize = (statresult==0) ? sb.st_size : 0;
    }
    grp->bcast(statresult);
    grp->bcast(statsize);
  }
  if (restart && statresult==0 && statsize) {
    BcastStateInBin si(grp,restartfile);
    if (keyval->exists("override")) {
      si.set_override(new PrefixKeyVal(keyval,"override"));
    }
    char *suf = strrchr(restartfile,'.');
    if (!strcmp(suf,".wfn")) {
      mole << SavableState::key_restore_state(si,"mole");
      ExEnv::out0() << endl
                   << indent << "Restored <" << mole->class_name()
                   << "> from " << restartfile << endl;

      opt << keyval->describedclassvalue("opt");
      if (opt.nonnull())
        opt->set_function(mole.pointer());
    }
    else {
      opt << SavableState::key_restore_state(si,"opt");
      if (opt.nonnull()) {
        mole << opt->function();
        ExEnv::out0() << endl << indent
             << "Restored <Optimize> from " << restartfile << endl;
      }
    }
  } else {
    mole << keyval->describedclassvalue("mole");
    opt << keyval->describedclassvalue("opt");
  }

  if (mole.nonnull()) {
    MolecularFormula mf(mole->molecule());
    ExEnv::out0() << endl << indent
         << "Molecular formula " << mf.formula() << endl;
    if (checkpoint) {
      mole->set_checkpoint();
      if (grp->me() == 0) mole->set_checkpoint_file(mole_ckpt_file);
      else mole->set_checkpoint_file(devnull);
      mole->set_checkpoint_freq(checkpoint_freq);
    }
  }

  if (checkpoint && opt.nonnull()) {
    opt->set_checkpoint();
    if (grp->me() == 0) opt->set_checkpoint_file(ckptfile);
    else opt->set_checkpoint_file(devnull);
  }

  // see if frequencies are wanted

  Ref<MolecularHessian> molhess;
  molhess << keyval->describedclassvalue("hess");
  Ref<MolecularFrequencies> molfreq;
  molfreq << keyval->describedclassvalue("freq");
  
  int check = (options.retrieve("c") != 0);
  int limit = atoi(options.retrieve("l"));
  if (limit) {
    Ref<Wavefunction> wfn; wfn << mole;
    if (wfn.nonnull() && wfn->ao_dimension()->n() > limit) {
      ExEnv::out0() << endl << indent
           << "The limit of " << limit << " basis functions has been exceeded."
           << endl;
      check = 1;
    }
  }

  if (check) {
    ExEnv::out0() << endl << indent
         << "Exiting since the check option is on." << endl;
    exit(0);
  }
  
  if (tim.nonnull()) tim->change("calc");

  int do_energy = keyval->booleanvalue("do_energy",truevalue);
  
  int do_grad = keyval->booleanvalue("do_gradient",falsevalue);

  int do_opt = keyval->booleanvalue("optimize",truevalue);
  
  int do_pdb = keyval->booleanvalue("write_pdb",falsevalue);
  
  int print_mole = keyval->booleanvalue("print_mole",truevalue);
  
  int print_timings = keyval->booleanvalue("print_timings",truevalue);
  
  // sanity checks for the benefit of reasonable looking output
  if (opt.null()) do_opt=0;
  
  ExEnv::out0() << endl << indent
       << "MPQC options:" << endl << incindent
       << indent << "matrixkit     = <"
       << SCMatrixKit::default_matrixkit()->class_name() << ">" << endl
       << indent << "filename      = " << molname << endl
       << indent << "restart_file  = " << restartfile << endl
       << indent << "restart       = " << (restart ? "yes" : "no") << endl
       << indent << "checkpoint    = " << (checkpoint ? "yes" : "no") << endl
       << indent << "savestate     = " << (savestate ? "yes" : "no") << endl
       << indent << "do_energy     = " << (do_energy ? "yes" : "no") << endl
       << indent << "do_gradient   = " << (do_grad ? "yes" : "no") << endl
       << indent << "optimize      = " << (do_opt ? "yes" : "no") << endl
       << indent << "write_pdb     = " << (do_pdb ? "yes" : "no") << endl
       << indent << "print_mole    = " << (print_mole ? "yes" : "no") << endl
       << indent << "print_timings = " << (print_timings ? "yes" : "no")
       << endl << decindent;

  delete[] restartfile;
  delete[] ckptfile;
  
  int ready_for_freq = 1;
  if (mole.nonnull()) {
    if (((do_opt && opt.nonnull()) || do_grad)
        && !mole->gradient_implemented()) {
      ExEnv::out0() << indent
           << "WARNING: optimization or gradient requested but the given"
           << endl
           << "         MolecularEnergy object cannot do gradients."
           << endl;
    }

    if (do_opt && opt.nonnull() && mole->gradient_implemented()) {
      int result = opt->optimize();
      if (result) {
        ExEnv::out0() << indent
             << "The optimization has converged." << endl << endl;
        ExEnv::out0() << indent
             << scprintf("Value of the MolecularEnergy: %15.10f",
                         mole->energy())
             << endl << endl;
      } else {
        ExEnv::out0() << indent
             << "The optimization has NOT converged." << endl << endl;
        ready_for_freq = 0;
      }
    } else if (do_grad && mole->gradient_implemented()) {
      mole->do_gradient(1);
      ExEnv::out0() << endl << indent
           << scprintf("Value of the MolecularEnergy: %15.10f",
                       mole->energy())
           << endl;
      if (mole->value_result().actual_accuracy()
          > mole->value_result().desired_accuracy()) {
        ExEnv::out0() << indent
             << "WARNING: desired accuracy not achieved in energy" << endl;
      }
      ExEnv::out0() << endl;
      // Use result_noupdate since the energy might not have converged
      // to the desired accuracy in which case grabbing the result will
      // start up the calculation again.  However the gradient might
      // not have been computed (if we are restarting and the gradient
      // isn't in the save file for example).
      RefSCVector grad;
      if (mole->gradient_result().computed()) {
        grad = mole->gradient_result().result_noupdate();
      }
      else {
        grad = mole->gradient();
      }
      if (grad.nonnull()) {
        grad.print("Gradient of the MolecularEnergy:");
        if (mole->gradient_result().actual_accuracy()
            > mole->gradient_result().desired_accuracy()) {
          ExEnv::out0() << indent
               << "WARNING: desired accuracy not achieved in gradient" << endl;
        }
      }
    } else if (do_energy && mole->value_implemented()) {
      ExEnv::out0() << endl << indent
           << scprintf("Value of the MolecularEnergy: %15.10f",
                       mole->energy())
           << endl << endl;
    }
  }

  if (tim.nonnull()) tim->exit("calc");

  // save this before doing the frequency stuff since that obsoletes the
  // function stuff
  if (savestate) {
    if (opt.nonnull()) {
      if (grp->me() == 0) {
        ckptfile = new char[strlen(molname)+6];
        sprintf(ckptfile,"%s.ckpt",molname);
      }
      else {
        ckptfile = new char[strlen(devnull)+1];
        strcpy(ckptfile, devnull);
      }

      StateOutBin so(ckptfile);
      SavableState::save_state(opt.pointer(),so);
      so.close();

      delete[] ckptfile;
    }

    if (mole.nonnull()) {
      if (grp->me() == 0) {
        if (wfn_file == 0) {
          wfn_file = new char[strlen(molname)+6];
          sprintf(wfn_file,"%s.wfn",molname);
        }
      }
      else {
        delete[] wfn_file;
        wfn_file = new char[strlen(devnull)+1];
        strcpy(wfn_file, devnull);
      }
  
      StateOutBin so(wfn_file);
      SavableState::save_state(mole.pointer(),so);
      so.close();

      delete[] wfn_file;
    }
  }

  // Frequency calculation.
  if (ready_for_freq && molfreq.nonnull()) {
    RefSymmSCMatrix xhessian;
    if (molhess.nonnull()) {
      // if "hess" input was given, use it to compute the hessian
      xhessian = molhess->cartesian_hessian();
    }
    else if (mole->hessian_implemented()) {
      // if mole can compute the hessian, use that hessian
      xhessian = mole->get_cartesian_hessian();
    }
    else if (mole->gradient_implemented()) {
      // if mole can compute gradients, use gradients at finite
      // displacements to compute the hessian
      molhess = new FinDispMolecularHessian(mole);
      xhessian = molhess->cartesian_hessian();
    }
    else {
      ExEnv::out0() << "mpqc: WARNING: Frequencies cannot be computed" << endl;
    }

    if (xhessian.nonnull()) {
      char *hessfile = SCFormIO::fileext_to_filename(".hess");
      MolecularHessian::write_cartesian_hessian(hessfile,
                                                mole->molecule(), xhessian);
      delete[] hessfile;

      molfreq->compute_frequencies(xhessian);
      // DEGENERACY IS NOT CORRECT FOR NON-SINGLET CASES:
      molfreq->thermochemistry(1);
    }
  }

  // see if any pictures are desired
  Ref<Render> renderer;
  renderer << keyval->describedclassvalue("renderer");
  Ref<RenderedObject> rendered;
  rendered << keyval->describedclassvalue("rendered");
  Ref<AnimatedObject> animated;
  animated << keyval->describedclassvalue("rendered");
  if (renderer.nonnull() && rendered.nonnull()) {
    if (tim.nonnull()) tim->enter("render");
    if (grp->me() == 0) renderer->render(rendered);
    if (tim.nonnull()) tim->exit("render");
  }
  else if (renderer.nonnull() && animated.nonnull()) {
    if (tim.nonnull()) tim->enter("render");
    if (grp->me() == 0) renderer->animate(animated);
    if (tim.nonnull()) tim->exit("render");
  }
  else if (renderer.nonnull()) {
    if (tim.nonnull()) tim->enter("render");
    int n = keyval->count("rendered");
    for (i=0; i<n; i++) {
      rendered << keyval->describedclassvalue("rendered",i);
      animated << keyval->describedclassvalue("rendered",i);
      if (rendered.nonnull()) {
        // make sure the object has a name so we don't overwrite its file
        if (rendered->name() == 0) {
          char ic[64];
          sprintf(ic,"%02d",i);
          rendered->set_name(ic);
        }
        if (grp->me() == 0) renderer->render(rendered);
      }
      else if (animated.nonnull()) {
        // make sure the object has a name so we don't overwrite its file
        if (animated->name() == 0) {
          char ic[64];
          sprintf(ic,"%02d",i);
          animated->set_name(ic);
        }
        if (grp->me() == 0) renderer->animate(animated);
      }
    }
    if (tim.nonnull()) tim->exit("render");
  }
  Ref<MolFreqAnimate> molfreqanim;
  molfreqanim << keyval->describedclassvalue("animate_modes");
  if (ready_for_freq && molfreq.nonnull()
      && molfreqanim.nonnull() && renderer.nonnull()) {
    if (tim.nonnull()) tim->enter("render");
    molfreq->animate(renderer, molfreqanim);
    if (tim.nonnull()) tim->exit("render");
  }

  if (mole.nonnull()) {
    if (print_mole)
      mole->print(ExEnv::out0());

    if (do_pdb && grp->me() == 0) {
      ckptfile = new char[strlen(molname)+5];
      sprintf(ckptfile, "%s.pdb", molname);
      ofstream pdbfile(ckptfile);
      mole->molecule()->print_pdb(pdbfile);
      delete[] ckptfile;
    }
    
  }
  else {
    ExEnv::out0() << "mpqc: The molecular energy object is null" << endl
         << " make sure \"mole\" specifies a MolecularEnergy derivative"
         << endl;
  }

  if (parsedkv->have_unseen()) {
    ExEnv::out0() << indent
         << "The following keywords in \"" << input << "\" were ignored:"
         << endl;
    ExEnv::out0() << incindent;
    parsedkv->print_unseen(ExEnv::out0());
    ExEnv::out0() << decindent;
    ExEnv::out0() << endl;
  }

  if (print_timings)
    if (tim.nonnull()) tim->print(ExEnv::out0());

  delete[] basename;
  delete[] molname;
  SCFormIO::set_default_basename(0);

  molfreqanim = 0;
  animated = 0;
  rendered = 0;
  renderer = 0;
  molfreq = 0;
  molhess = 0;
  opt = 0;
  mole = 0;
  integral = 0;
  debugger = 0;
  thread = 0;
  tim = 0;
  keyval = 0;
  parsedkv = 0;
  grp = 0;
  memory = 0;
  clean_up();

#if defined(HAVE_TIME) && defined(HAVE_CTIME)
  time(&t);
  tstr = ctime(&t);
#endif
  if (!tstr) {
    tstr = "UNKNOWN";
  }
  ExEnv::out0() << endl
               << indent << scprintf("End Time: %s", tstr) << endl;

  if (output != 0) {
    ExEnv::set_out(&cout);
    delete outstream;
  }

  return 0;
}

int
main(int argc, char *argv[])
{
  try {
    try_main(argc, argv);
  }
  catch (bad_alloc &e) {
    cout << argv[0] << ": ERROR: MEMORY ALLOCATION FAILED:" << endl
         << e.what()
         << endl;
    clean_up();
    throw;
  }
  catch (exception &e) {
    cout << argv[0] << ": ERROR: EXCEPTION RAISED:" << endl
         << e.what()
         << endl;
    clean_up();
    throw;
  }
  catch (...) {
    cout << argv[0] << ": ERROR: UNKNOWN EXCEPTION RAISED" << endl;
    clean_up();
    throw;
  }
  return 0;
}

/////////////////////////////////////////////////////////////////////////////

// Local Variables:
// mode: c++
// c-file-style: "ETS"
// End: